# Numerical conditioning

Charbel Tannous

► To cite this version:

Charbel Tannous. Numerical conditioning. Master. Analyse Numérique, UBO Brest, France. 2016, pp.23. hal-04943018

HAL Id: hal-04943018

https://hal.science/hal-04943018v1

Submitted on 12 Feb 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
# Numerical Conditioning

C. Tannous

<figure><img src="image_2.png" /></figure>

Université de Brest, Lab-STICC, CNRS-UMR 6285, F-29200 Brest, FRANCE

Conditioning is an important notion in Numerical Analysis and generally poorly described. In this work, many examples are given in order to illustrate this notion in polynomials, functions, non-linear equations, linear algebra, differential equations and solution of linear systems.

PACS numbers: 02.10.Yn, 02.10.Ud, 02.60.Cb

Keywords: Matrix theory, Linear Algebra, Numerical simulation, solution of equations

Version: February 11, 2025

## Contents

### I. Introduction

2

### II. Conditioning

3

A. Arithmetic operations conditioning

3

B. Recursion conditioning

3

C. Polynomial conditioning

4

D. Function conditioning

6

E. Non-linear equation conditioning

6

F. Differential equation conditioning

6

### III. Linear systems and Matrix conditioning

7

A. Impact of rounding on rank

7

B. Norms of Vectors and Matrices

8

C. Condition numbers

9

### IV. Singular value decomposition

9

A. Hilbert matrices

11

B. Wilkinson matrix

13

C. Pascal matrices

13

D. Solution of over-determined linear systems

14

### V. Iterative methods for solving linear systems

15

A. Jacobi Iterative Method

16

B. Gauss-Seidel Iterative Method

16

C. Relaxation Methods

17

D. Steepest Descent

17

E. Conjugate Gradient

18

### VI. Preconditioning techniques for iterative solvers of linear systems

19

A. Diagonal preconditioning

19

B. Binormalization

20

C. Singular linear systems preconditioning

21

### A. Extremum Problems in Higher dimensions

21

1. Powell Method

21

2. Broyden method

22

### References

23
# I. INTRODUCTION

Every problem in numerical analysis is described by a condition figure since the number representation in a computer is finite [1–3].

A numerical problem can be described as possessing an input $x$ and output $f(x)$, with their respective representation error $\delta x$ and $\delta f(x)$ (cf. Fig. 1).

<figure><img src="image_3.png" /><figcaption>**Fig.1:** Conditioning viewed as a non-linear filtering operation. Relative sensitivity is $|\delta f(x)/f(x)|$ whereas relative uncertainty is $|\delta x/x|$.</figcaption></figure>

The condition number measures how errors in the input (relative sensitivity) are reflected into output errors (relative uncertainty). Thus the condition number $\kappa$ is given by $|\delta f(x)/f(x)| = \kappa(|\delta x/x|)$, i.e. $\kappa = (|\delta f(x)/f(x)|)/(|\delta x/x|)$.

Thus a problem is:

* Well-conditioned if $\kappa$ is small,

* Ill-conditioned if $\kappa$ is large,

* Ill-posed when $\kappa = \infty$

Generally $\kappa = \infty$ in the case of polynomials or functions with multiple roots.

For a function, we have in the the discrete case: $\text{cond}=\kappa = |\frac{\delta f(x)}{f(x)} / (\frac{\delta x}{x})|$.

The condition number $\kappa$ may be also viewed as an amplifying factor of the input error.

**Note:** In the analytical continuous case: $\text{cond}=\kappa = |\frac{f(x+\delta x)-f(x)}{f(x)} / (\frac{\delta x}{x})| \approx |\frac{xf'(x)}{f(x)}|$

In the case of a function of several variables $f(x_1, x_2, \dots, x_n)$, there are partial condition coefficients $a_i$ given by:

$$ \frac{\delta f(x)}{f(x)} = \sum_{i=1}^{n} \left( \left( \frac{\partial f}{\partial x_i} \right) \right) \frac{\delta x_i}{f(x)} = \sum_{i=1}^{n} a_i \frac{\delta x_i}{x_i} \quad \text{with} \quad a_i = \frac{\delta x_i}{f(x)} \left( \frac{\partial f}{\partial x_i} \right) \qquad (1) $$

Partial condition numbers are $\text{cond}_i = |a_i|$, $i = 1...n$ and the function $f(x_1, x_2, \dots, x_n)$ is:
Well-conditioned when $\forall i \text{ cond}_i \le 1$,
Ill-conditioned when $\exists i$ such that $\text{cond}_i > 1$.

In linear systems ($Ax = b$) the matrix case, there are several definitions for the condition number depending on the matrix norm chosen for $A$.

Illustrating conditioning of linear systems, such as:

$$ Ax = b \qquad (2) $$

can be made with the following example:

$$ Ax = \begin{pmatrix} 1 & 1 \\ a & 1 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \end{pmatrix} = b = \begin{pmatrix} 1 \\ 2 \end{pmatrix} $$

The solution of this system is $x_1 = 10, x_2 = -9$ when $a = 1.1$.

Let us suppose the factor $a$ is in fact $a = 1.05$ then the solution becomes $x_1 = 20, x_2 = -19$. This means that a 5% change in a single coefficient of matrix $A$ induced 100% change in the solution $x = \begin{pmatrix} x_1 \\ x_2 \end{pmatrix}$ [4]. The error amplification $\kappa \approx 20$ might be explained further below.

**Note:** Alternatively, the number of significant digits in a numerical problem is off by $n$ digits when $\kappa \propto 10^n$.
## II. CONDITIONING

### A. Arithmetic operations conditioning

Arithmetic operations can be represented by functions of two variables:

* Addition is represented by $f(x_1, x_2) = x_1 + x_2$

* Subtraction is represented by $f(x_1, x_2) = x_1 - x_2$

* Multiplication is represented by $f(x_1, x_2) = x_1 \cdot x_2$

* Division is represented by $f(x_1, x_2) = x_1/x_2$

Applying the rules pertaining to functions of many variables, we get:

* Addition:

$$ \frac{\delta f(x)}{f(x)} = \frac{x_1}{(x_1 + x_2)} \frac{\delta x_1}{x_1} + \frac{x_2}{(x_1 + x_2)} \frac{\delta x_2}{x_2} \qquad (3) $$

Considering $x_1, x_2 > 0$ the partial conditions are $\le 1$, thus addition is well-conditioned.

* Subtraction:

$$ \frac{\delta f(x)}{f(x)} = \frac{x_1}{(x_1 - x_2)} \frac{\delta x_1}{x_1} + \frac{x_2}{(x_1 - x_2)} \frac{\delta x_2}{x_2} \qquad (4) $$

Considering $x_1, x_2 > 0$ the partial conditions are $\le 1$, thus subtraction is well-conditioned when $x_1 \ne x_2$ and ill-conditioned when $x_1 \approx x_2$.

* Multiplication:

$$ \frac{\delta f(x)}{f(x)} = \frac{\delta x_1}{x_1} + \frac{\delta x_2}{x_2} \qquad (5) $$

Considering $x_1, x_2 > 0$ the partial conditions are 1, thus multiplication is well-conditioned.

* Division:

$$ \frac{\delta f(x)}{f(x)} = \frac{\delta x_1}{x_1} - \frac{\delta x_2}{x_2} \qquad (6) $$

Considering $x_1 > 0, x_2 \ne 0$ the partial conditions are 1, thus multiplication is well-conditioned.

### B. Recursion conditioning

Consider a recursion formula originating from the discretization of an Ordinary Differential Equation (ODE), a Partial Differential Equation (PDE) or some Finite Element Analysis:

$$ x_{i+2} = -(13/6)x_{i+1} + (5/2)x_i \qquad (7) $$

with initial conditions: $x_1 = 30, x_2 = 25$.

It can be solved with the substitution $x_i = r^i$ to yield its associated characteristic equation $r^2 = -(13/6)r + (5/2)$ whose roots are $r_1 = 5/6, r_2 = -3$.

The general solution is given by:

$$ x_i = C_1 r_1^i + C_2 r_2^i \qquad (8) $$

Accounting for the initial conditions $x_1 = 30, x_2 = 25$, the constants are evaluated $C_1 = 36, C_2 = 0$ and the solution obtained is $x_i = 36(5/6)^i$.
<figure><img src="image_4.png" /><figcaption>**Fig.2:** Values $|x(i)|$ given by recursion $x_{i+2} = -(13/6)x_{i+1} + (5/2)x_i$ in simple precision and double precision compared to the exact solution $x_i = 36(5/6)^i$. While they agree for $i < 15$, they diverge with even sign change above $i=15$, while the exact solution keep decreasing uniformly toward zero.</figcaption></figure>

Evaluation of the recursion in simple precision and double precision compared to the exact solution $x_i = 36(5/6)^i$ are displayed in Fig. 2.

Conditioning is due to a rounding error $\delta \neq 0$ impacting the solution given numerically by:

$$x_i = C_1 r_1^i \pm \delta r_2^i = 36(5/6)^i \pm \delta(3)^i \quad (9)$$

In float representation $\delta = 10^{-7}$ whereas in double it is $\delta = 10^{-16}$.

Validating the rounding effect is displayed in Fig. 3 where the recursion is compared to the rounded float representation $x_i = 36(5/6)^i \pm \delta(3)^i$ with the $\pm$ sign chosen randomly.

Another example is the Muller [5, 6] recursion:

$$x_{i+1} = 111 - (1130 - 3000/x_{i-1})/x_i \quad (10)$$

with initial conditions: $x_1 = 11/2, x_2 = 61/11$.

In exact arithmetic, this sequence increases monotonically toward 6 which is found by solving the fixed point relationship $r = \lim_{i \to \infty} x_i$:

$$r^3 - 111r^2 + 1130r - 3000 = 0, \quad \text{with solutions: } r_1 = 5, r_2 = 6, r_3 = 100 \quad (11)$$

Given the initial conditions, the sequence is forced to proceed toward the second fixed point $r_2 = 6$. In exact arithmetic, this is observed indeed as displayed in Fig. 4. Float computer arithmetic with a truncated representation induces the sequence to switch, after a large jump, to the third fixed point $r_3 = 100$ as displayed in Fig. 4. In the double precision case, the sequence proceeds smoothly to the third fixed point $r_3 = 100$.

Conditioning of the Muller recursion is proven by the extreme sensitivity to its initial conditions. Using exact arithmetic and starting it with $x_1 = 11/2, x_2 = 62/11$, it progresses to the third fixed point $r_3 = 100$ whereas it progresses to the second fixed point $r_2 = 6$ when $x_1 = 11/2, x_2 = 61/11$ as displayed in Fig. 4. The condition number is about $10^3$ since $x_2$ change from 61/11 to 62/11 induces fixed point drift from 6 to 100.

### C. Polynomial conditioning

Polynomial conditioning is best illustrated with the Wilkinson polynomial [3].
<figure><img src="image_5.png" /><figcaption>Fig.3: Values $|x(i)|$ given by recursion $x_{i+2} = -(13/6)x_{i+1} + (5/2)x_i$ in simple precision and compared to the rounded float representation $x_i = 36(5/6)^i \pm \delta(3)^i$ with $\pm$ value chosen randomly.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig.4: Exact and numerical (float and double) values $x(i)$ given by the Muller recursion $x_{i+1} = 111 - (1130 - 3000/x_{i-1})/x_i$ with initial conditions: $x_1 = 11/2, x_2 = 61/11$. (Left) In exact arithmetic, the sequence proceeds toward the second fixed point $r_2 = 6$. (Right) In float representation the sequence displays a double glitch before switching, after a large jump, to the third fixed point $r_3 = 100$ whereas in the double precision case, the sequence proceeds smoothly to the third fixed point $r_3 = 100$. Using exact arithmetic and starting the sequence with $x_1 = 11/2, x_2 = 62/11$, it proceeds toward $r_3 = 100$ (curve labeled 62/11) like the float and double cases.</figcaption></figure>

A Wilkinson polynomial is made from factors containing consecutive integers 1, 2, ...

$$p_n(x) = (x-1)(x-2)\cdots(x-n) \qquad (12)$$

For example, the 20-th order Wilkinson [3] polynomial is

$$p_{20}(x) = (x-1)(x-2)\cdots(x-20) \qquad (13)$$

The constant coefficient of $p(x)$, after expanding it, is 20!. This number can be evaluated with Stirling approxima-
tion [7]:

$$ \Gamma(n) = (n-1)! \sim \left(\frac{2\pi}{n}\right)^{1/2} \left(\frac{n}{e}\right)^n \left(1 + \frac{A_1}{n} + \frac{A_2}{n^2} + \frac{A_3}{n^3} + \dots \right), \quad n \to \infty \quad (14) $$

Thus $20! \approx 2.432 \times 10^{18}$.

### D. Function conditioning

Function conditioning is illustrated with the following example around $x = 1$:

* $f(x) = \ln(x)$ for which the condition number is given by: $|\frac{x f'(x)}{f(x)}| = |1/\ln(x)|$ implying that when $x \approx 1$, error of evaluation of $f(x)$ becomes very large.

Supposing $x = 1 - \delta$ with $\delta = 10^{-2}, 10^{-4}, 10^{-6}$ successively, the condition number is obtained from $|1/\ln(x)|$ as 99.5, 9999.5, 999999.5 respectively.

* Another example around $x = 1$ is given by:

$$ f(x) = 1 + \sqrt{|x - 1|} $$

The condition number of this function is examined about point $x = 1$.

In the continuous case, $x = 1$ is a singularity since $f'(x) = \frac{1}{2\sqrt{|x-1|}}$

Supposing $x = 1 - \delta$ with $\delta = 10^{-2}, 10^{-4}, 10^{-6}$ successively, the condition number is obtained as 4.59, 49.51 et 511.50 respectively.

* Another example around $x = 0$ is [5]:

$$ f(x) = \frac{e^x - 1}{x}, \quad \lim_{x \to 0} f(x) = 1 $$

Evaluating $f(x)$ at $x = 9.10^{-8}$ yields 1.324455766 in simple precision whereas rewriting $f(x) = \frac{e^x - 1}{\ln(e^x)}$ yields 1.00000005.

It is better to transform it into $\frac{e^x + 1}{x} \tanh(x/2)$ in order to avoid the subtraction operation known to be ill-conditioned.

### E. Non-linear equation conditioning

The simplest example is finding numerically the solution of a quadratic equation. Given $ax^2 + bx + c = 0$, the solution is well-known (since 17th century French mathematician François Viète) to be: $x_{1,2} = \frac{-b \pm \sqrt{b^2 - 4ac}}{2a}$ for $a \neq 0$.

Numerical errors occur when $b^2 \approx 4ac$ because of the ill-conditioned subtraction operation or when either $b^2, ac$ are large.

An example showing the effect of increasing $b$ on the resulting solutions is displayed in Table 1

### F. Differential equation conditioning

The simplest way to describe rounding error is through the description of the first-order ordinary differential equation (ODE) [4]:

$$ \frac{dy}{dx} = 4y(x) - 5e^{-x}, \quad y(0) = 1 \quad (15) $$

The solution is written as a superposition of the general solution: $y_G(x) = C_1 e^{4x}$ and the particular solution: $y_P(x) = C_2 e^{-x}$. Since the initial condition: $y(0) = 1$ the analytical solution is $y(x) = e^{-x}$, thus:

$$ y(x) = e^{-x} + [y(0) - 1]e^{4x} \quad (16) $$
<table><caption>**Table 1:** Numerical and exact solutions of the quadratic $ax^2 + bx + c = 0$ with $a = 1, c = 2$, $b$ variable. As $b$ increases, the values of the roots $x_1, x_2$ are checked with their sum that should be $-b$ and their product that should be $c = 2$. The results show that while the sums respect the checking test, the products do not, signalling the effect of finite arithmetic on the roots accuracy.</caption><thead><tr><th>b</th><th>x<sub>1</sub></th><th>x<sub>2</sub></th><th>(x<sub>1</sub> + x<sub>2</sub>)</th><th>x<sub>1</sub> × x<sub>2</sub></th><th>x<sub>1</sub>(exact)</th><th>x<sub>2</sub>(exact)</th></tr></thead><tbody><tr><td>300.</td><td>-0.006668</td><td>-299.99335</td><td>-300.</td><td>2.00038</td><td>-0.00666</td><td>-299.9933</td></tr><tr><td>600.</td><td>-0.003326</td><td>-599.9967</td><td>-600.</td><td>1.99583</td><td>-0.00333</td><td>-599.9966</td></tr><tr><td>900.</td><td>-0.002227</td><td>-899.9978</td><td>-900.</td><td>2.00499</td><td>-0.00222</td><td>-899.9977</td></tr><tr><td>1200.</td><td>-0.001647</td><td>-1199.9983</td><td>-1200.</td><td>1.97753</td><td>-0.00166</td><td>-1199.9983</td></tr><tr><td>1500.</td><td>-0.001342</td><td>-1499.9987</td><td>-1500.</td><td>2.01415</td><td>-0.00133</td><td>-1499.9986</td></tr><tr><td>1800.</td><td>-0.001098</td><td>-1799.9989</td><td>-1800.</td><td>1.97753</td><td>-0.00111</td><td>-1799.9988</td></tr><tr><td>2100.</td><td>-0.0009765</td><td>-2099.999</td><td>-2100.</td><td>2.05078</td><td>-9.5238E-4</td><td>-2099.9990</td></tr><tr><td>2400.</td><td>-0.0008544</td><td>-2399.999</td><td>-2400.</td><td>2.05078</td><td>-8.3333E-4</td><td>-2399.9991</td></tr><tr><td>2700.</td><td>-0.0007324</td><td>-2699.9993</td><td>-2700.</td><td>1.97753</td><td>-7.4074E-4</td><td>-2699.9992</td></tr><tr><td>3000.</td><td>-0.0006103</td><td>-2999.9995</td><td>-3000.</td><td>1.83105</td><td>-6.6666E-4</td><td>-2999.9993</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>**Fig.5:** (Color on line) Analytical solution of ODE [4] $y(x) = e^{-x}$ flanked by rounding impacted numerical solutions $y_\delta(x) = e^{-x} \pm \delta e^{4x}$ with $\delta = 10^{-7}$</figcaption></figure>

Numerically, the solution is impacted by a rounding error linked to a chosen precision of the initial value $y(0) = 1$ which is handled by the integration algorithm as $y(0) = 1 \pm \delta$ where $\delta = \pm 10^{-7}$ is the rounding error in single precision computation with a float (32 bit) representation of a real number [2]. This means that the computed solution is in fact $y(t) = e^{-x} \pm \delta e^{4x}$.

Fig. 5 shows the analytical solution $y(x) = e^{-x}$ surrounded by the rounding impacted solutions $y_\delta(t) = e^{-x} \pm \delta e^{4x}$.

### III. LINEAR SYSTEMS AND MATRIX CONDITIONING

#### A. Impact of rounding on rank

Consider the (3x3) matrix

$$A = \begin{bmatrix} 1 & 1 & -1 \\ 2 & 2 & -2 \\ 3 & 3 & -3 \end{bmatrix}$$

$\text{rank}(A)=1$ (maximum size of the non zero principal sub-determinant).
Let us suppose that $A$ is perturbed diagonally by an uncertainty $\delta = 10^{-3}$ such that:

$$ \tilde{A} = \begin{bmatrix} 1 + \delta & 1 & -1 \\ 2 & 2 + \delta & -2 \\ 3 & 3 & -3 + \delta \end{bmatrix} $$

Then $\mathrm{rank}(\tilde{A})=3$.

In order to diagnose $A$ and $\tilde{A}$ in order to find the origin of this discrepancy, consider the SVD (cf. section IV) of each matrix:

* The singular eigenvalues of $A$ are $\sigma_1 = \sqrt{42} \approx 6.48$, $\sigma_2 = \sigma_3 = 0$

* The singular eigenvalues of $\tilde{A}$ are $\sigma_1 = 6.48075$, $\sigma_2 = \sigma_3 = 10^{-5}$

For instance, the condition number is obtained from SVD as $\sigma_{max}/\sigma_{min}$ where $\sigma_{max}, \sigma_{min}$ are the largest, smallest singular eigenvalues of the matrix.

This implies $\mathrm{cond}(A) = \infty$ whereas $\mathrm{cond}(\tilde{A}) = 6.48075 \times 10^5$.

## B. Norms of Vectors and Matrices

A vector norm on $\mathbb{R}^n$ is an application (injection), $\| \cdot \|$, from $\mathbb{R}^n$ to $\mathbb{R}$ with the following properties [8, 9]:

1. $\|x\| \ge 0$ for all $x \in \mathbb{R}^n$.

2. $\|x\| = 0$ if and only if $x = 0$.

3. $\|\alpha x\| = |\alpha| \|x\|$ for all $\alpha \in \mathbb{R}$ and $x \in \mathbb{R}^n$.

4. $\|x + y\| \le \|x\| + \|y\|$ for all $x, y \in \mathbb{R}^n$.

The $\mathcal{L}_1, \mathcal{L}_2, \mathcal{L}_\infty$ norms for the vector $x = (x_1, x_2, \dots, x_n)^T$ are defined by

* $\|x\|_1 = \sum_{i=1}^n |x_i|$

* $\|x\|_2 = [\sum_{i=1}^n x_i^2]^{1/2}$

* $\|x\|_\infty = \max_{1 \le i \le n} |x_i|$

The sequence of vectors $x^{(k)}$ converges to $x$ in $\mathbb{R}^n$ with respect to the $\mathcal{L}_\infty$ norm if and only if $\lim_{k \to +\infty} x_i^{(k)} = x_i$, for each $i = 1, 2, \dots, n$.

A matrix norm on the set of all $n \times n$ matrices is a real-valued function, $\| \cdot \|$, satisfying for all $n \times n$ matrices $A$ and $B$ and all real numbers $\alpha$.

1. $\|A\| \ge 0$.

2. $\|A\| = 0$ if and only if $A$ is $0$, the matrix with all $0$ entries.

3. $\|\alpha A\| = |\alpha| \|A\|$.

4. $\|A + B\| \le \|A\| + \|B\|$.

5. $\|AB\| \le \|A\|\|B\|$.

There are other norms [8, 9] such as:

1. Frobenius: $\|A\|_F = \sqrt{\sum_{i,j=1}^n a_{ij}^2} = \sqrt{\mathrm{Tr}(A^T A)}$

2. $\mathcal{L}_\infty$ norm: Considering $\| \cdot \|$ as a vector norm on $\mathbb{R}$, then

$$ \|A\| = \max_{\|x\|=1} \|Ax\| $$

might be considered a matrix norm. If $A = (a_{ij})$ is an $n \times n$ matrix, then

$$ \|A\|_{\infty} = \max_{1 \le i \le n} \sum_{j=1}^{n} |a_{ij}|. $$
### C. Condition numbers

In the case of linear systems ($Ax = b$) the condition number is given by

$$ \frac{\|\delta x\|}{\|x\|} = \mathrm{cond}(A) \frac{\|\delta A\|}{\|A\|} \qquad (17) $$

Regardless of uncertainty affecting matrix $A$ elements such that $A \to A + \delta A$ or vector $b$ elements such that $b \to b + \delta b$, the condition number $\mathrm{cond}(A) = \kappa(A) = \|A\|\|A^{-1}\|$. Then for any natural norm

$$ \|x - \tilde{x}\| = \|A^{-1}r\| \le \|r\|\|A^{-1}\| $$

and if $x \ne 0$ and $b \ne 0$

$$ \frac{\|x - \tilde{x}\|}{\|x\|} \le \frac{\|A\|\|A^{-1}\|}{\|b\|\|r\|} \quad \text{with} \quad \|b\| \le \|A\|\|x\| $$

The condition number of the non-singular matrix $A$ is

$$ 1 \le \|I\| \le \|AA^{-1}\| \le \|A\|\|A^{-1}\| = \kappa(A). $$

A matrix $A$ is well-conditioned, if $\kappa(A)$ is close to 1, and ill-conditioned, when $\kappa(A)$ is significantly greater than 1.

Thus, approximate solutions of the linear system $Ax = b$ can be obtained as:

1. In the non-singular square matrix case: Supposing that $\tilde{x}$ is an approximation to the solution of the linear system $Ax = b$. $A$ is non-singular matrix and $r$ is the residual vector for $\tilde{x}$ is $r = A\tilde{x} - b$. The best approximation is obtained when the minimum residual is obtained, i.e. $\min\|r\| = \min\|A\tilde{x} - b\|$.

2. In the rectangular matrix case, when iterative methods for solving the symmetric linear system $Ax = b$ are used:

$$ \min \|x\|_2 \text{ such that } \min_x \|Ax - b\|_2, \qquad (18) $$

where $A \in \mathbb{R}^{m \times n}$ represents an over-determined system ($m > n$) with more equations than unknowns, $b \in \mathbb{R}^n$, $A \ne 0$, and $b \ne 0$ (cf. section IV D).

Note that there are other condition numbers [8, 9] such as:

* Von Neumann: $\mathrm{cond}_V(A) = |\lambda_{max}/\lambda_{min}|$ where $\lambda_{max}, \lambda_{min}$ are matrix $A$ largest and smallest eigenvalues.

* Turing: $\mathrm{cond}_T(A) = n(\max_{1 \le i,j \le n} |A_{ij}|)(\max_{1 \le i,j \le n} |A_{ij}^{-1}|)$ [10]

* Singular Eigenvalues: $\mathrm{cond}_S(A) = |\sigma_{max}/\sigma_{min}|$ where $\sigma_{max}, \sigma_{min}$ are matrix $A$ largest and smallest singular eigenvalues obtained with SVD (Singular Value Decomposition described in next section).

## IV. SINGULAR VALUE DECOMPOSITION

SVD is a generalization of the similarity transformation given by: $B = PDP^{-1}$ where $B$ is $n \times n$ square matrix to any $m \times n$ rectangular matrix $A$, i.e. $A \in \mathbb{R}^{m \times n}$. $D$ is the diagonalized $A$ matrix and $P$ is the matrix of base change from the original basis to the diagonal one. SVD transformation is given by: $A = U\Sigma V^T$ where $U \in \mathbb{R}^{m \times n}$, $\Sigma \in \mathbb{R}^{n \times m}$, $V \in \mathbb{R}^{m \times n}$. Thus $\Sigma$ is the counterpart of $D$ whereas $U, V$ are the counterparts of $P, P^{-1}$ matrices. The orthogonal matrices $U$ and $V$ are called the left and right singular vectors, respectively.

In practice, the singular eigenvalues and eigenvectors of a matrix $A$ evaluated by a singular value decomposition are such that:

$$ A = U\Sigma V^T \qquad (19) $$
$\Sigma$ is a diagonal matrix of singular values:

$$ \Sigma = \begin{pmatrix} \sigma_1 & 0 & 0 & \dots & 0 \\ 0 & \sigma_2 & 0 & \dots & \vdots \\ 0 & 0 & \ddots & \ddots & 0 \\ \vdots & \ddots & \ddots & \ddots & \vdots \\ 0 & \dots & 0 & 0 & \sigma_n \end{pmatrix}. \qquad (20) $$

Since $\Sigma$ is a diagonal matrix it can be simply expressed as: $\Sigma = \text{diag}(\sigma_1, \sigma_2, \dots, \sigma_n)$.

The singular eigenvalues are sorted in descending order such that:

$$ \sigma_1 \ge \sigma_2 \ge \sigma_3 \dots \ge \sigma_n \ge 0 $$

$\sigma_1$ is the Euclidean $\mathcal{L}_2$ norm of $A$. Note that $\|U\| = \|V\| = 1$ and $UU^T = VV^T = I$.

SVD solves the "Generalized eigenvalue problem" by supposing two sets of singular eigenvectors for $A$ such that:

* $Av_j = \sigma_j u_j$

* $A^T u_j = \sigma_j v_j$

This implies that $A^T A = V\Sigma U^T U\Sigma V^T = V\Sigma^2 V^T$ leading to the eigenvalue relation: $A^T A v_j = \sigma_j^2 v_j$.

Similarly: $AA^T = U\Sigma V^T V\Sigma U^T = U\Sigma^2 U^T$ leading to the eigenvalue relation: $AA^T u_j = \sigma_j^2 u_j$.

The SVD algorithm is first applied to the 2x2 matrix: $A = \begin{bmatrix} 3.000 & 5.000 \\ 4.000 & 0.000 \end{bmatrix}$

The eigenvalues of $A^T A == \begin{bmatrix} 25.000 & 15.000 \\ 15.000 & 25.000 \end{bmatrix}$

are given by $\lambda_1 = 40, \lambda_2 = 10$, implying $\sigma_1 = \sqrt{40} \approx 6.3, \sigma_2 = \sqrt{10} \approx 3.1$.

Thus matrix $\Sigma = \begin{bmatrix} \sigma_1 & 0 \\ 0 & \sigma_2 \end{bmatrix}$, is expressed as $\Sigma = \text{diag}(\sigma_1, \sigma_2)$.

The singular eigenvectors of $A^T A$ are: $v_1 = \begin{bmatrix} \frac{1}{\sqrt{2}} \\ \frac{1}{\sqrt{2}} \end{bmatrix}$ and $v_2 = \begin{bmatrix} -\frac{1}{\sqrt{2}} \\ \frac{1}{\sqrt{2}} \end{bmatrix}$ yielding matrix $V$ built column-wise from $v_1, v_2$ columns as: $V = \begin{bmatrix} \frac{1}{\sqrt{2}} & -\frac{1}{\sqrt{2}} \\ \frac{1}{\sqrt{2}} & \frac{1}{\sqrt{2}} \end{bmatrix}$.

Similarly, the singular eigenvectors of $AA^T$ are determined from the relations: $u_j = Av_j/\sigma_j, j = 1, 2$ yielding:

$u_1 = \begin{bmatrix} \frac{2}{\sqrt{5}} \\ \frac{1}{\sqrt{5}} \end{bmatrix}$ and $u_2 = \begin{bmatrix} \frac{1}{\sqrt{5}} \\ -\frac{1}{\sqrt{5}} \end{bmatrix}$

Thus matrix $U$ built column-wise from $u_1, u_2$ columns as: $U = \begin{bmatrix} \frac{2}{\sqrt{5}} & \frac{1}{\sqrt{5}} \\ \frac{1}{\sqrt{5}} & -\frac{2}{\sqrt{5}} \end{bmatrix}$.

**Note:** The non-singular eigenvalues of $A$ are $\lambda_1(A) = (3 + \sqrt{89})/2 \approx 6.2, \lambda_2(A) = (3 - \sqrt{89})/2 \approx -3.2$, implying $\sigma_1 \approx \lambda_1(A), \sigma_2 \approx |\lambda_2(A)|$.

To summarize, the steps for SVD decomposition of $A = U\Sigma V^T \in \mathbb{R}^{m \times n}$ where $U \in \mathbb{R}^{m \times n}, \Sigma \in \mathbb{R}^{n \times n}, V \in \mathbb{R}^{n \times n}$ are given by:

* Determine $A^T A$ eigenvalues: $\lambda_i, i = 1, n$ and singular eigenvalues: $\sigma_i = \sqrt{\lambda_i}, i = 1, n$

* Build $(n \times n)$ matrix $\Sigma = \text{diag}(\sigma_1, \sigma_2, \sigma_3, \dots, \sigma_n)$

* Determine $A^T A$ eigenvectors $v_j, j = 1, \dots, n$ as columns making the $(n \times n)$ $V$ matrix

* Determine $AA^T$ eigenvectors $u_j, j = 1, \dots, n$ from the relations $u_j = Av_j/\sigma_j, j = 1, \dots, n$ and build $(m \times n)$ matrix $U$ from $u_j$ vectors as columns.
Let us apply this procedure to the 2x2 matrix mentioned in the Introduction section:

$$A\boldsymbol{x} = \begin{pmatrix} 1 & 1 \\ a & 1 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \end{pmatrix} = \boldsymbol{b} = \begin{pmatrix} 1 \\ 2 \end{pmatrix}$$

For $a=1.1$, the singular values are: $\sigma_1 = 2.051249, \sigma_2 = 0.048751$ yielding $\kappa = 42.7$, and for $a=1.05$, the singular values are: $\sigma_1 = 2.025312, \sigma_2 = 0.024688$ yielding $\kappa = 84.2$ explaining the 20 amplification factor evaluated previously.

Another example is given by a square 3x3 matrix $A$: $A = \begin{bmatrix} 2.000 & 1.000 & 0.000 \\ -4.000 & -3.000 & 3.000 \\ 4.000 & 1.000 & -2.000 \end{bmatrix}$

$\mathcal{L}_2$ norm requires SVD decomposition of matrix $A$ such that $A = U\Sigma V^T$.

The SVD algorithm [2, 11] yields:

$$\Sigma = \begin{bmatrix} 7.563 & 0.000 & 0.000 \\ 0.000 & 1.365 & 0.000 \\ 0.000 & 0.000 & 0.969 \end{bmatrix}, \text{ expressed as } \Sigma = \text{diag}(7.563, 1.365, 0.969).$$

$$V^T = \begin{bmatrix} -0.786 & -0.415 & 0.459 \\ -0.613 & 0.622 & -0.486 \\ -0.084 & -0.663 & -0.744 \end{bmatrix} \quad \text{and} \quad U = \begin{bmatrix} -0.263 & -0.442 & -0.858 \\ 0.762 & -0.640 & 0.097 \\ -0.592 & -0.628 & 0.505 \end{bmatrix},$$

Matrix $\Sigma$ yields $\sigma_{max} = 7.563$, i.e. $A$ matrix $\mathcal{L}_2$ norm.

Another example is for a (3x2) rectangular matrix:

$$A = \begin{bmatrix} 1.000000 & 4.000000 \\ 2.000000 & 5.000000 \\ 3.000000 & 0.000000 \end{bmatrix}.$$

$$\Sigma = \begin{bmatrix} 6.851909 & 0 \\ 0 & 2.837490 \end{bmatrix}, \text{ i.e. } \Sigma = \text{diag}(6.851909, 2.837490).$$

$$U = \begin{bmatrix} -0.226926 & 0.594362 \\ -0.040388 & 0.785758 \\ 0.973074 & 0.171222 \end{bmatrix} \quad \text{and} \quad V^T = \begin{bmatrix} 0.920363 & -0.391065 \\ 0.391065 & 0.920363 \end{bmatrix}.$$

SVD performs exceptional analysis on any $n \times m$ matrix $A$ by ordering its eigenvalues and paving the way to diagnose its properties and cure them. Basically SVD reveals the internal structure of the matrix and allows to redesign it at will. Some of the analysis, structuring, repair and applications of SVD are:

1. Evaluation of the condition number from $\sigma_1/\sigma_n$

2. Solving linear systems that might be singular, over-determined or underdetermined by finding $\min \|Ax - b\|$.

3. Regularization of $A$ by approximating it with a reduced set of eigenvalues larger than a chosen cutoff with the transformation: $\Sigma = \text{diag}(\sigma_1, \sigma_2, \sigma_3...\sigma_n) \rightarrow \tilde{\Sigma} = \text{diag}(\sigma_1, \sigma_2, \sigma_3...\sigma_p)$ where $p < n$ corresponds to the cutoff singular eigenvalue. This implies a reduction of the condition number since it went from $\sigma_1/\sigma_n$ to $\sigma_1/\sigma_p$ where $\sigma_p > \sigma_n$.

4. In Digital Signal Processing, SVD can help filtering noisy data [11].

5. In Image processing, SVD can help restore a blurred image [12].

6. SVD can also intervene in image compression [13].

## A. Hilbert matrices

A general ($n \times n$) Hilbert matrix is defined by: $A_n = [a_{ij}^{(n)}]$ such that $a_{ij}^{(n)} = \frac{1}{i+j-1}$, $i, j = 1, 2, \dots, n$.

Thus:
$$A_n = \begin{bmatrix} 1 & 1/2 & 1/3 & \cdots & 1/n \\ 1/2 & 1/3 & 1/4 & \cdots & 1/(n+1) \\ 1/3 & 1/4 & 1/5 & \cdots & 1/(n+2) \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1/n & 1/(n+1) & 1/(n+2) & \cdots & 1/(2n+1) \end{bmatrix}$$

The Von Neumann condition number is $cond_V(A_n) = e^{3.5n}$ whereas the Turing value is $cond_T(A_n) = \text{constant} \times e^{3.525n}$ [14].

Consider a Hilbert matrix of order 4 [3]:

$$A_4 = \begin{bmatrix} 1 & 1/2 & 1/3 & 1/4 \\ 1/2 & 1/3 & 1/4 & 1/5 \\ 1/3 & 1/4 & 1/5 & 1/6 \\ 1/4 & 1/5 & 1/6 & 1/7 \end{bmatrix}$$

LU decomposition yields:

$$L = \begin{bmatrix} 1 & 0 & 0 & 0 \\ 1/2 & 1 & 0 & 0 \\ 1/3 & 1 & 1 & 0 \\ 1/4 & 9/10 & 3/2 & 1 \end{bmatrix}, \quad U = \begin{bmatrix} 1 & 1/2 & 1/3 & 1/4 \\ 0 & 1/12 & 1/12 & 3/40 \\ 0 & 0 & 1/180 & 1/120 \\ 0 & 0 & 0 & 1/2800 \end{bmatrix}$$

The inverse matrix is given column by column [2] by: $Ax_i = e_i$ avec $i = 1, 2, 3, 4$ et $[e_i]_j = \delta_{ij}$. Writing $Ly_i = e_i$ avec $Ux_i = y_i$, the inverse matrix writes:

$$A_4^{-1} = \begin{bmatrix} 16 & -120 & 240 & -140 \\ -120 & 1200 & -2700 & 1680 \\ 240 & -2700 & 6480 & -4200 \\ -140 & 1680 & -4200 & 2800 \end{bmatrix}$$

Hilbert matrices are strongly ill-conditioned. In effect, $cond(A)=||A||||A^{-1}||$. The $L_2$ A norm is obtained by SVD as $||A||_2 = 1.5$ since:

$$\Sigma = diag(1.500214, 0.169141, 0.006738, 0.000097)$$

with $A = U\Sigma V^T$ where:

$$U = \begin{bmatrix} -0.792 & 0.582 & -0.179 & 0.029 \\ -0.451 & -0.370 & 0.741 & -0.328 \\ -0.322 & -0.509 & -0.100 & 0.791 \\ -0.252 & -0.514 & -0.638 & -0.514 \end{bmatrix}, \quad V^T = \begin{bmatrix} -0.792 & -0.451 & -0.322 & -0.252 \\ 0.582 & -0.370 & -0.509 & -0.514 \\ -0.179 & 0.741 & -0.100 & -0.638 \\ 0.029 & -0.328 & 0.791 & -0.514 \end{bmatrix},$$

The condition number is obtained from matrix $\Sigma$ as $\sigma_{max}/\sigma_{min} = 15000$ by evaluating the $L_2$ norm of $A_4^{-1}$ by SVD that yields $||A_4^{-1}||_2 = 10341$.

Thus the condition number is: $||A_4||||A_4^{-1}|| = 1.5 \times 10341 = 15511$. Analytically, the $L_2$ norm is: $cond(A) \approx O((1+\sqrt{2})^{4n}/\sqrt{n})$ that becomes $\exp(3.5n)$ for large $n$ ($n$ is matrix size). A norm $L_1$ is $25/12$ whereas $A^{-1}$ norm is 13620 thus $cond_V(A_4) = 28375$.

**Note:** Consequently the number of significant digits of the linear system $Ax = b$ is off by 4 since $cond_V(A_4) = 28375 \propto 10^4$
## B. Wilkinson matrix

The 4x4 Wilkinson matrix is lower triangular [4]. It is a strongly ill-conditioned and given by:

$$A = \begin{bmatrix} 0.9143 \times 10^{-4} & 0.0000 & 0.0000 & 0.0000 \\ 0.8762 & 0.7156 \times 10^{-4} & 0.0000 & 0.0000 \\ 0.7943 & 0.8143 & 0.9504 \times 10^{-4} & 0.0000 \\ 0.8017 & 0.6123 & 0.7165 & 0.7123 \times 10^{-4} \end{bmatrix}$$

$$A^{-1} = \begin{bmatrix} 10937.32 & -1.33919 \times 10^{8} & 1.14732 \times 10^{12} & -1.15397 \times 10^{16} \\ 0 & 13974.28 & -1.197312 \times 10^{8} & 1.20425 \times 10^{12} \\ 0 & 0 & 10521.88 & -1.0583 \times 10^{8} \\ 0 & 0 & 0 & 14039.02 \end{bmatrix}$$

$\Sigma$ is given by: $\Sigma = diag(1.7299, 0.5931, 0.4981, 8.6656 \times 10^{-17})$

$$U = \begin{bmatrix} -0.000042 & 0.000084 & -0.000045 & 1.000000 \\ -0.405866 & 0.806814 & -0.429329 & -0.000104 \\ -0.616957 & 0.104702 & 0.780001 & 0.000000 \\ -0.674267 & -0.581453 & -0.455275 & -0.000000 \end{bmatrix} \quad \text{and} \quad V^T = \begin{bmatrix} -0.801301 & -0.529067 & -0.279294 & -0.000028 \\ 0.546179 & -0.456421 & -0.702402 & -0.000070 \\ -0.244142 & 0.715380 & -0.654695 & -0.000065 \\ -0.000000 & 0.000000 & -0.000099 & 1.000000 \end{bmatrix}$$

The norms are given by: $\|A\|_1 = 2.13057$, $\|A\|_{\infty} = 2.4722$, $\|A\|_F = 1.8954$.

On the other hand: $\|A^{-1}\|_1 = 1.1540 \times 10^{16}$, $\|A^{-1}\|_{\infty} = 1.1540 \times 10^{16}$, $\|A^{-1}\|_F = 1.1539 \times 10^{16}$.

**Conclusion:** All condition numbers $\|A\|\|A^{-1}\|$ are around $10^{16}$.

## C. Pascal matrices

An $(n \times n)$ Pascal matrix $P_n$ is given by combinatorial factors:

$$P_{ij} = \frac{(i+j-2)!}{(i-1)!(j-1)!} = \binom{i+j-2}{j-1} \quad (21)$$

It has a condition number given by: $\kappa(P_n) \sim 16^n/(n\pi)$

The 4x4 Pascal matrix is:

$$P_4 = \begin{bmatrix} 1 & 1 & 1 & 1 \\ 1 & 2 & 3 & 4 \\ 1 & 3 & 6 & 10 \\ 1 & 4 & 10 & 20 \end{bmatrix}$$

SVD yields the following: $\Sigma = diag(26.3047, 2.20345, 0.4538, 0.03802)$,

$$U = \begin{bmatrix} -0.06019 & -0.5304 & 0.7873 & -0.3087 \\ -0.2012 & -0.6403 & -0.1632 & 0.7231 \\ -0.4581 & -0.3918 & -0.5321 & -0.5946 \\ -0.8638 & 0.3939 & 0.2654 & 0.1684 \end{bmatrix} \quad \text{and} \quad V^T = \begin{bmatrix} -0.06019 & -0.2012 & -0.4581 & -0.8638 \\ -0.5304 & -0.6403 & -0.3918 & 0.3939 \\ 0.7873 & -0.1632 & -0.5321 & 0.2654 \\ -0.3087 & 0.7231 & -0.5946 & 0.1684 \end{bmatrix}$$

The condition number is given by $\sigma_{max}/\sigma_{min} = 876$ while $16^n/(n\pi) \sim 5215$.
### D. Solution of over-determined linear systems

In over-determined linear systems $Ax = b$, $A$ is a rectangular $m \times n$ matrix with $m > n$, meaning there are more equations than unknowns, it is possible to find a solution by multiplying $A$ with its transpose $A^T$ such that: $A^T A x = A^T b$, then the solution is found from: $x = (A^T A)^{-1} A^T b$.

A difficulty might arise when $A$ is singular and $(A^T A)^{-1}$ does not exist, thus it is possible to use SVD to overcome the singularity and find a minimum-length solution (also called pseudo-inverse matrix solution) that differs from the residual approach defined by: $min||r|| = min||A\tilde{x} - b||$.

The residual method might lead to a non-unique solution, thus it is better to turn to a minimum-length approach that leads to a unique solution: $min ||x||_2$ such that $min_x ||Ax - b||_2$.

A pseudo-inverse matrix $A^+$ has been defined by Moore and Penrose [1]. Using SVD, $A^+$ denotes the $(n \times m)$ pseudo-inverse of $A$ defined as $A^+ = V\Sigma^{-1}U^T$ where $V, \Sigma, U$ are the SVD components of $A = U\Sigma V^T$.

As an illustration, consider the following example:

$$A = \begin{bmatrix} -11 & 2 \\ 2 & 3 \\ 2 & -1 \end{bmatrix}$$

SVD decomposition yields: $\Sigma = diag(11.4783, 3.35387)$ revealing the singularity originating from a single zero eigenvalue.

$$U = \begin{bmatrix} -0.9737 & -0.09387 & 0.2078 \\ 0.1327 & -0.9743 & 0.1818 \\ 0.1854 & 0.2046 & 0.9611 \end{bmatrix} \quad \text{and} \quad V^T = \begin{bmatrix} 0.9885 & -0.1511 \\ -0.1511 & -0.9885 \end{bmatrix}$$

Note that $V^T$ is symmetric thus $V = V^T$ and matrix size compatibility requires: $\Sigma$ of size 2x2, $U$ of size 3x3 and $V$ of size 2x2.

Evaluation of $A^+ = V\Sigma^{-1}U^T$ requires inverting $\Sigma$ given by: $\Sigma^{-1} = diag(1/\sigma_1, 1/\sigma_2, \dots)$, however matrix size compatibility requires that $\Sigma^{-1}$ to be of size 2x3, thus:

$$\Sigma^{-1} = \begin{bmatrix} 1/11.4783 & 0 & 0 \\ 0 & 1/3.35387 & 0 \\ 0 & 0 & 0 \end{bmatrix}$$

The pseudo-inverse matrix $A^+ = V\Sigma^{-1}U^T$ is 2x3 obtained as:

$$A^+ = \begin{bmatrix} -0.0796 & 0.05532 & 0.0067 \\ 0.0404 & 0.28541 & 0.0627 \end{bmatrix}$$

Validation of the result is given by Moore-Penrose conditions: $A^+AA^+ = A^+$ and $AA^+A = A$. Indeed we get:

$$AA^+A = \begin{bmatrix} -11.0002 & 2.00024 \\ 2.00007 & 2.99979 \\ 2.00003 & -1.00001 \end{bmatrix} \equiv A \quad \text{and} \quad A^+AA^+ = \begin{bmatrix} -0.07963 & 0.05531 & 0.006751 \\ 0.04048 & 0.2854 & -0.06274 \end{bmatrix} \equiv A^+$$

Note incidentally that $A^+A$ is given by:

$$A^+A = \begin{bmatrix} 1.000017 & -3.103 \times 10^{-5} \\ 9.862 \times 10^{-6} & 0.99994 \end{bmatrix}$$
implying $A^+$ indeed plays the role of an inverse matrix.

Consequently, the minimum-length solution of the over-determined system is: $Ax = b$ is $x^+ = A^+b$.

Selecting $b = \begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix}$ the minimum-length solution is $x^+ = \begin{pmatrix} -0.079600 \\ 0.040400 \end{pmatrix}$ with residual $\|Ax^+ - b\|_2 = 0.207810$.

In order to validate the solution, we use Powell method [2] (cf. Appendix) to minimize the residual written in the $L_2$ norm (cf. Fig. 6), as:

$$ \|Ax - b\|_2 = \sqrt{(-11x_1 + 2x_2 - 1)^2 + (2x_1 + 3x_2)^2 + (2x_1 - x_2)^2} \text{ with } x = \begin{pmatrix} x_1 \\ x_2 \end{pmatrix}. \text{ The Powell minimum solution is} $$
$$ x_P = \begin{pmatrix} -0.079622 \\ 0.040486 \end{pmatrix} \text{ with norms: } \|x_P\|_1 = 0.12010705, \|x_P\|_2 = 0.089323275 \text{ and residual } \|Ax_P - b\|_2 = 0.207810 $$

agree with the Pseudo-inverse method.

However, if we use the $L_1$ norm, the residual (cf. Fig. 6) is written as:

$$ \|Ax - b\|_1 = |-11x_1 + 2x_2 - 1| + |2x_1 + 3x_2| + |2x_1 - x_2| \text{ with we get another solution } x_P = \begin{pmatrix} -0.081079 \\ 0.054053 \end{pmatrix} \text{ with norms:} $$

$\|x_P\|_1 = 0.1351316, \|x_P\|_2 = 0.09744477$ and residual $\|Ax_P - b\|_1 = 0.216237$ agreeing with the Pseudo-inverse method.

<figure><img src="image_8.png" /><figcaption>Fig.6: $L_1$ (at left) and $L_2$ (at right) norm residuals $\|Ax - b\|_1$ and $\|Ax - b\|_2$ as function of variables $x_1, x_2$. $L_1$ surface shows two minima close to the $x_1, x_2$ plane origin whereas $L_2$ surface shows a line of minima close to the $x_1, x_2$ plane origin.</figcaption></figure>

In order to reveal the nature of the minimum, Broyden method (cf. Appendix) is the preferred tool that allows to address and minimize both residuals simultaneously:

$$ \|Ax - b\|_2 = \sqrt{(-11x_1 + 2x_2 - 1)^2 + (2x_1 + 3x_2)^2 + (2x_1 - x_2)^2} $$

$$ \text{and } \|Ax - b\|_1 = |-11x_1 + 2x_2 - 1| + |2x_1 + 3x_2| + |2x_1 - x_2|. $$

The solution is $x_B = \begin{pmatrix} -0.079811 \\ 0.047653 \end{pmatrix}$ with norms: $\|x_B\|_1 = 0.12746358, \|x_B\|_2 = 0.09295451$ and residuals $\|Ax_B - b\|_1 = 0.250714, \|Ax_B - b\|_2 = 0.209660$ that both agree with the Pseudo-inverse method.

## V. ITERATIVE METHODS FOR SOLVING LINEAR SYSTEMS

Iterative solving methods of linear systems are now treated since Preconditioning techniques for them are tackled in the next section.
## A. Jacobi Iterative Method

Given a vector $x_0 \in \mathbb{R}^n$, the goal is to generate a sequence $\{x_k\}_{k=1}^\infty$ that converges to the solution [15] of (2). Suppose that $a_{i,i} \neq 0$ for all $i$, then we can rewrite (2) as

$$x_i = \frac{1}{a_{i,i}} \left( b_i - \sum_{\substack{j=1 \\ j \neq i}}^{n} a_{i,j} x_j \right)$$

for $i = 1, 2, \dots, n$. This formula motivates the following algorithm.

### Algorithm

1. Choose a vector $x_0$ close (if possible) to the solution of $Ax = b$.

2. Given the vector $x_k$, compute the vector $x_{k+1}$ as follows:

$$x_{k+1,i} = \frac{1}{a_{i,i}} \left( b_i - \sum_{\substack{j=1 \\ j \neq i}}^{n} a_{i,j} x_{k,j} \right) \quad (22)$$

for $i = 1, 2, \dots, n$.

3. Repeat (2) until $\|x_{k+1} - x_k\| < \epsilon$, where $\epsilon$ is given.

However, conditions are required for matrix $A$ to ensure that the sequence $\{x_k\}_{k=1}^\infty$ converges to a solution of $Ax = b$.

## B. Gauss-Seidel Iterative Method

As for Jacobi iterative method [15], given a vector $x_0 \in \mathbb{R}^n$, the goal is to generate a sequence $\{x_k\}_{k=1}^\infty$ that converges to the solution of (2).

If we use $x_{k+1,1}, x_{k+1,2}, \dots, x_{k+1,i-1}$ instead of $x_{k,1}, x_{k,2}, \dots, x_{k,i-1}$ in the formula (22) to compute $x_{k+1,i}$, hoping that $x_{k+1,1}, x_{k+1,2}, \dots, x_{k+1,i-1}$ are better approximations of the first $(i-1)$ coordinates of the solution of (2) than $x_{k,1}, x_{k,2}, \dots, x_{k,i-1}$, then perhaps that will get a new sequence $\{x_k\}_{k=0}^\infty$ that converges faster to the solution of (2). This motivates the following algorithm.

### Algorithm

1. Choose a vector $x_0$ closed to the solution of $Ax = b$ (if possible).

2. Given the vector $x_k$, compute the vector $x_{k+1}$ as follows:

$$x_{k+1,i} = \frac{1}{a_{i,i}} \left( b_i - \sum_{j=1}^{i-1} a_{i,j} x_{k+1,j} - \sum_{j=i+1}^{n} a_{i,j} x_{k,j} \right) \quad (23)$$

for $i = 1, 2, \dots, n$.

3. Repeat (2) until $\|x_{k+1} - x_k\| < \epsilon$, where $\epsilon$ is given.

As for the Jacobi iterative method, we need conditions on the matrix $A$ to ensure that the sequence $\{x_k\}_{k=1}^\infty$ converges to a solution of $Ax = b$.

We consider systems of the form $Ax = b$, where $A$ is a strictly positive definite matrix (this also means that $A$ is symmetric).
### C. Relaxation Methods

As for Jacobi and Gauss-Seidel iterative methods, given a vector $x_0 \in \mathbb{R}^n$, the goal is to generate a sequence $\{x_k\}_{k=1}^\infty$ that converges to the solution of (2). The classical relaxation methods [15] are given in the following algorithm [15].

**Algorithm**

1. Choose a real number $\omega$ between 0 and 2. $\omega$ is a relaxation parameter and its choice will be justified later.

2. Choose a vector $x_0$ closed to the solution of $Ax = b$ (if possible).

3. Given the vector $x_k$, compute the vector $x_{k+1}$ as follows:

$$x_{k+1,i} = x_{k,i} + \frac{\omega}{a_{i,i}} \left( b_i - \sum_{j=1}^{i-1} a_{i,j}x_{k+1,j} - \sum_{j=i}^{n} a_{i,j}x_{k,j} \right) \quad (24)$$

for $i = 1, 2, \dots, n$.

4. Repeat (3) until $\|x_{k+1} - x_k\| < \epsilon$, where $\epsilon$ is given.

The previous algorithm is called an under-relaxation method for $0 < \omega < 1$, and an over-relaxation method or a successive over-relaxation (SOR) method for $1 < \omega < 2$.

We can write $Ax = b$ with $A = D - L - U$ as:

$$-Lx = -Dx + Ux + b.$$

Multiplying both sides by a non-zero factor $\omega$ and adding $Dx$ on both sides yield

$$(D - \omega L)x = (1 - \omega)Dx + \omega Ux + \omega b.$$

Finally, multiplying by $(D - \omega L)^{-1}$ from the left on both sides of the equality above gives

$$x = (D - \omega L)^{-1} ((1 - \omega)D + \omega U)x + \omega(D - \omega L)^{-1}b. \quad (25)$$

With $T = (D - \omega L)^{-1}((1 - \omega)D + \omega U)$ and $c = \omega(D - \omega L)^{-1}b$, the equation (25) becomes $x = Tx + c$, which is equivalent to $Ax = b$.

In a similar manner, defining the backward SOR method:

$$(D - \omega U)x_{k+1} = [\omega L + (1 - \omega)D]x_k + \omega b,$$

and combining it with a step of SOR method, the symmetric SOR method (SSOR), is obtained:

$$x_{k+1} = B_s(\omega)x_k + b_{\omega}, \quad (26)$$

where

$$B_s(\omega) = (D - \omega U)^{-1}(\omega L + (1 - \omega)D)(D - \omega L)^{-1}(\omega U + (1 - \omega)D), \\ b_{\omega} = \omega(2 - \omega)(D - \omega U)^{-1}D(D - \omega L)^{-1}b \quad (27)$$

### D. Steepest Descent

The basis for the method of steepest descent is the result of the following Proposition.

Let $A$ be a strictly positive definite matrix and $b \in \mathbb{R}^n$. The solution of the linear system $Ax = b$ is the point where the quadratic function

$$g: \mathbb{R}^n \to \mathbb{R}, \quad x \mapsto \langle x, Ax \rangle - 2 \langle x, b \rangle \quad (28)$$
reaches its absolute minimum.

The previous proposition suggests the following algorithm.

**Algorithm**

1. Choose a vector $\mathbf{x}_0$ closed to the solution of $Ax = \mathbf{b}$ (if possible).

2. Given the vector $\mathbf{x}_k$, choose a direction $\mathbf{u}_k$ such that $\langle \mathbf{u}_k, \mathbf{b} - A\mathbf{x}_k \rangle \neq 0$. If no such vector exists, then $A\mathbf{x}_k = \mathbf{b}$ and the solution has been found.

3. Compute

$$t_k = \frac{\langle \mathbf{u}_k, \mathbf{b} - A\mathbf{x}_k \rangle}{\langle \mathbf{u}_k, A\mathbf{u}_k \rangle}$$

and let $\mathbf{x}_{k+1} = \mathbf{x}_k + t_k \mathbf{u}_k$.

4. Repeat (2) to (3) until $\|\mathbf{x}_{k+1} - \mathbf{x}_k\| < \epsilon$, where $\epsilon$ is given.

The third step of the steepest descent algorithm is inferred as in the proof of Proposition V D by considering $q(t) = g(\mathbf{x}_k + t\mathbf{u}_k)$ where $g$ is defined in eq. 28.

The vectors $\mathbf{u}_k$ are selected to speed up algorithm convergence. In particular, the size of $t_k$ should be controlled in order to avoid overshooting the solution.

The next proposition shows that, in theory, the steepest descent algorithm ends after a finite number of iterations. However, this does not generally the case in the computer implementation of this algorithm due to roundoff errors, ill-conditioning, ...

Let $A$ be a strictly positive definite matrix and $\mathbf{b} \in \mathbb{R}^n$. Suppose that the vectors $\mathbf{u}_1, \mathbf{u}_2, \dots, \mathbf{u}_n$ are $A$-orthogonal vectors in $\mathbb{R}^n$; namely, $\langle \mathbf{u}_i, A\mathbf{u}_j \rangle = 0$ for $i \neq j$. Then the steepest descent algorithm produces the solution of $Ax = \mathbf{b}$ after $n$ steps.

**E. Conjugate Gradient**

The conjugate gradient algorithm is a special case of the steepest descent algorithm, where the vectors $\mathbf{u}_j$ are chosen such that the vectors $\mathbf{r}_j = \mathbf{b} - A\mathbf{x}_j$ are mutually orthogonal.

**Algorithm**

1. Choose a vector $\mathbf{x}_0$ closed to the solution of $Ax = \mathbf{b}$ (if possible).

2. Let $\mathbf{r}_0 = \mathbf{b} - A\mathbf{x}_0$ and $\mathbf{u}_0 = \mathbf{r}_0$.

3. Given the vectors $\mathbf{u}_k \neq \mathbf{0}$ and $\mathbf{r}_k$, compute

$$t_k = \frac{\langle \mathbf{r}_k, \mathbf{r}_k \rangle}{\langle \mathbf{u}_k, A\mathbf{u}_k \rangle}.$$

4. Given the vectors $\mathbf{x}_k$ and $\mathbf{r}_k$, let $\mathbf{x}_{k+1} = \mathbf{x}_k + t_k \mathbf{u}_k$ and $\mathbf{r}_{k+1} = \mathbf{r}_k - t_k A\mathbf{u}_k$.

5. Stop if $\|\mathbf{r}_{k+1}\|_2^2 < \epsilon$, where $\epsilon$ is given.

6. Compute

$$s_k = \frac{\langle \mathbf{r}_{k+1}, \mathbf{r}_{k+1} \rangle}{\langle \mathbf{r}_k, \mathbf{r}_k \rangle}.$$

7. Let $\mathbf{u}_{k+1} = \mathbf{r}_{k+1} + s_k \mathbf{u}_k$.
8. Repeat (3) to (7) until the condition in (5) is satisfied.

The next theorem shows that the $t_k$'s used in the conjugate gradient algorithm are of the form

$$t_k = \frac{\langle b - A\mathbf{x}_k, \mathbf{u}_k \rangle}{\langle \mathbf{u}_k, A\mathbf{u}_k \rangle}$$

as required for the steepest descent method.

**Theorem:**

The vectors $\mathbf{x}_k$ and $\mathbf{r}_k$ of the conjugate gradient algorithm satisfy $\mathbf{r}_k = \mathbf{b} - A\mathbf{x}_k$ and $\langle \mathbf{r}_i, \mathbf{r}_j \rangle = 0$ for $i \neq j$ as long as $\mathbf{u}_k \neq 0$. ■

# VI. PRECONDITIONING TECHNIQUES FOR ITERATIVE SOLVERS OF LINEAR SYSTEMS

A preconditioner [16] for a linear system $A\mathbf{x} = \mathbf{b}$ solver allows to solve the same problem with fewer iterations.

Preconditioning transforms the linear system into $M^{-1}A\mathbf{x} = M^{-1}\mathbf{b}$ where $M^{-1}$ is an inverse conditioning matrix that allows to speed convergence and reduce the number of iterations leading to the solution.

The recipe for finding matrix $M$ is given below:

* Jacobi preconditioning is done with $M = D$

* Gauss-Seidel preconditioning is done with $M = D + L$

* SOR preconditioning is done with $M = \omega^{-1}D + L$

* SSOR (Symmetric SOR) preconditioning is done with $M = (2 - \omega)^{-1}(\omega^{-1}D + L)(\omega^{-1}D)^{-1}(\omega^{-1}D + U)$

If $A$ is symmetric and positive definite, the SSOR method is convergent if $0 < \omega < 2$. Typically, the SSOR method with an optimal choice of the relaxation parameter converges more slowly than the corresponding SOR method.

Generally the construction of preconditioners is in fact problem-dependent. If not much is known about the structure of $A$, we can only consider general methods such as diagonal preconditioning, binormalization... These methods require access to the nonzero elements of $A$. (They are not applicable if $A$ exists only as an operator for returning the product $A\mathbf{x}$.)

Benzi[17] is a comprehensive survey of preconditioning techniques. We discuss three methods for symmetric $A$ requiring access to the nonzero $A_{ij}$.

## A. Diagonal preconditioning

If $A$ has entries with very different magnitudes, diagonal preconditioning might improve its condition. When $A$ is diagonally dominant and non-singular, we define:

$$D = \operatorname{diag}(d_1, \dots, d_n), \quad \text{with} \quad d_j = 1/|A_{jj}|^{1/2} \qquad (29)$$

Instead of solving $A\mathbf{x} = \mathbf{b}$, we solve $DAD\mathbf{y} = Db\mathbf{y}$, with $\mathbf{x} = D\mathbf{y}$ and $DAD$ still diagonally dominant and non-singular with all entries $\le 1$ in magnitude.

A matrix $(n \times n)$ $A$ is diagonally dominant when:

$$|A_{ii}| > \sum_{\substack{j=1 \\ j \neq i}}^{n} |A_{ij}|, \forall i = 1, \dots, n \qquad (30)$$

More generally, if $A$ is not diagonally dominant and possibly singular, it is possible to avoid division-by-zero errors by choosing a parameter $\delta > 0$ to redefine a $\delta$ dependent $d_j(\delta)$ such that:

$$d_j(\delta) = 1 / \max\{\delta, \sqrt{|A_{jj}|}, \max_{i \neq j} |A_{ij}|\}, \quad j = 1, \dots, n. \qquad (31)$$

**Examples [16]:**
1. If $A = \begin{bmatrix} -1 & 10^{-8} \\ 10^{-8} & 1 & 10^4 \\ 10^4 & 0 & 0 \end{bmatrix}$, then $\kappa(A) \approx 10^4$.

Let $\delta = 1$, $D = \begin{bmatrix} 1 & 10^{-2} & 1 \\ & 10^{-2} & 1 \\ & & 1 \end{bmatrix}$, in eq. 31.

Then $DAD = \begin{bmatrix} -1 & 10^{-10} & 1 \\ 10^{-10} & 10^{-4} & 1 \\ 1 & 0 & 0 \end{bmatrix}$, and $\kappa(DAD) \approx 1$.

2. $A = \begin{bmatrix} 10^{-4} & 10^{-8} & 1 \\ 10^{-8} & 10^{-4} & 10^{-8} \\ 10^{-8} & 0 & 0 \end{bmatrix}$ contains mostly very small entries, and $\kappa(A) \approx 10^{10}$.

Let $\delta = 10^{-8}$ and $D = \begin{bmatrix} 10^2 & 10^2 & 1 \\ 10^2 & 10^8 & 10^8 \end{bmatrix}$.

Then $DAD = \begin{bmatrix} 1 & 10^{-4} & 1 \\ 10^{-4} & 1 & 10^2 \\ 10^2 & 0 & 0 \end{bmatrix}$ and $\kappa(DAD) \approx 10^2$.

**Note:** The choice of $\delta$ makes a critical difference in this case [16], since $\delta = 1$ yields $D = I$.

## B. Binormalization

Golub et al. [18] scale a symmetric matrix by a series of $k$ diagonal matrices on both sides until all rows and columns of the scaled matrix have unit 2-norm:

$$DAD = D_k \cdots D_1 AD_1 \cdots D_k \quad (32)$$

### Example [16]

If $A = \begin{bmatrix} 10^{-8} & 1 \\ 1 & 10^{-8} & 10^4 \\ 10^4 & 0 \end{bmatrix}$, then $\kappa(A) \approx 10^{12}$.

With a single sweep of binormalization, we obtain $D = diag(8.1 \times 10^{-3}, 6.6 \times 10^{-5}, 1.5)$, and

$$DAD \approx \begin{bmatrix} 6.5 \times 10^{-1} & 5.3 \times 10^{-1} & 0 \\ 5.3 \times 10^{-1} & 0 & 1 \\ 0 & 1 & 0 \end{bmatrix} \quad \text{and} \quad \kappa(DAD) \approx 2.6 \text{ even though the rows and columns have not converged to}$$

one in the two-norm.

In contrast, diagonal preconditioning (eq. 31) defined by $\delta = 1$ and $D = diag(1, 10^{-4}, 10^{-4})$ reduces the condition number to approximately $10^4$.
### C. Singular linear systems preconditioning

Since preconditioning transforms the linear system into $M^{-1}Ax = M^{-1}b$ where $M^{-1}$ is an inverse conditioning matrix, $M$ should be non-nonsingular.

The non preconditioned solution solves $Ax \approx b$, but is not necessarily a minimum-length solution.

Let us illustrate the procedure with the following example:

$$A = \begin{bmatrix} 1 & 1 & 0 & 0 \\ 1 & 1 & 1 & 0 \\ 0 & 1 & 0 & 1 \\ 0 & 0 & 1 & 0 \end{bmatrix} \quad \text{and} \quad b = \begin{bmatrix} 6 \\ 9 \\ 6 \\ 3 \end{bmatrix}.$$

$rank(A) = 3$ and $Ax = b$ is a singular compatible system.

The minimum-length [16] solution is obtained from SVD decomposition of $A$ to yield $x^+ = A^+b = [2 \ 4 \ 3 \ 2]^T$.

With binormalization [18], we construct the matrix $D = diag([0.84201 \ 0.81228 \ 0.30957 \ 3.2303])$.

The minimum-length solution of the diagonally preconditioned problem: $DADy = Db$ is:

$$y^+ = \begin{bmatrix} 3.5739 & 3.6819 & 9.6909 & 0.93156 \end{bmatrix}^T$$

Then $x = Dy^+ = \begin{bmatrix} 3.0092 & 2.9908 & 3.0000 & 3.0092 \end{bmatrix}^T$ is a solution of $Ax = b$, however $x \neq x^+$.

## APPENDIX A: EXTREMUM PROBLEMS IN HIGHER DIMENSIONS

Finding extrema [2, 19, 20] of a multivariate function $f(x) \equiv f(x_1, x_2, \dots, x_n)$ amounts to finding zeroes $x^*$ of its partial derivatives $(\frac{\partial f}{\partial x_i})$. The extremum is a maximum or minimum depending on whether the Hessian matrix $H_f(x)$ (or local convexity) whose components are defined by the second derivatives: $(\frac{\partial^2 f}{\partial x_i \partial x_j})_{x^*}$ is negative or positive at point $x^*$.

In other words, the Hessian $H_f(x)$ is the Jacobian of partial derivatives $\frac{\partial}{\partial x_1}, \frac{\partial}{\partial x_2}, \dots, \frac{\partial}{\partial x_n}$ of $f$ with respect to $x_1, x_2, \dots, x_n$ so that

$$H_f(\mathbf{x}) = \begin{vmatrix} \frac{\partial^2 f}{\partial x_1^2} & \frac{\partial^2 f}{\partial x_1 \partial x_2} & \frac{\partial^2 f}{\partial x_1 \partial x_3} & \dots & \frac{\partial^2 f}{\partial x_1 \partial x_n} \\ \frac{\partial^2 f}{\partial x_2 \partial x_1} & \frac{\partial^2 f}{\partial x_2^2} & \frac{\partial^2 f}{\partial x_2 \partial x_3} & \dots & \frac{\partial^2 f}{\partial x_2 \partial x_n} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \frac{\partial^2 f}{\partial x_n \partial x_1} & \frac{\partial^2 f}{\partial x_n \partial x_2} & \frac{\partial^2 f}{\partial x_n \partial x_3} & \dots & \frac{\partial^2 f}{\partial x_n^2} \end{vmatrix}$$

### 1. Powell Method

Powell method is a procedure [2, 19] for carrying out minimization in $N$-dimensional space with no derivatives. It performs repeated application of a 1D optimization procedure. Thus Powell method is considered as a direction-set method for minimization [2] which means that if we start at a point $P$ in $N$-dimensional space, and proceed from there along a vector direction $n$, then any function of $N$ variables $f(P)$ can be minimized along the line $n$ by the 1D optimization method. Powell method belongs to a type of procedures that optimize $f(x)$ through iterative updating of a candidate vector $x^*$, where updates are given by $x^* \leftarrow x^* + \lambda^* e_j$, with $\lambda^*$ defined as a 1D optimum of $f(x^* + \lambda^* e_j)$. Vectors $\{e_j\}$ make a basis for $\mathbb{R}^d$, and cycles are performed through basis $j = 1, \dots, d$ as $x^*$ is being updated.

A major drawback occurs when $\{e_j\}$ basis is not well-aligned with function $f$ contours. The algorithm converges very slowly making erratic walks on $f$ surface instead of converging steadily as illustrated in Fig. 7.
<figure><img src="image_9.png" /><figcaption>Fig.7: Powell iterations from an initial point progressing steadily toward an extremum.</figcaption></figure>

The optimal basis for a given starting point $x^*$ is obtained from the eigenvectors [2, 19] of $H_f(x^*)$, where $H_f(u)$ is the Hessian of $f$ evaluated at point $u$, i.e. the $d \times d$ matrix with value $(\frac{\partial^2 f}{\partial x_i \partial x_j})_u$.

The set of directions is optimal when $f$ is quadratic, the exact minimum being reached in at most $d$ iterations [2, 19]. Since Powell method [2, 19] does not use first nor second derivatives of $f$, the optimal basis cannot be used. Powell method uses instead an approximation to the optimal basis updated after every full cycle of the algorithm. Let $x_0^1$ be an initial guess. The $i^{th}$ cycle of Powell method is as follows:

For $j = 1, \dots, d$:
$$ \lambda^* = \max_{\lambda} f(x_{j-1}^i + \lambda e_j). $$
If $j \le d$, $x_{j+1}^i = x_{j-1}^i + \lambda^* e_j$.
If $j < d$, set $e_j = e_{j+1}$.
If $j = d$, set $x_0^{i+1} = x_d^i$, $e_d = x_d^i - x_0^i$.

The main step is that basis vectors $e_j$ are re-written with the increment in $x$ for every cycle through a basis. This vector points in a direction where the function is known to increase substantially, so it should be mostly aligned with, rather than perpendicular to the contours of $f$ [2, 19]. Note that after every set of $d$ iterations of the above procedure, $\{e_j\}$ basis will have been completely re-written (this occurs after $d^2$ 1D optimizations are achieved). Moreover, it is important to orthonormalize $\{e_d, \dots, e_1\}$ before the next set of $d$ iterations begins. Moreover, orthonormalizing in reverse order, induces the next set to search along the direction found in the final cycle of the preceding one.

## 2. Broyden method

Broyden method is considered as a multidimensional version of the secant method [2].

Broyden update for solving a system of nonlinear equations is given by:

$$ F(x) = 0, $$

where $F: \mathbb{R}^n \to \mathbb{R}^n$ is differentiable. We denote its Jacobian as $J(x) \in \mathbb{R}^{n \times n}$. We need to iteratively update the approximate Jacobian $\{B_k\}$ in order to approach the true Jacobian matrix, while updating parameters $\{x_k\}$ as:

$$ x_{k+1} = x_k - B_k^{-1} F(x_k). $$

The alternative matrix $B_k$ should satisfy the secant condition, that is

$$ B_k u_k = y_k, \text{ with } u_k = x_{k+1} - x_k, y_{k+1} = F(x_{k+1}) - F(x_k). \quad (A1) $$

However, we could not obtain a definite solution with only eq. (A1) because it is an under-determined linear system with an infinite number of solutions, since there are less constraints than intrinsic parameters.

The most common procedure is to select a Jacobian approximation $B_{k+1}$ as close as possible to $B_k$ in the absolute norm sense, i.e. $||B_{k+1} - B_k||$ should be as small as possible. Broyden update [1, 2] defined as follows

$$ B_{k+1} = B_k + \frac{(y_k - B_k u_k) u_k^\top}{u_k^\top u_k}, \quad (A2) $$
is obtained by solving the following problem:

$$ \min_{\mathbf{B}} \|\mathbf{B} - \mathbf{B}_k\|_F \Rightarrow B\mathbf{u}_k = \mathbf{y}_k. $$

Updating the inverse of matrix $G_k := B_k^{-1}$ is as follows [1, 2]:

$$ G_{k+1} = G_k - \frac{(G_k \mathbf{y}_k - \mathbf{u}_k) \mathbf{u}_k^\top G_k}{\mathbf{u}_k^\top G_k \mathbf{y}_k}. $$

Broyden updating operations are done in $O(n^2)$ time [1] (rather than $O(n^3)$). For large problems, the complexity of matrix $B_k$ gets more involved as one proceeds, imposing memory limit to keep only a few recent updates to the initial (often diagonal) Jacobian approximation, might resolve this issue.

It is interesting to find that Broyden method converges in a superlinear fashion [2]. In addition, Broyden method works in over-determined (more constraints than parameters) and under-determined problems [1, 4] (less constraints than parameters) finding always an solution $\mathbf{x}_{opt}$ that is close to the ideal one $\mathbf{x}_{ideal}$ in the absolute norm sense, i.e. $||\mathbf{x}_{opt} - \mathbf{x}_{ideal}||$ is minimized.

[1] J. Stoer R. and Bulirsch, *Introduction to Numerical Analysis*, Second Edition, Springer-Verlag, (New York, 1992).

[2] W. H. Press, W. T. Vetterling, S. A. Teukolsky and B. P. Flannery, *Numerical Recipes in C: The Art of Scientific Computing*, Second Edition, Cambridge University Press (New-York, 1992).

[3] J. H. Wilkinson, *Rounding errors in algebraic processes*, Dover Publications, Inc., New York, (1994).

[4] J. R. Rice, *Numerical Methods, Software and Analysis*, Second Edition, Academic Press, New-York (1993).

[5] N. J. Higham, *N.J. Higham, Accuracy and Stability of Numerical Algorithms*, SIAM, Philadelphia (1996)

[6] J-F Muller, *Arithmétique des ordinateurs*, Chapter 2 p. 48, Masson (1989).

[7] C. M. Bender, S. A. Orszag, *Advanced Mathematical Methods for Scientists and Engineers*, McGraw-Hill, New York, (1978).

[8] R. Fitzpatrick, *Lecture Notes on Numerical Analysis*, University of Texas at Austin (2008).

[9] J. Westlake, *A Handbook of Numerical Matrix Inversion and Solution of Linear Equations*, John Wiley and Sons, New-York (1968)

[10] A. M. Turing, *Rounding-off Errors in Matrix Processes*, Quart. Jour. Mech. Applied Math., Vol 1, 287 (1948)

[11] C. D. Meyer, *Matrix Analysis and Applied Linear Algebra*, SIAM, Philadelphia (2000).

[12] W. K. Pratt, *Digital Image Processing*, Wiley, New-York (1991)

[13] R. C Gonzalez and R. E. Woods, *Digital Image Processing*, Third Edition, Pearson Prentice Hall, Upper Saddle River N. J. (2008)

[14] R. T. Gregory and D. L. Karney, *A Collection of Matrices for Testing Computational Algorithms*, Krieger Publishing Company Huntington, NY (1978)

[15] B. Dionne, *MAT3380 and graduate courses in Numerical Analysis*, Department of Mathematics and Statistics, University of Ottawa (2023).

[16] S.-C. T. Choi, C. C. Paige and M. A. Saunders, *MINRES-QLP: A Krylov subspace method for indefinite or singular symmetric systems*, SIAM Journal on Scientific Computing, Vol. 33, Issue 4 (2011).

[17] M. Benzi, *Preconditioning techniques for large linear systems: a survey*, J. Comput. Phys., 182(2), 418 (2002).

[18] O. E. Livne and G. H. Golub, *Scaling by binormalization*, Numer. Algorithms, 35(1), 97 (2004).

[19] D.E. Potter, *Computational physics*, John Wiley & Sons, London, New York (1973).

[20] M. R. Garey and D. S Johnson, *Computers and Intractability: A guide to the theory of NP completeness*, Bell Telephone Laboratories (1979).