6 Generalized Jacobi Methods

In this section, we will develop a family of related methods for iteratively estimating a solution to a system of linear equations

Given a system like \[ Ax = b \] we prefer not to invert \(A\) to find \(x\), a difficult task. What our methods will do is take an initial guess \(x'\) and then transform it in a way to produce a new guess \(x''\) which is closer to \(x\). If we repeat the process over and over, then we will can construct a sequence \(x_0,x_1,...\) whose limit is \(x\).

Throughout this section we will state various results of the form “if \(||T||\) is small, then some limit converges”, where the norm is given. From our careful study of norms and spectra and the equivalence of norms, all of these results can be improved to take the form “if \(\rho(T)\) is small, then …” because if \(\rho(T)\) is small, there exists a norm for which the claim holds, and we know that a limit which converges for one norm converges to the same value for all norms.

Here, we take a more cautious form below because we only proved half of the norm equivalence theorem, and so that we can talk about specific errors or estimates (since that usually means we have a norm in mind).

6.1 Neumann Series

Our hope is to invert a matrix without using matrix inversion, just multiplication and addition. This sounds like an odd task, but it is actually something we know how to do from calculus.

Namely, the geometric series tells us that \[ \frac 1 {1-x} = \sum_{k=0}^\infty x^k,\] as long as \(|x| < 1\). So we can invert something \(1-x\) as long as \(x\) is small.

Remarkably, the geometric series still works for linear maps!

Theorem 6.1 Let \(T: V\to V\) be a linear map (or matrix) from a vector space to itself. If \(||T|| < 1\) for some induced norm, then the following limit converges: \[\sum_{k=1}^\infty T^k.\] Moreover, its limit is \((I-T)^{-1}\).

Proof. First, lt’s show that \(I-T\) is invertible. If not, there is some nonzero \(v\in V\) such that \((I-T)v = 0\). Rearranged, we get \(Tv = v\), and so \(|Tv|/|v| = 1\). But that ratio appears in the set whose supremum is \(||T||\), and so that would imply \(||T|| \geq 1\), a contradiction.

We can now prove the series converges with almost the same argument one uses in a calculus course (see HW1). Namely, the partial sums

\[S_n = \sum_{k=0}^n T^k\] satisfy the identity

\[(I-T)S_n = I - T^{n+1}.\]

By submultiplicativity, we know \(||T^{n+1}|| \leq ||T||^{n+1}\), and powers of a number less than \(1\) go to zero. So, the companion sequence \((I-T)S_n\) converges to the identity. By continuity and the invertibility of \((I-T)\), the sequence \(S_n\) converges, and its limit must be \((I-T)^{-1}\).

A series expression naturally gives us a way to estimate the inverse – just truncate the series at a large value of \(n\)! We usually apply this to \(I-T\): the result tells us that matrices which are close to the identity can be inverted by a series.

Corollary 6.1 Suppose \(T\) is such that \(||I-T|| < 1\). Then \(T\) is invertible, and \[T^{-1}= \sum_{k=0}^\infty (I-T)^k.\]

In fact, this can be improved slightly with some algebraic tricks. Namely, if \(||I-BA|| < 1\) then

\[\begin{align*} A^{-1} &= \sum_{k=0}^{\infty}(I - BA)^{k}B,\\ B^{-1} &= \sum_{k=0}^{\infty}A(I - BA)^{k}. \end{align*}\]

Proof. The first claims (about \(T\)) follow immediately from replacing \(T\) by \(I-T\) in Theorem 6.1. As for the second, note \(A^{-1}= (BA)^{-1}B\) and \(B^{-1}= A (BA)^{-1}\) and use the series to calculate the inverse of \(BA\). Since matrix multiplication is continuous, we can move it into the limit.

The second form is more natural than it appears. If we want to solve \(Ax=b\) without inverting \(A\), we may still carry out some partial simplifications row reductions. These can be represented by left multiplication, leading to a system like \(BA\). It stands to reason that \(BA\) is closer to the identity than \(A\), so \(||I-BA||\) could be small enough, even when \(||I-A||\) is not.

The method itself is straightforward to implement in code. A better implementation would save \((I-A)^k\) for reuse at the next step.

def NS(norm,A,b,steps):
    if norm(1-A) >= 1:
        raise ValueError('The norm of $I-A$ must be less than $1$ for the Neumann series to converge')
    k = 0
    x = b
    while k < steps:
        x = x + (1-A)^k*b
    return x

Now that we have a numerical method, we would like to undersand the error – how much accuracy do you gain at each step? Here, it is controlled in a fairly straightforward way.

Proposition 6.1 Suppose \(||I-A|| = \delta < 1\). Form the sequence \[x_n = \sum_{k=0}^n (I-A)^k b.\]

The error of the \(n\)th approximation is bounded: \[|x-x_n| < \delta^{n+1} |x|,\] and also by \[|x-x_n| < \delta^{n} ||A^{-1}- I|| |b|.\]

In particular, the error improves by a factor of \(\delta\) at each step.

Proof. Well, the error is \[\left|\sum_{k=n+1}^\infty (I-A)^k b\right| = \left|(I-A)^{n+1} \sum_{k=0}^\infty (I-A)^k b\right.\]

By submultiplicativity and upon replacing the series with \(x\), we obtain \[|x-x_n| < \delta^{n} |(I-A)x|.\]

We could apply submultiplicativity again to get the first error bound. Or we could note \(x=A^{-1}b\) and then bound \[|(I-A)x| = |(I-A)A^{-1}b| = |(A^{-1}- I) b| \leq ||A^{-1}- I|| |b|.\]

The second error estimate is slightly weaker in general, but has the advantage of being independent of the solution \(x\).

6.2 A Faster Method

The estimates from Neumann series seem straightforward to calculate – and they are. However, it is possible to speed up the calculation. This technique dates back to at least the Han dynasty around 200BC, though it is often credited to an English mathematician who rediscovered it about two millenia later. It was originally used to estimate roots of polynomials. In fact, the synthetic division you may have learned a long, long time ago is a way of carrying it out by hand.

Consider a polynomials like

\[f(x) = a_0 + a_1 x + a_2 x^2 +... + a_n x^n,\]

which we typically evaluate by calculating all the \(x_k\), multiplying by the \(a_i\), then summing. This takes about \(n(n+1)/2\) operations.

But we can rearrange it into a form which requires fewer operations, \[f(x) = a_0 + x (a_1 + x(a_2 + ... x(a_{n-1}+xa_n))).\] There are \(n\) multiplications and \(n\) additions, far fewer than the usual approach!! This works just as well for matrices.

In our situation, all the coefficients are \(1\). If we let \(f_n(x) = 1 + x + ... + x^n\), then the recursive shape relates them as \[f_{n+1}(x) = 1 + x f_n(x).\]

In our matrix situation, this turns into \[x_{n+1} = f_{n+1}(I-T)b = b + (I-T)f_n(I-T)b = b + (I-T) x_n.\]

This recursion is the standard form of the method.

In fact, the initial \(x_0\) does not matter for this form. Assuming that \(x_n \to x\) and take the limit of both sides: \[x = b + (I-T)x\] which is readily rearranged into \[Tx = b.\]

6.3 Splitting Matrices and Variations

We’ll conclude by using Corollary 6.1 to produce some more series for numerically solving \(Ax = b\), covering situations where \(||I-A||\) is large. These include the Richardson, Jacobi, and Gauss-Siedel methods, and their natural generalizations.

Given a linear system \(Ax = b\), and any other linear map or matrix \(B\). Combining our speed improvements with Corollary 6.1 leads us to the recursion \[x_{n+1} = Bb + (I-BA) x_n,\] which converges to \(x\) as \(n\to \infty\).

How should one choose \(B\)? Well, in practice one finds a matrix \(Q\) which is easy to invert, and sets \(B=Q^{-1}\). As discussed in the motivation before Corollary 6.1, one way to find \(B\) is to partially invert \(A\), using row operations, in which case it is naturally an inverse.

This tweak lets us rewrite the recursion as follows: \[Qx_{n+1} = b + (Q - A) x_n.\]

This has the advantage of not using an inverse at all. One can evaluate the right hand side directly. If \(Q\) is upper or lower triangular, it is easy to solve the resulting system without calculating an inverse.

The matrix \(Q\) is called a “splitting matrix”. The reason for the name is that \(Q\) is often taken to be a piece of the matrix \(A\), and so it appears from “splitting” the matrix \(A\) into \(A=Q+E\).

There are three main variations, each corresponding to a different choice of \(Q\).

Let \(Q=I\), recovering the original approach. This is sometimes called Richardson’s method.
Let \(Q\) be the diagonal of \(A\). This is called the Jacobi method. This is a natural choice for diagonally dominant matrices, for example.
Let \(Q\) be the upper (or lower) triangular part of \(A\). This is called the Gauss-Seidel method.

6.4 Exercises

Exercise 6.1 Apply iterative methods to estimate a solution to \[\begin{bmatrix} 4 & 5 \\ 3 & 5 \end{bmatrix} x = \begin{bmatrix} 2\\3 \end{bmatrix}.\]

Pick any admissible splitting matrix besides \(A\), and any initial guess besides the actual solution, and determine the result after two iterations.
Determine the quantity \(\delta = ||I-Q^{-1}A||\).
Determine the actual solution by inverting the matrix.
Compare the actual solution to your approximate solution from (a) using the \(\infty\) norm. Then, compare to the error estimate theorem from class. What do you notice?

Exercise 6.2 Suppose \(||I-Q^{-1}A|| = \delta < 1\), so \(Q\) is a good choice of splitting matrix for our iterative method. Show that the errors satisfy the following recursion

\[|x-x_k| \leq \frac{\delta}{1-\delta} |x_k - x_{k-1}|.\]

Exercise 6.3 Suppose \(||I-Q^{-1}A|| < 1\). Directly verify that the series given by \(x_0 = b\) and the recursion \[Qx_{n+1} = b + (Q - A) x_n\] converges to a solution to \(Ax = b\).

Exercise 6.4 Improve Exercise 6.3 by showing that

The recursion converges for any starting vector \(x_0\).
Its limit is \(b\) for any initial \(x_0\).

Exercise 6.5 Show that the Jacobi method always converges for diagonally dominant matrices by using the induced \(L_\infty\) norm.

Exercise 6.6 Prove that the set of invertible matrices is open: if \(A\) is invertible, then there is some \(\epsilon > 0\) such that \(B\) is invertible when \(||B-A|| < \varepsilon\).

Exercise 6.7 In contrast to Exercise 6.6, show that the set of diagonalizable matrices isn’t open: for all diagonalizable \(A\) and for all \(\varepsilon > 0\), there is some matrix \(B\) that \(||B-A|| < \varepsilon\) but \(B\) is not diagonalizable.

6.5 Programming Problems

Exercise 6.8 Implement the faster version of the Neumann series (Richardson) method.

Exercise 6.9 Implement our general iterative method for solving \(Ax = b\). Your function should take as input a norm, matrix \(A\), vector \(b\), guess \(x_0\), and splitting matrix \(Q\). It should raise an error if the norm condition is not satisfied.