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
-
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 to0.0
.
- Returns
Indicates if \(f\left(s_{\ast}\right) = r\) (where \(s_{\ast}\) is
s_val
and \(r\) isrhs_val
).- Return type