Numerical certification

The cornerstone of Smale's alpha theory is bounding the error for Newton's method. Smale's 1986 work[1] introduced the quantity ${\displaystyle \alpha }$, which quantifies the convergence of Newton's method. More precisely, let ${\displaystyle F}$ be a system of analytic functions in the variables ${\displaystyle x}$, ${\displaystyle D}$ the derivative operator, and ${\displaystyle N}$ the Newton operator. The quantities

$${\displaystyle \beta (f,x)=\|x-N(x)\|=\|Df(x)^{-1}f(x)\|}$$

$${\displaystyle \gamma (f,x)=\sup _{k\geq 2}\left\|{\frac {Df(x)^{-1}D^{k}f(x)}{k!}}\right\|^{\frac {1}{k-1}}}$$

and

$${\displaystyle \alpha (f,x)=\beta (f,x)\gamma (f,x)}$$

are used to certify a candidate solution. In particular, if

$${\displaystyle \alpha (f,x)<{\frac {13-3{\sqrt {17}}}{4}},}$$

then ${\displaystyle x}$ is an approximate solution for ${\displaystyle f}$, i.e., the candidate is in the domain of quadratic convergence for Newton's method. In other words, if this inequality holds, then there is a root ${\displaystyle x^{\ast }}$ of ${\displaystyle F}$ so that iterates of the Newton operator converge as

$${\displaystyle \left\|N^{k}(x)-x^{\ast }\right\|\leq {\frac {1}{2^{2^{k}-1}}}\|x-x^{\ast }\|.}$$

The software package alphaCertified provides an implementation of the alpha test for polynomials by estimating ${\displaystyle \beta }$ and ${\displaystyle \gamma }$.[2]

Suppose ${\displaystyle G:\mathbb {R} ^{n}\rightarrow \mathbb {R} ^{n}}$ is a function whose fixed points correspond to the roots of ${\displaystyle F}$. For example, the Newton operator has this property. Suppose that ${\displaystyle I}$ is a region, then,

1. If ${\displaystyle G}$ maps ${\displaystyle I}$ into itself, i.e., ${\displaystyle G(I)\subseteq I}$, then by Brouwer fixed-point theorem, ${\displaystyle G}$ has at least one fixed point in ${\displaystyle I}$, and, hence ${\displaystyle F}$ has at least one root in ${\displaystyle I}$.
2. If ${\displaystyle G}$ is contractive in a region containing ${\displaystyle I}$, then there is at most one root in ${\displaystyle I}$.

There are versions of the following methods over the complex numbers, but both the interval arithmetic and conditions must be adjusted to reflect this case.

In the univariate case, Newton's method can be directly generalized to certify a root over an interval. For an interval ${\displaystyle J}$, let ${\displaystyle m(J)}$ be the midpoint of ${\displaystyle J}$. Then, the interval Newton operator applied to ${\displaystyle J}$ is

${\displaystyle IN(J)=m(J)-F(m(J))/F'(J).}$

In practice, any interval containing ${\displaystyle F'(J)}$ can be used in this computation. If ${\displaystyle x}$ is a root of ${\displaystyle F}$, then by the mean value theorem, there is some ${\displaystyle c\in J}$ such that ${\displaystyle F(m(J))-F'(c)(m(J)-x)=F(x)=0}$. In other words, ${\displaystyle F(m(J))=F'(c)(m(J)-x)}$. Since ${\displaystyle F'(J)}$ contains the inverse of ${\displaystyle F}$ at all points of ${\displaystyle J}$, it follows that ${\displaystyle m(J)-x\in F(m(J))/F'(J)}$. Therefore, ${\displaystyle x=m(J)-(m(J)-x)\in IN(J)}$.

Furthermore, if ${\displaystyle 0\not \in F'(J)}$, then either ${\displaystyle m(J)}$ is a root of ${\displaystyle F}$ and ${\displaystyle IN(J)=\{m(J)\}}$ or ${\displaystyle m(J)\not \in IN(J)}$. Therefore, ${\displaystyle J\cap N(J)}$ is at most half the width of ${\displaystyle J}$. Therefore, if there is some root of ${\displaystyle F}$ in ${\displaystyle J}$, the iterative procedure of replacing ${\displaystyle J}$ by ${\displaystyle J\cap IN(J)}$ will converge to this root. If, on the other hand, there is no root of ${\displaystyle F}$ in ${\displaystyle J}$, this iterative procedure will eventually produce an empty interval, a witness to the nonexistence of roots.

See interval Newton method for higher dimensional analogues of this approach.

Let ${\displaystyle Y}$ be any ${\displaystyle n\times n}$ invertible matrix in ${\displaystyle GL(n,\mathbb {R} )}$. Typically, one takes ${\displaystyle Y}$ to be an approximation to ${\displaystyle F'(y)^{-1}}$. Then, define the function ${\displaystyle G(x)=x-YF(x).}$ We observe that ${\displaystyle x}$ is a fixed of ${\displaystyle G}$ if and only if ${\displaystyle x}$ is a root of ${\displaystyle F}$. Therefore the approach above can be used to identify roots of ${\displaystyle F}$. This approach is similar to a multivariate version of Newton's method, replacing the derivative with the fixed matrix ${\displaystyle Y}$.

We observe that if ${\displaystyle J}$ is a compact and convex region and ${\displaystyle y\in J}$, then, for any ${\displaystyle x\in J}$, there exist ${\displaystyle c_{1},\dots ,c_{n}\in J}$ such that

${\displaystyle G(y)-G(x)={\begin{bmatrix}\nabla g_{1}(c_{1})^{T}\\\vdots \\\nabla g_{n}(c_{n})^{T}\end{bmatrix}}(y-x).}$

Let ${\displaystyle G'(J)}$ be the Jacobian matrix of ${\displaystyle G}$ evaluated on ${\displaystyle J}$. In other words, the entry ${\displaystyle (G'(J))_{ij}}$ consists of the image of ${\displaystyle {\frac {\partial g_{i}}{\partial x_{j}}}}$ over ${\displaystyle J}$. It then follows that ${\displaystyle G(x)\in G(y)+\nabla G(J)(x-y),}$ where the matrix-vector product is computed using interval arithmetic. Then, allowing ${\displaystyle x}$ to vary in ${\displaystyle J}$, it follows that the image of ${\displaystyle G}$ on ${\displaystyle J}$ satisfies the following containment: ${\displaystyle G(J)\subset G(y)+G'(J)(J-y),}$ where the calculations are, once again, computed using interval arithmetic. Combining this with the formula for ${\displaystyle G}$, the result is the Krawczyck operator

${\displaystyle K_{y,Y}(J)=y-YF(y)+(I-F'(J))(J-y),}$

where ${\displaystyle I}$ is the identity matrix.

If ${\displaystyle K_{y,Y}(J)\subset J}$, then ${\displaystyle G}$ has a fixed point in ${\displaystyle J}$, i.e., ${\displaystyle F}$ has a root in ${\displaystyle J}$. On the other hand, if the maximum matrix norm using the supremum norm for vectors of all matrices in ${\displaystyle I-F'(J)}$ is less than ${\displaystyle 1}$, then ${\displaystyle G}$ is contractive within ${\displaystyle J}$, so ${\displaystyle G}$ has a unique fixed point.

A simpler test, when ${\displaystyle J}$ is an axis-aligned parallelepiped, uses ${\displaystyle y=m(J)}$, i.e., the midpoint of ${\displaystyle J}$. In this case, there is a unique root of ${\displaystyle F}$ if

${\displaystyle \|K(X)-m(X)\|<{\frac {w(X)}{2}},}$

where ${\displaystyle w(X)}$ is the length of the longest side of ${\displaystyle J}$.

• Miranda test (Yap, Vegter, Sharma)

Interval arithmetic can be used to provide an a priori numerical certificate by computing intervals containing unique solutions. By using intervals instead of plain numeric types during path tracking, resulting candidates are represented by intervals. The candidate solution-interval is itself the certificate, in the sense that the solution is guaranteed to be inside the interval.

Numerical algebraic geometry solves polynomial systems using homotopy continuation and path tracking methods. By monitoring the condition number for a tracked homotopy at every step, and ensuring that no two solution paths ever intersect, one can compute a numerical certificate along with a solution. This scheme is called a priori path tracking.[3]

Non-certified numerical path tracking relies on heuristic methods for controlling time step size and precision.[4] In contrast, a priori certified path tracking goes beyond heuristics to provide step size control that guarantees that for every step along the path, the current point is within the domain of quadratic convergence for the current path.