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

Tuple [ numpy.ndarray, float ]

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