(a) Multiple-valued functions.

We use this to refer to situations where given input data may have many possible (and equally valid) outputs. For instance, given a complex univariate polynomial $P$, compute a zero ${\it\zeta}$ of $P$. If $P$ has degree $d$ then it has at most $d$ zeros (and, generically, it has exactly $d$ zeros) and any of them is a possible output of a zero-finding algorithm. If ${\it\zeta}_{1},\ldots ,{\it\zeta}_{d}$ are these zeros, then definition (1.2) (or (1.3)) yields a measure of condition $\mathsf{acond}(P,{\it\zeta}_{i})$ for each one of them. But how does one define the condition of $P$? There is no canonical way to do so. But three possibilities arise by taking the best, worst, or average condition for $P$ over its zeros.

More generally, for an input $a$, let $\mathsf{Sol}(a)$ be its associated set of solutions (i.e. all the possible outputs for $a$). If for each $y\in \mathsf{Sol}(a)$ we have a number ${\it\xi}(a,y)$ quantifying the condition of $a$ with respect to $y$, we may define the condition ${\it\xi}(a)$ of $a$ by taking any of the following options:

$$\begin{eqnarray}{\it\xi}(a):=\inf _{y\in \mathsf{Sol}(a)}{\it\xi}(a,y),\quad {\it\xi}(a):=\underset{y\in \mathsf{Sol}(a)}{\mathbb{E}}{\it\xi}(a,y),\quad {\it\xi}(a):=\sup _{y\in \mathsf{Sol}(a)}{\it\xi}(a,y),\end{eqnarray}$$

where the expectation in the middle expression is for some distribution on $\mathsf{Sol}(A)$. Which of the three we choose will depend on the problem at hand, and sometimes on the algorithmic approach we use to solve it.

(b) Condition on manifolds.

Another way in which the setting (1.1) may be too restrictive has to do with the spaces where data and outputs live. As well as the obvious fact that ${\it\varphi}$ may be only partially defined, for a number of problems these spaces are manifolds more general than Euclidean spaces. For instance, if $({\it\lambda},v)$ is an eigenpair of a matrix $A\in \mathbb{C}^{n\times n}$, the eigenvector $v$ naturally lives in the projective space $\mathbb{P}(\mathbb{C}^{n})$, which motivates the following general framework for condition numbers.

Let $X$ be a manifold of inputs, $Y$ a manifold of outputs, and let $V\subseteq X\times Y$ be a submanifold of ‘solutions’ to some computational problem. We assume that $X$ and $V$ have the same dimension $n$ to guarantee a local uniqueness of solutions. Consider the projection ${\it\pi}_{1}:V\rightarrow X,(x,y)\mapsto x$, and its derivative $D{\it\pi}_{1}(x,y):T_{(x,y)}V\rightarrow T_{x}X$, which is the restriction of the projection

$$\begin{eqnarray}T_{x}X\times T_{y}Y\rightarrow T_{x}X,\quad ({\dot{x}},{\dot{y}})\mapsto {\dot{x}},\end{eqnarray}$$

to the subspace $T_{(x,y)}V$. We define the set of ill-posed solutions as

(1.4)$$\begin{eqnarray}{\rm\Sigma}^{\prime }:=\big\{(x,y)\in V\mid \mathsf{rank}\,D{\it\pi}_{1}(x,y)<\dim X\big\}.\end{eqnarray}$$

If $(x_{0},y_{0})\in V\setminus {\rm\Sigma}^{\prime }$, then $D{\it\pi}_{1}(x_{0},y_{0}):T_{(x_{0},y_{0})}V\rightarrow T_{x_{0}}X$ is a linear isomorphism. The implicit function theorem shows that ${\it\pi}_{1}:V\rightarrow X$ can be locally inverted around $(x_{0},y_{0})$. Its inverse $x\mapsto (x,G(x))$ is given by the solution map $G:X\supseteq U\rightarrow Y$ defined on an open neighbourhood $U$ of $x_{0}$. So we have $(x,G(x))\in V$ for all $x\in U$. The derivative

$$\begin{eqnarray}DG(x_{0}):T_{x_{0}}X\rightarrow T_{y_{0}}Y\end{eqnarray}$$

will be called the condition map. The inverse of $D{\it\pi}_{1}(x_{0},y_{0})$ is given by

(1.5)$$\begin{eqnarray}T_{x_{0}}X\rightarrow T_{(x_{0},y_{0})}V,\quad {\dot{x}}\mapsto ({\dot{x}},DG(x_{0})({\dot{x}})).\end{eqnarray}$$

A common situation is that in which $V$ is given as the zero set of a smooth map $F:X\times Y\rightarrow \mathbb{R}^{n}$. In this case we have the following characterization of ${\rm\Sigma}^{\prime }$, in terms of the function $F$ itself,

(1.6)$$\begin{eqnarray}{\rm\Sigma}^{\prime }=\big\{(x,y)\in V\mid \mathsf{rank}\,\partial F/\partial y(x,y)<n\big\}.\end{eqnarray}$$

Here the partial derivative $\partial F/\partial y(x,y):T_{y}Y\rightarrow \mathbb{R}^{n}$ is defined as the restriction of $DF(x,y)$ to $T_{y}Y$.

Now suppose that $X$ and $Y$ are Riemannian manifolds. That is, we have an inner product on each tangent space $T_{x}X$ and $T_{y}Y$ that varies smoothly with $x$ and $y$, respectively. In particular, $T_{x}X$ and $T_{y}Y$ are normed vector spaces. In this case, we may define the (absolute) normwise condition number

$$\begin{eqnarray}\mathsf{acond}^{G}(x_{0}):=\Vert DG(x_{0})\Vert =\max _{\Vert {\dot{x}}\Vert =1}\Vert DG(x_{0})({\dot{x}})\Vert\end{eqnarray}$$

as the operator norm of $DG(x_{0})$. We note that in this general framework, it does not make sense to define relative condition numbers. However, implicitly, we can model relative notions of condition numbers by choosing the manifolds appropriately. For instance, working with projective spaces means studying ratios, which accounts for a relative notion.

This framework relates to our previous discussion. For an element $x_{0}$ in the manifold of inputs $X$ we have finitely many points $(x_{0},y_{1}),\ldots ,(x_{0},y_{s})$ in the fibre ${\it\pi}_{1}^{-1}(x_{0})$. Each of them has a corresponding solution map $G_{j}$, $j=1,\ldots ,s$, and a condition number $\mathsf{acond}^{G_{j}}(x_{0})$. The condition of $x_{0}$ will depend on the computational problem we are considering associated with the geometric situation $V\subseteq X\times Y$. As described in (a) above, the three typical choices are

$$\begin{eqnarray}\displaystyle \mathsf{acond}(x_{0}) & := & \displaystyle \inf _{j\leq s}\mathsf{acond}^{G_{j}}(x_{0}),\nonumber\\ \displaystyle \mathsf{acond}(x_{0}) & := & \displaystyle \underset{j\leq s}{\mathbb{E}}\mathsf{acond}^{G_{j}}(x_{0}),\nonumber\\ \displaystyle \mathsf{acond}(x_{0}) & := & \displaystyle \sup _{j\leq s}\mathsf{acond}^{G_{j}}(x_{0}).\nonumber\end{eqnarray}$$

Note that the last two choices force us to define the set of ill-posed inputs to be the set ${\rm\Sigma}:={\it\pi}_{1}({\rm\Sigma}^{\prime })$, whereas the first choice requires

$$\begin{eqnarray}{\rm\Sigma}:=\big\{x\in X\mid (x,y)\in {\rm\Sigma}^{\prime }\text{for all }(x,y)\in {\it\pi}^{-1}(x)\big\}.\end{eqnarray}$$

(c) Finite-valued functions.

Here we consider a function ${\it\varphi}:\mathbb{R}^{m}\rightarrow Y$ where $Y$ is a finite set (say $Y=\{y_{1},\ldots ,y_{k}\}$). In this situation the quantity $\mathsf{acond}^{{\it\varphi}}(a)$ in (1.2) is of little use. It is immediate to check that $\mathsf{acond}^{{\it\varphi}}(a)=\infty$ when $a$ is in the boundary of some $S_{j}:=\{a\in \mathbb{R}^{m}\mid {\it\varphi}(a)=y_{j}\}$, for $j=1,\ldots ,k$, and that $\mathsf{acond}^{{\it\varphi}}(a)=0$ otherwise.

The union ${\rm\Sigma}$ of the boundaries between the sets $S_{j}$ is composed of data $a$ for which there is no hope that a finite-precision computation with input $a$ will yield a reliable output.

An idea championed by Jim Renegar (Reference Renegar1994b, Reference Renegar1995a, Reference Renegar1995b) is to define the condition for finite-valued problems to be the inverse of the distance to ill-posedness, that is,

$$\begin{eqnarray}\mathsf{cond}^{{\it\varphi}}(a):=\frac{1}{\mathsf{dist}(a,{\rm\Sigma})},\end{eqnarray}$$

or a relative inverse distance

(1.7)$$\begin{eqnarray}\mathsf{cond}^{{\it\varphi}}(a):=\frac{\Vert a\Vert }{\mathsf{dist}(a,{\rm\Sigma})}\end{eqnarray}$$

if the scaling properties of the problem suggest that this is better, which is often the case. Other approaches to defining conditions for these problems have been used as well.

Notable within the class of finite-valued problems is the subclass of decision problems. These consist of a question whose answer is either Yes or No (e.g. Is a symmetric matrix $A$ positive semidefinite? Does a given linear programme have feasible points?) Although less relevant in the context of numerical algorithms (as opposed to their relevance to discrete computations, where they form the backbone of complexity theory) they nonetheless form a class deserving analysis.

(a) Loss of precision.

Assume we are solving the system $Ax=b$ (or inverting the matrix $A$) with a real-number algorithm ${\mathcal{A}}$ executed using floating-point arithmetic, and let ${\it\epsilon}_{\mathsf{mach}}$ be the epsilon-machine (or round-off unit) at hand (for a primer on finite-precision computations see Higham Reference Higham1996, Chapters 2 and 3). This means that all floating-point numbers occurring during the execution of ${\mathcal{A}}$ are rounded to $|\!\log {\it\epsilon}_{\mathsf{mach}}|$ significant digits (here $\log$ refers to logarithm in base 10). In particular, this is the precision of the entries in the data $(A,b)$. On the other hand, the number of correct digits in the $i$th component of the computed solution $\tilde{x}$ is approximately $-\log \mathsf{RelError}(\tilde{x}_{i})$, where

$$\begin{eqnarray}\mathsf{RelError}(\tilde{x}_{i}):=\frac{\mathsf{Error}(\tilde{x}_{i})}{|x_{i}|}=\frac{|x_{i}-\tilde{x}_{i}|}{|x_{i}|}.\end{eqnarray}$$

Therefore, the value

$$\begin{eqnarray}\mathsf{LoP}(x_{i}):=\log \frac{\mathsf{RelError}(\tilde{x}_{i})}{{\it\epsilon}_{\mathsf{mach}}}=|\!\log {\it\epsilon}_{\mathsf{mach}}|-|\!\log \mathsf{RelError}(\tilde{x}_{i})|\end{eqnarray}$$

quantifies the loss of precision in the computation of $x_{i}$. If we take the precision of $\tilde{x}$ to be

$$\begin{eqnarray}-\log \mathsf{RelError}(\tilde{x}):=-\log \frac{\Vert x-\tilde{x}\Vert }{\Vert x\Vert },\end{eqnarray}$$

then the loss of precision in the computation of $x$ is

(3.1)$$\begin{eqnarray}\mathsf{LoP}(x):=\log \frac{\mathsf{RelError}(\tilde{x})}{{\it\epsilon}_{\mathsf{mach}}}=|\!\log {\it\epsilon}_{\mathsf{mach}}|+\log \mathsf{RelError}(\tilde{x}).\end{eqnarray}$$

Next assume that algorithm ${\mathcal{A}}$ corresponds to computing $x$ via the Householder QR factorization method. The following result is a well-known backward error analysis for this algorithm (see Higham Reference Higham1996, Theorem 19.3).

In all that follows, for any square matrix $M$, $\Vert M\Vert _{F}$ and $\Vert M\Vert$ denote the Frobenius and spectral norms of $M$, respectively.

Theorem 3.1. Let $A\in \mathbb{R}^{n\times n}$ be invertible and $b\in \mathbb{R}^{n}$. If the system $Ax=b$ is solved using the Householder QR factorization method, then the computed solution $\tilde{x}$ satisfies

$$\begin{eqnarray}\tilde{A}\tilde{x}=\tilde{b},\end{eqnarray}$$

where $\tilde{A}$ and $\tilde{b}$ satisfy the relative error bounds

$$\begin{eqnarray}\Vert \tilde{A}-A\Vert _{F}\leq n{\it\gamma}_{sn}\Vert A\Vert _{F}\quad \text{and}\quad \Vert \tilde{b}-b\Vert \leq n{\it\gamma}_{sn}\Vert b\Vert\end{eqnarray}$$

for a small constant $s$ and with ${\it\gamma}_{k}:=k{\it\epsilon}_{\mathsf{mach}}/(1-k{\it\epsilon}_{\mathsf{mach}}).$ ◻

We can tie this result to the general discussion on condition numbers in Section 1.1. To do so, we define the set of ill-posed matrices

$$\begin{eqnarray}{\rm\Sigma}:=\{A\in \mathbb{R}^{n\times n}\mid \det (A)=0\},\end{eqnarray}$$

write ${\mathcal{D}}:=\mathbb{R}^{n\times n}\setminus {\rm\Sigma}$ and let ${\it\varphi}:{\mathcal{D}}\times \mathbb{R}^{n}\rightarrow \mathbb{R}^{n}$ be given by ${\it\varphi}(A,b)=A^{-1}b$. Taking

(3.2)$$\begin{eqnarray}\mathsf{RelError}(\tilde{A},\tilde{b}):=\max \biggl\{\frac{\Vert \tilde{A}-A\Vert }{\Vert A\Vert },\frac{\Vert \tilde{b}-b\Vert }{\Vert b\Vert }\biggr\},\end{eqnarray}$$

Theorem 3.1 (and the bound $\Vert \cdot \Vert _{F}\leq \sqrt{n}\Vert \cdot \Vert$) shows that $\tilde{x}={\it\varphi}(\tilde{A},\tilde{b})$ with

$$\begin{eqnarray}\mathsf{RelError}(\tilde{A},\tilde{b})\leq n^{3/2}{\it\gamma}_{sn}.\end{eqnarray}$$

We now fix a normwise measure for the relative error of $\tilde{x}$ (i.e. $\mathsf{RelError}(\tilde{x}):=\Vert \tilde{x}-x\Vert /\Vert x\Vert$), consider the relative condition number $\mathsf{cond}^{{\it\varphi}}(A,b)$ with respect to the relative errors in both the input $(A,b)$ and the solution $x$, and set ${\it\epsilon}_{\mathsf{mach}}\rightarrow 0$. Then, ${\it\gamma}_{k}=k{\it\epsilon}_{\mathsf{mach}}+o({\it\epsilon}_{\mathsf{mach}})$ and (recall definition (1.3))

$$\begin{eqnarray}\mathsf{RelError}(\tilde{x})\leq \mathsf{RelError}(\tilde{A},\tilde{b})\,\mathsf{cond}^{{\it\varphi}}(A,b)+o({\it\epsilon}_{\mathsf{mach}}).\end{eqnarray}$$

We conclude that this relative error satisfies

(3.3)$$\begin{eqnarray}\frac{\Vert \tilde{x}-x\Vert }{\Vert x\Vert }\leq sn^{5/2}{\it\epsilon}_{\mathsf{mach}}\mathsf{cond}^{{\it\varphi}}(A,b)+o({\it\epsilon}_{\mathsf{mach}}),\end{eqnarray}$$

and consequently, that the loss of precision is bounded by

(3.4)$$\begin{eqnarray}\mathsf{LoP}(A^{-1}b)\leq \frac{5}{2}\log n+\log \mathsf{cond}(A,b)+\log s+o(1).\end{eqnarray}$$

(b) Characterizations of$\mathsf{cond}^{{\it\varphi}}(A,b)$.

We begin by proving an assertion we claimed early in Section 1.2. Define the map ${\it\kappa}:{\mathcal{D}}\rightarrow \mathbb{R}$ by

(3.5)$$\begin{eqnarray}{\it\kappa}(A):=\Vert A\Vert \Vert A^{-1}\Vert .\end{eqnarray}$$

Then ${\it\kappa}(A)\geq 1$, since $1=\Vert \text{Id}\Vert \leq \Vert A\Vert \Vert A^{-1}\Vert ={\it\kappa}(A)$.

Theorem 3.2. Let $\mathsf{cond}^{{\it\varphi}}(A,b)$ be as above. Then,

$$\begin{eqnarray}\mathsf{cond}^{{\it\varphi}}(A,b)={\it\kappa}(A)+\frac{\Vert A^{-1}\Vert \Vert b\Vert }{\Vert A^{-1}b\Vert }.\end{eqnarray}$$

In particular, we have

$$\begin{eqnarray}{\it\kappa}(A)\leq \mathsf{cond}^{{\it\varphi}}(A,b)\leq 2{\it\kappa}(A).\end{eqnarray}$$

Proof. Let $\tilde{A}=A-E$ and $\tilde{b}=b+f$. By definition, $\Vert E\Vert \leq \mathsf{R}\Vert A\Vert$ and $\Vert f\Vert \leq \mathsf{R}\Vert b\Vert$, where for simplicity, $\mathsf{R}=\mathsf{RelError}(A,b)$. We have, for $\mathsf{R}\rightarrow 0$,

$$\begin{eqnarray}\displaystyle (A-E)^{-1} & = & \displaystyle A^{-1}(\text{Id}-EA^{-1})^{-1}=A^{-1}(\text{Id}+EA^{-1}+o(\mathsf{R}))\nonumber\\ \displaystyle & = & \displaystyle A^{-1}+A^{-1}EA^{-1}+o(\mathsf{R}).\nonumber\end{eqnarray}$$

This implies, writing $x:=A^{-1}b$ and $\tilde{x}:=\tilde{A}^{-1}\tilde{b}$,

(3.6)$$\begin{eqnarray}\tilde{x}-x=(A-E)^{-1}(b+f)-x=A^{-1}Ex+A^{-1}f+o(\mathsf{R}).\end{eqnarray}$$

Taking norms we conclude that

$$\begin{eqnarray}\displaystyle \Vert \tilde{x}-x\Vert & \leq & \displaystyle \Vert A^{-1}E\Vert \Vert x\Vert +\Vert A^{-1}\Vert \Vert f\Vert +o(\mathsf{R})\nonumber\\ \displaystyle & \leq & \displaystyle \Vert A^{-1}\Vert \Vert E\Vert \Vert x\Vert +\Vert A^{-1}\Vert \Vert f\Vert +o(\mathsf{R})\nonumber\\ \displaystyle & \leq & \displaystyle \Vert A^{-1}\Vert \Vert A\Vert \Vert x\Vert \mathsf{R}+\Vert A^{-1}\Vert \Vert b\Vert \mathsf{R}+o(\mathsf{R}),\nonumber\end{eqnarray}$$

and hence

$$\begin{eqnarray}\frac{\Vert \tilde{x}-x\Vert }{\mathsf{R}\Vert x\Vert }\leq {\it\kappa}(A)+\frac{\Vert A^{-1}\Vert \Vert b\Vert }{\Vert x\Vert },\end{eqnarray}$$

which shows the upper bound in the claimed equality.

For the corresponding lower bound we choose $y\in \mathbb{R}^{n}$ such that $\Vert y\Vert =1$ and $\Vert A^{-1}y\Vert =\Vert A^{-1}\Vert$. Further, we choose $v=x/\Vert x\Vert \in \mathbb{R}^{n}$ so that $\Vert v\Vert =1$ and $v^{\text{T}}x=\Vert x\Vert$. Now we put

(3.7)$$\begin{eqnarray}E:=\mathsf{R}\Vert A\Vert \,yv^{\text{T}},\quad f:=\pm \mathsf{R}\Vert b\Vert \,y.\end{eqnarray}$$

We note that

$$\begin{eqnarray}\Vert E\Vert =\mathsf{R}\Vert A\Vert ,\quad \Vert f\Vert =\mathsf{R}\Vert b\Vert ,\end{eqnarray}$$

the first equality holding since $\Vert yv^{\text{T}}\Vert =\Vert y\Vert \Vert v\Vert =1$. We have

$$\begin{eqnarray}A^{-1}Ex=\mathsf{R}\Vert A\Vert \,A^{-1}y\,v^{\text{T}}x=\mathsf{R}\Vert A\Vert \Vert x\Vert A^{-1}y\end{eqnarray}$$

and hence $\Vert A^{-1}Ex\Vert ={\it\kappa}(A)\Vert x\Vert \mathsf{R}$. Similarly, $A^{-1}f=\pm R\Vert b\Vert A^{-1}y$ and $\Vert A^{-1}f\Vert =\Vert A^{-1}\Vert \Vert b\Vert \mathsf{R}$. Since $A^{-1}Ex$ and $A^{-1}f$ are both proportional to $A^{-1}y$, we obtain from (3.6)

$$\begin{eqnarray}\Vert \tilde{x}-x\Vert ={\it\kappa}(A)\Vert x\Vert \mathsf{R}+\Vert A^{-1}\Vert \Vert b\Vert \mathsf{R},\end{eqnarray}$$

if we choose the sign for $f$ in (3.7) appropriately. This proves the claimed lower bound.

The next result shows that ${\it\kappa}$ actually coincides with the condition number for the problem of matrix inversion.

Theorem 3.3. Let ${\it\psi}:{\mathcal{D}}\rightarrow \mathbb{R}^{n\times n}$ be given by ${\it\psi}(A)=A^{-1}$. We measure the relative error on the data space and solution space with respect to spectral and Euclidean norms, respectively. Then we have

$$\begin{eqnarray}\mathsf{cond}^{{\it\psi}}(A)={\it\kappa}(A).\end{eqnarray}$$

Proof. Let $E\in \mathbb{R}^{n\times n}$ be such that $\tilde{A}=A-E$. Then $\mathsf{RelError}(A)=\Vert E\Vert /\Vert A\Vert$. As in the proof of Theorem 3.2, we have for $\Vert E\Vert \rightarrow 0$,

(3.8)$$\begin{eqnarray}\Vert \tilde{A}^{-1}-A^{-1}\Vert =\Vert A^{-1}EA^{-1}\Vert +o(\Vert E\Vert ).\end{eqnarray}$$

Hence, $\Vert A^{-1}EA^{-1}\Vert \leq \Vert A^{-1}\Vert \Vert E\Vert A^{-1}\Vert .$ Consequently, we obtain

$$\begin{eqnarray}\mathsf{RelError}(A^{-1})=\frac{\Vert \tilde{A}^{-1}-A^{-1}\Vert }{\Vert A^{-1}\Vert }\leq \Vert A^{-1}\Vert \Vert E\Vert +o(\Vert E\Vert )).\end{eqnarray}$$

We conclude that

$$\begin{eqnarray}\frac{\mathsf{RelError}(A^{-1})}{\mathsf{RelError}(A)}\leq \Vert A\Vert \Vert A^{-1}\Vert +o(1)\end{eqnarray}$$

and hence $\mathsf{cond}^{{\it\psi}}(A)\leq {\it\kappa}(A)$.

To prove the reverse inequality it is enough to find arbitrarily small matrices $E$ such that $\Vert A^{-1}EA^{-1}\Vert =\Vert A^{-1}\Vert ^{2}\Vert E\Vert$, since then we can proceed from (3.8) as we did in Theorem 3.2 from (3.6).

To do so, let $y\in \mathbb{R}^{n}$ be such that $\Vert y\Vert =1$ and $\Vert A^{-1}y\Vert =\Vert A^{-1}\Vert$. Define $x:=\Vert A^{-1}\Vert ^{-1}A^{-1}y$, so that $A^{-1}y=\Vert A^{-1}\Vert x$ and $\Vert x\Vert =\Vert y\Vert =1$. For any $B\in \mathbb{R}^{n\times n}$ we have

$$\begin{eqnarray}\Vert A^{-1}BA^{-1}\Vert \geq \Vert A^{-1}BA^{-1}y\Vert =\Vert A^{-1}\Vert \cdot \Vert A^{-1}Bx\Vert .\end{eqnarray}$$

Since $\Vert x\Vert =\Vert y\Vert =1$ there exists an orthogonal matrix $B\in \mathbb{R}^{n\times n}$ such that $Bx=y$. Therefore,

$$\begin{eqnarray}\Vert A^{-1}BA^{-1}\Vert \geq \Vert A^{-1}\Vert \cdot \Vert A^{-1}y\Vert =\Vert A^{-1}\Vert ^{2}.\end{eqnarray}$$

Taking $E={\it\delta}B$ with arbitrarily small ${\it\delta}$ finishes the proof.

(c) The condition number theorem.

We mentioned in Section 2 that Demmel (Reference Demmel1987) showed that many problems admit conic condition numbers. The first result of this kind, which preceded the expression ‘condition number’, was proved by Eckart and Young (Reference Eckart and Young1936). It states that for $A\in \mathbb{R}^{n\times n}\setminus {\rm\Sigma}$, we have

(3.9)$$\begin{eqnarray}\Vert A^{-1}\Vert =\frac{1}{\mathsf{dist}(A,{\rm\Sigma})}\end{eqnarray}$$

where $\mathsf{dist}$ is the distance induced by either the spectral or the Frobenius norm. It follows that

$$\begin{eqnarray}{\it\kappa}(A)=\frac{\Vert A\Vert }{\mathsf{dist}(A,{\rm\Sigma})}.\end{eqnarray}$$

(d) Probabilistic analyses.

Now we prove an average analysis for $\log {\it\kappa}(A)$. The bounds we obtain are not the sharpest possible (see Theorem 3.12) but we can use them to derive sufficiently good bounds for the average loss of precision.

For an $n\times n$ matrix $A$, let $\hat{A}_{i}$ refer to $A$ with the $i$th column removed and, for an $n\times (n-1)$ matrix $A$ and $b\in \mathbb{R}^{n}$, let $(A,b)_{i}$ refer to the $n\times n$ matrix obtained from $A$ by inserting $b$ as the $i$th column. Also, let us denote by $A_{i}$ the $i$th column of $A$.

The next lemma relates the norm of the inverse of an $n\times n$ matrix to a geometric condition on the columns of the matrix.

Proof. That $\Vert A^{-1}\Vert =x$ means that there is a unit vector $b\in \mathbb{R}^{n}$ such that $v=A^{-1}b$ and $\Vert v\Vert =x$. Choose $i$ such that $|v_{i}|\geq x/\sqrt{n}$. By Cramer’s rule

$$\begin{eqnarray}v_{i}=\frac{\det (\hat{A}_{i},b)_{i}}{\det A}.\end{eqnarray}$$

Hence $\det (\hat{A}_{i},{v_{i}}^{-1}b)_{i}=\det A$ and $\det (\hat{A}_{i},A_{i}-{v_{i}}^{-1}b)_{i}=0$. As the columns of $A$ are independent, this implies that $A_{i}-{v_{i}}^{-1}b$ is spanned by the other columns. Also, $\Vert {v_{i}}^{-1}b\Vert \leq \sqrt{n}/x$ so the statement holds for $B_{i}:=-{v_{i}}^{-1}b$.

We shall let ${\mathcal{M}}_{(n,m)}$ denote the space of real $n\times m$ matrices, and if $m=n$ we shall denote this space simply by ${\mathcal{M}}_{n}$. Let $\mathbb{S}^{n^{2}-1}\subset {\mathcal{M}}_{n}$ be the unit sphere with respect to the Frobenius norm. Also, let $\mathbb{B}^{n^{2}-n}$ be the unit ball in the space ${\mathcal{M}}_{(n,n-1)}$ and for $i=1,\ldots ,n$ the map

For any ${\it\varepsilon}>0$ we consider $N_{{\it\varepsilon}}\subset \mathbb{S}^{n^{2}-1}$ the set of $A\in \mathbb{S}^{n^{2}-1}$ such that $A$ is not invertible or $\Vert A^{-1}\Vert >1/{\it\varepsilon}$. Equation (3.9) provides a geometric characterization of $N_{{\it\varepsilon}}$, namely

$$\begin{eqnarray}N_{{\it\varepsilon}}=\{A\in \mathbb{S}^{n^{2}-1}\mid \mathsf{dist}_{F}(A,{\rm\Sigma})<{\it\varepsilon}\}\end{eqnarray}$$

where $\mathsf{dist}_{F}$ is the distance in ${\mathcal{M}}_{n}$ induced by the Frobenius norm.

Our goal is to estimate the volume of $N_{{\it\varepsilon}}$. To do so, consider for ${\it\varepsilon}>0$ and $i=1,\ldots ,n$ the set

$$\begin{eqnarray}\displaystyle W_{{\it\varepsilon},i} & = & \displaystyle \big\{A\in \mathbb{S}^{n^{2}-1}\mid A_{i}\text{ is within }\sqrt{n}\,{\it\varepsilon}\text{ of the}\nonumber\\ \displaystyle & & \displaystyle \text{space spanned by the columns of }\hat{A}_{i}\big\}.\nonumber\end{eqnarray}$$

Lemma 3.5. For all ${\it\varepsilon}>0$,

$$\begin{eqnarray}\mathsf{vol}(N_{{\it\varepsilon}})\leq \mathop{\sum }_{i=1}^{n}\mathsf{vol}(W_{{\it\varepsilon},i}).\end{eqnarray}$$

Now we turn to $W_{{\it\varepsilon},i}$.

Proposition 3.6. For all ${\it\varepsilon}>0$, $B\in {\mathcal{M}}_{(n,n-1)}$, and $i=1,\ldots ,n$, we have

$$\begin{eqnarray}\mathsf{vol}(W_{{\it\varepsilon},i}\cap {\it\pi}_{i}^{-1}B)\leq 2\sqrt{n}\,{\it\varepsilon}\mathsf{vol}(\mathbb{S}_{(1-\Vert B\Vert _{F}^{2})^{1/2}}^{n-2}),\end{eqnarray}$$

where $\mathbb{S}_{r}^{m-1}$ is the sphere of radius $r$ in $\mathbb{R}^{m}$.

Proof. If $A\in W_{{\it\varepsilon},i}\cap {\it\pi}_{i}^{-1}B$ then $\Vert A_{i}\Vert =(1-\Vert B\Vert _{F}^{2})^{1/2}$, that is,

$$\begin{eqnarray}A_{i}\in \mathbb{S}_{(1-\Vert B\Vert _{F}^{2})^{1/2}}^{n-1},\end{eqnarray}$$

and $A_{i}$ is within $\sqrt{n}{\it\varepsilon}$ of the space spanned by the columns of $B$.

The volume of $W_{{\it\varepsilon},i}\cap {\it\pi}_{i}^{-1}B$ is thus at most the volume of the cylinder of base $\mathbb{S}_{(1-\Vert B\Vert _{F}^{2})^{1/2}}^{n-2}$ and height $2\sqrt{n}\,{\it\varepsilon}$.

Lemma 3.7. Let $\mathbb{B}^{m-\ell }\subset \mathbb{R}^{m-\ell }$ be the unit ball and $k\geq 0$. Then

$$\begin{eqnarray}\frac{1}{\mathsf{vol}(\mathbb{S}^{m-1})}\int _{x\in \mathbb{B}^{m-\ell }}\frac{\mathsf{vol}(\mathbb{S}_{(1-\Vert x\Vert ^{2})^{1/2}}^{k})}{(1-\Vert x\Vert ^{2})^{1/2}}\,\text{d}x_{1}\ldots \,\text{d}x_{m-\ell }=\frac{{\it\pi}^{\frac{k-\ell +1}{2}}{\it\Gamma}(\frac{m}{2})}{{\it\Gamma}(\frac{m+k-\ell +1}{2})}.\end{eqnarray}$$

Proof. The left-hand side in the statement equals

$$\begin{eqnarray}\frac{\mathsf{vol}(\mathbb{S}^{k})}{\mathsf{vol}(\mathbb{S}^{m-1})}\int _{x\in \mathbb{B}^{m-\ell }}(1-\Vert x\Vert ^{2})^{\frac{k-1}{2}}\,\text{d}x_{1}\ldots \,\text{d}x_{m-\ell },\end{eqnarray}$$

which using spherical coordinates equals

$$\begin{eqnarray}\frac{\mathsf{vol}(\mathbb{S}^{k})\mathsf{vol}(\mathbb{S}^{m-\ell -1})}{\mathsf{vol}(\mathbb{S}^{m-1})}\int _{0}^{1}(1-r^{2})^{\frac{k-1}{2}}\cdot r^{m-\ell -1}\,\text{d}r.\end{eqnarray}$$

Finally, using that for $p,q\in \mathbb{R}$,

$$\begin{eqnarray}\int _{0}^{1}t^{p}(1-t^{2})^{q}\,\text{d}t=\frac{1}{2}\frac{{\it\Gamma}(\frac{p+1}{2}){\it\Gamma}(q+1)}{{\it\Gamma}(\frac{p+1}{2}+q+1)},\end{eqnarray}$$

the fact that

$$\begin{eqnarray}\mathsf{vol}(\mathbb{S}^{p-1})=\frac{2{\it\pi}^{p/2}}{{\it\Gamma}(\frac{p}{2})},\end{eqnarray}$$

and some elementary computations, we obtain the right-hand side.

We can now estimate the volume of $N_{{\it\varepsilon}}$.

Proposition 3.8. For any ${\it\varepsilon}>0$ and any $n\geq 2$,

$$\begin{eqnarray}\frac{\mathsf{vol}(N_{{\it\varepsilon}})}{\mathsf{vol}(\mathbb{S}^{n^{2}-1})}\leq {\it\varepsilon}n^{5/2}.\end{eqnarray}$$

Proof. By Lemma 3.5 it suffices to show that

$$\begin{eqnarray}\frac{\mathsf{vol}(W_{{\it\varepsilon},i})}{\mathsf{vol}(\mathbb{S}^{n^{2}-1})}\leq {\it\varepsilon}n^{3/2}.\end{eqnarray}$$

Let ${\it\chi}_{W_{{\it\varepsilon},i}}$ be the characteristic function of $W_{{\it\varepsilon},i}$ and let ${\it\pi}:\mathbb{S}^{n^{2}-1}\rightarrow \mathbb{B}^{n^{2}-n}$ be the restriction of the orthogonal projection. Then

$$\begin{eqnarray}\displaystyle & & \displaystyle \frac{\mathsf{vol}(W_{{\it\varepsilon},i})}{\mathsf{vol}(\mathbb{S}^{n^{2}-1})}\nonumber\\ \displaystyle & & \displaystyle \qquad =\frac{1}{\mathsf{vol}(\mathbb{S}^{n^{2}-1})}\int _{A\in \mathbb{S}^{n^{2}-1}}{\it\chi}_{W_{{\it\varepsilon},i}}(A)\,\text{d}A\nonumber\\ \displaystyle & & \displaystyle \qquad =\frac{1}{\mathsf{vol}(\mathbb{S}^{n^{2}-1})}\int _{B\in \mathbb{B}^{n^{2}-n}}\int _{Q\in {\it\pi}^{-1}(B)}{\it\chi}_{W_{{\it\varepsilon},i}}(Q)(1-\Vert B\Vert _{F}^{2})^{-1/2}\,\text{d}Q\,\text{d}B\nonumber\\ \displaystyle & & \displaystyle \qquad =\frac{1}{\mathsf{vol}(\mathbb{S}^{n^{2}-1})}\int _{B\in \mathbb{B}^{n^{2}-n}}\mathsf{vol}(W_{{\it\varepsilon},i}\cap {\it\pi}_{i}^{-1}(B))(1-\Vert B\Vert _{F}^{2})^{-1/2}\,\text{d}B\nonumber\\ \displaystyle & & \displaystyle \qquad \leq \frac{2\sqrt{n}{\it\varepsilon}}{\mathsf{vol}(\mathbb{S}^{n^{2}-1})}\int _{B\in \mathbb{B}^{n^{2}-n}}\mathsf{vol}(\mathbb{S}_{(1-\Vert B\Vert _{F}^{2})^{1/2}}^{n-2})(1-\Vert B\Vert _{F}^{2})^{-1/2}\,\text{d}B,\nonumber\end{eqnarray}$$

the second equality by the transformation formula (using that the normal Jacobian of ${\it\pi}$ at $A$ is $(1-\Vert {\it\pi}(A)\Vert _{F}^{2})^{-1/2}$), and the inequality by Proposition 3.6. Using Lemma 3.7 (with $m=n^{2}$, $k=n-2$ and $\ell =n$) we deduce that the last line is at most

$$\begin{eqnarray}\frac{2\sqrt{n}{\it\varepsilon}}{{\it\pi}^{1/2}}\frac{{\it\Gamma}(\frac{n^{2}}{2})}{{\it\Gamma}(\frac{n^{2}-1}{2})}\leq \frac{2^{1/2}}{{\it\pi}^{1/2}}{\it\varepsilon}n^{3/2}\leq {\it\varepsilon}n^{3/2},\end{eqnarray}$$

the first inequality since

$$\begin{eqnarray}\frac{{\it\Gamma}(x+1/2)}{\sqrt{x}{\it\Gamma}(x)}<1,\quad \text{for all}~x>0.\square\end{eqnarray}$$

Theorem 3.9. The probability that ${\it\kappa}(A)>1/{\it\varepsilon}$ for the standard Gaussian probability distribution $N(0,\text{Id})$ on ${\mathcal{M}}_{n}$ is less than or equal to $n^{5/2}{\it\varepsilon}$. Also,

$$\begin{eqnarray}\underset{A\sim N(0,\text{Id})}{\mathbb{E}}\log {\it\kappa}(A)\leq \frac{5}{2}\log n+O(1).\end{eqnarray}$$

Proof. For the first statement in the theorem, define a second condition number ${\it\kappa}_{F}$ of $n\times n$ matrices by

$$\begin{eqnarray}{\it\kappa}_{F}(A)=\Vert A^{-1}\Vert \,\Vert A\Vert _{F}.\end{eqnarray}$$

Clearly, ${\it\kappa}(A)\leq {\it\kappa}_{F}(A)$. Hence, it suffices to prove the statement for ${\it\kappa}_{F}$. Also, since the latter is homogeneous of degree 0, it suffices to prove the theorem for the uniform distribution on the sphere. And finally, since $N_{{\it\varepsilon}}$ is the set of $A\in \mathbb{S}^{n^{2}-1}$ where ${\it\mu}(A)\geq 1/{\it\varepsilon}$, Proposition 3.8 gives the desired bound.

The bound for the expectation follows from integrating the bound for the probability tail.

The same ideas allow one to derive bounds for the tail and expectation of the condition number of complex matrices. The set $N_{{\it\varepsilon}}^{\mathbb{C}}$ is defined as in the real case, and one can prove that

$$\begin{eqnarray}\frac{\mathsf{vol}(N_{{\it\varepsilon}}^{\mathbb{C}})}{\mathsf{vol}(\mathbb{S}^{2n^{2}-1})}\leq n^{4}{\it\varepsilon}^{2}.\end{eqnarray}$$

The following result for the set ${\mathcal{M}}_{n}^{\mathbb{C}}$ of complex $n\times n$ matrices is an easy corollary of this bound.

Theorem 3.10. The probability that ${\it\mu}(A)>1/{\it\varepsilon}$ or that ${\it\kappa}(A)>1/{\it\varepsilon}$ for either the standard Gaussian probability distribution on ${\mathcal{M}}_{n}^{\mathbb{C}}$ or the uniform distribution on $\mathbb{S}^{2n^{2}-1}$ (the unit sphere in ${\mathcal{M}}_{n}^{\mathbb{C}}$) is less than or equal to $n^{4}{\it\varepsilon}^{2}$. Also,

$$\begin{eqnarray}\underset{A\sim N_{\mathbb{C}}(0,\text{Id})}{\mathbb{E}}\log {\it\kappa}(A)\leq 2\log n+O(1).\end{eqnarray}$$

(e) Average loss of precision in linear equation solving.

The bound in (3.4) and Theorem 3.9 immediately yield the following result on the expected loss of precision in linear equation solving.

Theorem 3.11. For all $b\in \mathbb{R}^{n}$,

$$\begin{eqnarray}\mathbb{E}_{A\sim N(0,\text{Id})}\mathsf{LoP}(A,b)\leq 5\log n+O(1).\square\end{eqnarray}$$

(f) Additional results.

The condition number ${\it\kappa}$ is doubtless better studied than any other condition number. Because of this, the content of the previous pages is just a small fraction of what is known about them. Without intending to fill this gap entirely, we next state two results that either sharpen or add to this content.

The first, by Alan Edelman (Reference Edelman1988), shows a sharp bound for $\mathbb{E}\log {\it\kappa}(A)$.

Theorem 3.12. For both real and complex standard Gaussian $n\times n$ matrices $A$ we have

$$\begin{eqnarray}\mathbb{E}\log {\it\kappa}(A)=\log n+C+O(1),\end{eqnarray}$$

where $C=1.537$ in the real case and $C=0.982$ in the complex case. ◻

The second one is the best smoothed analysis, as of today, for ${\it\kappa}(A)$. It is due to Mario Wschebor (Reference Wschebor2004).

Theorem 3.13. There exists a universal constant $K\geq 1$ such that for all $\overline{A}\in \mathbb{R}^{n\times n}$ with $\Vert \overline{A}\Vert \leq 1$, all $0<{\it\sigma}\leq 1$, and all $t\geq 1$, we have

$$\begin{eqnarray}\underset{A\sim N(\overline{A},{\it\sigma}^{2}\text{Id})}{\mathsf{Prob}}\{{\it\kappa}(A)\geq t\}\leq Kn\,\frac{1}{{\it\sigma}t}.\end{eqnarray}$$

Consequently, for all ${\it\sigma}\in (0,1]$ and all $\overline{A}\in \mathbb{R}^{n\times n}$ with $\Vert \overline{A}\Vert \leq 1$, we have

$$\begin{eqnarray}\mathbb{E}_{A\sim N(\overline{A},{\it\sigma}^{2}\text{Id})}(\log {\it\kappa}(A))\leq \log n+\log \frac{1}{{\it\sigma}}+2+\log K.\square\end{eqnarray}$$

(a) The algorithm.

For $x\in \mathbb{R}^{n}$ we define its $A$-norm by

(3.13)$$\begin{eqnarray}\Vert x\Vert _{A}:=(x^{\text{T}}Ax)^{1/2}\end{eqnarray}$$

coming from the scalar product

$$\begin{eqnarray}\left\langle x,y\right\rangle _{A}:=x^{\text{T}}Ay.\end{eqnarray}$$

We say that $x,y$ are $A$-orthogonal (or conjugated) if $\left\langle x,y\right\rangle _{A}=0$. Let $(d_{0},\ldots ,d_{n-1})$ be an $A$-orthogonal basis of $\mathbb{R}^{n}$, that is, $\left\langle d_{i},d_{j}\right\rangle _{A}=0$ for $i\neq j$. Moreover, let $\bar{x}$ be the solution of $Ax=b$ and $x_{0}\in \mathbb{R}^{n}$ be an initial point. Expand

$$\begin{eqnarray}\bar{x}-x_{0}=\mathop{\sum }_{i=0}^{n-1}{\it\alpha}_{i}d_{i},\end{eqnarray}$$

with ${\it\alpha}_{i}\in \mathbb{R}$. Then, by $A$-orthogonality,

$$\begin{eqnarray}{\it\alpha}_{k}=\frac{\left\langle \bar{x}-x_{0},d_{k}\right\rangle _{A}}{\Vert d_{k}\Vert _{A}^{2}}.\end{eqnarray}$$

Define

$$\begin{eqnarray}x_{k}:=\mathop{\sum }_{i=0}^{k-1}{\it\alpha}_{i}d_{i}+x_{0}.\end{eqnarray}$$

Then we have $\left\langle x_{k}-x_{0},d_{\ell }\right\rangle _{A}=0$, for all $\ell =k,\ldots ,n-1$. In particular, taking $\ell =k$,

$$\begin{eqnarray}\left\langle \bar{x}-x_{0},d_{k}\right\rangle _{A}=\left\langle \bar{x}-x_{k},d_{k}\right\rangle _{A}=d_{k}^{\text{T}}A(\bar{x}-x_{k})=d_{k}^{\text{T}}r_{k},\end{eqnarray}$$

with the residue $r_{k}:=b-Ax_{k}=A\bar{x}-Ax_{k}$. We obtain that

$$\begin{eqnarray}{\it\alpha}_{k}=\frac{\left\langle d_{k},r_{k}\right\rangle }{\Vert d_{k}\Vert _{A}^{2}}.\end{eqnarray}$$

Note that ${\it\alpha}_{k}$ depends only on $d_{k}$ and $r_{k}$.

Proof. For part (a) note that the point $x_{k}-x_{0}$ is the $A$-orthogonal projection of $\bar{x}-x_{0}$ onto $\mathsf{span}\{d_{0},\ldots ,d_{k-1}\}$. Therefore, $x_{k}-x_{0}$ minimizes the $A$-distance to $\bar{x}$,

$$\begin{eqnarray}x\;\mapsto \;\Vert x-\bar{x}\Vert _{A}=\Vert (x-x_{0})-(\bar{x}-x_{0})\Vert _{A},\end{eqnarray}$$

on $x_{0}+\mathsf{span}\{d_{0},\ldots ,d_{k-1}\}$.

For part (b), consider the function

$$\begin{eqnarray}f(x):=\frac{1}{2}x^{\text{T}}Ax-b^{\text{T}}x.\end{eqnarray}$$

The gradient of $f$ at $x$ is $\mathsf{grad}\,f(x)=Ax-b$, and it is easy to see (Bürgisser and Cucker Reference Bürgisser and Cucker2013, § 5.2) that

$$\begin{eqnarray}f(x)=\frac{1}{2}\Vert x-\bar{x}\Vert _{A}^{2}+f(\bar{x}).\end{eqnarray}$$

Part (a) then implies that $x_{k}$ minimizes $f$ on $x_{0}+\mathsf{span}\{d_{0},\ldots ,d_{k-1}\}$. Hence, $\mathsf{grad}\,f(x_{k})$ is orthogonal to $\mathsf{span}\{d_{0},\ldots ,d_{k-1}\}$. But $\mathsf{grad}\,f(x_{k})=-r_{k}$, which proves the first statement in (b). The second statement follows from the fact that $r_{n}$ is then orthogonal to $d_{0},\ldots ,d_{n-1}$ (and hence it must be 0). This, in turn, implies $x_{n}=\bar{x}$.

Thus far, we have assumed an $A$-orthogonal basis $\{d_{0},\ldots ,d_{n-1}\}$ at hand. The conjugate gradient method actually computes one, along with the computation of the $x_{k}$, by sequentially choosing the directions $d_{k}$ as conjugates of the gradients $-r_{k}$. Algorithm 1 shows how this is done (a proof of the $A$-orthogonality of the resulting set can be found in Bürgisser and Cucker Reference Bürgisser and Cucker2013, Theorem 5.12).

(b) First complexity considerations.

We begin by looking at the cost (i.e. number of arithmetic operations performed) of each iteration of the algorithm. The cost of computing an inner product, such as $\langle r_{k},d_{k}\rangle$, is $2n-1$. Consequently, the cost of a matrix–vector multiplication, such as $Ax_{k}$, is $2n^{2}-n$. It follows that computing an $A$-inner product costs $2n^{2}+n-1$. At each iteration of Conj_Grad the computation of each of $x_{k+1}$ and $d_{k+1}$ thus takes $O(n)$ arithmetic operations, and those of ${\it\alpha}_{k}$, $r_{k+1}$, and ${\it\beta}_{k}$ take $2n^{2}+O(n)$ each (note that for the latter we use the already computed $\Vert d_{k}\Vert _{A}^{2}$). That is, the cost of an iteration of Conj_Grad is $6n^{2}+O(n)$.

It follows from this bound that the cost of computing the solution $\bar{x}$ of $Ax=b$ is $6n^{3}+O(n^{2})$. This bound is larger than the $\frac{2}{3}n^{3}$ bound for the cost of Gaussian elimination, making the conjugate gradient a poor choice for the computation of $\bar{x}$.

We note at this juncture, however, that unlike Gaussian elimination, the conjugate gradient produces a sequence $x_{k}$ of approximations of $\bar{x}$. If we are only interested in a good approximation, as opposed to the exact value of $\bar{x}$, then early termination of the conjugate gradient might outperform Gaussian elimination. To understand this issue we need to estimate the speed of convergence of the former.

(c) Convergence analysis.

Let ${\it\epsilon}_{k}:=x_{k}-\bar{x}$ be the $k$th error.

Theorem 3.21. The error at the $k$th step of the conjugate gradient satisfies

$$\begin{eqnarray}\Vert {\it\epsilon}_{k}\Vert _{A}\leq 2\biggl(\frac{\sqrt{{\it\kappa}(A)}-1}{\sqrt{{\it\kappa}(A)}+1}\biggr)^{k}\Vert {\it\epsilon}_{0}\Vert _{A}.\end{eqnarray}$$

Corollary 3.22. For all ${\it\delta}>0$ we have $\Vert {\it\epsilon}_{k}\Vert _{A}\leq {\it\delta}\Vert {\it\epsilon}_{0}\Vert _{A}$ whenever

$$\begin{eqnarray}k\geq \ln \biggl(\frac{\sqrt{{\it\kappa}(A)}+1}{\sqrt{{\it\kappa}(A)}-1}\biggr)^{-1}\ln \biggl(\frac{2}{{\it\delta}}\biggr)\approx \frac{1}{2}\sqrt{{\it\kappa}(A)}\ln \biggl(\frac{1}{{\it\delta}}\biggr).\end{eqnarray}$$

With the proof of Theorem 3.21 in mind, we introduce some notation. We ${\mathcal{P}}_{k}$ denote the linear space of all real polynomials in one variable $X$ with degree at most $k$. We also write ${\mathcal{Q}}_{k}$ for the subset of ${\mathcal{P}}_{k}$ of polynomials with constant coefficient 1.

Theorem 3.23. The error ${\it\epsilon}_{k}$ at the $k$th step of the conjugate gradient method satisfies

$$\begin{eqnarray}\Vert {\it\epsilon}_{k}\Vert _{A}^{2}\leq \min _{q\in {\mathcal{Q}}_{k}}\max _{j\leq n}q({\it\lambda}_{j})^{2}\cdot \Vert {\it\epsilon}_{0}\Vert _{A}^{2},\end{eqnarray}$$

where ${\it\lambda}_{1},\ldots ,{\it\lambda}_{n}$ are the eigenvalues of $A$.

Proof. By Proposition 3.20(a) we know that $x_{k}$ minimizes the $A$-distance of $\bar{x}$ to the affine space

$$\begin{eqnarray}x_{0}+S_{k-1}=x_{0}+\mathsf{span}\{d_{0},\ldots ,d_{k-1}\}=x_{0}+\mathsf{span}\{r_{0},Ar_{0},\ldots ,A^{k-1}r_{0}\}.\end{eqnarray}$$

An element $x$ of $x_{0}+S_{k-1}$ can therefore be written as

$$\begin{eqnarray}x=x_{0}-p(A)r_{0},\end{eqnarray}$$

with $p\in {\mathcal{P}}_{k-1}$, and conversely, for any such polynomial we obtain $x\in x_{0}+Sk-1$. Using $r_{0}=b-Ax_{0}=A(\bar{x}-x_{0})$ we get

$$\begin{eqnarray}\displaystyle x-\bar{x} & = & \displaystyle x_{0}-\bar{x}-p(A)A(\bar{x}-x_{0})=x_{0}-\bar{x}+Ap(A)(x_{0}-\bar{x})\nonumber\\ \displaystyle & = & \displaystyle (\text{Id}+Ap(A))(x_{0}-\bar{x}).\nonumber\end{eqnarray}$$

It follows that the error ${\it\epsilon}_{k}=x_{k}-\bar{x}$ at the $k$th step of the conjugate gradient method satisfies

(3.14)$$\begin{eqnarray}\Vert {\it\epsilon}_{k}\Vert _{A}=\min _{x\in x_{0}+S_{k-1}}\Vert x-\bar{x}\Vert _{A}=\min _{p\in {\mathcal{P}}_{k-1}}\Vert (\text{Id}+Ap(A)){\it\epsilon}_{0}\Vert _{A}.\end{eqnarray}$$

Suppose that $v_{1},\ldots ,v_{n}$ is an orthonormal basis of eigenvectors of $A$ corresponding to the eigenvalues ${\it\lambda}_{1}\geq \cdots \geq {\it\lambda}_{n}>0$. Write

$$\begin{eqnarray}{\it\epsilon}_{0}=x_{0}-\bar{x}=\mathop{\sum }_{j=1}^{n}{\it\xi}_{j}v_{j}\end{eqnarray}$$

for some ${\it\xi}_{1},\ldots ,{\it\xi}_{n}\in \mathbb{R}$. Then

$$\begin{eqnarray}\Vert {\it\epsilon}_{0}\Vert _{A}^{2}=\mathop{\sum }_{j=1}^{n}\mathop{\sum }_{k=1}^{n}{\it\xi}_{j}{\it\xi}_{k}\;v_{j}^{\text{T}}Av_{k}=\mathop{\sum }_{j=1}^{n}{\it\lambda}_{j}{\it\xi}_{j}^{2}.\end{eqnarray}$$

Moreover, for any polynomial $p\in {\mathcal{P}}_{k-1}$,

$$\begin{eqnarray}(\text{Id}+Ap(A)){\it\epsilon}_{0}=\mathop{\sum }_{j=1}^{n}(1+{\it\lambda}_{j}p({\it\lambda}_{j})){\it\xi}_{j}v_{j}.\end{eqnarray}$$

Therefore

$$\begin{eqnarray}\Vert (\text{Id}+Ap(A)){\it\epsilon}_{0}\Vert _{A}^{2}=\mathop{\sum }_{j=1}^{n}(1+{\it\lambda}_{j}p({\it\lambda}_{j}))^{2}{\it\lambda}_{j}{\it\xi}_{j}^{2}\leq \max _{j\leq n}(1+{\it\lambda}_{j}p({\it\lambda}_{j}))^{2}\mathop{\sum }_{j=1}^{n}{\it\lambda}_{j}{\it\xi}_{j}^{2},\end{eqnarray}$$

and using (3.14),

$$\begin{eqnarray}\Vert {\it\epsilon}_{k}\Vert _{A}^{2}\leq \min _{p\in {\mathcal{P}}_{k-1}}\max _{j\leq n}(1+{\it\lambda}_{j}p({\it\lambda}_{j}))^{2}\mathop{\sum }_{j=1}^{n}{\it\lambda}_{j}{\it\xi}_{j}^{2}.\end{eqnarray}$$

The result now follows by observing that ${\mathcal{Q}}_{k}=1+X{\mathcal{P}}_{k-1}$.

Theorem 3.23 is hard to apply in concrete situations. It depends on all the eigenvalues of $A$ and it optimizes a function of them over the space ${\mathcal{Q}}_{k}$. It is nevertheless a building block in the proof of Theorem 3.21. We proceed to see why.

Suppose that ${\it\lambda}_{n}=a$ and ${\it\lambda}_{1}=b$, for $0<a<b$. What are the optimal polynomials $q\in {\mathcal{Q}}_{k}$ if nothing is known about the location of the eigenvalues ${\it\lambda}$ except that ${\it\lambda}\in [a,b]$? In this case we have to minimize the quantity

$$\begin{eqnarray}\max _{a\leq {\it\lambda}\leq b}q({\it\lambda})^{2}\end{eqnarray}$$

over all real polynomials $q\in {\mathcal{Q}}_{k}$. This minimization problem can also be considered for arbitrary $a,b\in \mathbb{R}$ with $a<b$. In the particular case $a=-1$, $b=1$ it turns out that its solution is given by the Chebyshev polynomials $T_{k}$ (defined by the equation $\cos (k{\it\phi})=T_{k}(\cos ({\it\phi}))$).

Proposition 3.24. For all $0<a<b$ we have

$$\begin{eqnarray}\min _{\substack{ \mathsf{deg}q\leq k \\ q(0)=1}}\;\max _{a\leq {\it\lambda}\leq b}q({\it\lambda})^{2}=\frac{1}{T_{k}(\frac{b+a}{b-a})^{2}}.\end{eqnarray}$$

Proof. The polynomial

$$\begin{eqnarray}p({\it\lambda}):=\frac{T_{k}(\frac{b+a-2{\it\lambda}}{b-a})}{T_{k}(\frac{b+a}{b-a})}\end{eqnarray}$$

has degree $k$, satisfies $p(0)=1$, and

$$\begin{eqnarray}\max _{a\leq {\it\lambda}\leq b}p({\it\lambda})^{2}=\frac{1}{T_{k}(\frac{b+a}{b-a})^{2}}=:c^{2}.\end{eqnarray}$$

(Note that ${\it\lambda}\mapsto (b+a-2{\it\lambda})/(b-a)$ maps $[a,b]$ onto $[-1,1]$ and $T_{k}$ has maximal value 1 on this interval.)

In order to show the optimality of $p$, suppose that for some polynomial $q\in {\mathcal{Q}}_{k}$,

$$\begin{eqnarray}\max _{a\leq {\it\lambda}\leq b}|q({\it\lambda})|<\max _{a\leq {\it\lambda}\leq b}|p({\it\lambda})|.\end{eqnarray}$$

Since $p$ alternately takes the values $c,-c,c,-c,\ldots$ exactly $k+1$ times in the interval $[a,b]$, we conclude that the polynomial $p-q$ has at least $k$ roots in this interval (intermediate value theorem). Moreover, $p(0)-q(0)=1-1=0$. Hence $p-q$ has degree at most $k$ and $k+1$ distinct roots. It follows that $p=q$.

Proof of Theorem 3.21.

To estimate the optimal value in Proposition 3.24, note that for $z\in \mathbb{C}$, $z=\cos ({\it\phi})+i\sin ({\it\phi})=x+i\sqrt{1-x^{2}}$, we have

$$\begin{eqnarray}T_{k}(x)=\mathsf{Re}(z^{k})=\frac{1}{2}(z^{k}+\overline{z}^{k})=\frac{1}{2}((x+\sqrt{x^{2}-1})^{k}+(x-\sqrt{x^{2}-1})^{k}).\end{eqnarray}$$

Again take $a={\it\lambda}_{n}$, $b={\it\lambda}_{1}$, the smallest and largest eigenvalues of $A$. Then ${\it\kappa}={\it\kappa}(A)=b/a$. For

$$\begin{eqnarray}x=\frac{b+a}{b-a}=\frac{{\it\kappa}+1}{{\it\kappa}-1}\end{eqnarray}$$

we get

$$\begin{eqnarray}x^{2}-1=\frac{({\it\kappa}+1)^{2}-({\it\kappa}-1)^{2}}{({\it\kappa}-1)^{2}}=\frac{4{\it\kappa}}{({\it\kappa}-1)^{2}}\end{eqnarray}$$

and therefore

$$\begin{eqnarray}x+\sqrt{x^{2}-1}=\frac{{\it\kappa}+1+2\sqrt{{\it\kappa}}}{{\it\kappa}-1}=\frac{(\sqrt{{\it\kappa}}+1)^{2}}{{\it\kappa}-1}=\frac{\sqrt{{\it\kappa}}+1}{\sqrt{{\it\kappa}}-1}.\end{eqnarray}$$

Hence

$$\begin{eqnarray}T_{k}(x)\geq \frac{1}{2}\biggl(\frac{\sqrt{{\it\kappa}}+1}{\sqrt{{\it\kappa}}-1}\biggr)^{k}.\square\end{eqnarray}$$

(d) Conjugate gradient on random data.

Theorem 3.21 gives bounds for the rate of convergence of the conjugate gradient in terms of the condition ${\it\kappa}(A)$ of the matrix $A$. An average (or smoothed) analysis of ${\it\kappa}(A)$ would therefore result in an average (resp. smoothed) analysis for this rate of convergence. The question that emerges is which probability distribution we draw $A$ from.

In the case of unstructured or triangular matrices $A$, we have endowed the corresponding class of matrices with a Gaussian distribution. In our present case, with $A$ positive definite, this choice is no longer available. A look at our original computational problem may, however, shed some light. Matrix $A$ is obtained as $A=R^{\text{T}}R$. It then makes sense to consider $R$ as our primary random data – and for $R$ we can assume Gaussianity – and endow $A$ with the distribution inherited from that of $R$. In the particular case that $R$ is standard Gaussian, this inherited distribution of $A$ is known as the Wishart distribution with $q$degrees of freedom. It extends the chi-square distribution, since it coincides with a ${\it\chi}_{q}^{2}$ when $n=1$.

By the Singular Value Decomposition, there exist orthogonal matrices $U$ and $V$ and positive reals ${\it\sigma}_{1}\geq {\it\sigma}_{2}\geq \cdots \geq {\it\sigma}_{n}>0$ such that $R=U\mathsf{diag}({\it\sigma}_{1},{\it\sigma}_{2},\ldots ,{\it\sigma}_{n})V^{\text{T}}$. Then

$$\begin{eqnarray}A=R^{\text{T}}R=V\mathsf{diag}({\it\sigma}_{1}^{2},{\it\sigma}_{2}^{2},\ldots ,{\it\sigma}_{n}^{2})V^{\text{T}}\end{eqnarray}$$

and, consequently, ${\it\kappa}(A)={\it\kappa}_{\dagger }^{2}(R)$.

It follows from the above that the analysis of ${\it\kappa}(A)$ for a Wishart distribution reduces to the analysis of ${\it\kappa}_{\dagger }(R)$ when $R$ is standard Gaussian (and similarly for ‘non-centred Wishart distributions’).

Corollary 3.22 shows that the number $k$ of iterations needed to decrease $\Vert x_{0}-\bar{x}\Vert _{A}$ by a factor of ${\it\delta}$ is proportional to $\sqrt{{\it\kappa}(A)}$, that is, proportional to ${\it\kappa}_{\dagger }(R)$. Theorem 3.18 shows that for $q$ large enough and ${\it\lambda}=(n-1)/q$ small enough we have, for all $0<{\it\sigma}\leq 1$,

$$\begin{eqnarray}\sup _{\Vert \overline{R}\Vert \leq 1}\mathbb{E}_{R\sim N(\overline{R},{\it\sigma}^{2}\text{Id})}({\it\kappa}_{\dagger }(R))\leq \frac{20.1}{1-{\it\lambda}}.\end{eqnarray}$$

Thus, if $A$ is obtained as $R^{\text{T}}R$ for a large elongated rectangular matrix $R$, then we should expect to effect the decrease mentioned above with about

$$\begin{eqnarray}\frac{1}{2}\cdot \frac{20.1}{1-{\it\lambda}}\ln \biggl(\frac{1}{{\it\delta}}\biggr)\end{eqnarray}$$

iterations (where ${\it\lambda}=q/n$ is the elongation). Since each iteration takes $6n^{2}+O(n)$ arithmetic operations, the expected cost of the computation is

$$\begin{eqnarray}3n^{2}\frac{20.1}{1-{\it\lambda}}\ln \frac{1}{{\it\delta}}+O(n)=\frac{60.3n^{2}}{1-{\it\lambda}}\ln \frac{1}{{\it\delta}}+O(n).\end{eqnarray}$$

The leading term in this expression is smaller than the $\frac{2}{3}n^{3}$ operations performed by Gaussian elimination as long as

$$\begin{eqnarray}{\it\delta}\geq e^{-\frac{n(1-{\it\lambda})}{91}}.\end{eqnarray}$$

For large $n$ (and ${\it\lambda}$ not too close to 1) this bound produces very small values of ${\it\delta}$, and therefore Conj_Grad yields, on average (both for a Gaussian distribution of data and for Gaussian perturbations of arbitrary data), remarkably good approximations of $\bar{x}$.

We conclude that for reasonable ${\it\delta}>0$ and large $n$, computing an approximation of the solution $\bar{x}$ of $Ax=b$ using Algorithm Conj_Grad will typically be faster than computing $\bar{x}$ with Gaussian elimination.

(a) The feasibility problem and its condition.

A homogenization of the linear programming feasibility programme can be stated as follows. Given $A\in \mathbb{R}^{m\times n}$, $n>m$, consider the primal feasibility problem

and the dual feasibility problem

We say that $A$ is primal feasible or dual feasible when (PF) or (DF), respectively, is satisfied. In both cases we talk about strict feasibility when the satisfied inequality is strict.

Let ${\mathcal{F}}_{P}$ and ${\mathcal{F}}_{D}$ denote the sets of matrices $A$ where (PF) and (DF) are satisfied, respectively, and let ${\mathcal{F}}_{P}^{\circ }$ and ${\mathcal{F}}_{D}^{\circ }$ be the sets of strictly primal and strictly dual feasible matrices. Finally, let

$$\begin{eqnarray}{\rm\Sigma}:={\mathcal{F}}_{P}\cap {\mathcal{F}}_{D}.\end{eqnarray}$$

It is well known that, modulo a set of lower dimension, primal (dual) feasibility is equivalent to dual (resp. primal) strict infeasibility. The following result (see Bürgisser and Cucker Reference Bürgisser and Cucker2013, Theorem 6.14) is more precise. To state it, we let $\mathsf{int}(M)$, $\overline{M}$, and $\partial M=\overline{M}\setminus \mathsf{int}(M)$ denote the interior, closure, and boundary of a subset $M$ of Euclidean space.

Theorem 3.25. Both ${\mathcal{F}}_{P}$ and ${\mathcal{F}}_{D}$ are closed subsets of $\mathbb{R}^{m\times n}$. In addition, this space is partitioned as

$$\begin{eqnarray}\mathbb{R}^{m\times n}=\mathsf{int}({\mathcal{F}}_{P}^{\circ })\cup \mathsf{int}({\mathcal{F}}_{D}^{\circ })\cup {\rm\Sigma},\end{eqnarray}$$

and we have

$$\begin{eqnarray}{\rm\Sigma}=\partial {\mathcal{F}}_{P}=\partial {\mathcal{F}}_{D}.\end{eqnarray}$$

Furthermore,

$$\begin{eqnarray}\displaystyle {\mathcal{F}}_{P}^{\circ } & \supseteq & \displaystyle \mathsf{int}({\mathcal{F}}_{P}),\nonumber\\ \displaystyle {\mathcal{F}}_{P}^{\circ }\setminus \mathsf{int}({\mathcal{F}}_{P}) & = & \displaystyle {\mathcal{F}}_{P}^{\circ }\cap {\rm\Sigma}=\{A\in {\rm\Sigma}\mid \mathsf{rank}\,A<m\},\nonumber\\ \displaystyle \text{and}\quad {\mathcal{F}}_{D}^{\circ } & = & \displaystyle \mathsf{int}({\mathcal{F}}_{D}).\square\nonumber\end{eqnarray}$$

We see that for matrices in ${\rm\Sigma}$, arbitrarily small perturbations can lead to a change with respect to feasibility. In contrast, in the set ${\mathcal{D}}=\mathbb{R}^{m\times n}\setminus {\rm\Sigma}$ the following problem is well-defined:

$$\begin{eqnarray}\text{Given }A\in {\mathcal{D}},\text{decide whether }A\in {\mathcal{F}}_{P}^{\circ }\text{ or }A\in {\mathcal{F}}_{D}^{\circ }.\end{eqnarray}$$

We call this the polyhedral cone feasibility problem. For all $A\in {\rm\Sigma}$ the problem is ill-posed.

The polyhedral cone feasibility problem is a decision problem (as described in Section 1.1(c)). We will define for it a measure of condition in terms of a relative distance to ill-posedness. We have just seen that the set of ill-posed instances is the boundary between feasible and infeasible instances. This motivates the following definition.

Definition 3.26. Let $A\in \mathbb{R}^{m\times n}$ be given with non-zero columns $a_{i}$. Suppose $A\not \in {\rm\Sigma}$. Then, $A\in {\mathcal{F}}_{S}^{\circ }$ for $S\in \{P,D\}$. We define

$$\begin{eqnarray}{\rm\Delta}(A):=\sup \biggl\{{\it\delta}>0~\bigg|~\forall A^{\prime }\in \mathbb{R}^{m\times n}\;\biggl(\max _{i\leq n}\frac{\Vert a_{i}^{\prime }-a_{i}\Vert }{\Vert a_{i}\Vert }<{\it\delta}\;\Rightarrow \;A^{\prime }\in {\mathcal{F}}_{S}^{\circ }\biggr)\biggr\},\end{eqnarray}$$

where $a_{i}^{\prime }$ stands for the $i$th column of $A^{\prime }$. The GCC condition number of $A$ is defined as

$$\begin{eqnarray}\mathscr{C}(A):=\frac{1}{{\rm\Delta}(A)}.\end{eqnarray}$$

If $A\in {\rm\Sigma}$, we set ${\rm\Delta}(A)=0$ and $\mathscr{C}(A)=\infty$.

We have written the definition in such a way that it becomes clear that we measure the relative size of the perturbation for each row $a_{i}$, where the relativization is with respect to the norm of $a_{i}$. Also, it is clear from the definition that ${\rm\Delta}(A)$ is scale-invariant in the sense that

$$\begin{eqnarray}{\rm\Delta}([{\it\lambda}_{1}a_{1},\ldots ,{\it\lambda}_{n}a_{n}])={\rm\Delta}([a_{1},\ldots ,a_{n}])\quad \text{for}~{\it\lambda}_{i}>0.\end{eqnarray}$$

For the analysis of ${\rm\Delta}$ we may therefore assume, without loss of generality, that $\Vert a_{i}\Vert =1$ for all $i$. Hence we can see the matrix $A$ with columns $a_{1},\ldots ,a_{n}$ as an element in the product $(\mathbb{S}^{m-1})^{n}$ of spheres. The scale invariance of $\mathscr{C}(A)$, together with the characterization

$$\begin{eqnarray}\Vert A\Vert _{12}=\max _{i\leq n}\Vert a_{i}\Vert\end{eqnarray}$$

immediately yields the following result.

Proposition 3.27.

$$\begin{eqnarray}\mathscr{C}(A)=\frac{\Vert A\Vert _{12}}{d_{12}(A,{\rm\Sigma})}.\square\end{eqnarray}$$

A number of different characterizations of $\mathscr{C}(A)$ are shown in Bürgisser and Cucker (Reference Bürgisser and Cucker2013, Chapter 6).

(b) Algorithms and complexity.

During the late 1970s and early 1980s two algorithms proved that linear programming problems (with integer coefficients) can be solved in polynomial time: the ellipsoid algorithm and interior-point methods. These algorithms are naturally described as real number algorithms. But when considered over real (as opposed to integer) data their halting time is no longer bounded by a function of the input size only. Instead, the condition of the data at hand appears to play a role.

Of course, one may argue that the condition of the data is now replacing the bit-size of the largest component in this data. But the fact remains that this bit-size is directly related to the size of the input whereas the condition of the input is not (and it is not even known a priori).

Of these two methods, the most efficient (both in theory and in practice) is the interior-point method. There is a vast literature on this method (the book by Renegar Reference Renegar2000 being a very good introduction to them for general conic programming). For our purposes, it should suffice to state a result with the complexity of an interior-point algorithm solving the polyhedral cone feasibility problem. The proof of this result (and a self-contained introduction to interior-point methods for linear programming) can be found in Chapter 9 of Bürgisser and Cucker (Reference Bürgisser and Cucker2013).

Theorem 3.28. We exhibit an interior-point algorithm which, on input $A\not \in {\rm\Sigma}$, returns either a strictly feasible primal solution, i.e., a point $x\in \mathbb{R}^{n}$ such that $Ax=0$, $x>0$, or a strictly feasible dual solution, i.e., a point $y\in \mathbb{R}^{m}$ such that $A^{\text{T}}y<0$. In both cases, the algorithm halts after at most

$$\begin{eqnarray}O(\sqrt{n}(\log n+\log \mathscr{C}(A)))\end{eqnarray}$$

iterations. The total number of arithmetic operations is bounded by

$$\begin{eqnarray}O(n^{3.5}(\log n+\log \mathscr{C}(A))).\square\end{eqnarray}$$

A version of this result which takes into consideration the effects of finite precision, and shows complexity bounds of the same order, can be found in Cucker and Peña (Reference Cucker and Peña2002).

(c) Condition of random data.

Theorem 3.28 is an instance of a common phenomenon in the analysis of iterative algorithms (which we already met regarding the conjugate gradient method). The cost of each iteration is bounded by a function of the size or dimension of the input (in our case, it is $O(n^{3})$ where, we recall, $A\in \mathbb{R}^{m\times n}$ with $n>m$) but the number of iterations performed by the algorithm depends on the condition of the input at hand. To get an idea of how efficient the algorithm is ‘in general’, a probabilistic analysis is a must.

There is a sequence of papers on probabilistic analyses of $\mathscr{C}(A)$. The sharpest bound for the expectation of $\log \mathscr{C}(A)$ appears in Bürgisser, Cucker and Lotz (Reference Bürgisser, Cucker and Lotz2010). It shows that

$$\begin{eqnarray}\mathbb{E}_{A\sim N(0,\text{Id})}\ln \mathscr{C}(A)\leq 2\ln (m+1)+3.31.\end{eqnarray}$$

As for a smoothed analysis, Amelunxen and Bürgisser (Reference Amelunxen and Bürgisser2012) showed that

$$\begin{eqnarray}\sup _{\overline{A}\in (\mathbb{S}^{m-1})^{n}}\mathbb{E}_{A\sim N(\overline{A},{\it\sigma}^{2}\text{Id})}\log \mathscr{C}(A)=O\biggl(\log \frac{n}{{\it\sigma}}\biggr).\end{eqnarray}$$

These two bounds carry a very clear message: the cost of solving a typical polyhedral feasibility problem (or a typical perturbation of an arbitrary one) is $O(n^{3.5}\log n)$. The contribution of the condition of the input to this cost is very minor.

(a) The problem.

The linear programming problem in standard form is to find $x\in \mathbb{R}^{n}$ solving

where $A\in \mathbb{R}^{m\times n}$, $b\in \mathbb{R}^{m}$ and $c\in \mathbb{R}^{n}$ are given. The well-known Duality Theorem states that if a minimizer $x^{\ast }$ exists then a maximizer $y^{\ast }$ exists as well for the dual problem

and we have $c^{\text{T}}x^{\ast }=b^{\text{T}}y^{\ast }$. It additionally ensures that this is the case if and only if both the primal and dual are feasible. That is, there exists a pair $(x,y)\in \mathbb{R}^{n+m}$ such that $Ax=b$, $x\geq 0$ and $A^{\text{T}}y\leq c$. In this case, we say that the triple $d=(A,b,c)$ is feasible.

For any subset $B$ of $\{1,2,\ldots ,n\}$, we let $A_{B}$ denote the submatrix of $A$ obtained by removing from $A$ all the columns with index not in $B$. If $x\in \mathbb{R}^{n}$, $x_{B}$ is defined analogously. Also, let $N:=\{1,2,\ldots ,n\}\setminus B$. Then $A_{N}$ and $x_{N}$ are defined similarly. We say that $B$ is a basis for $d=(A,b,c)$ when $|B|=m$ and $A_{B}$ is invertible.

In this case, we may solve $A_{B}x=b$. Consider the point $x^{\ast }\in \mathbb{R}^{n}$ defined by $x_{N}^{\ast }=0$ and $x_{B}^{\ast }=x$. Clearly, $Ax^{\ast }=b$. We say that $x^{\ast }$ is the primal basic solution of (OP) associated with $B$. If, in addition, $x^{\ast }\geq 0$, which is equivalent to $x_{B}^{\ast }\geq 0$, then we say that $x^{\ast }$ is a primal basic feasible solution.

Similarly, we may solve $A_{B}^{\text{T}}y=c_{B}$. The point $y^{\ast }$ thus obtained is said to be the dual basic solution of  associated with $B$. If, in addition, $A^{\text{T}}y^{\ast }\leq c$, then $y^{\ast }$ is said to be a dual basic feasible solution.

There is a rationale for the word ‘optimal’ in this definition. If $B$ is an optimal basis and $x^{\ast },y^{\ast }$ are its associated basic optimal solutions, then $x^{\ast }$ and $y^{\ast }$ are optimizers for the problems (OP) and (OD), respectively. Furthermore, optimal bases exist (Bürgisser and Cucker Reference Bürgisser and Cucker2013, Theorem 8.13).

Theorem 3.30 suggests a strategy to compute optimizers for feasible data with full-rank matrices, namely, to compute an optimal basis $B$ and its associated basic solutions $x^{\ast }$ and $y^{\ast }$. These points are vertices in the polyhedra of feasible points for (OP) and (OD), respectively. The simplex algorithm searches for such a solution by examining a cleverly constructed sequence in this set of vertices. Interior-point methods instead produce a sequence of points in the interior of these polyhedra that leads to the desired vertex.

(b) Ill-posedness and condition.

The problems of computing optimizers and optimal bases are, strictly speaking, different. But both their sets of ill-posed data and their condition numbers are, as one can expect, very closely related.

Even though we will only be interested in feasible data (in this example we compute optimizers for data of this type) we still need to mention infeasible data at the boundary of the set of feasible data. A more detailed discussion of ill-posedness, with connections to the notion of degeneracy, can be found in Bürgisser and Cucker (Reference Bürgisser and Cucker2013, § 8.4).

We observe that infeasible triples should be considered as ill-posed if and only if arbitrarily small perturbations can make these triples feasible (and, consequently, create optimizers for them). We thus define

$$\begin{eqnarray}{\rm\Sigma}_{{\mathcal{I}}}:=\{d\mid d\text{ is infeasible and }\forall {\it\varepsilon}>0\;\exists d^{\prime }\text{feasible with }\Vert d-d^{\prime }\Vert \leq {\it\varepsilon}\}\end{eqnarray}$$

and call this the set of infeasible ill-posed triples.

For the problem of computing optimizers one can show that the property of having more than one pair $(x^{\ast },y^{\ast })$ of optimizers can be destroyed by arbitrarily small perturbations of the triple $d=(A,b,c)$. Also, for problems having a unique pair of optimizers, small perturbations of the data will preserve this character and the resulting pair of optimizers will be close to the original one. Consequently, we define the set of ill-posed triples for this problem to be

$$\begin{eqnarray}{\rm\Sigma}_{\mathsf{opt}}:=\{d\mid d\text{ has at least two optimizers }(x,y)\}\cup {\rm\Sigma}_{{\mathcal{I}}}.\end{eqnarray}$$

Let us focus next on the problem of computing an optimal basis. In contrast with the preceding problem, this one is discrete-valued: the output for any given datum $d$ is an $m$-element subset in $[n]$. Therefore, we can take the set ${\rm\Sigma}_{\mathsf{OB}}$ of ill-posed triples as the boundary between these possible outputs. In other words, if $\mathsf{OB}(d)$ denotes the set of optimal bases for a triple $d$, we define

$$\begin{eqnarray}{\rm\Sigma}_{\mathsf{OB}}:=\{d\mid \forall {\it\varepsilon}>0\;\exists d^{\prime }~\text{such that}~\mathsf{OB}(d)\neq \mathsf{OB}(d^{\prime })~\text{and}~\Vert d-d^{\prime }\Vert \leq {\it\varepsilon}\}.\end{eqnarray}$$

Note that ${\rm\Sigma}_{{\mathcal{I}}}$ corresponds to triples $d$ in ${\rm\Sigma}_{\mathsf{OB}}$ with $\mathsf{OB}(d)=\emptyset$.

Following the above we say that $d$ is feasible well-posed when there exists a unique optimal solution $(x^{\ast },y^{\ast })$ for $d$. In this case, we write $d\in {\mathcal{W}}$. If more than one optimal solution exists, we say that $d$ is feasible ill-posed.

Let ${\mathcal{B}}:=\{B\subseteq [n]\mid |B|=m\}$. Theorem 3.30(c) implies that if $d\in {\mathcal{W}}$ is feasible well-posed with optimal solution $(x^{\ast },y^{\ast })$, then there is a unique $B\in {\mathcal{B}}$ such that $(x^{\ast },y^{\ast })$ is the basic optimal solution associated with $B$. We may then define, for any $B\in {\mathcal{B}}$, the set

$$\begin{eqnarray}{\mathcal{W}}_{B}:=\{d\in {\mathcal{W}}\mid B\text{ is the optimal basis for }d\}.\end{eqnarray}$$

The set ${\mathcal{W}}$ is thus partitioned by the family $\{{\mathcal{W}}_{B}\mid B\in {\mathcal{B}}\}$. Figure 3.1 schematically summarizes the situation in the space of data.

Figure 3.1.

The situation in the space of data.

The left of the diagram is not relevant to us: it represents the set ${\mathcal{I}}$ of infeasible triples. At the right of the diagram we have the set ${\mathcal{W}}$ of well-posed feasible triples, which is partitioned into the subsets ${\mathcal{W}}_{B}$. The dotted lines separating these subsets correspond to data with more than one optimal basis. The two continuous segments are the boundary between feasible and infeasible triples. The union of these lines – the dotted lines and the continuous segments – forms the set ${\rm\Sigma}_{\mathsf{opt}}={\rm\Sigma}_{\mathsf{OB}}$ of ill-posed triples for both the optimal solution and the optimal basis problems.

The picture we just described allows for a very natural definition of condition number. We only need a norm in the space of data triples, and for this we take

$$\begin{eqnarray}\Vert d\Vert =\biggl\|\left(\begin{array}{@{}cc@{}}A & b\\ c^{\text{T}} & 0\end{array}\right)\biggr\|.\end{eqnarray}$$

Definition 3.32. Let $d$ be feasible well-posed. We define the distance to ill-posedness to be

$$\begin{eqnarray}{\it\varrho}(d)=\inf \{\Vert {\rm\Delta}d\Vert \mid d+{\rm\Delta}d\in {\rm\Sigma}_{\mathsf{opt}}\}.\end{eqnarray}$$

Moreover, we define the RCC condition number $\mathscr{K}(d)$ of $d$ by

$$\begin{eqnarray}\mathscr{K}(d)=\frac{\Vert d\Vert }{{\it\varrho}(d)}.\end{eqnarray}$$

If $d$ is feasible but not well-posed, we let ${\it\varrho}(d)=0$ and $\mathscr{K}(d)=\infty$.

(c) Algorithms and complexity.

There is no shortage of algorithms for the linear programming optimization problem. A distinctive feature of the one described in Cheung and Cucker (Reference Cheung and Cucker2006) is that it works with finite precision. Another feature is that both complexity and accuracy are analysed in terms of $\mathscr{K}(d)$ (which, somewhat conveniently, was introduced and studied in Cheung and Cucker Reference Cheung and Cucker2004).

Theorem 3.33. There exists a finite-precision algorithm which, given a feasible well-posed triple $d=(A,b,c)$, finds an optimal basis $B$ for $d$. The machine precision varies during the execution of the algorithm (allowing for rereading the input $d$ with higher accuracy) and the finest required precision satisfies

$$\begin{eqnarray}{\it\epsilon}_{\mathsf{mach}}^{\ast }=\frac{1}{O((m+n+1)^{26}\mathscr{K}(d)^{16})}.\end{eqnarray}$$

The total number of arithmetic operations is bounded by

$$\begin{eqnarray}O(n^{3.5}(\log n+\log \mathscr{K}(d))\log \log (\mathscr{K}(d)+1)).\square\end{eqnarray}$$

The feasibility of the triple $d$ needs to be ascertained before the execution of the algorithm in Theorem 3.33, but it is shown in Cheung and Cucker (Reference Cheung and Cucker2006) that this can be done within the accuracy and complexity bounds in the theorem.

The bound for ${\it\epsilon}_{\mathsf{mach}}^{\ast }$ shows that the largest number of digits needed for the floating-point numbers occurring during the computation is bounded as

(3.15)$$\begin{eqnarray}O(\log (n)+\log \mathscr{K}(d)).\end{eqnarray}$$

In Bürgisser and Cucker (Reference Bürgisser and Cucker2013, Chapter 11) there is a version of Theorem 3.33 for infinite precision. Needless to say, the corresponding proof is much simplified. There is also a description of a different algorithm with a slightly smaller complexity bound (without the $\log \log \mathscr{K}(d)$ factor) but requiring an initial pair $(x_{0},y_{0})$ satisfying some properties.

One final remark is essential. Once we have computed the optimal basis $B$ for the triple $d$ we still want to compute the optimizers. For this, recall that we need to compute

$$\begin{eqnarray}x_{B}^{\ast }=A_{B}^{-1}b\quad \text{and}\quad y^{\ast }=A_{B}^{-\text{T}}c_{B}.\end{eqnarray}$$

Whereas the complexity of this computation is dominated by that of computing $B$, its accuracy analysis (such as the one in Section 3.1) depends on ${\it\kappa}(A_{B})={\it\kappa}(A_{B}^{\text{T}})$. However, the inequality

$$\begin{eqnarray}{\it\kappa}(A_{B})\leq \mathscr{K}(d)\end{eqnarray}$$

(see Bürgisser and Cucker Reference Bürgisser and Cucker2013, Proposition 11.22) reduces all our needs to one: understanding how large $\mathscr{K}(d)$ is in general.

(d) Probabilistic analyses.

An average analysis of $\mathscr{K}(d)$ was obtained in Cheung and Cucker (Reference Cheung and Cucker2013).

Theorem 3.34. We have

$$\begin{eqnarray}\mathbb{E}_{d\sim N(0,\text{Id})}(\log \mathscr{K}(d)\mid d\in {\mathcal{W}})\leq \frac{5}{2}\log (n+1)+\log (m+1)+\log 24e,\end{eqnarray}$$

as well as

$$\begin{eqnarray}\mathbb{E}_{d\sim N(0,\text{Id})}(\log \log \mathscr{K}(d)\mid d\in {\mathcal{W}})=\log \log n+O(1)\end{eqnarray}$$

and

$$\begin{eqnarray}\mathbb{E}_{d\sim N(0,\text{Id})}(\log \mathscr{K}(d)\log \log \mathscr{K}(d)\mid d\in {\mathcal{W}})=\log n\log \log n+O(1).\square\end{eqnarray}$$

Theorems 3.33 and 3.34 combine to yield the following.

Corollary 3.35. There exists a finite-precision algorithm that, given a triple $d=(A,b,c)$, finds an optimal basis $B$ for $d$. The average cost of the algorithm, on Gaussian data $d$, conditioned to $d$ being feasible well-posed, is bounded by

$$\begin{eqnarray}O(n^{3.5}\log n\log \log n).\end{eqnarray}$$

The machine precision varies during the execution of the algorithm. The largest number of digits needed for the floating-point numbers occurring in the computation is bounded on average as

$$\begin{eqnarray}\mathbb{E}_{d\sim N(0,\text{Id})}\log {\it\epsilon}_{\mathsf{mach}}^{\ast }(d)=O(\log n).\square\end{eqnarray}$$

Remarkably, as of today, the problem of showing a smoothed analysis version of Theorem 3.34 remains open.

(a) The problem.

Let $\mathbf{d}=(d_{1},\ldots ,d_{n})$ with $d_{i}\in \mathbb{N}$, $d_{i}\geq 2$. We let ${\mathcal{H}}_{\mathbf{d}}$ denote the linear space of homogeneous polynomial systems $f=(f_{1},\ldots ,f_{n})$, where $f_{i}\in \mathbb{C}[X_{0},\ldots ,X_{n}]$ is homogeneous of degree $d_{i}$.

Steve Smale (Reference Smale, Arnold, Atiyah, Lax and Mazur2000) proposed a list of 18 problems for the mathematicians of the twenty-first century. The 17th problem asks whether there exists an algorithm that, given $f\in {\mathcal{H}}_{\mathbf{d}}$, returns an approximate zero of $f$ in time polynomial (on average) in the size of $f$. Here ‘approximate zero’ is in the sense of Definition 3.36 (with respect to a projective Newton’s method tailored to this problem by Mike Shub Reference Shub, Hirsch, Marsden and Shub1993) and ‘average’ is with respect to the standard normal distribution for the coefficients of the $f_{i}$ when these polynomials are written in the Weyl basis (a simple modification to the usual basis that makes the associated Hermitian product unitarily invariant: for details see Bürgisser and Cucker Reference Bürgisser and Cucker2013, Section 16.1). The ‘size’ of $f$ is the total number of these coefficients or, equivalently, the dimension $N$ of ${\mathcal{H}}_{\mathbf{d}}$ over $\mathbb{C}$. This yields

$$\begin{eqnarray}N=\mathop{\sum }_{i=1}^{n}\binom{n+d_{i}}{d_{i}}.\end{eqnarray}$$

Attempts to give an affirmative answer to Smale’s 17th problem have revolved around homotopy algorithms. The works by Beltrán and Pardo (Reference Beltrán and Pardo2009, Reference Beltrán and Pardo2011) provided a partial answer in the form of a randomized algorithm whose expected running time (with respect to both the randomizations and the input $f$) is polynomial in $N$. A further development was given in Bürgisser and Cucker (Reference Bürgisser and Cucker2011), where a deterministic algorithm working in nearly polynomial time was described. All these works heavily relied on the analysis of a condition number. In what follows we define this condition number, specify the choice of step-length in the homotopy algorithm, and show (sketching the proofs) how a smoothed analysis leads to the results mentioned above.

(b) Condition numbers.

To define the condition for Smale’s 17th problem we use the setting described in Section 1.1(b). Here the manifold of inputs is ${\mathcal{H}}_{\mathbf{d}}$, the manifold of outputs is $\mathbb{P}^{n}:=\mathbb{P}(\mathbb{C}^{n+1})$, and the submanifold $V\subset {\mathcal{H}}_{\mathbf{d}}\times \mathbb{P}(\mathbb{C}^{n+1})$ of solutions is the set of pairs $(f,{\it\zeta})$ with $f({\it\zeta})=0$.

It turns out that the set of ill-posed solutions is then (recall (1.6))

(3.16)$$\begin{eqnarray}{\rm\Sigma}^{\prime }=\big\{(f,{\it\zeta})\in V\mid \mathsf{rank}\,Df({\it\zeta})<n\big\}\end{eqnarray}$$

and that, for $(f,{\it\zeta})\in V$, $(f,{\it\zeta})\in {\rm\Sigma}^{\prime }$ if and only if $z$ is a multiple zero of $f$.

Recall the projections ${\it\pi}_{1}$ and ${\it\pi}_{2}$:

We also define the set ${\rm\Sigma}:={\it\pi}_{1}({\rm\Sigma}^{\prime })$ consisting of those systems $f\in {\mathcal{H}}_{\mathbf{d}}$ that have some multiple zero ${\it\zeta}$.

For a pair $(f_{0},{\it\zeta}_{0})\in V\setminus {\rm\Sigma}^{\prime }$ we have defined the absolute condition number $\mathsf{acond}(f_{0},{\it\zeta}_{0})$ to be the operator norm of the solution map $G:U\rightarrow \mathbb{P}^{n},G(f)={\it\zeta}$ that locally inverts ${\it\pi}_{1}:V\rightarrow {\mathcal{H}}_{\mathbf{d}}$ around $(f_{0},{\it\zeta}_{0})$ on an open neighbourhood $U$ of $f_{0}$ in ${\mathcal{H}}_{\mathbf{d}}$. To make the condition independent of the scaling of $f$ we take, for $(f,{\it\zeta})\in V\setminus {\rm\Sigma}^{\prime }$,

(3.17)$$\begin{eqnarray}{\it\mu}(f,{\it\zeta}):=\Vert f\Vert \cdot \Vert DG(f)\Vert ,\end{eqnarray}$$

where the operator norm is with respect to the standard Hermitian product (with respect to the Weyl’s basis) on ${\mathcal{H}}_{\mathbf{d}}$ and the Hermitian product on $T_{{\it\zeta}}\mathbb{P}^{n}$ coming from the standard one on $\mathbb{C}^{n+1}$.

The condition number ${\it\mu}(f,{\it\zeta})$ can also be expressed in terms of the derivative $Df({\it\zeta}):\mathbb{C}^{n+1}\rightarrow \mathbb{C}^{n}$ of $f:\mathbb{C}^{n+1}\rightarrow \mathbb{C}^{n}$. Euler’s identity implies that ${\it\zeta}$ lies in the kernel of $Df({\it\zeta})$. So, if $Df({\it\zeta})$ is of full rank, its kernel is $\mathbb{C}{\it\zeta}$ and $T_{{\it\zeta}}$ equals its orthogonal complement. The inverse of the restriction $Df({\it\zeta})_{|T_{{\it\zeta}}}$ is described by the Moore–Penrose inverse $Df({\it\zeta})^{\dagger }$, and we have

$$\begin{eqnarray}{\it\mu}(f,{\it\zeta})=\Vert f\Vert \cdot \Vert Df({\it\zeta})^{\dagger }\mathsf{diag}(\Vert {\it\zeta}\Vert ^{d_{i}-1})\Vert .\end{eqnarray}$$

The function ${\it\mu}$ is a condition number in the classical sense: it bounds the speed of change of the solutions when the data change. Indeed, if $[0,1]\rightarrow V,\,t\mapsto (f_{t},{\it\zeta}_{t})\in V$ is a smooth curve such that $f_{t}\in \mathbb{S}({\mathcal{H}}_{\mathbf{d}})$ for all $t$, then

$$\begin{eqnarray}\Vert \dot{{\it\zeta}}_{t}\Vert \leq {\it\mu}(f_{t},{\it\zeta}_{t})\Vert \dot{f_{t}}\Vert .\end{eqnarray}$$

For technical reasons, it has proved to be convenient to work with a slight modification to ${\it\mu}$. We define the normalized condition number ${\it\mu}_{\text{norm}}(f,{\it\zeta})$ for $(f,{\it\zeta})\in V$ by

(3.18)$$\begin{eqnarray}{\it\mu}_{\text{norm}}(f,{\it\zeta}):=\Vert f\Vert \cdot \Vert Df({\it\zeta})^{\dagger }\mathsf{diag}(\sqrt{d_{i}}\Vert {\it\zeta}\Vert ^{d_{i}-1})\Vert .\end{eqnarray}$$

Note that the introduction of the $\sqrt{d_{i}}$ factors is the only change compared with ${\it\mu}(f,{\it\zeta})$. In particular, we have

(3.19)$$\begin{eqnarray}{\it\mu}(f,{\it\zeta})\leq {\it\mu}_{\text{norm}}(f,{\it\zeta})\leq \sqrt{D}\,{\it\mu}(f,{\it\zeta})\end{eqnarray}$$

where $D:=\max _{i}d_{i}$. Furthermore, since our algorithms work with approximations of zeros only, it will be convenient to extend ${\it\mu}_{\text{norm}}(f,z)$ to the case that $z$ is not a zero of $f$.

Definition 3.37. For $f\in {\mathcal{H}}_{\mathbf{d}}$ and $z\in (\mathbb{C}^{n+1})_{\ast }$ we define the normalized condition number ${\it\mu}_{\text{norm}}(f,z)$ as

$$\begin{eqnarray}{\it\mu}_{\text{norm}}(f,z):=\Vert f\Vert \cdot \Vert Df(z)|_{T_{z}}^{-1}\mathsf{diag}(\sqrt{d_{i}}\,\Vert z\Vert ^{d_{i}-1})\Vert\end{eqnarray}$$

if $Df(z)|_{T_{z}}:T_{z}\rightarrow \mathbb{C}^{n}$ is invertible. Otherwise, we set ${\it\mu}_{\text{norm}}(f,z):=\infty$.

By Euler’s identity, if $f(z)=0$, this is consistent with the original definition (3.18):

$$\begin{eqnarray}{\it\mu}_{\text{norm}}(f,z)=\Vert f\Vert \cdot \Vert Df(z)^{\dagger }\mathsf{diag}(\sqrt{d_{i}}\Vert z\Vert ^{d_{i}-1})\Vert .\end{eqnarray}$$

We note, however, that using this formula in the case $f(z)\neq 0$ would lead to a different notion of ${\it\mu}_{\text{norm}}(f,z)$ (although the difference is small if $z$ is close to ${\it\zeta}$).

The condition number ${\it\mu}_{\text{norm}}$ has several properties. It is invariant under scaling: that is, for ${\it\lambda}_{1},{\it\lambda}_{2}\in \mathbb{C}_{\ast }$, we have ${\it\mu}_{\text{norm}}({\it\lambda}_{1}f,{\it\lambda}_{2}z)={\it\mu}_{\text{norm}}(f,z)$. It is bounded below: ${\it\mu}_{\text{norm}}(f,z)\geq \sqrt{n}$ for all $(f,z)\in {\mathcal{H}}_{\mathbf{d}}\times (\mathbb{C}^{n+1})_{\ast }$. And, crucially, it controls the size of the basin of quadratic attraction of $f$ around a zero ${\it\zeta}$ (see Bürgisser and Cucker Reference Bürgisser and Cucker2011, Theorem 2.2). Let $d_{\mathbb{P}}$ denote the Riemannian distance in projective space $\mathbb{P}^{n}$ (i.e. the angle).

Theorem 3.38. Let $(f,{\it\zeta})\in V$ and $x\in \mathbb{P}^{n}$ such that

$$\begin{eqnarray}d_{\mathbb{P}}(x,{\it\zeta})\leq \frac{3-\sqrt{7}}{D^{3/2}{\it\mu}_{\text{norm}}(f,{\it\zeta})}.\end{eqnarray}$$

Then $x$ is an approximate zero of $f$ with associated zero ${\it\zeta}$.◻

(c) The homotopy algorithm.

It is this last theorem (along with a bound on the Lipschitz constant for ${\it\mu}_{\text{norm}}$) that allows for a very simple choice of step-size in the homotopy, namely,

$$\begin{eqnarray}\displaystyle {\rm\Delta}{\it\tau} & := & \displaystyle \frac{6.67\times 10^{-3}}{d_{\mathbb{S}}(f,g)D^{3/2}{\it\mu}_{\text{norm}}^{2}(q_{i},z_{i})},\nonumber\\ \displaystyle {\it\tau}_{i+1} & := & \displaystyle \min \{{\it\tau}_{i}+{\rm\Delta}{\it\tau},1\},\nonumber\end{eqnarray}$$

where $d_{\mathbb{S}}(f,g)$ denotes the angle between $f$ and $g$. Algorithm 2 – given a system $f$ and an initial pair $(g,{\it\zeta})$ – endowed with this choice of step-size returns an approximate zero of $f$. The following result (see Bürgisser and Cucker Reference Bürgisser and Cucker2011, Theorem 3.1, or Bürgisser and Cucker Reference Bürgisser and Cucker2013, Theorem 17.3) shows that this is so, and gives a first estimate of the cost of the homotopy.

Theorem 3.39. Suppose that the segment $[f,g]$ does not intersect the discriminant variety ${\rm\Sigma}$. Then the homotopy stops after at most $K$ iterations with

$$\begin{eqnarray}K:=K(f,g,{\it\zeta}):=\leq C\,D^{3/2}\,d_{\mathbb{ S}}(f,g)\int _{0}^{1}{\it\mu}_{\text{norm}}^{2}(q_{{\it\tau}},{\it\zeta}_{{\it\tau}})\,\text{d}{\it\tau}\end{eqnarray}$$

for some constant $C$. The returned point $z$ is an approximate zero of $f$ with associated zero ${\it\zeta}_{1}$.◻

The number $K$ of iterations performed by the algorithm is closely related to the total cost since it is easy to see (Bürgisser and Cucker Reference Bürgisser and Cucker2013, Proposition 16.32) that each iteration is performed with $O(N)$ arithmetic operations.

We note that the ingredients in the integral in Theorem 3.39, the angular distances and the condition number ${\it\mu}_{\text{norm}}$ are scale-invariant. This is natural since the number of iterations of the homotopy algorithm is invariant under scaling of $f$ and/or $g$.

(d) Random data and/or initial pairs.

The bound in Theorem 3.39 depends on an integral which, in turn, depends on $f$ and the pair $(g,{\it\zeta})$. In order to obtain more concrete bounds one needs to eliminate this integral, a computation that requires either specific choices for $f$ or $(g,{\it\zeta})$ or assumes them drawn from a probability distribution. In the case of $f$, the latter amounts to working in the average complexity framework of Smale’s 17th problem. In the case of $(g,{\it\zeta})$, it allows for the use of randomized algorithms. For now, we will consider randomizing both $f$ and the pair $(g,{\it\zeta})$.

Note that to do the latter we need two additional ingredients: a probability measure on the solution manifold $V$ and an algorithmic procedure to draw pairs from this measure.

For the first, we adopt the measure ${\it\rho}_{\mathsf{st}}$ (we call it the standard distribution on $V$) that is informally described by the following procedure:

$$\begin{eqnarray}\displaystyle & & \displaystyle \mathtt{draw}~g~\mathtt{from}~N(0,\text{Id})\nonumber\\ \displaystyle & & \displaystyle \mathtt{draw}~{\it\zeta}~\mathtt{from~the~uniform~distribution~on~the}~{\mathcal{D}}~\mathtt{zeros~of}~g\nonumber\end{eqnarray}$$

Here ${\mathcal{D}}:=\prod _{i=1}^{n}d_{i}$ is the Bézout number. Systems in ${\mathcal{H}}_{\mathbf{d}}$ have, with probability 1, ${\mathcal{D}}$ different zeros. While this description has the advantage of clarity, it does not suggest an algorithmic procedure to actually draw pairs from ${\it\rho}_{\mathsf{st}}$: for this we should have to have at hand the ${\mathcal{D}}$ zeros of $g$!

Such a procedure (along with the randomization approach for Smale’s 17th problem) was devised by Beltrán and Pardo. It is described in Beltrán and Pardo (Reference Beltrán and Pardo2011) (a more efficient variant can be found in Bürgisser and Cucker Reference Bürgisser and Cucker2013, Section 17.6). It allows for a Las Vegas type algorithm – denoted by $\mathsf{LV}$ – that, with input $f\in {\mathcal{H}}_{\mathbf{d}}$, first draws $(g,{\it\zeta})$ from $V$ and then computes an approximate zero of $f$ by passing $f$ and the initial pair $(g,{\it\zeta})$ to Algorithm 2.

The central question is: What is the expected complexity of $\mathsf{LV}$? Here the expectation is with respect to both the input $f$ and the pair $(g,{\it\zeta})$. The expected number of iterations of the call to Algorithm 2 is then

$$\begin{eqnarray}\mathbb{E}_{f\sim N(0,\text{Id})}\;\mathbb{E}_{(g,{\it\zeta})\sim {\it\rho}_{\mathsf{st}}}K(f,g,{\it\zeta}),\end{eqnarray}$$

a quantity that, using Theorem 3.39 and the description of ${\it\rho}_{\mathsf{st}}$, can be bounded by

(3.20)$$\begin{eqnarray}CD^{3/2}d_{\mathbb{ S}}(f,g)\mathbb{E}_{f\sim N(0,\text{Id})}\;\mathbb{E}_{g\sim N(0,\text{Id})}\frac{1}{{\mathcal{D}}}\mathop{\sum }_{j=1}^{{\mathcal{D}}}\int _{0}^{1}{\it\mu}_{\text{norm}}^{2}(q_{{\it\tau}},{\it\zeta}_{{\it\tau}}^{(j)})\,\text{d}{\it\tau}.\end{eqnarray}$$

Here $\{(q_{{\it\tau}},{\it\zeta}_{{\it\tau}}^{(j)})\}_{{\it\tau}\in [0,1]}$ is the lifting of the segment $[g,f]$ corresponding to the $j$th zero of $g$. To further analyse (3.20), we note that neither $f$ nor $g$ explicitly appears in $q_{{\it\tau}}$. They both appear appear implicitly through a complicated function of ${\it\tau}$. This function can be simplified by changing the parametrization of the segment $[g,f]$ from the ratio ${\it\tau}$ of angles to a ratio $t$ of Euclidean distances. Elementary geometry shows that this change of parametrization is given by

$$\begin{eqnarray}t({\it\tau})=\frac{\Vert g\Vert }{\Vert f\Vert (\sin {\it\alpha}\cot ({\it\tau}{\it\alpha})-\cos {\it\alpha})+\Vert g\Vert },\end{eqnarray}$$

where ${\it\alpha}=d_{\mathbb{S}}(f,g)$; elementary calculus shows that (3.20) is bounded by

$$\begin{eqnarray}CD^{3/2}d_{\mathbb{ S}}(f,g)\mathbb{E}_{f\sim N(0,\text{Id})}\;\mathbb{E}_{g\sim N(0,\text{Id})}\Vert f\Vert \Vert g\Vert \frac{1}{{\mathcal{D}}}\mathop{\sum }_{j=1}^{{\mathcal{D}}}\int _{0}^{1}\frac{{\it\mu}_{\text{norm}}^{2}(q_{t},{\it\zeta}_{t}^{(j)})}{\Vert q_{t}\Vert ^{2}}\,\text{d}t.\end{eqnarray}$$

A further argument (involving Gaussian distributions truncated at $\sqrt{2N}$ and scale invariance: see Bürgisser and Cucker Reference Bürgisser and Cucker2011 for details) shows that this last expression is bounded by

$$\begin{eqnarray}2CD^{3/2}d_{\mathbb{ S}}(f,g)N\mathbb{E}_{f\sim N(0,\text{Id})}\;\mathbb{E}_{g\sim N(0,\text{Id})}\frac{1}{{\mathcal{D}}}\mathop{\sum }_{j=1}^{{\mathcal{D}}}\int _{0}^{1}\frac{{\it\mu}_{\text{norm}}^{2}(q_{t},{\it\zeta}_{t}^{(j)})}{\Vert q_{t}\Vert ^{2}}\,\text{d}t.\end{eqnarray}$$

Further, Tonelli’s theorem allows us to move the integral to the left of the expectations and obtain

(3.21)$$\begin{eqnarray}2CD^{3/2}d_{\mathbb{ S}}(f,g)N\int _{0}^{1}\mathbb{E}_{f\sim N(0,\text{Id})}\;\mathbb{E}_{g\sim N(0,\text{Id})}\frac{{\it\mu}_{\mathsf{av}}^{2}(q_{t})}{\Vert q_{t}\Vert ^{2}}\,\text{d}t,\end{eqnarray}$$

where

$$\begin{eqnarray}{\it\mu}_{\mathsf{av}}^{2}(q_{t}):=\frac{1}{{\mathcal{D}}}\mathop{\sum }_{j=1}^{{\mathcal{D}}}{\it\mu}_{\text{norm}}^{2}(q_{t},{\it\zeta}_{t}^{(j)})\end{eqnarray}$$

is the mean square condition number of $q_{t}$. The advantage of (3.21) is that, for a fixed $t\in [0,1]$, the distribution of $q_{t}=tf+(1-t)g$ is Gaussian with centre $0$ and covariance matrix $(t^{2}+(1-t)^{2})\text{Id}$. That is, we can replace (3.21) with the simpler expression

(3.22)$$\begin{eqnarray}2CD^{3/2}d_{\mathbb{ S}}(f,g)N\int _{0}^{1}\mathbb{E}_{q_{t}\sim N(0,(t^{2}+(1-t)^{2})\text{Id})}\;\frac{{\it\mu}_{\mathsf{av}}^{2}(q_{t})}{\Vert q_{t}\Vert ^{2}}\,\text{d}t.\end{eqnarray}$$

We are therefore left with a need to bound the expectation in this expression (for a fixed $t\in [0,1]$). Before proceeding with this task, though, it will be useful to discuss the more involved framework of a deterministic choice of the initial system $g$.

For the time being, we will still assume such a $g$ in the sphere $\mathbb{S}({\mathcal{H}}_{\mathbf{d}})$ and a zero ${\it\zeta}$ of $g$ which is randomly chosen (as is the input $f$). In this case, the reasoning above, replacing the starting bound (3.20) with

$$\begin{eqnarray}CD^{3/2}d_{\mathbb{ S}}(f,g)\mathbb{E}_{f\sim N(0,\text{Id})}\;\frac{1}{{\mathcal{D}}}\mathop{\sum }_{j=1}^{{\mathcal{D}}}\int _{0}^{1}{\it\mu}_{\text{norm}}^{2}(q_{{\it\tau}},{\it\zeta}_{{\it\tau}}^{(j)})\,\text{d}{\it\tau},\end{eqnarray}$$

ends with the bound

(3.23)$$\begin{eqnarray}CD^{3/2}d_{\mathbb{ S}}(f,g)\sqrt{2N}\Vert g\Vert \int _{0}^{1}\mathbb{E}_{q_{t}\sim N((1-t)g,t^{2}\text{Id})}\;\frac{{\it\mu}_{\mathsf{av}}^{2}(q_{t})}{\Vert q_{t}\Vert ^{2}}\,\text{d}t.\end{eqnarray}$$

Again, we are left with a need to bound the expectation in this expression, but this time the expectation leads to a smoothed (as opposed to average) analysis.

(e) Smoothed analysis of${\it\mu}_{\mathsf{av}}^{2}(q)/\Vert q\Vert ^{2}$.

This is the most technically involved part of the example. The result (see Bürgisser and Cucker Reference Bürgisser and Cucker2011 § 8, for details of the proof) is as follows.

Theorem 3.40. For $\overline{q}\in {\mathcal{H}}_{\mathbf{d}}$ and ${\it\sigma}>0$ we have

$$\begin{eqnarray}\mathbb{E}_{q\sim N(\overline{q},{\it\sigma}^{2}\text{Id})}\biggl(\frac{{\it\mu}_{\mathsf{av}}^{2}(q)}{\Vert q\Vert ^{2}}\biggr)\leq \frac{e(n+1)}{2{\it\sigma}^{2}}.\end{eqnarray}$$

Sketch of the proof.

For ${\it\zeta}\in (\mathbb{C}^{n+1})_{\ast }$ we consider the subspace $R_{{\it\zeta}}$ of ${\mathcal{H}}_{\mathbf{d}}$ consisting of all systems $h$ that vanish at ${\it\zeta}$ to higher order:

$$\begin{eqnarray}R_{{\it\zeta}}:=\big\{h\in {\mathcal{H}}_{\mathbf{d}}\mid h({\it\zeta})=0,Dh({\it\zeta})=0\big\}.\end{eqnarray}$$

We further decompose the orthogonal complement $R_{{\it\zeta}}^{\bot }$ of $R_{{\it\zeta}}$ in ${\mathcal{H}}_{\mathbf{d}}$. Let $L_{{\it\zeta}}$ denote the subspace of $R_{{\it\zeta}}^{\bot }$ consisting of the systems vanishing at ${\it\zeta}$, and let $C_{{\it\zeta}}$ denote its orthogonal complement in $R_{{\it\zeta}}^{\bot }$. Then we have an orthogonal decomposition

(3.24)$$\begin{eqnarray}{\mathcal{H}}_{\mathbf{d}}=C_{{\it\zeta}}\oplus L_{{\it\zeta}}\oplus R_{{\it\zeta}}\end{eqnarray}$$

parametrized by ${\it\zeta}\in (\mathbb{C}^{n+1})_{\ast }$. Note that the spaces $C_{{\it\zeta}}$, $L_{{\it\zeta}}$, and $R_{{\it\zeta}}$ depend only on $[{\it\zeta}]\in \mathbb{P}^{n}$.

Let ${\it\varphi}_{{\mathcal{H}}_{\mathbf{d}}}$ denote the density of the Gaussian $N(\overline{q},{\it\sigma}^{2}\text{Id})$ on ${\mathcal{H}}_{\mathbf{d}}$, where $\overline{q}\in {\mathcal{H}}_{\mathbf{d}}$ and ${\it\sigma}>0$. For fixed ${\it\zeta}\in \mathbb{S}(\mathbb{C}^{n+1})$ we decompose the mean $\overline{q}$ as

$$\begin{eqnarray}\overline{q}=\overline{k}_{{\it\zeta}}+\overline{g}_{{\it\zeta}}+\overline{h}_{{\it\zeta}}\in C_{{\it\zeta}}\oplus L_{{\it\zeta}}\oplus R_{{\it\zeta}}\end{eqnarray}$$

according to (3.24). If we denote by ${\it\varphi}_{C_{{\it\zeta}}}$, ${\it\varphi}_{L_{{\it\zeta}}}$, and ${\it\varphi}_{R_{{\it\zeta}}}$ the densities of the Gaussian distributions in the spaces $C_{{\it\zeta}}$, $L_{{\it\zeta}}$, and $R_{{\it\zeta}}$ with covariance matrices ${\it\sigma}^{2}\text{Id}$ and means $\overline{k}_{{\it\zeta}},\overline{M}_{{\it\zeta}}$, and $\overline{h}_{{\it\zeta}}$, respectively, then the density ${\it\varphi}_{{\mathcal{H}}_{\mathbf{d}}}$ factors as

(3.25)$$\begin{eqnarray}{\it\varphi}_{{\mathcal{H}}_{\mathbf{d}}}(k+g+h)={\it\varphi}_{C_{{\it\zeta}}}(k)\cdot {\it\varphi}_{L_{{\it\zeta}}}(g)\cdot {\it\varphi}_{R_{{\it\zeta}}}(h).\end{eqnarray}$$

Let $\mathscr{M}:=\mathbb{C}^{n\times (n+1)}$ and consider the manifold

$$\begin{eqnarray}{\hat{W}}:=\big\{(M,{\it\zeta})\in \mathscr{M}\times \mathbb{S}(\mathbb{C}^{n+1})\mid M{\it\zeta}=0\big\},\end{eqnarray}$$

whose projection $p_{2}:{\hat{W}}\rightarrow \mathbb{S}(\mathbb{C}^{n+1}),(M,{\it\zeta})\mapsto {\it\zeta}$, has the fibres

$$\begin{eqnarray}W_{{\it\zeta}}:=\{M\in \mathscr{M}\mid M{\it\zeta}=0\}.\end{eqnarray}$$

One can prove (see Bürgisser and Cucker Reference Bürgisser and Cucker2013, Proposition 16.16) that we have isometric linear maps

(3.26)$$\begin{eqnarray}W_{{\it\zeta}}\rightarrow L_{{\it\zeta}},\quad M\mapsto g_{M,{\it\zeta}}:=(\sqrt{d_{i}}\,\langle X,{\it\zeta}\rangle ^{d_{i}-1}\mathop{\sum }_{j}m_{ij}X_{j}).\end{eqnarray}$$

Thus the Gaussian density ${\it\varphi}_{L_{{\it\zeta}}}$ on $L_{{\it\zeta}}$ induces a Gaussian density ${\it\varphi}_{W_{{\it\zeta}}}$ on the fibre $W_{{\it\zeta}}$ with covariance matrix ${\it\sigma}^{2}\text{Id}$ via this map, so that ${\it\varphi}_{W_{{\it\zeta}}}(M)={\it\varphi}_{L_{{\it\zeta}}}(g_{M,{\it\zeta}})$.

Now consider the lifting $\hat{V}\subseteq {\mathcal{H}}_{\mathbf{d}}\times \mathbb{S}(\mathbb{C}^{n+1})$ of the solution variety $V\subseteq {\mathcal{H}}_{\mathbf{d}}\times \mathbb{P}^{n}$ induced by the natural projection $\mathbb{S}(\mathbb{C}^{n+1})\rightarrow \mathbb{P}^{n}$.

We can derive from the Gaussian distribution ${\it\varphi}_{{\mathcal{H}}_{\mathbf{d}}}$ on ${\mathcal{H}}_{\mathbf{d}}$ a probability distribution ${\it\rho}_{\hat{V}}$ on  $\hat{V}$ (naturally extending ${\it\rho}_{\mathsf{st}}$) as follows. One draws $(q,{\it\zeta})$ at random from $\hat{V}$ by first choosing $q\in {\mathcal{H}}_{\mathbf{d}}$ from $N(\overline{q},{\it\sigma}^{2}\text{Id})$, then choosing one of its ${\mathcal{D}}$ zeros $[{\it\zeta}]\in \mathbb{P}^{n}$ at random from the uniform distribution on $\{1,\ldots ,{\mathcal{D}}\}$, and finally choosing a representative ${\it\zeta}$ in the unit circle $[{\it\zeta}]\cap \mathbb{S}^{n}$ uniformly at random.

One can then prove that

(3.27)$$\begin{eqnarray}\mathbb{E}_{q\sim N(\overline{q},{\it\sigma}^{2}\text{Id})}\biggl(\frac{{\it\mu}_{\mathsf{av}}^{2}(q)}{\Vert q\Vert ^{2}}\biggr)=\mathbb{E}_{(q,{\it\zeta})\sim {\it\rho}_{\hat{V}}}\biggl(\frac{{\it\mu}_{\text{norm}}^{2}(q,{\it\zeta})}{\Vert q\Vert ^{2}}\biggr).\end{eqnarray}$$

Put ${\rm\Delta}:=\mathsf{diag}(\sqrt{d_{i}})$. Then the condition number ${\it\mu}_{\text{norm}}(q,{\it\zeta})$ can be described in terms of the linearization map

$$\begin{eqnarray}{\rm\Psi}:\hat{V}\rightarrow {\hat{W}},\quad (q,{\it\zeta})\mapsto (M,{\it\zeta}),\quad \text{where }M:={\rm\Delta}^{-1}Dq({\it\zeta}),\end{eqnarray}$$

as follows:

$$\begin{eqnarray}\frac{{\it\mu}_{\text{norm}}(q,{\it\zeta})}{\Vert q\Vert }=\Vert M^{\dagger }\Vert ,\quad \text{where }(M,{\it\zeta})={\rm\Psi}(q,{\it\zeta}).\end{eqnarray}$$

Hence

(3.28)$$\begin{eqnarray}\mathbb{E}_{(q,{\it\zeta})\sim {\it\rho}_{\hat{V}}}\biggl(\frac{{\it\mu}_{\text{norm}}^{2}(q,{\it\zeta})}{\Vert q\Vert ^{2}}\biggr)=\mathbb{E}_{M\sim {\it\rho}_{\mathscr{M}}}(\Vert M^{\dagger }\Vert ^{2}),\end{eqnarray}$$

where ${\it\rho}_{\mathscr{M}}$ denotes the pushforward density of ${\it\rho}_{\hat{V}}$ with respect to the map $p_{1}\circ {\rm\Psi}:\hat{V}\rightarrow \mathscr{M}$.

A general identity in probability theory tells us that

(3.29)$$\begin{eqnarray}\mathbb{E}_{M\sim {\it\rho}_{\mathscr{M}}}(\Vert M^{\dagger }\Vert ^{2})=\mathbb{E}_{{\it\zeta}\sim {\it\rho}_{\mathbb{S}(\mathbb{C}^{n+1})}}(\mathbb{E}_{M\sim \tilde{{\it\rho}}_{W_{{\it\zeta}}}}(\Vert M^{\dagger }\Vert ^{2})),\end{eqnarray}$$

where ${\it\rho}_{\mathbb{S}(\mathbb{C}^{n+1})}$ is the pushforward density of ${\it\rho}_{\hat{V}}$ with respect to $p_{2}\circ {\rm\Psi}:\hat{V}\rightarrow \mathbb{S}(\mathbb{C}^{n+1})$ and $\tilde{{\it\rho}}_{W_{{\it\zeta}}}$ denotes the conditional density on the fibre $W_{{\it\zeta}}$. This conditional density turns out to be (and this is intricate) of the form

(3.30)$$\begin{eqnarray}\tilde{{\it\rho}}_{W_{{\it\zeta}}}(M)=c_{{\it\zeta}}^{-1}\cdot \det (MM^{\ast })\,{\it\varphi}_{W_{{\it\zeta}}}(M),\end{eqnarray}$$

with $c_{{\it\zeta}}$ denoting a normalization factor. At this stage one can use invariance under the unitary group acting on $\mathbb{C}^{n+1}$ to reduce the study of the inner expectation in (3.29) to the case ${\it\zeta}=(1,0,\ldots ,0)$. In this case we can identify $W_{{\it\zeta}}$ with $\mathbb{C}^{n\times n}$, and $\tilde{{\it\rho}}_{W_{{\it\zeta}}}$ takes the form

$$\begin{eqnarray}\tilde{{\it\rho}}(A)=c^{-1}|\det A|^{2}{\it\varphi}(A),\end{eqnarray}$$

where $c:=\mathbb{E}_{A\sim N(0,\text{Id})}|\det A|^{2}$ and ${\it\varphi}$ is the standard Gaussian density in $\mathbb{C}^{n\times n}$. The density $\tilde{{\it\rho}}$ is thus closely related to a standard Gaussian, and one can prove that

(3.31)$$\begin{eqnarray}\mathbb{E}_{M\sim \tilde{{\it\rho}}_{W_{{\it\zeta}}}}(\Vert M^{\dagger }\Vert ^{2})=\mathbb{E}_{A\sim \tilde{{\it\rho}}}(\Vert A^{-1}\Vert ^{2})\leq \frac{e(n+1)}{2{\it\sigma}^{2}}.\end{eqnarray}$$

This implies by (3.29) that

$$\begin{eqnarray}\mathbb{E}_{M\sim {\it\rho}_{\mathscr{M}}}(\Vert M^{\dagger }\Vert ^{2})\leq \frac{e(n+1)}{2{\it\sigma}^{2}},\end{eqnarray}$$

and this bound, substituted in (3.28) and hence in (3.27), completes this sketch of the proof.◻