Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It only takes a minute to sign up.
Sign up to join this community
Descent Methods, like Gradient and Conjugate Gradient ones, have a nice geometric interpretation and I really love them. What about Jacobi, Gauss-Seidel or other matrix splitting methods? I can't see any intuition behind the formulas.
A posteriori I know that these methods works because of Banach Fixed Point theorem but I still need a more stimulant picture.
Any ideas?
I'm not really sure what do you mean by intuition behind the methods, but all classical iterations for solving $Ax=b$ are based on splitting $$ A=M-N, $$ where
The equation $Ax=b$ is then equivalent to $$Mx=Nx+b=Mx+(b-Ax),$$ which "naturally" defines a fixed-point iteration $$ x_k=M^{-1}Nx_{k-1}+M^{-1}b=x_{k-1}+M^{-1}(b-Ax_{k-1}). $$ Note that the second expression is more suitable in practice as it does not involve the "residual" matrix $N=M-A$.
There are many ways how to choose $M$ (assume $A=D-L-U$ is a splitting of $A$ to diagonal and lower and upper triangular parts):
Splitting methods have a property that the errors are transformed in the same way independently of the iteration: $$ x-x_k=(I-M^{-1}A)(x-x_{k-1})=(I-M^{-1}A)^k(x-x_0). $$ From the spectral point of view, this has an advantage that (in contrast with CG or steepest descent) the spectral components of the initial error are damped with (asymptotically) the same rate from iteration to iteration. This has an important application, e.g., in multigrid methods, where methods like Jacobi and Gauss-Seidel are often used as smoothers since they efficiently damp (for PDE-related problems) the oscilatory error components associated with large eigenvalues.
Another particular motivation for Jacobi and Gauss-Seidel methods is that they zero out certain components of the residual vector $b-Ax$. Given an approximation $x_k$, the $i$th component of $x_{k+1}$ of the Jacobi iterate is computed simply by solving the $i$th equation of $Ax=b$: $$\tag{$*$} (x_{k+1})_i=\frac{1}{a_{ii}}\left(b_i-\sum_{j\neq i}a_{ij}(x_k)_j\right). $$ If $z_i=[(x_k)_1,\ldots,(x_k)_{i-1},(x_{k+1})_i,(x_k)_{i+1},\ldots,(x_k)_n]^T$, we have $(b-Az_i)_i=0$. If the Jacobi update is executed in a loop for $i=1,\ldots,n$, the components $(x_{k+1})_j$, $j=1,\ldots,i-1$, can be already used on the right-hand side of ($*$), which leads to the update $$\tag{$**$} (x_{k+1})_i=\frac{1}{a_{ii}}\left(b_i-\sum_{j<i}a_{ij}(x_{k+1})_j-\sum_{j>i}a_{ij}(x_k)_j\right). $$ This leads to the (forward) Gauss-Seidel iteration. Switching the loop order and the indices $k$ and $k+1$ on the right-hand side in ($**$) gives the backward Gauss-Seidel iteration. Again, we can interpret it as zeroing out a certain component of the residual vector.
These methods are highly intuitive and natural, but the idea isn't always explained clearly. Jacobi and Gauss-Seidel are the simplest things you would try if you were going to invent your own iterative method for solving $Ax = b$.
The first equation in $Ax = b$ is \begin{equation} A_{11} x_1 + A_{12} x_2 + \cdots + A_{1n} x_n = b_1. \end{equation} This equation suggests a way to update our estimate of $x_1$: \begin{equation} x_1^{k+1} = \frac{b_1 - A_{12} x_2^k - \cdots - A_{1n} x_n^k}{A_{11}}. \end{equation} Our estimates of $x_2,\ldots,x_n$ are updated in the same way, all updates performed in parallel. Then we've got the full vector $x^{k+1}$, and we're ready for the next iteration. That's the Jacobi method.
The Gauss-Seidel method is just like the Jacobi method, except that you update the variables one at a time (rather than in parallel), and during each update you use the most recent value for each variable.
As an afterthought you can express these methods as splitting methods using matrix notation. But that obscures how obvious they are to invent.
In order to see the idea of how the Jacobi Method works, the straightforward approach is to use the very useful Contraction Mapping principle, which is roughly, if there is a constant $0\le c<1$ such that $\|F(x_{k+1}) − F(x_{k}) \|\le c\|x_{k+1} – x_{k}\|$, then you can see with a little thought, that the iteration $x_{k+1} = F(x_{k})$ will converge to a fixed point.
Now would be a good time to review the Contraction Mapping principle, since iterative methods are an easy downhill slide from there.
Thus we need to convert $Ax=b$ into a form $x_{k+1} = F(x_{k})$. The most apparent conversion is to extract the diagonal elements of $A$, because a diagonal matrix has a simple inverse, \begin{align*}
Ax &= b, \\
(A - D + D)x &= b, \\
(A - D)x + Dx &= b, \\
Dx &= b - (A - D)x, \\
x &= D^{-1}b - D^{-1} (A - D)x,
\end{align*}
and iterate as,
$$x_{k+1} = D^{-1}b - D^{-1}(A - D) x_k$$
This is the Jacobi iterative method.
Importantly, notice above, that if you have a value for x that satisfies the relation Dx = b - (A - D)x, then x also satisfies Ax = b.
Now there is no guarantee that this will be a contraction mapping, but often it is. And it will contract when the diagonal elements are dominant.
What would motivate you to do this? An easy-peasy motivation starts with the rudimentary relaxation method. First assume that the diagonal elements of A are much larger than the off-diagonal elements of A. In this case the inverse of A is approximately the inverse of D, which is simple to calculate.
Compute a residual:
$$r = b - Ax_k$$
Compute an approximate correction:
$$Correction = D^{-1}r$$
Then,
\begin{align*} x_{k+1} &= correction + x_k, \\ x_{k+1} &= D^{-1}r + x_k, \\ x_{k+1} &= D^{-1}(b - Ax_k) + x_k \\ x_{k+1} &= D^{-1}b - D^{-1}(A - D)x_k, \\ \end{align*}
And this is the Jacobi method from the relaxation method. There you have it.
The idea is that you assign to each equation a dominant variable and solve the equation for this variable. The assignment should result in a bijective relation between equations and variables
Then use these partial solutions in some kind of fixed point scheme, the dominant ones are the Gauß-Seidel and Jacobi schemes. This works especially well if the original system matrix is a small perturbation of a diagonal matrix.
Closely related are sparse reordering methods for sparse matrices, these often also start with the formulation of an assignment problem equations to variables, and solving it with, e.g., the Hungarian algorithm.
Let me assume that $A$ is symmetric and invertible. Solving the linear system $Ax = b$ is equivalent to minimizing quadratic function $\frac{1}{2}x^T A x - b x$. Classical splitting methods such as Jacobi, Gauss-Seidel etc. define an invertible matrix $P$ and use the trivial identity $A = P + (A - P)$ to rewrite the linear system $Ax = b$ like this: $$ Px = (P-A)x + b \\ x = x + P^{-1}(b - Ax) \\ x^{k+1} := x^{k} + P^{-1}(b - Ax^{k}) \\ $$ where the last line is an update rule for the iterative method. This iterative method becomes gradient descent if you set $P$ to be the identity (or scaled identity for if you want to adjust the step size). For Jacobi, you set $P$ to the diagonal of $A$ and so on for other methods. Any setting of $P$ can be seen as a Quasi-Newton method -- the idea of Quasi-Newton methods being that we try to transform the gradient in attempt to get better convergence. Intuitively, $P$ should be an approximation of $A$ that is easily invertible. The gradient descent's $P = I$ is one extreme. The opposite extreme is full (i.e., non-Quasi) Newton's method, which corresponds to setting $P$ to $A$. The update rule then ends up saying: "if you can directly compute $A^{-1} b$ then you don't need to iterate at all."
A good approach for an intuitive understanding of the Jacobi method is to pretend that you are presented with a system of linear equations to solve, and you need to invent a way to solve them. The exercise below shows how you would naturally invent the general idea of using iterative solutions of linear equations, including the Jacobi method.
P.S. First recall the general concept of any iterative relaxation:
Step 1: Calculate a residual value that measures a difference between a current approximate solution and the true solution.
Step2: Perform an operation on the input variables that will reduce the residual value of a new approximate solution.
Step 3: Perform these steps until the residual is small enough.
A system of linear equations L1, L2, L3… has variables x1, x2, x3…, and is expressed in matrix form as Ax=b. The core idea developed below is solving a system of linear equations by inventing an iterative relaxation method of shrinking the length of a residual error vector toward zero (residual vector = Ax-b). Then the method is modified to perform a group of iterations in parallel in order to produce essentially the Jacobi method. The Jacobi method to solve linear equations is evolved by starting from a well-known, hill-climbing, multidimensional, gradient-free minimization method that shrinks the length of a residual error vector (residual vector = Ax-b) to zero, at which point the x vector then solves the linear equations.
For teaching purposes, we assume that the linear equations have strict diagonal dominance in order to ensure convergence below. Each linear equation L1, L2, L3… defines a “plane” (in 2D, 3D, 4D…) and is associated with the corresponding variable x1, x2, x3…. First, a trial point x1, x2, x3… is guessed. At each iteration, the method then cyclically adjusts each x1, x2, x3… in turn along the corresponding axis direction until the length of the residual vector is minimized in that direction, which needs numerous calculations of the length of the residual vector. Next, and importantly, you notice that when each x1, x2, x3… is adjusted, that the minimum value of the length of the residual vector always occurs when the trial point is adjusted to lie in the “plane” of the linear equation corresponding to the variable x1, x2, x3…. Because the equations are linear, this adjustment can be directly calculated. Thus the residual vector no longer needs to be calculated, and the calculation is removed, leaving only the iterative algorithm. This method of cyclically adjusting the trial point is seen to converge because it reduces the length of the residual vector at each iteration.
At this point, you essentially have the Gauss-Seidel method. Finally you notice that the adjustments above to x1, x2, x3… can be made independently and simultaneously, which is essentially the Jacobi method. When the algorithm is expressed algebraically in an iterative matrix form, it becomes $x^{k+1} = x^{k} - D^{-1}(Ax^k - b)$, where D is the diagonal matrix.
Backing up a bit, you see that changing x1 may reduce the residual element of the first equation L1, but also changes the residuals of the other equations L2, L3, L4…, sometimes increasing the individual residual element length and sometimes decreasing the individual residual element length, but we easily notice that the total decreases are always more than the total increases: the sum of the lengths (i.e., absolute value) of all the residual elements (the L1 metric) always decreases, and we realize that this means that eventually the residual vector will be reduced to zero. This is ensured because of diagonal dominance, as assumed above. Thus we will use the L1 metric for the length of the residual vector. In other words, you are adjusting x1 in order to reduce the total residual vector, not specifically to reduce the residual of the first equation L1, but because of diagonal dominance, it happens that the residual of the first linear equation will shrink the most. Similarly, adjusting x2 will reduce the residual of the second liner equation the most. It is old and well-known to use the L1 metric in minimization methods, as discussed in many references on minimization.
We see that $D^{-1}$ is used as an approximate inverse matrix by using the diagonal terms of A as an approximation of the original matrix. You can see that this is not a bad idea when the diagonal terms are strongly dominant. In other words, when the diagonal terms are dominant, then the system of linear equations is approximated by setting the off-diagonal terms to zero; that is, each linear equation is approximated by the single largest term. The idea of an approximate inverse matrix is another valuable idea from the Jacobi method.
This iterative process motivates the use of general iterative methods to solve linear equations.
For a more detailed discussion of an intuitive explanation of the Jacobi method, see the document at https://doi.org/10.6084/m9.figshare.5969263.v1
