Appendix: Linear Solvers

In Section 3.6, we present an overview of solution methods for discretized systems of equations, followed by an introduction to Newton-based methods for solving nonlinear equations. Here, we review the solvers for linear systems required to solve for each step of Newton-based methods.^[1]

B.1 Systems of Linear Equations¶

If the equations are linear, they can be written as

where $A$ is a square ($n \times n$) matrix, and $b$ is a vector, and neither of these depends on $u$. If this system of equations has a unique solution, then the system and the matrix $A$ are nonsingular. This is equivalent to saying that $A$ has an inverse, $A^{-1}$. If $A^{-1}$ does not exist, the matrix and the system are singular.

A matrix $A$ is singular if its rows (or equivalently, its columns) are linearly dependent (i.e., if one of the rows can be written as a linear combination of the others).

If the matrix $A$ is nonsingular and we know its inverse $A^{-1}$, the solution of the linear system (Equation B.1) can be written as $x=A^{-1}b$. However, the numerical methods described here do not form $A^{-1}$. The main reason for this is that forming $A^{-1}$ is expensive: the computational cost is proportional to $n^3$.

For practical problems with large $n$, it is typical for the matrix $A$ to be sparse, that is, for most of its entries to be zeros. An entry $A_{ij}$ represents the interaction between variables $i$ and $j$. When solving differential equations on a discretized grid, for example, a given variable $i$ only interacts with variables $j$ in its vicinity in the grid. These interactions correspond to nonzero entries, whereas all other entries are zero. Sparse linear systems tend to have a number of nonzero terms that is proportional to $n$. This is in contrast with a dense matrix, which has $n^2$ nonzero entries. Solvers should take advantage of sparsity to remain efficient for large $n$.

We rewrite the linear system (Equation B.1) as a set of residuals,

To solve this system of equations, we can use either a direct method or an iterative method. We explain these briefly in the rest of this appendix, but we do not cover more advanced techniques that take advantage of sparsity.

B.2 Conditioning¶

The distinction between singular and nonsingular systems blurs once we have to deal with finite-precision arithmetic. Systems that are singular in the exact sense are ill-conditioned when a small change in the data (entries of $A$ or $b$) results in a large change in the solution. This large sensitivity of the solution to the problem parameters is an issue because the parameters themselves have finite precision. Then, any imprecision in these parameters can lead to significant errors in the solution, even if no errors are introduced in the numerical solution of the linear system.

The conditioning of a linear system can be quantified by the condition number of the matrix, which is defined as the scalar

where any matrix norm can be used.

Because $\|A\| \cdot \|A^{-1}\| \ge \|A A^{-1}\|$, we have

for all matrices. A matrix $A$ is well-conditioned if $\mathrm{cond}(A)$ is small and ill-conditioned if $\mathrm{cond}(A)$ is large.

B.3 Direct Methods¶

The standard way to solve linear systems of equations with a computer is Gaussian elimination, which in matrix form is equivalent to LU factorization. This is a factorization (or decomposition) of $A$, such as $A = L U$, where $L$ is a unit lower triangular matrix, and $U$ is an upper triangular matrix, as shown in Figure B.1.

The factorization transforms the matrix $A$ into an upper triangular matrix $U$ by introducing zeros below the diagonal, one column at a time, starting with the first one and progressing from left to right. This is done by subtracting multiples of each row from subsequent rows. These operations can be expressed as a sequence of multiplications with lower triangular matrices $L_i$,

After completing these operations, we have $U$, and we can find $L$ by computing $L=L_1^{-1} L_2^{-1} \cdots L_{n-1}^{-1}$.

Once we have this factorization, we have $LUu=b$. Setting $Uu$ to $y$, we can solve $Ly=b$ for $y$ by forward substitution. Now we have $Uu=y$, which we can solve by back substitution for $u$.

This process is not stable in general because roundoff errors are magnified in the backward substitution when diagonal elements of $A$ have a small magnitude. This issue is resolved by partial pivoting, which interchanges rows to obtain more favorable diagonal elements.

Cholesky factorization is an LU factorization specialized for the case where the matrix $A$ is symmetric and positive definite. In this case, pivoting is unnecessary because the Gaussian elimination is always stable for symmetric positive-definite matrices. The factorization can be written as

where $D = \text{diag} [U_{11}, \ldots, U_{nn}]$.

This can be expressed as the matrix product

where $G = L D^{1/2}$.

B.4 Iterative Methods¶

Although direct methods are usually more efficient and robust, iterative methods have several advantages:

• Iterative methods make it possible to trade between computational cost and precision because they can be stopped at any point and still yield an approximation of $u$. On the other hand, direct methods only get the solution at the end of the process with the final precision.

• Iterative methods have the advantage when a good guess for $u$ exists. This is often the case in optimization, where the $u$ from the previous optimization iteration can be used as the guess for the new evaluations (called a warm start).

• Iterative methods do not require forming and manipulating the matrix $A$, which can be computationally costly in terms of both time and memory. Instead, iterative methods require the computation of the residuals $r(u)=Au-b$ and, in the case of Krylov subspace methods, products of $A$ with a given vector. Therefore, iterative methods can be more efficient than direct methods for cases where $A$ is large and sparse. All that is needed is an efficient process to get the product of $A$ with a given vector, as shown in Figure B.2.

Iterative methods are divided into stationary methods (also known as fixed-point iteration methods) and Krylov subspace methods.

B.4.1 Jacobi, Gauss–Seidel, and SOR¶

Fixed-point methods generate a sequence of iterates $u_{1}, \ldots, u_{k}, \ldots$ using a function

starting from an initial guess $u_0$. The function $G(u)$ is devised such that the iterates converge to the solution $u^*$, which satisfies $r(u^*)=0$.

Many stationary methods can be derived by splitting the matrix such that $A=M-N$. Then, $Au=b$ leads to $Mu = Nu + b$, and substituting this into the linear system yields

Because $Nu = Mu - Au$, substituting this into the previous equation results in the iteration

Defining the residual at iteration $k$ as

we can write

The splitting matrix $M$ is fixed and constructed so that it is easy to invert. The closer $M^{-1}$ is to the inverse of $A$, the better the iterations work.

We now introduce three stationary methods corresponding to three different splitting matrices.

The Jacobi method consists of setting $M$ to be a diagonal matrix $D$, where the diagonal entries are those of $A$. Then,

In component form, this can be written as

Using this method, each component in $u_{k+1}$ is independent of each other at a given iteration; they only depend on the previous iteration values, $u_k$, and can therefore be done in parallel.

The Gauss–Seidel method is obtained by setting $M$ to be the lower triangular portion of $A$ and can be written as

where $E$ is the lower triangular matrix. Because of the triangular matrix structure, each component in $u_{k+1}$ is dependent on the previous elements in the vector, but the iteration can be performed in a single forward sweep. Writing this in component form yields

Unlike the Jacobi iterations, a Gauss–Seidel iteration cannot be performed in parallel because of the terms where $j < i$, which require the latest values. Instead, the states must be updated sequentially. However, the advantage of Gauss–Seidel is that it generally converges faster than Jacobi iterations.

The successive over-relaxation (SOR) method uses an update that is a weighted average of the Gauss–Seidel update and the previous iteration,

where $\omega$, the relaxation factor, is a scalar between 1 and 2. Setting $\omega=1$ yields the Gauss–Seidel method. SOR in component form is as follows:

With the correct value of $\omega$, SOR converges faster than Gauss–Seidel.

B.4.2 Conjugate Gradient Method¶

The conjugate gradient method applies to linear systems where $A$ is symmetric and positive definite. This method can be adapted to solve nonlinear minimization problems (see Section 4.4.2).

We want to solve a linear system (Equation B.2) iteratively. This means that at a given iteration $u_{k}$, the residual is not necessarily zero and can be written as

Solving this linear system is equivalent to minimizing the quadratic function

This is because the gradient of this function is

Thus, the gradient of the quadratic is the residual of the linear system,

We can express the path from any starting point to a solution $u^*$ as a sequence of $n$ steps with directions $p_{k}$ and length $\alpha_{k}$:

Substituting this into the quadratic (Equation B.20), we get

The conjugate gradient method uses a set of $n$ vectors $p_{k}$ that are conjugate with respect to matrix $A$. Such vectors have the following property:

Using this conjugacy property, the double-sum term can be simplified to a single sum,

Then, Equation B.24 simplifies to

Because each term in this sum involves only one direction $p_{k}$, we have reduced the original problem to a series of one-dimensional quadratic functions that can be minimized one at a time. Each one-dimensional problem corresponds to minimizing the quadratic with respect to the step length $\alpha_{k}$. Differentiating each term and setting it to zero yields the following:

Now, the question is: How do we find this set of directions? There are many sets of directions that satisfy conjugacy. For example, the eigenvectors of $A$ satisfy Equation B.25.^[2]However, it is costly to compute the eigenvectors of a matrix. We want a more convenient way to compute a sequence of conjugate vectors.

The conjugate gradient method sets the first direction to the steepest-descent direction of the quadratic at the first point. Because the gradient of the function is the residual of the linear system (Equation B.22), this first direction is obtained from the residual at the starting point,

Each subsequent direction is set to a new conjugate direction using the update

where $\beta$ is set such that $p_{k+1}$ and $p_{k}$ are conjugate with respect to $A$.

We can find the expression for $\beta$ by starting with the conjugacy property that we want to achieve,

Substituting the new direction $p_{k+1}$ with the update (Equation B.30), we get

Expanding the terms and solving for $\beta$, we get

For each search direction $p_{k}$, we can perform an exact line search by minimizing the quadratic analytically. The directional derivative of the quadratic at a point $x$ along the search direction $p$ is as follows:

By setting this derivative to zero, we can get the step size that minimizes the quadratic along the line to be

The numerator can be written as a function of the residual alone. Replacing $p_{k}$ with the conjugate direction update (Equation B.30), we get

Here we have used the property of the conjugate directions stating that the residual vector is orthogonal to all previous conjugate directions, so that ${r_{i}}^\intercal p_{i}$ for $i=0,1,\ldots,k-1$.^[3] Thus, we can now write,

The numerator of the expression for $\beta$ (Equation B.33) can also be written in terms of the residual alone. Using the expression for the residual (Equation B.19) and taking the difference between two subsequent residuals, we get

Using this result in the numerator of $\beta$ in Equation B.33, we can write

where we have used the property that the residual at any conjugate residual iteration is orthogonal to the residuals at all previous iterations, so ${r_{k+1}}^\intercal r_{k} = 0$.^[4]

Now, using this new numerator and using Equation B.37 to write the denominator as a function of the previous residual, we obtain

We use this result in the nonlinear conjugate gradient method for function minimization in Section 4.4.2.

The linear conjugate gradient steps are listed in Algorithm 2. The advantage of this method relative to the direct method is that $A$ does not need to be stored or given explicitly. Instead, we only need to provide a function that computes matrix-vector products with $A$. These products are required to compute residuals ($r = Au-b$) and the $Ap$ term in the computation of $\alpha$. Assuming a well-conditioned problem with good enough arithmetic precision, the algorithm should converge to the solution in $n$ steps.^[5]

B.4.3 Krylov Subspace Methods¶

Krylov subspace methods are a more general class of iterative methods.^[6] The conjugate gradient is a special case of a Krylov subspace method that applies only to symmetric positive-definite matrices. However, more general Krylov subspace methods, such as the generalized minimum residual (GMRES) method, do not have such restrictions on the matrix. Compared with stationary methods of Section B.4.1, Krylov methods have the advantage that they use information gathered throughout the iterations. Instead of using a fixed splitting matrix, Krylov methods effectively vary the splitting so that $M$ is changed at each iteration according to some criteria that use the information gathered so far. For this reason, Krylov methods are usually more efficient than stationary methods.

Like stationary iteration methods, Krylov methods do not require forming or storing $A$. Instead, the iterations require only matrix-vector products of the form $A v$, where $v$ is some vector given by the Krylov algorithm. The matrix-vector product could be given by a black box, as shown in Figure B.2.

For the linear conjugate gradient method (Section B.4.2), we found conjugate directions and minimized the residual of the linear system in a sequence of these directions.

Krylov subspace methods minimize the residual in a space,

where $x_0$ is the initial guess, and $\mathcal{K}_k$ is the Krylov subspace,

In other words, a Krylov subspace method seeks a solution that is a linear combination of the vectors $r_0, A r_0, \ldots, A^{k-1} r_0$. The definition of this particular sequence is convenient because these terms can be computed recursively with the matrix-vector product black box as $r_0, A(r_0), A(A(r_0)), A(A(A(r_0))), \ldots$ .

Under certain conditions, it can be shown that the solution of the linear system of size $n$ is contained in the subspace $\mathcal{K}_n$.

Krylov subspace methods (including the conjugate gradient method) converge much faster when using preconditioning. Instead of solving $Ax=b$, we solve

where $M$ is the preconditioning matrix (or simply preconditioner). The matrix $M$ should be similar to $A$ and correspond to a linear system that is easier to solve.

The inverse, $M^{-1}$, should be available explicitly, and we do not need an explicit form for $M$. The matrix resulting from the product $M^{-1} A$ should have a smaller condition number so that the new linear system is better conditioned.

In the extreme case where $M=A$, that means we have computed the inverse of $A$, and we can get $x$ explicitly. In another extreme, $M$ could be a diagonal matrix with the diagonal elements of $A$, which would scale $A$ such that the diagonal elements are 1.^[7]Krylov subspace solvers require three main components: (1) an orthogonal basis for the Krylov subspace, (2) an optimal property that determines the solution within the subspace, and (3) an effective preconditioner. Various Krylov subspace methods are possible, depending on the choice for each of these three components. One of the most popular Krylov subspace methods is the GMRES.4^[8]“‘