Iterative methods for solving linear systems
Contents
10. Iterative methods for solving linear systems#
Consider the system of consistent linear equations:
Three iterative methods that we will use for solving the linear systems are the Jacobi method, the Gauss-Seidel Methods, and the Preconditioned Conjugate Gradient (PCCG) method.
These methods can be used to solve normal equations shown in the previous section and Mixed Model Equations (MME) we will covered later.
Consider MME, the left-hand-side (LHS) of the MME is represented by \(\mathbf{A}\) and the right-hand-side (RHS) by \(\mathbf{b}\).
The LHS of the MME is often too large to store in memory as a “fully-stored” matrix. However, \(\mathbf{A}\) is often very sparse. Thus, it is may be possible to store only the non-zero elements of \(\mathbf{A}\) and compute \(\mathbf{Ax_n}\), using sparse matrix methods.
10.1. The Jacobi method#
In the Jacobi method, the solution at iteration \(t+1\) can be written as:
where \(\mathbf{D}\) is the diagonal component of \(\mathbf{A}\).
10.1.1. Weighted Jacobi algorithm#
Convergence can often be improved by using a weighted Jacobi approach:
for \(0 < \alpha < 1\).
Straightforward application of the Jacobi method to solve the normal equations or mixed model equations would require first computing the left-hand-side (LHS) and right-hand-side (RHS) of the normal equations or mixed model equations, and then using equation (3) until convergence.
10.1.2. Numerical Example#
using LinearAlgebra, Distributions, DataFrames, Random
data = DataFrame(x=[1,1,2,2,2,2,3,3,4,1],y=[1.1,1.2,1.9,1.2,2.0,1.7,1.0,1.7,1.1,1.7])
n = size(data,1)
p = length(unique(data[:,:x]))
X = zeros(n,p);
for i = 1:n
j = data[:,:x][i]
X[i,j] = 1.0
end
X = [ones(n) X];
y = data[:,:y];
lhs = X'X;
rhs = X'y;
D = Diagonal(lhs)
A = lhs + I*0.0
b = rhs
x = ones(5)
for i = 1:20
x=inv(D)*(b-A*x)+x
println(x)
end
[0.45999999999999996, 0.33333333333333337, 0.7, 0.3500000000000001, 0.10000000000000009]
[1.0000000000000002, 0.8733333333333334, 1.24, 0.8900000000000001, 0.6400000000000001]
[0.45999999999999985, 0.33333333333333315, 0.6999999999999997, 0.34999999999999987, 0.09999999999999987]
[1.0000000000000002, 0.8733333333333334, 1.24, 0.8900000000000002, 0.6400000000000002]
[0.45999999999999985, 0.33333333333333315, 0.6999999999999997, 0.34999999999999976, 0.09999999999999976]
[1.0000000000000002, 0.8733333333333334, 1.24, 0.8900000000000002, 0.6400000000000002]
[0.45999999999999985, 0.33333333333333315, 0.6999999999999997, 0.34999999999999976, 0.09999999999999976]
[1.0000000000000002, 0.8733333333333334, 1.24, 0.8900000000000002, 0.6400000000000002]
[0.45999999999999985, 0.33333333333333315, 0.6999999999999997, 0.34999999999999976, 0.09999999999999976]
[1.0000000000000002, 0.8733333333333334, 1.24, 0.8900000000000002, 0.6400000000000002]
[0.45999999999999985, 0.33333333333333315, 0.6999999999999997, 0.34999999999999976, 0.09999999999999976]
[1.0000000000000002, 0.8733333333333334, 1.24, 0.8900000000000002, 0.6400000000000002]
[0.45999999999999985, 0.33333333333333315, 0.6999999999999997, 0.34999999999999976, 0.09999999999999976]
[1.0000000000000002, 0.8733333333333334, 1.24, 0.8900000000000002, 0.6400000000000002]
[0.45999999999999985, 0.33333333333333315, 0.6999999999999997, 0.34999999999999976, 0.09999999999999976]
[1.0000000000000002, 0.8733333333333334, 1.24, 0.8900000000000002, 0.6400000000000002]
[0.45999999999999985, 0.33333333333333315, 0.6999999999999997, 0.34999999999999976, 0.09999999999999976]
[1.0000000000000002, 0.8733333333333334, 1.24, 0.8900000000000002, 0.6400000000000002]
[0.45999999999999985, 0.33333333333333315, 0.6999999999999997, 0.34999999999999976, 0.09999999999999976]
[1.0000000000000002, 0.8733333333333334, 1.24, 0.8900000000000002, 0.6400000000000002]
[A*x b]
5×2 Matrix{Float64}:
20.0 14.6
5.62 4.0
8.96 6.8
3.78 2.7
1.64 1.1
10.1.3. Function for Weighed Jacobi Iteration#
function WJacobi(A,x,b,p=0.7;tolerance=0.000001,printout_frequency=10,maxiter=1000) #optional arguments keyword arguments
n = size(A,1) #number of linear equations
D = diag(A)
error = b - A*x
diff = sum(error.^2)/n
iter = 0
while (diff > tolerance) && (iter < maxiter)
iter += 1
error = b - A*x
x_temp = error./D + x
x = p*x_temp + (1-p)*x
diff = sum(error.^2)/n
if iter%printout_frequency == 0
println(iter," ",diff)
end
end
return x
end
WJacobi (generic function with 2 methods)
starting_value = zeros(size(lhs,2))
sol=WJacobi(lhs,starting_value,rhs,0.7)
[A*sol A*(lhs\rhs)]
10 3.8029144498098243e-6
5×2 Matrix{Float64}:
14.6006 14.6
4.00018 4.0
6.80024 6.8
2.70012 2.7
1.10006 1.1
10.1.4. Exercise 1#
Use the WJocobi function given above to get a solution to the system equation (1) for the Matrix \(\bo{A}\) and the vector \(\bo{b}\) of this numerical example.
10.1.5. Why does weighted Jacobi work better?#
Equation (1) can be rearranged as follows to get
where \(\bo{B} = -\bo{D}^{-1}\bo{E}\) and \(\bo{z} = \bo{D}^{-1}\bo{b}\). Functional iteration using Equation (5) results in the Jacobi method:
Let \(\bo{x}\) be a solution for equation (1). Then, equation (6) can be written as
where \(\bo{e}^{t} = \bo{x}^t - \bo{x}\). But, from equation (5), \(\bo{x} = \bo{B}\bo{x} + \bo{z}\), and therefore equation (7) can be written as
or
which implies
For the matrix \(\bo{B}\) in equation (9), let \(\bs{\Lambda}\) be a diagonal matrix with its eigen values and \(\bo{R}\) a matrix with the corresponding eigen vectors. Provided that the matrix \(\bo{R}\) is non-singular, it can be shown that \(\bo{B}^t = \bo{R}\bo{\Lambda}^{t}\bo{R}^{-1}\), and thus
10.1.5.1. Eigen values and Eigen vectors of \(\mathbf{B}\)#
Can see below that \(\bo{B}\) has two non-zero eigen values equal to 1 and -1. Thus, equation (11) does not converge.
A
5×5 Matrix{Float64}:
10.0 3.0 4.0 2.0 1.0
3.0 3.0 0.0 0.0 0.0
4.0 0.0 4.0 0.0 0.0
2.0 0.0 0.0 2.0 0.0
1.0 0.0 0.0 0.0 1.0
E = A - D
5×5 Matrix{Float64}:
0.0 3.0 4.0 2.0 1.0
3.0 0.0 0.0 0.0 0.0
4.0 0.0 0.0 0.0 0.0
2.0 0.0 0.0 0.0 0.0
1.0 0.0 0.0 0.0 0.0
E + D == A
true
B = -inv(D)*E
5×5 Matrix{Float64}:
-0.0 -0.3 -0.4 -0.2 -0.1
-1.0 -0.0 -0.0 -0.0 -0.0
-1.0 -0.0 -0.0 -0.0 -0.0
-1.0 -0.0 -0.0 -0.0 -0.0
-1.0 -0.0 -0.0 -0.0 -0.0
d,R = eigen(B)
Eigen{Float64, Float64, Matrix{Float64}, Vector{Float64}}
values:
5-element Vector{Float64}:
-1.0000000000000002
-0.0
-0.0
-0.0
1.0
vectors:
5×5 Matrix{Float64}:
-0.447214 -1.58794e-17 2.61582e-17 -4.38609e-18 0.447214
-0.447214 -0.621287 -0.440536 -0.171607 -0.447214
-0.447214 0.676061 -0.0910583 -0.0773749 -0.447214
-0.447214 -0.280128 0.888451 -0.0773749 -0.447214
-0.447214 -0.280128 -0.0910583 0.97907 -0.447214
Λ = diagm(0=>d)
5×5 Matrix{Float64}:
-1.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 1.0
Ri = inv(R)
5×5 Matrix{Float64}:
-1.11803 -0.33541 -0.447214 -0.223607 -0.111803
5.83939e-16 -0.580115 0.838839 -0.206979 -0.0517447
2.55472e-16 -0.566304 -0.202051 0.818869 -0.0505128
1.58186e-16 -0.525063 -0.187337 -0.187337 0.899737
1.11803 -0.33541 -0.447214 -0.223607 -0.111803
e0 = [1,0,0,0,0]
Ri*e0
5-element Vector{Float64}:
-1.1180339887498945
5.839389629374014e-16
2.5547195007681714e-16
1.5818570868256753e-16
1.1180339887498947
R*Λ
5×5 Matrix{Float64}:
0.447214 0.0 0.0 0.0 0.447214
0.447214 0.0 0.0 0.0 -0.447214
0.447214 0.0 0.0 0.0 -0.447214
0.447214 0.0 0.0 0.0 -0.447214
0.447214 0.0 0.0 0.0 -0.447214
R*Λ^2
5×5 Matrix{Float64}:
-0.447214 0.0 0.0 0.0 0.447214
-0.447214 0.0 0.0 0.0 -0.447214
-0.447214 0.0 0.0 0.0 -0.447214
-0.447214 0.0 0.0 0.0 -0.447214
-0.447214 0.0 0.0 0.0 -0.447214
R*Λ^100*Ri*e0
5-element Vector{Float64}:
1.000000000000011
1.0968138704326957e-14
1.0857116401864442e-14
1.0857116401864442e-14
1.0857116401864442e-14
R*Λ^101*Ri*e0
5-element Vector{Float64}:
-1.1301205611714504e-14
-1.0000000000000109
-1.0000000000000109
-1.0000000000000109
-1.0000000000000109
R*Λ^102*Ri*e0
5-element Vector{Float64}:
1.0000000000000113
1.1190183309251989e-14
1.1079161006789473e-14
1.1079161006789473e-14
1.1079161006789473e-14
10.1.5.2. Weighted Jacobi#
The weighted Jacobi algorithm can be obtained by writing \(\bo{x}\) as
and substituting equation (5) for the first term on right to get
Using this equation for functional iteration, the weighted Jacobi algorithm can be written as:
where \(\bo{C} = (1-\alpha)\bo{I} + \alpha\bo{B}\).
As reasoned before, equation (12) can be written as
or
where \(\bs{\Lambda}\) and \(\bo{R}\) contain the eigen values and eigen vectors of \(\bo{C}\).
Can see below that the largest eigen value of \(\bo{C}\) is 1.0, and the rest, in absolute value, are smaller than 1. This implies that the system of equations given by equation (1) is singular and has an infinite number of solutions.
10.1.6. Exercise#
Using the solution obtained above in equation 13, obtain two other solutions to equation (1).
10.1.6.1. Eigen values and Eigen vectors of \(\bo{C}\)#
α = 0.7
C = I*(1-α) + B*α
5×5 Matrix{Float64}:
0.3 -0.21 -0.28 -0.14 -0.07
-0.7 0.3 -0.0 -0.0 -0.0
-0.7 -0.0 0.3 -0.0 -0.0
-0.7 -0.0 -0.0 0.3 -0.0
-0.7 -0.0 -0.0 -0.0 0.3
d,R = eigen(C)
Λ = diagm(0=>d)
5×5 Matrix{Float64}:
-0.4 0.0 0.0 0.0 0.0
0.0 0.3 0.0 0.0 0.0
0.0 0.0 0.3 0.0 0.0
0.0 0.0 0.0 0.3 0.0
0.0 0.0 0.0 0.0 1.0
Ri = inv(R)
5×5 Matrix{Float64}:
-1.11803 -0.33541 -0.447214 -0.223607 -0.111803
-0.0 0.87178 -0.290593 -0.290593 -0.290593
0.0 0.0 0.881917 -0.440959 -0.440959
-0.0 0.0 0.0 -0.714143 0.714143
1.11803 -0.33541 -0.447214 -0.223607 -0.111803
e0 = [1,0,0,0,0]
h = Ri*e0
5-element Vector{Float64}:
-1.118033988749895
0.0
0.0
0.0
1.118033988749895
R*Λ^100
5×5 Matrix{Float64}:
-7.18645e-41 -2.00184e-69 9.98453e-70 1.05185e-70 0.447214
-7.18645e-41 4.13825e-53 -9.73972e-54 3.60836e-54 -0.447214
-7.18645e-41 -1.77354e-53 2.92192e-53 3.60836e-54 -0.447214
-7.18645e-41 -1.77354e-53 -2.92192e-53 -3.24753e-53 -0.447214
-7.18645e-41 -1.77354e-53 -2.92192e-53 3.9692e-53 -0.447214
R*Λ^100*Ri*e0
5-element Vector{Float64}:
0.4999999999999945
-0.4999999999999943
-0.4999999999999945
-0.4999999999999945
-0.4999999999999945
R*Λ^101*Ri*e0
5-element Vector{Float64}:
0.4999999999999945
-0.4999999999999943
-0.4999999999999945
-0.4999999999999945
-0.4999999999999945
R[:,1]*h[1]
5-element Vector{Float64}:
0.5
0.4999999999999999
0.5000000000000001
0.5000000000000001
0.5000000000000001
sol
5-element Vector{Float64}:
0.7300306184191997
0.6033641761387336
0.9700301932664003
0.6200308132809003
0.37003125614840027
e0 = [1,1,0,0,0]
sol1 = sol + R*Λ^101*Ri*e0
sol2 = sol1 + R*Λ^101*Ri*e0
[sol sol1 sol2]
5×3 Matrix{Float64}:
0.730031 1.08003 1.43003
0.603364 0.253364 -0.0966358
0.97003 0.62003 0.27003
0.620031 0.270031 -0.0799692
0.370031 0.0200313 -0.329969
[A*sol A*sol1 A*sol2]
5×3 Matrix{Float64}:
14.6006 14.6006 14.6006
4.00018 4.00018 4.00018
6.80024 6.80024 6.80024
2.70012 2.70012 2.70012
1.10006 1.10006 1.10006
10.2. The Gauss–Seidel method#
The Gauss–Seidel method can be obtained by first writing equation (1) as
where \(\mathbf{D}\) is a diagonal matrix with the diagonal elements of \(\mathbf{A}\), and \(\mathbf{L}\) and \(\mathbf{U}\) are matrices with the strictly lower and upper triangular components of \(\mathbf{A}\).
Functional iteration using equation (20) as:
is the Gauss-Seidel algorithm, which is equivalent to
The Gauss-Seidel algorithm given by equation (20) only requires the inverse of a diagonal matrix, and thus it is used for iteration.
10.2.1. Function for Gauss-Seidel Iteration#
function GaussSeidel(A,x,b;tolerance=0.000001,printout_frequency=10,maxiter=1000)
n = size(A,1)
for i = 1:n
x[i] = (b[i] - A[:,i]'x)/A[i,i] + x[i]
end
error = b - A*x
diff = sum(error.^2)/n
iter = 0
while (diff > tolerance) & (iter < maxiter)
iter += 1
for i = 1:n
x[i] = (b[i] - A[:,i]'x)/A[i,i] + x[i]
end
error = b - A*x
diff = sum(error.^2)/n
if iter%printout_frequency == 0
println(iter," ",diff)
end
end
return x
end
GaussSeidel (generic function with 1 method)
Random.seed!(31415)
n,p = 10_000, 1_000
M = randn(n,p)
A = M'M
b = M'randn(n);
x = zeros(p)
sol = GaussSeidel(A,x,b,printout_frequency=1);
1 19.959515013708828
2 0.8863762408459797
3 0.04852706759111071
4 0.0030786785399093204
5 0.00020239643667934267
6 1.3516227962181506e-5
7 9.224218116730581e-7
[A*sol b]
1000×2 Matrix{Float64}:
-2.35089 -2.34826
-62.8561 -62.8591
44.9825 44.9797
-47.9957 -47.9973
77.6082 77.6065
127.588 127.588
-4.43427 -4.43347
18.1429 18.1412
253.432 253.431
-159.59 -159.588
10.3759 10.3761
67.1881 67.1872
123.459 123.46
⋮
70.0611 70.0611
3.5686 3.5685
-145.141 -145.141
125.552 125.552
68.3659 68.3661
13.6079 13.6079
56.9186 56.9186
43.1212 43.1212
-14.3976 -14.3976
31.7366 31.7366
210.281 210.281
-58.7355 -58.7355
10.2.2. Weighted Gauss-Seidel#
function WGaussSeidel(A,x,b,p=0.7;tolerance=0.000001,printout_frequency=10,maxiter=1000)
n = size(A,1)
for i = 1:n
x[i] = (b[i] - A[:,i]'x)/A[i,i] + x[i]
end
error = b - A*x
diff = sum(error.^2)/n
iter = 0
while (diff > tolerance) & (iter < maxiter)
iter += 1
oldx = copy(x)
for i = 1:n
x[i] = (b[i] - A[:,i]'x)/A[i,i] + x[i]
end
x = p.*x + (1-p).*oldx
error = b - A*x
diff = sum(error.^2)/n
if iter%printout_frequency == 0
println(iter," ",diff)
end
end
return x
end
WGaussSeidel (generic function with 2 methods)
x = zeros(p)
sol = WGaussSeidel(A,x,b,1.1,printout_frequency=1);
1 18.684321462215486
2 0.5419988686885832
3 0.017502284939803422
4 0.0006220700405498619
5 1.6854226657005534e-5
6 6.264257528033032e-7
10.2.3. Convergence#
As we will see below, the second form given by equation (21) is useful to study the convergence properties of the algorithm. To study the convergence of Gauss-Seidel, equation (21) can be written as
for \(\bo{B} = -(\bo{D + L})^{-1}\bo{U}\) and \(\bo{z} = (\bo{D + L})^{-1}\bo{b}\). Given \(\bo{x}\) is a solution to equation (1) and \(\bo{e}^t = \bo{x}^t - \bo{x}\), equation (22) can be written as:
Can see below that the matrix \(\bo{B}\) for Gauss-Seidel in this example has only one non-zero eigen value, which is equal to 1.0. Thus, in this example, Gauss-Seidel will converge in one iteration.
D = Diagonal(A)
L = LowerTriangular(A-D)
U = UpperTriangular(A-D)
B = -inv(D+L)*U
d,R = eigen(B)
Λ = Diagonal(d)
1000×1000 Diagonal{ComplexF64, Vector{ComplexF64}}:
-0.079898-0.00141444im ⋅ … ⋅
⋅ -0.079898+0.00141444im ⋅
⋅ ⋅ ⋅
⋅ ⋅ ⋅
⋅ ⋅ ⋅
⋅ ⋅ … ⋅
⋅ ⋅ ⋅
⋅ ⋅ ⋅
⋅ ⋅ ⋅
⋅ ⋅ ⋅
⋅ ⋅ … ⋅
⋅ ⋅ ⋅
⋅ ⋅ ⋅
⋮ ⋱
⋅ ⋅ ⋅
⋅ ⋅ ⋅
⋅ ⋅ … ⋅
⋅ ⋅ ⋅
⋅ ⋅ ⋅
⋅ ⋅ ⋅
⋅ ⋅ ⋅
⋅ ⋅ … ⋅
⋅ ⋅ ⋅
⋅ ⋅ ⋅
⋅ ⋅ ⋅
⋅ ⋅ 0.269199+0.0im
10.2.4. Exercise 2#
Modify the Gauss-Seidel function given above for a weighted Gauss-Seidel algorithm. Compare convergence of the new algorithm with that given above for a range of weights between 0 and 1.