We want to solve the linear system given by \begin{equation} Ax = b, \end{equation} where $A$ is a non-singular $n\times n$ matrix.
You have already solved this problem using Gaussian Elimination (and its partial pivoted version) which has an assymptotic cost $\mathcal{O}(n^3)$.
In this homework you will try to solve the system in lower complexity using an iterative method. In particular, you will implement the Jacobi and Gauss-Seidel iterations and you will study its limitations.
Both Jacobi and Gauss-Iterations can be seen as fixed point methods, used by a matrix splitting.
Given a square matrix $A$ you will define a splitting in two matrices $N$ and $M$ such that $A = N+M$. If you suppose that $N$ is invertible you can write the system \begin{equation} Ax = b, \end{equation} as \begin{equation} Nx = b - Mx, \end{equation} and you can define the fixed point iteration as \begin{equation} Nx^{n+1} = b - Mx^{n}, \end{equation} or equivalently \begin{equation} x^{n+1} = N^{-1}\left ( b - Mx^{n} \right), \end{equation} in which the inverse of $N$ is never computed explicitily.
In this case you can show that the convergence speed is can be determined from the spectral radius of the iteration matrix: \begin{equation} T = N^{-1}M. \end{equation}
If we define \begin{equation} f(x) := N^{-1}\left ( b - Mx^{n} \right), \end{equation} then we are solving the fixed point problem \begin{equation} x = f(x). \end{equation}
We know that in this case, the fixed point iteration given by $x^{n+1} = f(x^n)$ converges if and only if there exist $\kappa <1$ such that \begin{equation} \| f(x) - f(y) \| \leq \kappa \|x - y\|, \qquad \forall x,y \in R^n, \end{equation} where we are using the Euclidean metric given by \begin{equation} \| x \| = \sqrt{\sum_{i = 1}^n |x_i|^2 }. \end{equation}
In such case, the error is given by \begin{equation} \| x^{n+1} - x\| \leq \frac{\kappa^{n+1}}{1 - \kappa} \| x^{0} - x\|, \end{equation} in other words, the convergence and its ratio ar given by $\kappa$.
Using the expression for $f$ we have that \begin{align} \| f(x) - f(y) \| & = \| N^{-1}(b - Mx) - N^{-1}(b - My) \|, \\ & = \| N^{-1} M(x -y) \|, \\ & \leq \|N^{-1} M \| \|x - y \|; \end{align} or \begin{equation} \kappa = \|N^{-1} M \| = \|T \|. \end{equation}
When we are using the Euclidean norm, we have that \begin{equation} \|T \| = \max_{z \neq 0} \frac{\|Tz \|}{\| z\|} = \rho(T). \end{equation} where the spectral radius $\rho(T)$ is given by \begin{equation} \rho(T) = \max_{i = 1,..,n} |\lambda_i|, \qquad \mbox{ where } \lambda_i \mbox{ are the eigenvalues of } T \end{equation}
The Jacobi iteration correspond to a fixed point method, in which the Matrix splitting is given by \begin{equation} N = D, \qquad M= L+U, \end{equation} where $D$ is the diagonal of $A$ and $L$ and $U$ are the (strictly) upper and lower triangular parts of $A$.
In order to extract the diagonal, upper triangular and lower triangular matrices from $A$ we will use the built-in functions in Julia.
In [3]:
# size of the matrices
m = 10
# generate a random matrix
A = rand(m,m) + m*eye(m)
# get the digonal part
D = diagm(diag(A),0)
# to get the upper triangular part of the matrix
U = triu(A,1)
#and the lower part
L = tril(A,-1)
# we check that this is indeed a matrix splitting
println("Error in the splitting = ",norm(A -(L+D+U)))
Q1.a Implement the Jacobi iteration with input:
Your function will have the following optional parameters
The ouput of your function is the final approximation $x$ of your vector, of in the case that history is true, it will output a matrix with the all the intermediate approximations of size $ n \times \mbox{number of iterations to converge}$.
Your function will raise an error if it didn't converge in the $Nmax$ iterations.
Hint:
In [4]:
function JacobiIt(A,b; Nmax = 30, ϵ=1e-6, history = false)
(n,m) = size(A)
n!=m && error("Dimension mistmach")
# initial guess
x0 = zeros(b)
history && (xHist = x0)
absB = norm(b)
# D
D = diagm(diag(A),0)
# L+U
M = triu(A,1) + tril(A,-1)
for i = 1:Nmax
x = D\(b - M*x0)
history && (xHist = hcat(xHist,x))
if norm(x-x0)/absB < ϵ
history ? (return xHist[:,2:end]) : (return x)
end
x0 = x
end
error("Tolerance not achieved within the maximum number of iterations")
end
Out[4]:
The Gauss-Seidel iteration is analogous to the Jacobi iteration; however, it uses a different splitting, namely \begin{equation} N = D+U, \qquad M= L, \end{equation}
Q1.a Implement the Gauss-Seidel iteration with input:
Your function will have the following optional parameters
The ouput of your function is the final approximation $x$ of your vector, of in the case that history is true, it will output a matrix with the all the intermediate approximations of size $ n \times \mbox{number of iterations to converge}$
Your function will raise an error if it didn't converge in the $Nmax$ iterations, and if the dimension of $A$ and $b$ are not consistent.
Hint:
In [5]:
function GaussSeidelIt(A,b; Nmax = 30, ϵ=1e-6, history = false)
x0 = zeros(b)
history && (xHist = x0)
absB = norm(b)
N = triu(A,0)
# to get the upper triangular part of the matrix
M = tril(A,-1)
for i = 1:Nmax
x = N\(b - M*x0)
history && (xHist = hcat(xHist,x))
if norm(x-x0)/absB < ϵ
history ? (return xHist[:,2:end]) : (return x)
end
x0 = x
end
error("Tolerance not achieved within the maximum number of iterations")
end
Out[5]:
In [6]:
m = 100
A = rand(m,m) + m*eye(m)
b = rand(m,1)
@time XJac = JacobiIt(A,b, history = true)
println("The residual is : ", norm(A*XJac[:,end]- b)/norm(b) )
println("number of iterations is : ",size(XJac,2))
@time XGS = GaussSeidelIt(A,b, history = true)
println("The residual is : ", norm(A*XGS[:,end]- b)/norm(b) )
println("number of iterations is : ",size(XGS,2))
Q3.a How can you explain the one converges faster than the other? Write a small script that computes the spectral radius of the iteration matrix for both algorithms and use it on your explanation.
In [7]:
# we compute the spectral radius of the iteration matrix for both method
# Jacobi
N = diagm(diag(A),0)
M = triu(A,1)+ tril(A,-1)
TJac = N\M;
λ = eigvals(TJac)
println("Spectral radius of the iteration of the Jacobi iteration is ", maximum(abs(λ)));
# Gauss Seidel
N = triu(A,0)
M = tril(A,-1)
TGS = N\M;
λ = eigvals(TGS)
println("Spectral radius of the iteration of the Jacobi iteration is ", maximum(abs(λ)));
Answer: The number of iteration to converge is determined by the spectral radius of the iteration matrix. A smaller radius will mean a faster convergence, as you can see above, the spectral radius in the case of the Gauss Seidel iteration is much smaller than the one for Jacobi.
Now you will use your algorithms to solve a slighly different system
In [8]:
m = 100;
A = rand(m,m) + sqrt(m)*eye(m);
b = rand(m,1);
With the Jacobi iteration,
In [9]:
@time XJac = JacobiIt(A,b, history = true)
println("number of iterations is : ",size(XJac,2))
and with the Gauss Seidel iteration
In [10]:
@time XGS = GaussSeidelIt(A,b,Nmax = 60, history = true)
println("number of iterations is : ",size(XGS,2))
Q3.b How can you explain the one converges and the other just fails? Write a small script that computes the spectral radius of the iteration matrix for both algorithms and use it on your explanation.
In [11]:
# we compute the spectral radius of the iteration matrix for both method
# Jacobi
N = diagm(diag(A),0)
M = triu(A,1)+ tril(A,-1)
TJac = N\M;
λ = eigvals(TJac)
println("Spectral radius of the iteration of the Jacobi iteration is ", maximum(abs(λ)));
# Gauss Seidel
N = triu(A,0)
M = tril(A,-1)
TGS = N\M;
λ = eigvals(TGS)
println("Spectral radius of the iteration of the Jacobi iteration is ", maximum(abs(λ)));
Answer: We can observe that Jacobi iteration diverges, this is because the spectral radius of the iteration matrix is greater than 1, meaning that it has a instable mode. On the other hand, the spectral radius of the iteration matrix for Gauss-Seidel is well below one, meaning that the method converges.
Q4.a what is the complexity of one iteration of the Gauss-Seidel method? and about Jacobi?
Answer:
Q4.b What would be the condition on the number of iterations in order such that the Jacobi and Gauss-Seidel iteration would have a better complexity that Gaussian Elimination
Answer:
Q4.c Run a benchmark with the randomly generated systems above. Can you achieve quadratic convergence?
In [136]:
# write your whole script here