## Thursday, April 22, 2021

### Growth of the Gaussian pivoting algorithm

Gaussian elimination” is the standard method for solving systems of linear equations that runs by choosing one pivot variable in one of the equations and eliminating it from the other equations by a suitable linear combination. These other equations have one less variable and one may iterate the process, leading to a system that has a triangular form and can be solved in a relatively straightforward way. In the so-called Gauss-Jordan method, the pivot variables are eliminated from all of the equations, and this leads to the row reduced echelon form of the initial system, a new system in which the pivot variables are explicitly solved in terms of the remaining variables; it also has the merit of being independent of the intermediate steps, but this would be a story for another night.

How the name of Gauss is attributed to this method is also another story, that is recounted in great detail by J. Grcar. As he explains, systems of linear equations and their solution already appear on Babylonian tablets (2000 BC), on the Egyptian Rhind papyrus (1550 BC), in the Chinese Nine chapters on the mathematical art (200 AD), the Arithmetica of the Greek mathematician Diophantus (250 AD), the Āryabhaṭīya of the Hindu mathematician Āryabhaṭa (499 AD). Such systems were essentially absent of the mathematics treatises during the Renaissance and it is to Newton (17th century) that we owe its revival in Western Europe. At the beginning of the 19th century, in relation with the problem in celestial mechanics/geodesy of fitting together multiple imprecise measurements, Gauss and Legendre invented the least square methods. This involved a system of linear equations which Gauss solved by what he called “common elimination”. In the 20th century, the development of computing machines and numerical analysis led to further work, from Cholesky (in relation with geodesy), Crout, to von Neumann and Goldstine and their $LDU$ decomposition.

Whoever had to perform elimination by hand knows that the computations are rapidly tedious and often lead to more and more complicated fractions.

When computer calculations are done with floating point algebra, the difficulty of rounding errors appears. If in a linear system, say $Ax=b$,  the matrices $A$ and $b$ are only known up to an error, so that the system that is actually solved would rather be $(A+E)x=b+\delta b$, and it is of an obvious importance to compare its solution with the solution of the initial system. One way to make this comparison involves the inverse of $A$, which is unknown at this stage. The product of the norms $\|A\| \,\|A^{-1}\|$ is the conditioning number of $A$, and one can not avoid the problem of ill-conditioned matrices, which will inherently lead to lack of precision.

But when floating points numbers are used in the computation, a new problem appears, even when one restricts to well-conditioned systems. Floating points numbers are of the form $\pm a\cdot 10^e$ (in base 10, say), where $a$ is a real number between $1$ and $10$ which is known up to a fixed number of decimal places. In other words, floating points numbers are known up to a relative error. Let us examine the consequence for the subtraction of floating point numbers. Consider two such numbers, say $x=a\cdot 10^e$ and $x'=a'\cdot 10^e$, with the same exponent $e$, and such that $a$ and $a'$ are close. Their difference is given by $x'-x=(a'-a)\cdot 10^e$, but since $a$ and $a'$ are close, their difference $y=x'-x$ is no more between $1$ and $10$, but may be, say of size $10^{-5}$, so that the floating point expression of $x'-x$ is $(10^5(a'-a))\cdot 10^{e-5}=b\cdot 10^{e-5}$; the problem is that the last 5 decimals of $b$ are absolutely unknown, which leads to a relative error for $y$ which is $10^5$ as big as expected!

Now, by its very essence, the elimination method features a lot of such subtractions, hence it is inherently not very well suited with floating points numbers.

In the 1960s, the American mathematician Wilkinson analysed the situation very precisely.  He showed that for the Gaussian elimination, the main parameter is the relative size of the matrices that intervene in the process. To set up some notation, imagine that the initial system $Ax=b$ is transformed, step by step, into a series of equivalent systems $A^{(r)}x=b^{(r)}$, where $A^{(r)}=(a_{ij}^{(r)})$ is a matrix whose first $r$ lines are in triangular form, and the remaining $n-r$ lines still need to be worked on. To get the matrix $A^{(r+1)}$, one subtract a multiple $m_{ir}$ of the $r$th row from the $i$th row he multipliers, for $i$ ranging from $r+1$ to $n$, where the multipliers $m_{ir}$ are defined by
$m_{ir}=a_{ir}^{(r)}/a_{rr}^{(r)}.$
Large multipliers lead to lack of precision, but if complete pivoting method is used, one has $|m_{ir}|\leq 1$. In this case, one observes that a bound $|a_{ij}^{(r)}|\leq M$ for the coefficients of $A^{(r)}$ leads to the bound $|a_{ij}^{(r+1)}|\leq 2M$ at the next step. At the last step, one gets $|a_{ij}^{(n)}|\leq 2^{n-1}M$, where $M=\sup(|a_{ij}|)$. Consequently, the relevant constant to be estimated,
$R(A) = \sup_r \frac{\sup_{i,j}|a_{ij}^{(r)}|}{\sup_{i,j}|a_{ij}|},$
satisfies $R(A)\leq 2^{n-1}$.

In fact, Wilkinson (1961) gave a much better bound. Let $B^{(r)}$ be the square matrix of size $n-r$ that has to be worked on after the $r$th step and let $b_{r}$ be the maximum size of its coefficients, the size of the chosen pivot since one does complete pivoting. One has the following fomula for its determinant:
$\det(B^{(r)})= b_{r+1}\cdots b_n.$
Moreover, the Euclidean norms of the the columns of $B^{(r)}$ are bounded above by $\sqrt{n-r} b_{r+1}$ and Hadamard inequality (“the volume of a parallelepiped is smaller than the product of the sizes of its edges”) implies that
$| \det(B^{(r)})| \leq (n-r)^{(n-r)/2} b_{r+1}^{n-r}.$
Together, these relations lead to an upper bound
$R(A) \leq n^{1/2} \left( 2\cdot 3^{1/2}\cdot 4^{1/3}\cdots n^{1/(n-1)}\right)^{1/2},$
roughly $n^{\log(n)/4}$, a bound that Wilkinson considers to be a “severe overestimate”, and in Wilkinson (Rounding Errors in Algebraic Processes, Dover, 1963, bottom of p. 97), he even notes that “No matrix has yet been discovered for which $R(A)>n$.

This statement remained known as Wilkinson's conjecture, although Wilkinson himself did not state it as such. Tornheim (1964) proved that Hadamard matrices — matrices with entries $\pm1$ and pairwise orthogonal columns and rows — satisfy $R(A)\geq n$ and this led Cryer (1968) to formally state the conjecture and to suspect that $R(A)<n$ unless $A$ is a Hadamard matrix. Some matrices have been shown such that $R(A)\geq (n+1)/2$ (Higham & Higham, 1989) and random matrices seems to have an $R$-factor roughly $n^{1/2}$ (Trefethen & Schreiber, 1990). In fact, Hadamard matrices of size $n=12$ (Edelman & Mascarenhas, 1995) or $n=16$ (Kravvaritis & Mtrouli, 2009) satisfy $R(A)=n$, but this is definitely nontrivial.

However, Gould (1991) could exhibit a matrix $A$ of size $n=13$ with  growth factor $R(A)=13.0205$, thus providing a counterexample to Wilkinson's “conjecture”. To that aim, he reduced the question to finding a solution of a nonlinear programming problem with roughly $n^3/3\approx 700$ variables, fortunately a sparse one, a computation he did using a programming package he had developed with Conn and Toint. The matrix he gives has integer coefficients of size up to 20 digits!

But the story does not end here!

Edelman tried to replicate Gould's computations using the computer algebra softwares Mathematica and Maple — what he found is the growth factor of Gould's matrix is around $7.355$, consistently in both softwares! As he writes, one of these softwares could have had a bug, but it is rather unlikely that both of them had had the same one.

What happened, and you will agree that there is much irony in this story, is that Gould had performed his computations using floating point algebra. A near tie in the 6th pivot lead to an incorrect choice of pivot and to an erroneous computation, even within the double precision that he has used.

Fortunately, Edelman (1992) showed that changing one coefficient by $10^{-7}$ in Gould's matrix yields a growth of $13.02$, so that Wilkinson's “conjecture” is definitely incorrect.

## Friday, April 2, 2021

### On the Hadamard-Lévy theorem, or is it Banach-Mazur?

During the preparation of an agrégation lecture on connectedness, I came across the following theorem, attributed to Hadamard–Lévy:

Theorem. — Let $f\colon \mathbf R^n\to\mathbf R^n$ be a $\mathscr C^1$-map which is proper and a local diffeomorphism. Then $f$ is a global diffeomorphism.

In this context, that $f$ is proper means that $\| f(x)\| \to+\infty$ when $\| x\|\to+\infty$, while, by the inverse function theorem, the condition that $f$ is a local diffeomorphism is equivalent to the property that its differential $f'(x)$ is invertible, for every $x\in\mathbf R^n$. The conclusion is that $f$ is a diffeomorphism from $\mathbf R^n$ to itself; in particular, $f$ is bijective and its inverse is continuous.

This theorem is not stated in this form neither by Hadamard (1906), nor by Lévy (1920), but is essentially due to Banach & Mazur (1934) and it is the purpose of this note to clarify the history, explain a few proofs, as well as more recent consequences for partial differential equations.

A proper map is closed: the image $f(A)$ of a closed subset $A$ of $\mathbf R^n$ is closed in $\mathbf R^n$. Indeed, let $(a_m)$ be a sequence in $A$ whose image $(f(a_m))$ converges in $\mathbf R^n$ to an element $b$; let us show that there exists $a\in A$ such that $b=f(a)$. The properness assumption on $f$ implies that $(a_m)$ is bounded. Consequently, it has a limit point $a$, and $a\in A$ because $A$ is closed. Necessarily, $f(a)$ is a limit point of the sequence $(f(a_m))$, hence $b=f(a)$.

In this respect, let us note the following reinforcement of the previous theorem, due to Browder (1954):
Theorem (Browder). — Let $f\colon \mathbf R^n\to\mathbf R^n$ be a local homeomorphism. If $f$ is closed, then $f$ is a global homeomorphism.

A surprising aspect of these results and their descendents is that they are based on two really different ideas. Banach & Mazur and Browder are based on the notion of covering, with ideas of homotopy theory and, ultimately, the fact that $\mathbf R^n$ is simply connected. On the other hand, the motivation of Hadamard was to generalize to dimension $n$ the following elementary discussion in the one-dimensional case: Let $f\colon\mathbf R\to\mathbf R$ be a $\mathscr C^1$-function whose derivative is $>0$ everywhere (so that $f$ is strictly increasing); give a condition for $f$ to be surjective. In this case, the condition is easy to find: the indefinite integral $\int f'(x)\,dx$ has to be divergent both at $-\infty$ and $+\infty$. In the $n$-dimensional case, the theorems of Hadamard is the following:

Theorem.Let $f\colon\mathbf R^n\to\mathbf R^n$ be a $\mathscr C^1$-map. For $r\in\mathbf R_+$, let $\omega(r)$ be the infimum, for $x\in\mathbf R^n$ such that $\|x\|=r$, of the norm of the linear map $f'(x)^{-1}$; if $\int_0^\infty dr/\omega(r)=+\infty$, then $f$ is a global diffeomorphism.

In Hadamard's paper, the quantity $\omega(r)$ is described geometrically as the minor axis of the ellipsoid defined by $f'(x)$, and Hadamard insists that using the volume of this ellipsoid only, essentially given by the determinant of $f'(x)$, would not suffice to characterize global diffeomorphisms. (Examples are furnished by maps of the form $f(x_1,x_2)=(f_1(x_1),f_2(x_2))$. The determinant condition considers $f_1'(x_1)f_2'(x_2)$, while one needs individual conditions on $f'_1(x_1)$ and $f'_2(x_2)$.)

In fact, as explained in Plastock (1974), both versions (closedness hypothesis or quantitative assumptions on the differential) imply that the map $f$ is a topological covering of $\mathbf R^n$. Since the target $\mathbf R^n$ is simply connected and the source $\mathbf R^n$ is connceted, $f$ has to be a homeomorphism. I will explain this proof below, but I would first like to explain another one, due to Zuily & Queffelec (1995) propose an alternating proof which is quite interesting.

## A dynamical system approach

The goal is to prove that $f$ is bijective and, to that aim, we will prove that every preimage set $f^{-1}(b)$ is reduced to one element. Replacing $f$ by $f-b$, it suffices to treat the case of $b=0$. In other words, we wish to solve that the equation $f(x)=0$ has exactly one solution. For that, it is natural to try to start from some point $\xi\in\mathbf R^n$ and to force $f$ to decrease. This can be done by following the flow of the vector field given by $v(x)=-f'(x)^{-1}(f(x))$. This is a vector field on $\mathbf R^n$ and we can consider its flow: a map $\Phi$ defined on an open subset of $\mathbf R\times\mathbf R^n$ such that $\partial_t \Phi(t,x)=v(\Phi(t,x))$ for all $(t,x)$ and $\Phi(0,x)=x$ for all $x$. In fact, the Cauchy–Lipschitz theorem guarantees the existence of such a flow only if the vector field $v$ is locally Lipschitz, which happens if, for example, $f$ is assumed to be $\mathscr C^2$. In this case, there is even uniqueness of a maximal flow, and we will make this assumption, for safety. (In fact, the paper of De Marco, Gorni & Zampieri (1994) constructs the flow directly thanks to the hypothesis that the vector field is pulled back from the Euler vector field on $\mathbf R^n$.)

What are we doing here? Note that in $\mathbf R^n$, the opposite of the Euler vector field, defined by $u(y)=-y$, has a very simple solution: the flow lines are straight lines going to $0$. The formula above just pulls back this vector field $u$ via the local diffeomorphism $f$, and the flow lines of the vector field $v$ will just be the ones given by pull back by $f$, which will explain the behaviour described below.

In particular, let $a\in\mathbf R^n$ be such that $f(a)=0$ and let $U$ be a neighborhood of $a$ such that $f$ induces a diffeomorphism from $U$ to a ball around $0$. Pulling back the solution of the minus-Euler vector field by $f$, we see that once a flow line enters the open set $U$, it converges to $a$. The goal is now to prove that it will indeed enter such a neighborhood (and, in particular, that such a point $a$ exists).

We consider a flow line starting from a point $x$, that is, $\phi(t)=\Phi(t,x)$ for all times $t$. Let $g(t)= f(\phi(t))$; observe that $g$ satisfies $g'(t)=f'(\phi(t))(\phi'(t))=-g(t)$, hence $g(t)=g(0)e^{-t}$. Assume that the line flow is defined on $[0;t_1\mathopen[$, with $t_1<+\infty$. by what precedes, $g$ is bounded in the neighborhood of $t_1$; since $f$ is assumed to be proper, this implies that $\phi(t)$ is bounded as well. The continuity of the vector field $v$ implies that $\phi$ is uniformly continuous, hence it has a limit at $t_1$. We may then extend the line flow a bit right of $t_1$. As a consequence, the line flow is defined for all times, and $g(t)\to0$ when $t\to+\infty$. By the same properness argument, this implies that $\phi(t)$ is bounded when $t\to+\infty$, hence it has limit points $a$ which satisfy $f(a)=0$. Once $\phi$ enters an appropriate neighborhood of such a point, we have seen that the line flow automatically converges to some point $a\in f^{-1}(0)$.

Let us now consider the map $\lambda\colon\mathbf R^n\to f^{-1}(0)$ that associates with a point $\xi$ the limit of the line flow $t\mapsto \Phi(t,\xi)$ starting from the initial condition $\xi$. By continuity of the flow of a vector field depending on the initial condition, the map $\lambda$ is continuous. On the other hand, the hypothesis that $f$ is a local diffeomorphism implies that $f^{-1}(0)$ is a closed discrete subset of $\mathbf R^n$. Since $\mathbf R^n$ is connected, the map $\lambda$ is constant. Since one has $\lambda(\xi)=\xi$ for every $\xi\in f^{-1}(0)$, this establishes that $f^{-1}(0)$ is reduced to one element, as claimed.

Once $f$ is shown to be bijective, the fact that it is proper (closed would suffice) implies that its inverse bijection $f^{-1}$ is continuous. This concludes the proof.

## The theorem of Banach and Mazur

The paper of Banach and Mazur is written in a bigger generality. They consider multivalued continuous maps $F\colon X\to Y$ ($k$-deutige stetige Abbildungen) by which they mean that for every $x$, a subset $F(x)$ of $Y$ is given, of cardinality $k$, the continuity being expressed by sequences: if $x_n\to x$, one can order, for every $n$, the elements of $F(x_n)=\{y_{n,1},\dots,y_{n,k}\}$, as well as the elements of $F(x)=\{y_1,\dots,y_k\}$, in such a way that $y_{n,j}\to y_n$ for all $j$. (In their framework, $X$ and $Y$ are metric spaces, but one could transpose their definition to topological spaces if needed.) They say that such a map is decomposed (zerfällt) if there are continuous functions $f_1,\dots,f_k$ from $X$ to $Y$ such that $F(x)=\{f_1(x),\dots,f_k(x)\}$ for all $x\in X$.

In essence, the definition that Banach and Mazur are proposing contains as a particular case the finite coverings. Namely, if $p\colon Y\to X$ is a finite covering of degree $k$, then the map $x\mapsto p^{-1}(x)$ is a continuous $k$-valued map from $X$ to $Y$. Conversely, let us consider the graph $Z$ of $F$, namely the set of all points $(x,y)\in X\times Y$ such that $y\in F(x)$. Then the first projection $p\colon Z\to X$ is a covering map of degree $k$, but it is not clear that it has local sections.

It would however not be so surprising to 21st-century mathematicians that if one makes a suitable assumption of simple connectedness on $X$, then every such $F$ should be decomposed. Banach and Mazur assume that $X$ satisfies two properties:

1. The space $X$ is semilocally arcwise connected: for every point $x\in X$ and every neighborhood $U$ of $x$, there exists an open neighborhood $U'$ contained in $U$ such that for every point $x'\in U'$, there exists a path $c\colon[0;1]\to U$ such that $c(0)=x$ and $c(1)=x'$. (Semilocally means that the path is not necessarily in $U'$ but in $U$.)
2. The space $X$ is arcwise simply connected: two paths $c_0,c_1\colon[0;1]\to X$ with the same endpoints ($c_0(0)=c_1(0)$ and $c_0(1)=c_1(1)$) are strictly homotopic — there exists a continuous map $h\colon[0;1]\to X$ such that $h(0,t)=c_0(t)$ and $h(1,t)=c_1(t)$ for all $t$, and $h(s,0)=c_0(0)$ and $h(s,1)=c_0(1)$ for all $s$.

Consider a $k$-valued continuous map $F$ from $X$ to $Y$, where $X$ is connected. Banach and Mazur first prove that for every path $c\colon [0;1]\to X$ and every point $y_0\in F(c(0))$, there exists a continuous function $f\colon[0;1]\to Y$ such that $f(t)\in F(c(t))$ for all $t$. To that aim, the consider disjoint neighborhoods $V_1,\dots,V_k$ of the elements of $F(c(0))$, with $y_0\in V_1$, say, and observe that for $t$ small enough, there is a unique element in $F(c(t))\cap V_1$. This defines a bit of the path $c$, and one can go on. Now, given two paths $c,c'$ such that $c(0)=c'(0)$ and $c(1)=c'(1)$, and two maps $f,f'$ as above, they consider a homotopy $h\colon[0;1]\times[0;1]\to X$ linking $c$ to $c'$. Subdividing this square in small enough subsquares, one see by induction that $f(1)=f'(1)$. (This is analogous to the proof that a topological covering of the square is trivial.) Fixing a point $x_0\in X$ and a point $y_0\in F(x_0)$, one gets in this way a map from $X$ to $Y$ such that $F(x)$ is equal to $f(1)$, for every path $c\colon[0;1]\to X$ such that $c(0)=x_0$ and $c(1)=x$, and every continuous map $f\colon [0;1]\to Y$ such that $f(t)\in F(c(t))$ for all $t$ and $f(0)=y_0$. This furnishes a map from $X$ to $Y$, and one proves that it is continuous. If one considers all such maps, for all points in $F(x_0)$, one obtains the decomposition of the multivalued map $F$.

To prove their version of the Hadamard–Lévy theorem, Banach and Mazur observe that if $f\colon Y\to X$ is a local homeomorphism which is proper, then setting $F(x)=f^{-1}(y)$ gives a multivalued continuous map. It is not obvious that the cardinalities $k(x)$ of the sets $F(x)$ are constant, but this follows (if $X$ is connected) from the fact that $f$ is both a local homeomorphism and proper. Then $F$ is decomposed, so that there exist continuous maps $g_1,\dots,g_k\colon X\to Y$ such that $f^{-1}(x)=\{g_1(x),\dots,g_k(x)\}$ for all $x\in X$. This implies that $Y$ is the disjoint union of the $k$ connected subsets $g_j(X)$. If $Y$ is connected, then $f$ is a homeomorphism.

## The versions of Hadamard and Lévy, after Plastock

Hadamard considered the finite dimensional case, and Lévy extended it to the case of Hilbert spaces.

Plastock considers a Banach-space version of the theorem above: $f\colon E\to F$ is a $\mathscr C^1$-map between Banach spaces with invertible differentials and such that, setting $\omega(r)=\inf_{\|x\| = r}\|f'(x)^{-1}\|$, one has $\int_0^\infty \omega(r)\,dr=+\infty$. Of course, under these hypotheses, the Banach spaces $E$ and $F$ are isomorphic, but it may be useful that they are not identical. Note that $f(E)$ is open in $F$, and the proposition that will insure that $f$ is a global diffeomorphism is the following one, in the spirit of covering theory.

Proposition.(Assuming that $f$ is a local diffeomorphism.) It suffices to prove that the map $f$ satisfies the path lifting property: for every point $x\in E$ and every $\mathscr C^1$ map $c\colon[0;1]\to f(E)$ such that $c(0)=f(x)$, there exists a $\mathscr C^1$ map $d\colon[0;1]\to E$ such that $c(t)=f(d(t))$ for all $t$ and $d(0)=c$.

The goal is now to prove that $f$ satisfies this path lifting property. Using that $f$ is a local homeomorphism, one sees that lifts are unique, and are defined on a maximal subinterval of $[0;1]$ which is either $[0;1]$ itself, or of the form $[0;s\mathclose[$. To prevent the latter case, one needs to impose conditions on the norm $\| f'(x)^{-1}\|$ such as the one phrased in terms of $\omega(r)$ as in the Hadamard–Lévy theorem. In fact, Plastock starts with a simpler case.

Proposition.The path lifting property follows from the following additional hypotheses:

1. One has $\|f(x)\|\to+\infty$ when $\|x\|\to+\infty$;