4.4. Iterative Methods¶
The general idea for iterative methods proceeds as follows:
We write a (matrix) equation:
\[x = Tx + c\]in such a way that this equation is equivalent to solving \(Ax=b\).
We start with an initial guess \(x^{(0)}\) for the solution of \(Ax=b\).
We iterate
\[x^{(k+1)} = Tx^{(k)} + c\]If properly designed the sequence \(x^{(0)},x^{(1)},\ldots,x^{(k)},\ldots\) converges to \(x^\star\), which satisfies \(x^\star = Tx^\star + c\) and consequently \(Ax^\star = b\).
4.4.1. The Jacobi Method¶
We decompose
where \(D\) is the diagonal part, \(L\) is the lower triangular part, and \(U\) is the upper triangular part. Note that all diagonal elements in \(L\) and \(U\) are zero. Substituting this expression in \(Ax=b\) gives us the following equation:
The above equation can be cast into the iteration \(\boxed{x^{(k+1)}=D^{-1}(L+U)x^{(k)}+D^{-1}b}\) or \(\boxed{Dx^{(k+1)}=(L+U)x^{(k)}+b}\)
Solution Easy, since we need to solve a linear system of equations with diagonal coefficient matrix
\[\begin{split}\left[ \begin{array}{cccc} d_1 & & & \\ & d_2 & & \\ & & \ddots & \\ & & & d_n \end{array} \right] \left[ \begin{array}{c} x_1^{(k+1)} \\ x_2^{(k+1)} \\ \vdots \\ x_n^{(k+1)} \end{array} \right]=\left[ \begin{array}{c} c_1 \\ c_2 \\ \vdots \\ c_n \end{array} \right] \Rightarrow x_i^{(k+1)}=\frac{c_i}{d_i}\end{split}\]Convergence The Jacobi method is guaranteed to converge when \(A\) is diagonally dominant by rows.
Complexity Each iteration has a cost associated with:
- Solving \(Dx^{(k+1)}=c\) which requires \(n\) divisions.
- Computing \(x=(L+U)x^{(k)}+b\) which requires as many additions and multiplications as non-zero entries in \(A\) (worst case \(O(n^2)\), but could be \(O(n)\) for sparse matrices).
Stopping criteria \(||b-Ax^{(k)}||<\varepsilon\) or \(||x^{(k+1)}-x^{(k)}||<\varepsilon\)
There are three forms of this algorithm we will see, for different purposes:
Matrix form (for proofs) \(Dx^{(k+1)}=(L+U)x^{(k)}+b\).
Algorithm form (without in-place update). Each row of \(Dx^{(k+1)}=(L+U)x^{(k)}+b\) can be written as:
\[\begin{split}a_{ii}x_i^{(k+1)} &=& \enspace b_i-\sum_{j\neq i} a_{ij}x_j^{(k)} \\ \Rightarrow x_i^{(k+1)} &=& \enspace \frac{1}{a_{ii}}\left(b_i-\sum_{j\neq i}a_{ij}x_j^{(k)}\right)\end{split}\]Note that memory requirements in this version are high, as every iteration \(k\) stores a new copy of the solution vector \(x^{(k)}\).
In-place algorithm, maintains two copies \(x, x^\textsf{new}\) of the solution vector, and uses the following update rule:
\[x_i^{\textsf{new}} \leftarrow \frac{1}{a_{ii}}\left(b_i-\sum_{j\neq i}a_{ij}x_j\right)\]
The code below implements the in-place algorithm mentioned in step 3) above.
import numpy as np
def Jacobi(A,b):
m = A.shape[0]
n = A.shape[1]
if(m!=n):
print 'Matrix is not square!'
return
# initialize x and x_new
x = np.zeros(m)
x_n = np.zeros(m)
# counter for number of iterations
iterations = 0
# perform Jacobi iterations until convergence
while True:
for i in range(0,m):
x_n[i] = b[i]/A[i,i]
sum = 0
for j in range(0,m):
if(j!=i): sum+=A[i,j]*x[j]
x_n[i] -=sum/A[i,i]
# stopping criterion
if(np.linalg.norm(x-x_n,2)<.000001): break
# copy x_new into x
for i in range(0,m):
x[i]=x_n[i]
# display updated solution
print 'x: ',
print x
iterations+=1
print 'Iterations: %d'%iterations
def main():
A = np.matrix([[4.,2.,-2.],[4.,9.,-3.],[-2.,-3.,7.]])
b = np.array([2.,8.,10.])
Jacobi(A,b)
if __name__ == "__main__":
main()
4.4.2. The Gauss-Seidel Method¶
We again employ the decomposition \(A = D-L-U\) for solving \(Ax=b\):
At this point, we place \(x^{(k+1)}\) on the left hand side and \(x^{(k)}\) on the right hand side
The benefit of the Gauss-Seidel method (1) over Jacobi is the improved convergence. In terms of complexity, each iteration of equation (1) amounts to solving a lower triangular system via forward substitution, i.e., it incurs a cost of \(O(k)\), where \(k\) is the number of non-zero entries in \(A\). In terms of update, equation (1) takes the following form:
Without in-place update.
\[\begin{split}x_i^{(k+1)} \leftarrow \frac{1}{a_{ii}}\left(b_i-\sum_{j<i}a_{ij}x_j^{(k+1)}-\sum_{j>i}a_{ij}x_j^{(k)}\right)\end{split}\]Once again, this version incurs a high memory overhead because a new solution vector \(x^{(k)}\) is created for every iteration \(k\).
In place update, which as before, maintains only two copies \(x,x^\textsf{new}\) of the solution vector.
\[\begin{split}x_i^{\textsf{new}} \leftarrow \frac{1}{a_{ii}}\left(b_i-\sum_{j<i}a_{ij}x_j^{\textsf{new}}-\sum_{j>i}a_{ij}x_j\right)\end{split}\]Compare the above in-place update with that for Jacobi.
The real difference in performance is that Gauss-Seidel is generally serial in nature (although parallel variants do exist), while Jacobi is highly parallel. The code below implements the in-place Gauss-Seidel update mentioned above.
import numpy as np
def Gauss_Seidel(A,b):
m = A.shape[0]
n = A.shape[1]
if(m!=n):
print 'Matrix is not square!'
return
# initialize x and x_new
x = np.zeros(m)
x_n = np.zeros(m)
# counter for number of iterations
iterations = 0
# perform Gauss-Seidel iterations until convergence
while True:
for i in range(0,m):
x_n[i] = b[i]/A[i,i]
sum = 0
for j in range(0,m):
if(j<i): sum+=A[i,j]*x_n[j]
if(j>i): sum+=A[i,j]*x[j]
x_n[i] -=sum/A[i,i]
# stopping criterion
if(np.linalg.norm(x-x_n,2)<.000001): break
# copy x_new into x
for i in range(0,m):
x[i]=x_n[i]
# display updated x
print 'x: ',
print x
iterations+=1
print 'Iterations: %d'%iterations
def main():
A = np.matrix([[4.,2.,-2.],[4.,9.,-3.],[-2.,-3.,7.]])
b = np.array([2.,8.,10.])
Gauss_Seidel(A,b)
if __name__ == "__main__":
main()