Incremental Jacobi method

A well known iterative method for solving linear equations is easily derived by examining each of the $n$ equations in the linear system $A.w = b$ in isolation. If in the $i$'th equation $$ \sum_{j=1}^n a_{i,j}.w_j = b_i $$ we solve for the value of $w_i$, while assuming the other entries of $w$ remain fixed, we obtain: $$ w_i = ( b_i - \sum_{j \ne i} a_{i,j}.w_j ) / a_{i,i} $$ This suggests an iterative method defined by: $$ w_i^{(k)} = b_i - \sum_{j \ne i} a_{i,j}/a_{i,i} . w_j^{(k-1)} $$ Equivalently be written as: $$ w_i^{(k)} = w_i^{(k-1)} + b_i - \sum_{j=1}^n a_{i,j}/a_{i,i} . w_j^{(k-1)} $$ which is the Jacobi method. It will be assumed in the sequel that the equations system $A$ is always normed, which means (in Pascal):

   for i := 1 to N do
     d := 1/a[i,i];
   for j := 1 to N do
     a[i,j] := a[i,j] * d;

A necessary condition being that the main diagonal of $A$ is always non-zero. The main diagonal of the normed equations, hence, will be unity everywhere. The pseudo-code of the Jacobi method is then given by:

 Initialize:
   w := 0
 Iterations:
   w := w + (b - A.w)

Alternatively, an iterative method for solving the equations system $A.w = b$ may be devised by considering the fact that $A$ can be written as $I - M$ and therefore we can probably use the sum of the Geometric Series: $$ A^{-1} = \frac{I}{I - M} = I + M + M^2 + M^3 + M^4 + M^5 + ... $$ Herewith it is assumed that the matrices $M^n$ tend to become zero for large values of $n$. If such is the case, then the solution $w$ can be computed by: $$ w = A^{-1}.b = \frac{I}{I - M}.b = (I + M + M^2 + M^3 + M^4 + M^5 + ...).b $$ The (pseudo)code of the above method, substituting $M=(I-A)$, is given by:

 Initialize:
   r := b
   w := r
 Iterations:
   r := r - A.r
   w := w + r

The difference with the standard Jacobi method is that iterations are performed now upon the residual $r$ instead of the unknown solution $w$. Indeed $r$ is the solution of $A.r = b$ for $b = 0$, hence the iterations are according to standard Jacobi on $r$. Next, the solution is incremented with the new residual found, and the whole is assumed to converge. Hence it's quite sensible to refer to such an iterative process as an incremental Jacobi method.
Let's formulate the requirement that $M^n$ tends to become zero more precisely now. Any matrix $M$ can be written as: $$ M = U^{-1}.\Lambda.U \quad \mbox{where} \quad \Lambda = \left[ \begin{array}{ccccc} \lambda_1 & & & & \\ & \lambda_2 & & & \\ & & \lambda_3 & & \\ & & & . & \\ & & & & \lambda_n \end{array} \right] $$ It follows that: $$ M^k = (U^{-1}.\Lambda.U) (U^{-1}.\Lambda.U) (U^{-1}.\Lambda.U) ... (U^{-1}.\Lambda.U) = U^{-1}.\Lambda^k.U $$ Where: $$ \Lambda^k = \left[ \begin{array}{ccccc} \lambda_1^k & & & & \\ & \lambda_2^k & & & \\ & & \lambda_3^k & & \\ & & & . & \\ & & & & \lambda_n^k \end{array} \right] $$ If the absolute values of all eigenvalues are smaller than $1$, then: $$ \lim_{k \rightarrow \infty} \lambda_i^k = 0 $$ Then the transformation matrix $U$ ant its inverse are multiplied by the zero matrix. Consequently also the result $M^k$ must be zero. It thus seems that a necessary and sufficient condition for the iteration process to converge is that the absolute values of all eigenvalues of the matrix $M$ be smaller than $1$.
Convergence of the method can be judged further by Gersgorin's Circle Theorem, which states the following. Define the radius $R_i$ of the matrix-row $(i)$ by: $$ R_i = \sum_{j \ne i} | m_{i,j} | $$ Then each (complex) Eigenvalue $\lambda$ of the matrix $M$ is in at least one of the following disks in the complex plane: $$ \{ \lambda : | \lambda - m_{i,i} | \leq R_i \} $$ Now the diagonal elements of the iteration matrix $M$ are all equal to zero, because we assumed the matrix $A$ to be normalized and $M=I-A$ by definition. Furthermore, we find that all off-diagonal elements of $M$ are equal to $a_{i,j}/a_{i,i}$. This means that a sufficient condition for all eigenvalues of the iteration matrix be less than $1$ is: $$ | \lambda | \leq \sum_{j \ne i} | a_{i,j} | / | a_{i,i} | < 1 \qquad \Longleftrightarrow \qquad | a_{i,i} | > \sum_{j \ne i} | a_{i,j} | $$ Meaning that the original matrix $A$ should better be diagonally dominant.
With the incremental Jacobi method in mind, several variations on the theme can be easily thought of, such as an incremental Gauss-Seidel and an incremental Successive OverRelaxation (SOR) Method. Or an incremental Jacobi method with Preconditioning. The latter two actually have been implemented by the author, for a 3-D problem which is related to the Solar Wind.