# Algebraic Helpers¶

bezier._algebraic_intersection.bezier_roots(coeffs)

Compute polynomial roots from a polynomial in the Bernstein basis.

Note

This assumes the caller passes in a 1D array but does not check.

This takes the polynomial

$f(s) = \sum_{j = 0}^n b_{j, n} \cdot C_j.$

and uses the variable $$\sigma = \frac{s}{1 - s}$$ to rewrite as

\begin{split}\begin{align*} f(s) &= (1 - s)^n \sum_{j = 0}^n \binom{n}{j} C_j \sigma^j \\ &= (1 - s)^n \sum_{j = 0}^n \widetilde{C_j} \sigma^j. \end{align*}\end{split}

Then it uses an eigenvalue solver to find the roots of

$g(\sigma) = \sum_{j = 0}^n \widetilde{C_j} \sigma^j$

and convert them back into roots of $$f(s)$$ via $$s = \frac{\sigma}{1 + \sigma}$$.

For example, consider

\begin{split}\begin{align*} f_0(s) &= 2 (2 - s)(3 + s) \\ &= 12(1 - s)^2 + 11 \cdot 2s(1 - s) + 8 s^2 \end{align*}\end{split}

First, we compute the companion matrix for

$g_0(\sigma) = 12 + 22 \sigma + 8 \sigma^2$
>>> coeffs0 = np.asfortranarray([12.0, 11.0, 8.0])
>>> companion0, _, _ = bernstein_companion(coeffs0)
>>> companion0
array([[-2.75, -1.5 ],
[ 1.  ,  0.  ]])


then take the eigenvalues of the companion matrix:

>>> sigma_values0 = np.linalg.eigvals(companion0)
>>> sigma_values0
array([-2.  , -0.75])


after transforming them, we have the roots of $$f(s)$$:

>>> sigma_values0 / (1.0 + sigma_values0)
array([ 2., -3.])
>>> bezier_roots(coeffs0)
array([ 2., -3.])


In cases where $$s = 1$$ is a root, the lead coefficient of $$g$$ would be $$0$$, so there is a reduction in the companion matrix.

\begin{split}\begin{align*} f_1(s) &= 6 (s - 1)^2 (s - 3) (s - 5) \\ &= 90 (1 - s)^4 + 33 \cdot 4s(1 - s)^3 + 8 \cdot 6s^2(1 - s)^2 \end{align*}\end{split}
>>> coeffs1 = np.asfortranarray([90.0, 33.0, 8.0, 0.0, 0.0])
>>> companion1, degree1, effective_degree1 = bernstein_companion(
...     coeffs1)
>>> companion1
array([[-2.75 , -1.875],
[ 1.   ,  0.   ]])
>>> degree1
4
>>> effective_degree1
2


so the roots are a combination of the roots determined from $$s = \frac{\sigma}{1 + \sigma}$$ and the number of factors of $$(1 - s)$$ (i.e. the difference between the degree and the effective degree):

>>> bezier_roots(coeffs1)
array([3., 5., 1., 1.])


In some cases, a polynomial is represented with an “elevated” degree:

\begin{split}\begin{align*} f_2(s) &= 3 (s^2 + 1) \\ &= 3 (1 - s)^3 + 3 \cdot 3s(1 - s)^2 + 4 \cdot 3s^2(1 - s) + 6 s^3 \end{align*}\end{split}

This results in a “point at infinity” $$\sigma = -1 \Longleftrightarrow s = \infty$$:

>>> coeffs2 = np.asfortranarray([3.0, 3.0, 4.0, 6.0])
>>> companion2, _, _ = bernstein_companion(coeffs2)
>>> companion2
array([[-2. , -1.5, -0.5],
[ 1. ,  0. ,  0. ],
[ 0. ,  1. ,  0. ]])
>>> sigma_values2 = np.linalg.eigvals(companion2)
>>> sigma_values2
array([-1. +0.j , -0.5+0.5j, -0.5-0.5j])


so we drop any values $$\sigma$$ that are sufficiently close to $$-1$$:

>>> expected2 = np.asfortranarray([1.0j, -1.0j])
>>> roots2 = bezier_roots(coeffs2)
>>> np.allclose(expected2, roots2, rtol=8 * machine_eps, atol=0.0)
True

Parameters

coeffs (numpy.ndarray) – A 1D array of coefficients in the Bernstein basis.

Returns

A 1D array containing the roots.

Return type

numpy.ndarray

bezier._algebraic_intersection.lu_companion(top_row, value)

Compute an LU-factored $$C - t I$$ and its 1-norm.

Note

The output of this function is intended to be used with dgecon from LAPACK. dgecon expects both the 1-norm of the matrix and expects the matrix to be passed in an already LU-factored form (via dgetrf).

The companion matrix $$C$$ is given by the top_row, for example, the polynomial $$t^3 + 3 t^2 - t + 2$$ has a top row of -3, 1, -2 and the corresponding companion matrix is:

$\begin{split}\left[\begin{array}{c c c} -3 & 1 & -2 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \end{array}\right]\end{split}$

After doing a full cycle of the rows (shifting the first to the last and moving all other rows up), row reduction of $$C - t I$$ yields

$\begin{split}\left[\begin{array}{c c c} 1 & -t & 0 \\ 0 & 1 & -t \\ -3 - t & 1 & -2 \end{array}\right] = \left[\begin{array}{c c c} 1 & 0 & 0 \\ 0 & 1 & 0 \\ -3 - t & 1 + t(-3 - t) & 1 \end{array}\right] \left[\begin{array}{c c c} 1 & -t & 0 \\ 0 & 1 & -t \\ 0 & 0 & -2 + t(1 + t(-3 - t)) \end{array}\right]\end{split}$

and in general, the terms in the bottom row correspond to the intermediate values involved in evaluating the polynomial via Horner’s method.

>>> top_row = np.asfortranarray([-3.0, 1.0, -2.0])
>>> t_val = 0.5
>>> lu_mat, one_norm = lu_companion(top_row, t_val)
>>> lu_mat
array([[ 1.   , -0.5  ,  0.   ],
[ 0.   ,  1.   , -0.5  ],
[-3.5  , -0.75 , -2.375]])
>>> one_norm
4.5
>>> l_mat = np.tril(lu_mat, k=-1) + np.eye(3)
>>> u_mat = np.triu(lu_mat)
>>> a_mat = l_mat.dot(u_mat)
>>> a_mat
array([[ 1. , -0.5,  0. ],
[ 0. ,  1. , -0.5],
[-3.5,  1. , -2. ]])
>>> np.linalg.norm(a_mat, ord=1)
4.5

Parameters
• top_row (numpy.ndarray) – 1D array, top row of companion matrix.

• value (float) – The $$t$$ value used to form $$C - t I$$.

Returns

Pair of

• 2D array of LU-factored form of $$C - t I$$, with the non-diagonal part of $$L$$ stored in the strictly lower triangle and $$U$$ stored in the upper triangle (we skip the permutation matrix, as it won’t impact the 1-norm)

• the 1-norm the matrix $$C - t I$$

As mentioned above, these two values are meant to be used with dgecon.

Return type
bezier._algebraic_intersection.bezier_value_check(coeffs, s_val, rhs_val=0.0)

Check if a polynomial in the Bernstein basis evaluates to a value.

This is intended to be used for root checking, i.e. for a polynomial $$f(s)$$ and a particular value $$s_{\ast}$$:

Is it true that $$f\left(s_{\ast}\right) = 0$$?

Does so by re-stating as a matrix rank problem. As in bezier_roots(), we can rewrite

$f(s) = (1 - s)^n g\left(\sigma\right)$

for $$\sigma = \frac{s}{1 - s}$$ and $$g\left(\sigma\right)$$ written in the power / monomial basis. Now, checking if $$g\left(\sigma\right) = 0$$ is a matter of checking that $$\det\left(C_g - \sigma I\right) = 0$$ (where $$C_g$$ is the companion matrix of $$g$$).

Due to issues of numerical stability, we’d rather ask if $$C_g - \sigma I$$ is singular to numerical precision. A typical approach to this using the singular values (assuming $$C_g$$ is $$m \times m$$) is that the matrix is singular if

$\sigma_m < m \varepsilon \sigma_1 \Longleftrightarrow \frac{1}{\kappa_2} < m \varepsilon$

(where $$\kappa_2$$ is the 2-norm condition number of the matrix). Since we also know that $$\kappa_2 < \kappa_1$$, a stronger requirement would be

$\frac{1}{\kappa_1} < \frac{1}{\kappa_2} < m \varepsilon.$

This is more useful since it is much easier to compute the 1-norm condition number / reciprocal condition number (and methods in LAPACK are provided for doing so).

Parameters
• coeffs (numpy.ndarray) – A 1D array of coefficients in the Bernstein basis representing a polynomial.

• s_val (float) – The value to check on the polynomial: $$f(s) = r$$.

• rhs_val (Optional [ float ]) – The value to check that the polynomial evaluates to. Defaults to 0.0.

Returns

Indicates if $$f\left(s_{\ast}\right) = r$$ (where $$s_{\ast}$$ is s_val and $$r$$ is rhs_val).

Return type

bool