Table of Contents¶
versioninfo()
Iterative Methods for Solving Linear Equations¶
Introduction¶
So far we have considered direct methods for solving linear equations.
- Direct methods (flops fixed a priori) vs iterative methods:
- Direct method (GE/LU, Cholesky, QR, SVD): good for dense, small to moderate sized A
- Iterative methods (Jacobi, Gauss-Seidal, SOR, conjugate-gradient, GMRES): good for large, sparse, structured linear system, parallel computing, warm start
PageRank problem¶
A∈{0,1}n×n the connectivity matrix of webpages with entries aij={1if page i links to page j0otherwise. n≈109 in May 2017.
ri=∑jaij is the out-degree of page i.
Larry Page imagines a random surfer wandering on internet according to following rules:
- From a page i with ri>0
- with probability p, (s)he randomly chooses a link on page i (uniformly) and follows that link to the next page
- with probability 1−p, (s)he randomly chooses one page from the set of all n pages (uniformly) and proceeds to that page
- From a page i with ri=0 (a dangling page), (s)he randomly chooses one page from the set of all n pages (uniformly) and proceeds to that page
- From a page i with ri>0
The process defines a Markov chain on the space of n pages. Stationary distribution of this Markov chain gives the ranks (probabilities) of each page.
Stationary distribution is the top left eigenvector of the transition matrix P corresponding to eigenvalue 1. Equivalently it can be cast as a linear equation. (I−PT)x=0.
Gene Golub: Largest matrix computation in world.
GE/LU will take 2×(109)3/3/1012≈6.66×1014 seconds ≈20 million years on a tera-flop supercomputer!
Iterative methods come to the rescue.
Jacobi method¶
Solve Ax=b.
Jacobi iteration: x(t+1)i=bi−∑i−1j=1aijx(t)j−∑nj=i+1aijx(t)jaii.
With splitting: A=L+D+U, Jacobi iteration can be written as Dx(t+1)=−(L+U)x(t)+b, i.e., x(t+1)=−D−1(L+U)x(t)+D−1b=−D−1Ax(t)+x(t)+D−1b.
One round costs 2n2 flops with an unstructured A. Gain over GE/LU if converges in o(n) iterations. Saving is huge for sparse or structured A. By structured, we mean the matrix-vector multiplication Av is fast (O(n) or O(nlogn)).
Gauss-Seidel method¶
Gauss-Seidel (GS) iteration: x(t+1)i=bi−∑i−1j=1aijx(t+1)j−∑nj=i+1aijx(t)jaii.
With splitting, (D+L)x(t+1)=−Ux(t)+b, i.e., x(t+1)=−(D+L)−1Ux(t)+(D+L)−1b.
GS converges for any x(0) for symmetric and pd A.
Convergence rate of Gauss-Seidel is the spectral radius of the (D+L)−1U.
Successive over-relaxation (SOR)¶
SOR: x(t+1)i=ω(bi−i−1∑j=1aijx(t+1)j−n∑j=i+1aijx(t)j)/aii+(1−ω)x(t)i, i.e., x(t+1)=(D+ωL)−1[(1−ω)D−ωU]x(t)+(D+ωL)−1(D+L)−1ωb.
Need to pick ω∈[0,1] beforehand, with the goal of improving convergence rate.
Conjugate gradient method¶
Conjugate gradient and its variants are the top-notch iterative methods for solving huge, structured linear systems.
A UCLA invention! Hestenes and Steifel in 60s.
Solving Ax=b is equivalent to minimizing the quadratic function 12xTAx−bTx.
Kershaw's results for a fusion problem.
| Method | Number of Iterations |
|---|---|
| Gauss Seidel | 208,000 |
| Block SOR methods | 765 |
| Incomplete Cholesky conjugate gradient | 25 |
MatrixDepot.jl¶
MatrixDepot.jl is an extensive collection of test matrices in Julia. After installation, we can check available test matrices by
using MatrixDepot
mdinfo()
# List matrices that are positive definite and sparse:
mdlist(:symmetric & :posdef & :sparse)
# Get help on Poisson matrix
mdinfo("poisson")
# Generate a Poisson matrix of dimension n=10
A = matrixdepot("poisson", 10)
using UnicodePlots
spy(A)
# Get help on Wathen matrix
mdinfo("wathen")
# Generate a Wathen matrix of dimension n=5
A = matrixdepot("wathen", 5)
spy(A)
Numerical examples¶
The IterativeSolvers.jl package implements most commonly used iterative solvers.
Generate test matrix¶
using BenchmarkTools, IterativeSolvers, LinearAlgebra, MatrixDepot, Random
Random.seed!(280)
n = 100
# Poisson matrix of dimension n^2=10000, pd and sparse
A = matrixdepot("poisson", n)
@show typeof(A)
# dense matrix representation of A
Afull = convert(Matrix, A)
@show typeof(Afull)
# sparsity level
count(!iszero, A) / length(A)
spy(A)
# storage difference
Base.summarysize(A), Base.summarysize(Afull)
Matrix-vector muliplication¶
# randomly generated vector of length n^2
b = randn(n^2)
# dense matrix-vector multiplication
@benchmark $Afull * $b
# sparse matrix-vector multiplication
@benchmark $A * $b
Dense solve via Cholesky¶
# record the Cholesky solution
xchol = cholesky(Afull) \ b;
# dense solve via Cholesky
@benchmark cholesky($Afull) \ $b
Jacobi solver¶
It seems that Jacobi solver doesn't give the correct answer.
xjacobi = jacobi(A, b)
@show norm(xjacobi - xchol)
Reading documentation we found that the default value of maxiter is 10. A couple trial runs shows that 30000 Jacobi iterations are enough to get an accurate solution.
xjacobi = jacobi(A, b, maxiter=30000)
@show norm(xjacobi - xchol)
@benchmark jacobi($A, $b, maxiter=30000)
Gauss-Seidal¶
# Gauss-Seidel solution is fairly close to Cholesky solution after 15000 iters
xgs = gauss_seidel(A, b, maxiter=15000)
@show norm(xgs - xchol)
@benchmark gauss_seidel($A, $b, maxiter=15000)
SOR¶
# symmetric SOR with ω=0.75
xsor = ssor(A, b, 0.75, maxiter=10000)
@show norm(xsor - xchol)
@benchmark sor($A, $b, 0.75, maxiter=10000)
Conjugate Gradient¶
# conjugate gradient
xcg = cg(A, b)
@show norm(xcg - xchol)
@benchmark cg($A, $b)