This chapter is devoted to several iterative methods for solving the large sparse linear systems

\begin{equation} \label{EqPart4.1} {\bf A}\,{\bf x} = {\bf b} , \qquad {\bf A} \in \mathbb{C}^{n\times n} , \quad {\bf x}, \ {\bf b} \in \mathbb{C}^{n\times 1} . \end{equation}
As usual, we denote by ℂm×n the group of m-by-n matrices with complex coefficients; a similar notation is used for real-valued matrices, ℝm×n. Linear equations of the form \eqref{EqPart4.1} frequently arise from the solution of boundary-value problems for partial differential equations. In particylar, we conme across these linear systems when solving elliptic equations with either finite-difference approximations or with finite element approximations. Large sparse linear systems become apparent from many other practical problems too, of course, and the methods discussed here are important more generally.

Iterative methods are suited for use with sparse matrices. If A is dense, your best course of action is probably to factor A and solve the equation by backsubstitution. The time spent factoring a dense A is roughly equivalent to the time spent solving the system iteratively; and once A is factored, the system can be backsolved quickly for multiple values of b. Compare this dense matrix with a sparse matrix of larger size that fills the same amount of memory. The triangular factors of a sparse A usually have many more nonzero elements than A itself. Factoring may be impossible due to limited memory, and will be time-consuming as well; even the backsolving step may be slower than iterative solution. On the other hand, most iterative methods are memory-efficient and run quickly with sparse matrices.

Except when the matrix A in Eq.\eqref{EqPart4.1} has very special structure and fast direct methods apply, iterative methods are usually the method of choice for large sparse linear systems. All iterative methods can be divided into two categories: stationary methods and nonstationary iterative methods. Since the linear matrix/vector equarion \eqref{EqPart4.1} is not suitable for stationary iterative schemes, it requires preconditioning operation that transforms the given problem into fixed point problem

\[ {\bf x} = f({\bf x}) , \]
This equation is used to generate a stationary iterative scheme
\begin{equation} \label{EqPart4.2} {\bf x}^{(k+1)} = g \left( {\bf x}^{(k)} \right) , \qquad k= 0,1,2,\ldots , \end{equation}
which is actually a one-step recurvive method because every next approximation x(k+1) depends solely on one previous step, x(k). On the other hand, a nonstationary iterative algoritm evaluates every nest approximation using all previously obtained iterative values:
\begin{equation} \label{EqPart4.3} {\bf x}^{(k+1)} = g_k \left( {\bf x}^{(k)} , {\bf x}^{(k-1)} , \cdots , {\bf x}^{(0)} \right) , \qquad k= 0,1,2,\ldots . \end{equation}
For transferring the given linear system A x = b into fixed point form, it is convenient to split the coefficient matrix A into a sum of two matrices, one of which must be invertible,
\[ {\bf A} ={\bf S} + \left( {\bf A} - {\bf S} \right) , \qquad \det \mathbf{S} \ne 0 . \]
This allows us to rewrite the given equation (1) as
\[ {\bf A}\,{\bf x} = {\bf b} \qquad \Longrightarrow \qquad {\bf S}\,{\bf x} = {\bf b} + \left( {\bf S} - {\bf A} \right) {\bf x} . \]
Applying the inverse of S, we arrive at the equivalent equation
\[ {\bf x} = {\bf S}^{-1} {\bf b} + {\bf S}^{-1} \left( {\bf S} - {\bf A} \right) {\bf x} = {\bf S}^{-1} {\bf b} + \left( {\bf I} - {\bf S}^{-1} {\bf A} \right) , \]
where I is the identity matrix. This equation yields a classical stationary iterative scheme:
\begin{equation} \label{EqPart4.4} {\bf x}^{(k+1)} = {\bf S}^{-1} {\bf b} + \left( {\bf I} - {\bf S}^{-1} {\bf A} \right) {\bf x}^{(k)}, \qquad k=0,1,2,\ldots . \end{equation}

Stationary Iterative Methods

 

Theorem 1: Given vector b and splitting A = S + (A − S) of matrix A, with A and S nonsingular, the iteration scheme (4) converges if and only if the spectral radius \[ \rho \left( {\bf S}^{-1} {\bf M} \right) < 1 . \] Here ρ(X) is the largest modulus of eigenvalues of matrix X.

Proof of thsi main statement can be found in Darve & Wootters book.
Its statement follows from an important theorem of matrix theory:

Theorem 2: If H is a square matrix, then Hkx → 0 as k → ∞ for every vector x ∈ ℂn×1 if and only if each eigenvalue λ of H is less than one in modulus.

See Forsythe, G.E. and Wasow, W.R., Finite Difference Methods for Partial Differential Equations. New York: John Wiley & Sons, Inc., 1960.

Every square matrix H may be changed to a blockwise diagonal matrix S-1H S (the Jordan canonical form of H), each of whose diagonal blocks has the form \[ \begin{bmatrix} \lambda_i & 1 & 0 & \cdots & 0 \\ 0&\lambda_i & 1 & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0&0&0& \cdots & \lambda_i\end{bmatrix} . \] Then we have only to prove the theorem for an arbitrary block of the form above.

 

  1. Brigs, W.L., Multigrid Tutorial,
  2. Burden, R.L. and Faires, J.D., Numerical Analysis, nineth edition, Broks/Cole, Cengage, Boston, MA,
  3. Darve, E. and Wootters, M., Numerical Linear Algebra with Julia, SIAM, Philadelphia, 2021.
  4. Forsythe, G.E. and Wasow, W.R., Finite Difference Methods for Partial Differential Equations. New York: John Wiley & Sons, Inc., 1960.
  5. Loehr, N., Advanced Linear Algebra, CRC Press, 2014.
  8. Watkins, D.S., Fundamentals of Matrix Computations, Wiley; 3rd edition, 2010.
  9. White, R.E., Computational Linear Algebra with Applications and MATLAB Computations, CRC Press Boca Raton, FL, 2023. doi: 10.1201/9781003304128