[SOLVED] CS algorithm Computational Methods

Computational Methods
for Physicists and Materials Engineers
4
Systems of Linear Equations II Iterative Method

xi are unknowns to be determined Aij and bi are coefficients (knowns)
Matrix form:
System of linear equations
N linear equations: AxAxAxb
11 1 12 2 1N N 1 AxAxAxb
21 1 22 2 2N N 2
AxAxAxb N1 1 N2 2 NN N N
AAAxb 11 12 1N 1 1
21 22 2N 2 2 AAAxb
A A Ax b
N1 N2 NN N N A xb

https://www.cise.ufl.edu/research/sparse/matrices/
Rothberg/gearbox.html
Size of K is
153,746 153,746 = 23,637,832,516 9,080,404 nonzero elements
Iterative method: guess the solution and iteratively modify the guess
Iterative method Ax b
By direct method, we obtain the exact answer (if we dont consider rounding error) in the last step. The computational complexity for LU and QR is N3 if the coefficient matrix A is N N
For most of problems of practical interest, A is a largeN, sparse matrix E.g. FEM Ku = f
For aircraft flap actuator

Guess a solution x0 (probably its bad)
Residual
Correct the guess by residual
bAx0
x1 x0 C(bAx0)f(x0) Convergent? check error
How to correct guess s.t. error converges to 0?
If x* exists s.t. x* = f(x*)? (Fixed point theorem) How to choose f?
Correct the guess by residual
x2 x1 C(bAx1)f(x1)
Idea of iterative method Ax b
x1 x0
Definition of magnitude? Vector norm
How to estimate error?

Vector norm: size of a vector Norm x has the properties:
Positivity:
Definiteness:
Homogeneity:
Triangle inequality:
x 0
x = 0 if and only if x = 0 x = || x
x + y x + y
y x+y
forallx,yXandall
When x , norm is just absolute value |x|
x
2
x y x2 y2 2 xy x2 y2 2xy x y 2 x y x y
For example, lp norm:
N
p1/p
xx
p
i
i1

p i1 l norm is called maximum norm
i1,,N
Vector norm: size of a vector
2
i
i1
N p 1/p x x
p
i1 N
i

x1 x i1
i
21/2
x x xTx
N
l2 norm is called Euclidean norm because it is just the length of a vector
N
p1/p
xlim x maxx ii

2 x 3
i1
3 2 1/2
Vector norm: size of a vector
Equivalence of the norms
Theorem: on a finitedimensional vector space, all norms are equivalent. When one norm goes to 0, all the other norms go to 0
x1 xx
3
x1xi x1 x2 x3
x x x2 x2 x2 (lengthofx)
2
i 123 i1
x maxxmaxx,x,x i1,2,3 i 1 2 3
Onlywhenx1,x2,x3 0atthesametime,x1 0.Andatthesame time x2 and x 0

We can prove that
N
A 1 max A
Frobenius norm
AF A
2
Matrix norm: size of a matrix
Definition of matrix norm is based on vector norm
A sup Ax x 1
sup: supremum, the least upper bound; not necessarily maximum. E.g. {x | x < 2} only has supremum 2k1,,N i1ikN i1,,N k1A maxN 1/2A2A2ik i,k1 MNik i1 k1AxAx xsupAx xA x x x1 xA ikE.g. in 3DAAA 11 12 13 N Recall x 1 xAAAA 21 22 23 Matrix norm: size of a matrix AAA i1i 31 32 33 A max A max A A A maxc ,c ,c 3 ik 1k 2k 3k 1 2 31 k1,2,3 i1 k1,2,3 1 1 13 ik i1 i2 i3 1 2 3A max A max A A A maxr ,r ,r i1,2,3 k1 i1,2,3 1 1 1 3 1/2A A 2 A2 A2 A2 A2 A2 A2 A2 A2 A22ik 11 12 13 21 22 23 31 32 33 i,k1 import numpy as np# norms of vectors# norms of matricesA = np.array([[-1,1,1],[ 2, -2,0]])NumPy: norms of vectors/matricesnumpy.linalg.norm()a = np.array([-2, 1, 0])v_L1_norm = np.linalg.norm(a, 1)v_L2_norm = np.linalg.norm(a, 2)v_Linf_norm = np.linalg.norm(a, np.inf) print(v_L1_norm); print(v_L2_norm); print(v_Linf_norm)M_Lfro_norm = np.linalg.norm(A, ‘fro’)M_L1_norm = np.linalg.norm(A, 1)M_L2_norm = np.linalg.norm(A, 2)M_Linf_norm = np.linalg.norm(A, np.inf) print(M_Lfro_norm); print(M_L1_norm); print(M_L2_norm); print(M_Linf_norm)# default# defaultSolve the equations: E.g.f(x)x f(x) xxInitial guess: Update: Update:f(x) x0f(x)x1 x0 is small enough? NOUpdate: Solutionisxn+1:= f(xn) x* xn+1×1 x2:= f(x0) := f(x1)x2 x1 is small enough? NO xn+1 xn is small enough? YESGeneral idea of iterative methodAxb(AI)xbx or (AI)1x(AI)1bx or If it converges to x*, f(x*) = x* and x* is the solution (fixed point of f)However, does {x0, x1, x2, } converge? or lim xn1 xn 0 nExact solution:f (x) 0.1x 10 x 0.1x*10x*x* 10 11.111Iterative solution:General idea of iterative methodInitial guess: x0Update: x1Update: x2Update: x3Update: x4Update: x5If we think it is good enough, x* 11.111|x1 x0|=10 |x2 x1|=1|x3 x2| = 0.1 |x4 x3| = 0.01 |x5 x4| = 0.001:=0:= f(0):= f(10):= f(11):= f(11.1) := f(11.11)=0= 10= 11= 11.1= 11.11 = 11.111Example0.9Exact solution:f (x) 10x 0.1 x10 x * 0.1 x* x* 0.1 0.0111Iterative solution:Initial guess: x0 Update: x1 Update: x2 Update: x3 Update: x4 Update: x5 Iteration diverges:=0:= f(0):= f(0.1):= f(1.1):= f(11.1) := f(111.1)=0= 0.1= 1.1= 11.1= 111.1 = 1111.1|x1 x0| = 0.1 |x2 x1| =1|x3 x2| = 10 |x4 x3| = 100 |x5 x4| = 1000General idea of iterative methodWhether {x0, x1, x2, } converges or not depends on the form of function f(x). Which kind of f(x) guarantees convergence?Example9 From 1equation/1variable to Nequation/Nvariable f (x) x f(x) xf(x,x,,x) x 1 1 2 n 1f(x) is linear xx1 xf(x,x,,x) x 212n2fN(x1,x2,,xN) xN121 x1x22 x2f(x,x,,x)x 112 N1f(x) is nonlinear (next lecture)f(x,x,,x)x 212N20.5cosx1 0.5sinx2 0.5sinx1 0.5cosx2 x1 x2 f(x,x,,x) xN 1 2 NNWhether {x0, x1, x2, } converges or not depends on the form of function f(x). Which kind of f(x) guarantees convergence? Convergent sequenceA sequence {xn} = {x0, x1, x2, } is called convergent if there exists anelement x* s.t. xnlim xn x* 0 or limxn x* n nn0 1 2 3 4 5 6 7 8 9 10 11 12 10000 E.g. 11 lim 2n n2n 0 x 21/n , limx lim21/n 1n nnn 1nsin1 1 limnsin n n n from sympy import *n = symbols(‘n’, integer=True, nonnegative=True) pprint(limit(1/2**n, n, oo)) pprint(limit(2**(1/n), n, oo)) pprint(limit(n*sin(1/n), n, oo)) ContractionAn operator f is called contraction if there exists a constant q [0, 1) s.t.f(x)f(y) q xy We may think of f as a functionfor all x, y D 1variable example:f(x)0.1x10f (x) f (y) 0.1 x ySo, this f(x) is contractionf(x)10x10f (x) f (y) 10 x ySo, this f(x) is NOT contractionWhen f(x) is contraction, iteration for f(x) = x will convergeBanach fixed point theoremExistence of a unique fixed pointLet f be a contraction. There exists a unique element x* s.t. f(x*)x*This x* is called fixed point of f, because f maps x* to itself1variable example:f 100 0.1 100 10 Fixed point is x* 100f(x)0.1×10 999How to find the fixed pointStart from an arbitrary element x0. Define a sequence {xn} by recursion: xn1 :f(xn)We can find the fixed point bylimxn x* n 1variable example:f (x) 0.1x 3 (contraction)Exact solution for fixed point:0.1x*3x*x* 3 0.9xn 0.1xn1 3xn1 : f(xn)0.1xn 3 3)30.12 x (0.10 0.11)30.1(0.1x 0.12(0.1xn23)(0.10 0.11)30.13 xn2 n3n3(0.10 0.11 0.12)3 10.1nBanach fixed point theorem0.1 x0 (0.1 0.1 0.1 )30.1 x0 10.1 3n 0 1 n1 n1rn Sum of a geometric series: ar0 ar1 ar2 arn1 a 1r limxn 3 n 0.9 Exact solution for fixed point:10 x * 3 x* x* 3 9xn 10xn1 3xn1 : f(xn)10xn 3 3)3102 x (100 101)310(10x 102(10xn23)(100 101)3103 xn2 n3n3(100 101 102)3 110nBanach fixed point theorem1variable example:f (x) 10x 3 (not a contraction)10n x0 (100 101 10n1)310n x0 110 3limxn n xn 0.1xn1 0.30.1(0.1x 0.12(0.1xn2 0.3)(0.10 0.11)0.30.13 xn2 n3n3(0.10 0.11 0.12)0.3 10.1nBanach fixed point theorem1variable example:f (x) 10x 3 (not a contraction)We can rearrange the equation to construct a convergent recursion10x 3 x 0.1x 0.3 xf (x) 0.1x 0.3 (this is a contraction)xn1 : f(xn)0.1xn 0.3 0.3)0.30.12 x (0.10 0.11)0.30.1 x0 (0.1 0.1 0.1 )0.30.1 x0 10.1 0.3n 0 1 n1 nlimxn 3 n 9Banach fixed point theorem: proof Starting from an arbitrary x0 D, construct {xn} by recursion:Is this sufficient to prove convergence? NO E.g.xn1 :f(xn) (n0,1,) f is contractionx x f ( x ) f ( x ) q x x n1 n n n1 n n12 n nq f(xn1)f(xn2) q xn1 xn2 q x1 x0 0limxn1xn 0 nxn nlimxn1xn lim n1 n0n nBut, obviously, {xn} does not converge with nStronger convergence condition:Banach fixed point theorem: prooflim xmxn0 n ,mIf {xn} satisfies this condition, it is called Cauchy sequence xn nm2n lim x x lim 2n n lim1 n n,m m n n nSo, {xn} is NOT Cauchy sequence. It does not converge xn1 : f(xn ) f is contractionProve that {xn} is Cauchy sequence, i.e.,lim xmxn 0 n ,m Banach fixed point theorem: proofxn1 : f(xn ) f is contraction Prove that {xn} is Cauchy sequence ( q q q q ) x 1 x 0 0Is this sufficient to prove convergence? NOE.g.x 11n nnx x qn x x n1n 10xm xn xm xm1 xm1 xm2 xn2 xn1 xn1 xnxx x x x x x xtriangle inequalityn1m m1 m1 m2 n2 n1 m m1 n2 n1 n ,mn 1nlimxn lim1 e (irrationalnumber)n n nBut xn must be a rational number for any n. So, whichever n, xn can never reach the limit eOne additional condition:Banach fixed point theorem: proofxn1 :f(xn) For this sequence, there exists x* s.t.lim f(xn)f(x*) qlim xn x* 0 n nor limf(xn)f(x*) nSo, x* is the fixed pointxn U and limxn U nIf in the domain U all Cauchy sequences satisfy this condition, this domain U is called complete; as below, we only consider all elements in the complete domainIf f(x) is contraction, we can generate sequence by recursion:limxn x* nx*limxn1 limf(xn)f(x*) n n Priori (earlier) and posteriori (later) errorsLet f be a contraction with constant q < 1. The sequence constructed byxn1 :f(xn)converges to the unique fixed point x* with arbitrary initial guess x0Priori error x1 x0 bounds the exact error at step (n) via qnxn x* 1q x1 x0Posteriori error xn xn1 bounds the exact error at step (n) viaProof is simplePickx asx n1 0xx*q xx n 1q 1 0xn xm lim xn xm(qn qn1 qm2 qm1) x1 x0 lim(qn qn1 qm2 qm1) x1 x0mm nxnx* q xnxn1 1q min n for convergenceIf the desired accuracy is s.t. xn x* , how many iterative steps are needed? (i.e., n = ?)qn ln(1q)1 xnx*1qx1x0 n lnq lnqlnxx0.0 0.2 0.4 So, small q leads to fast convergenceq100 80 60 40 2000.6 0.8 1.0Posteriori error is used to check accuracy in the process of computation and terminate iterations when the desired accuracy is reached10 Construct an operator:f(x)Bxzf(x)f(y) B(xy) B xyB plays the role of constant q. f is contraction when B < 1 Hence, we can construct sequence by recursionxn1 :f(xn)Bxn zWith arbitrary initial guess x0, it converges to a unique fixed point x* s.t.Priori error estimateBn xnx*1B x1x0x*f(x*)Bx*z or (IB)x*zPosteriori error estimate B xnx*1B xnxn1A 11A22DA N1,N1 ANN 0AAA0 121,N11N 0AAA0 21Jacobi iteration Ax b ADAL AU2,N1 2N AL AUAA00A N1,1 N1,2 N1,N A AA 0 0 N1N2 N,N1 A 11 A1A1 A 2222D1 D A A1 N1,N1 A N1,N1 1 NN 0 121,N11N 0AAA0 21Jacobi iteration Ax b ADAL AUInverse of astraightforward 11 NN 0AAA2,N1 2N AL AUAA00A N1,1 N1,2 N1,N A AA 0 0 N1N2 N,N1 Adiagonal matrix isAssume that all diagonal elements are nonzero xD1(AL AU)xD1b(DAL AU)xb Dx(AL AU)xbRecall Banach fixed point theorem: when B < 1, we can constructBzsequence by recursionWith arbitrary initial guess x0, it converges to a unique fixed point x* s.t.Jacobi iteration:xn1 :D1(AL AU)xn D1bxn1 :Bxn z x*Bx*zxn1 :xn D1(bAxn) residual xn1 :D1(AL AU)xn D1bik ix N A x b ( i 1 , , N )n1,iki ii iin,k k1A AJacobi iteration algorithmeps=1e-10 #convergencecriterionforposteriorierror x = zeros(N) # initial guessfor n = 1, , max_nxprev = x # store x obtained from previous iterative stepfor i = 1, , N sum = 0Method of simultaneous displacementsfor k = 1, , N if k == i continuesum = sum – A[i][k] * xprev[k]x[i] = (sum + b[i]) / A[i][i]if norm(x xprev) <= eps breakif x is None: x = np.zeros(N)for count in range(1, max_n+1):xprev = x.copy()for i in range(N):Can be simplified by elementwise operationsxn1 :D1(AL AU)xn D1bik ix N A x b ( i 1 , , N )n1,iki ii iisum = 0for k in range(N):n,k k1A AJacobi iteration codedef jacobi(A, b, x=None, eps=1e-10, max_n=1000000): N = A.shape[0]if k == i: continuesum -= A[i,k] * xprev[k]x[i] = (sum + b[i]) / A[i,i]if np.linalg.norm(x – xprev) < eps: break print(f”Number of iteractions is {count}”) return xxn1 :D1(AL AU)xn D1bik ix N A x b ( i 1 , , N )n1,iki ii iin,k k1A AJacobi iteration codedef jacobi(A, b, x=None, eps=1e-10, max_n=1000000): N = A.shape[0]if x is None: x = np.zeros(N)D = np.diag(A)R = A – np.diag(D)# .diag(x): if x is 2D, returns 1D list of diagonal elements # if x is 1D, returns 2D list with x as diagonal elements for count in range(1, max_n+1):xprev = x.copy()x = (-np.matmul(R, xprev) + b) / Dif np.linalg.norm(x – xprev) < eps: break print(f”Number of iteractions is {count}”) return x Convergence condition xD1(AL AU)xD1bRecap:A max ABzN i1,,N k1B D1(AL AU) 1ikD1(AA) maxA A1 ork 1 , , N i k i1L Ui ki iN i1,,N k1 kiN A1max AJacobi iteration1 N A2A2D (ALAU) max A A 1 or i,k1 ik ii ik1 k1,,N i1 1/2ND1(AA) AA2 1 ik iiLU 2 i,k1 ik General trend: if {|Aii|} are large, convergence condition is satisfiedN 1/2 ik GaussSeidel iteration Axb, ADAL AU(DAL AU)xb (D AL )x AUx bx (D AL )1 AUx (D AL )1 b BzRecall Banach fixed point theorem: when B < 1, we can constructsequence by recursionWith arbitrary initial guess x0, it converges to a unique fixed point x* s.t.GaussSeidel iteration:xn1 :Bxn zx*Bx*zxn1 :(DAL)1AUxn (DAL)1bxn1 :(DAL)1AUxn (DAL)1bIn practice, at Step (n), we solve the system of linear equations:(DAL)xn1 AUxn bSince D + AL is lower triangular matrix, it can be solved simply byforward substitution (recall LU decomposition)A 11xy xn1,1y1A11 n1,1 1 A A 21 22x y x yAx A n1,2 2 n1,2 2 21 n1,1 22 N1 AAAxyxyAxA N1 N2k1 NN n1,N N n1,N Ni1 N Nxn1,iyi Ax A, y A x b Ax bikn1,k ii iU,ik n,k i k1 k1 ki1ikn,k iik ik ix i 1 A x N A x b ( i 1 , , N )n1,in1,kn,kk1 A ki1 A A ii ii iiNk n1,k NN xikxikx i n,k(i1,,N)n1,in1,kii ii iii1A NA bk1 Aki1 A AStep (n) & (n+1) data can be stored in same array Method of successive displacementsalready updatednot updated yetGaussSeidel iteration algorithmeps=1e-10 #convergencecriterionforposteriorierror x = zeros(N) # initial guessfor n = 1, , max_nxprev = xfor i = 1, , Nsum = 0for k = 1, , Nif k == i continue sum = sum A[i][k] * x[k]x[i] = (sum + b[i]) / A[i][i]if norm(x xprev) <= eps breakxikxikx i n,k(i1,,N)n1,in1,kii ii iidefgauss_seidel(A, b, x=None, eps=1e-10, max_n=1000000): N = A.shape[0]if x is None: x = np.zeros(N)for count in range(1, max_n+1):i1A NA bk1 Aki1 A AStep (n) & (n+1) data can be stored in same array Method of successive displacementsalready updatednot updated yetxprev = x.copy()for i in range(N):Cannot be written as elementwise operationsGaussSeidel iteration codesum = 0for k in range(N):if k == i: continueThis part can besum -= A[i,k] * x[k]x[i] = (sum + b[i]) / A[i,i]elementwiseif np.linalg.norm(x – xprev) < eps: break print(f”Number of iteractions is {count}”) return x xikxikx i n,k(i1,,N)n1,in1,kii ii iidefgauss_seidel(A, b, x=None, eps=1e-10, max_n=1000000): N = A.shape[0]if x is None: x = np.zeros(N)D = np.diag(A)i1A NA bk1 Aki1 A AStep (n) & (n+1) data can be stored in same array Method of successive displacementsalready updatednot updated yetR = A – np.diag(D)for count in range(1, max_n+1):xprev = x.copy()for i in range(N):GaussSeidel iteration codex[i] = (-np.inner(R[i,:], x) + b[i]) / D[i]if np.linalg.norm(x – xprev) < eps: break print(f”Number of iteractions is {count}”) return xStart from Jacobi method:Jacobi method with relaxation:xn1 xn D1(bAxn)Small q or B fast convergencexn1 (ID1A)xn D1bRelaxation method xn1 xn D1(bAxn)Weight > 0 is called relaxation parameter
Problem: which value of achieves the largest convergence rate?
B reaches minimum, i.e., the largest convergence rate, when 2
residual
Bz
2 max min
where max and min are the maximum and minimum eigenvalues of the matrix D1(AL + AU), i.e., the B of the original Jacobi with = 1

Relaxation method (SOR method)
Start from GaussSeidel method:
xn1 (DAL)1AUxn (DAL)1b
xn1 xn D bALxn1 (DAU)xn 1
GaussSeidel method with relaxation or Successive overrelaxation (SOR) method:
xn1 xn D bALxn1 (DAU)xn 1
(DAL)xn1 [(1)DAU]xn b Problem: which value of guarantees convergence of SOR?
Kahan theorem: A necessary condition for convergence of SOR is that 0 <<2Ostrowski theorem: If A is Hermitian (A = A) and positive definite (vAv > 0), SOR converges for all 0 < < 2 Stationary iterative methods vs. Krylov subspace methodsJacobi iteration GaussSeidel iteration SORmethodMxn1 Nxn b M = DN = (AL + AU)N = AUN=(1 1)DAUStationary iterative methodsM = D + AL M=1D+ALJacobi iteration, GaussSeidel iteration, and SOR method have largely lost their interest with the advent of Krylov subspace method. So, they are not implemented in NumPy or SciPyKrylov subspace methodsKrylov subspace method is the modern iterative method one the top 10 algorithms in 20th centuryGeneralized minimal residual method (GMRES) Conjugate gradient method (CG) A 3D space is spanned by 3 orthonormal basis vectors e1, e2, e3e3e1,1 e2,1 e3,1 ee , ee , ee 1 1,2 e 2 2,2 e 3 3,2 e x 1,3eTe 1, 2,3 ik 3,3i k ik 0, ik e2Any vector in 3D space is a linear combination of basis vectorsxx1e1 x2e2 x3e3Krylov subspace method3D space = span{e1, e2, e3}e1 2D space = span{e1, e2}2D space is a subspace of 3D spacee32D space 3D space xx1e1 x2e2x e2is a vector in the subspacee1Krylov subspace methode1,1 e2,1 e3,1 ee , ee , ee e31 1,2 e 2 2,2 e 3 3,2 e 3,3 span{e1, e2, e3, e4} is impossible, 1,3ei contains 3 components. So,e4e2 2,3Krylov subspace methodbecause e4 can be represented bye4 ae1 be2 ce3 e1The number of basis vectors dimension of basis vector K (A,v) spanM v,Av,,AM1 vKrylov subspace methodThe mth Krylov subspace:K (A,v)span 2 m1 m v,Av,A v,,A v Basis vectors are v, Av, A2v, , Am1vIf v is N1 and A is NN, m N (number of basis vectors dimension) K (A,v)span v K2(A,v)spanv,Av1K (A,v)span 2 3 v,Av,A vK1 K2 KM KM1 KNIts possible that all Km with m M are same K (A,v) K (A,v) M M1MAK (A,v) span Av,A2v,,AMvSo, KM is invariant w.r.t. AKrylov subspace method Ax bK (A,x)span 2 M1 M x,Ax,A x,,A xSolution x is a vector in KM(A, x), i.e., x is some linear combination of the basis vectorsKM is invariant: AKM(A, x) = KM(A, x)K (A,x) AK (A,x) span Ax,A2x,A3x,,AMx span 2 M1 K (A,b)MMb,Ab,A b,,A b M KM (A,x) KM (A,b)So, x KM(A, b), i.e., x is some linear combination of its basis vectors: x c1bc2Abc3A2bcMAM1b We can construct the Krylov subspace:K (A,b)span 2M1 bKrylov subspace method Ax bM b,Ab,A b,,A Then, construct the solution:x c1bc2Abc3A2bcMAM1bTwo problems:(i) How to construct the Krylov subspace? (ii) How to construct the solution? Step (1):Construct K1(A, b) = span{b}; Construct x = c1b Good enough? NOStep (2):Krylov subspace method Ax bIterative improvementConstruct K2(A, b) = span{b, Ab}; Construct x = c1b + c2Ab Good enough? NOStep (3):Construct K3(A, b) = span{b, Ab, A2b}; Construct x = c1b + c2Ab + c3A2b Good enough? NOKrylov subspace method Generalized minimal residual method (GMRES)(i) Construction of Krylov subspaceAx bVectors {b, Ab, A2b, , Am1b} are in general NOT orthogonal GramSchmidt process: Given a set of linearly independent basis vectors {b, Ab, A2b, } (m N), generate an orthonormal set of basis vectors {q1, q2, q3, }(ii) Construction of solutionThe solution is constructed by the basis vectors of mth Krylov subspace xm =c1q1 +c2q2 ++cmqmFind {c1, , cm} s.t. b Axm2 is minimized GramSchmidt processGiven {v1, v2, v3}They are linearly independent Anyvectorin3Dspacex=av1 +bv2 +cv3 Find orthonormal basis set {q1, q2, q3}v2v3v1 GramSchmidt processGiven {v1, v2, v3}uv,q 111u1 u1 2q1v1 = u1GramSchmidt processGiven {v1, v2, v3} uv(vTq)qvq(qTv), qu2 u2 22221121122q1Proj v2 to q1v2q2u2 Proj v3 to q2 q2GramSchmidt processGiven {v1, v2, v3}u v (vTq )q (vTq )q v q (qTv )q (qTv ), q u3 u3 23 3311 3223113 2233q1Proj v3 to q1q3u3v3q2GramSchmidt processGiven {v1, v2, v3}Find orthonormal basis set {q1, q2, q3}Anyvectorin3Dspacex=av1 +bv2 +cv3 Equivalently, x = aq1 + bq2 + cq3q3q1v1v2v3 GramSchmidt processGiven {v1, v2, v3, }Find orthonormal basis set {q1, q2, q3, }uv, q 111u1 u1 2u v q (qTv ), q 22112 2u2 u2 2uvq(qTv)q(qTv), q 331132233u3 u3 2ku v q(qTv ),q k1uk1 uk1 2k1 k1i1i i k1 Modified GramSchmidt process (Arnoldi)Given {b, Ab, A2b, }Find orthonormal basis set {q1, q2, q3, }v b, u v , q 1111u1 u1 2v Aq , u v q (qTv ), q 2122112 2u2 u2 2v Aq , u v q (qTv )q (qTv ), q 32331132233u3 u3 2kv Aq, u v q(qTv ),q k1uk1 uk1 2k1 k k1 k1i1i i k1 Modified GramSchmidt process (Arnoldi)Given {b, Ab, A2b, }Find orthonormal basis set {q1, q2, q3, }q1 qb b22Aq q (qT Aq ) 1111Aq q (qT Aq ) 11112Aq q(qTAq)q(qTAq)2112222 Aq q(qTAq)q(qTAq)q3 21122222Aq kkq (qT Aq ) iiki1 qk1kkiikAq q (qT Aq )i12 A q 1 q 1 H 1 1 q 2 H 2 1Aq qH qH qH q Ti A q k , kjjkModified GramSchmidt process (Arnoldi)Given {b, Ab, A2b, }Find orthonormal basis set {q1, q2, q3, }k1kkC2 112 222 332H j1 2m1 Aqm qHi im i1k k1i1 i1i ik i1kk i1 i1Aq qk q (qT Aq ) Aq iikkq (qT Aq ) iikAq Aqq C q(qTAq) qHq (qT Aq ) iikkk1 2iiki kik Aqkq(qTAq),ik1 A A| | 11 1N q1 qm A A| |Modified GramSchmidt process (Arnoldi) N1 NN H11 H12H1m HHHH 21 22 2,m1 2m | ||0 H32 H3,m1 H3mqqq 1 m m1 | ||00HH m1,m1 m1,m H1,m1 0 0 Hm,m1 Hmm 000Hm1,m AQm = Qm+1HmHessenberg matrixq1 b b2H qTAq 11 1 1Modified GramSchmidt process (Arnoldi) How to obtain Hessenberg matrix Hm? q Ti A q k , ki k Aqk q(qTAq) , ik1H ikjjk j1 2q:Aqq(qTAq)H q1 1 1 1 21 2 q1 q2 b b2H qTAq 11 1 1Modified GramSchmidt process (Arnoldi) How to obtain Hessenberg matrix Hm? q Ti A q k , ki k Aqk q(qTAq) , ik1H ikjjk j1 2q:Aqq(qTAq)H q1 1 1 1 21 2H qTAq, H qTAq 12 1 2 22 2 2q qTT2 q:Aqq(qAq)q(qAq)H q2 1 1 2 2 2 2 32 2 q1 q2 b b2H qTAq 11 1 1q3 q qTTT2q:Aqq(qAq)q(qAq)H q2 1 1 2 2 2 2 32 2Modified GramSchmidt process (Arnoldi) How to obtain Hessenberg matrix Hm? q Ti A q k , ki k Aqk q(qTAq) , ik1H ikjjk j1 2q:Aqq(qTAq)H q1 1 1 1 21 2H qTAq, H qTAq 12 1 2 22 2 2q qTTH qTAq, H qTAq, H qTAq 13 1 3 23 2 3 33 3 32 q:Aqq(qAq)q(qAq)q(qAq)H q3 1 1 3 2 2 3 3 3 3 43 2q1 q2 b b2H qTAq 11 1 1q3 q4 q qTTTq q22q:Aqq(qAq)q(qAq)q(qAq)H q3 1 1 3 2 2 3 3 3 3 43 2Modified GramSchmidt process (Arnoldi) How to obtain Hessenberg matrix Hm? q Ti A q k , ki k Aqk q(qTAq) , ik1H ikjjk j1 2q:Aqq(qTAq)H q1 1 1 1 21 2H qTAq, H qTAq 12 1 2 22 2 2q qTT2q:Aqq(qAq)q(qAq)H q2 1 1 2 2 2 2 32 2H qTAq, H qTAq, H qTAq 13 1 3 23 2 3 33 3 3 Find the optimized solutionGMRES: generalized minimal residual method Construct the solution: | | c1 xKcqcqcqq qQcm m 11 22 mm1 mmm | |c Minimize the residual:bb2q1b2Qm1e1, AxmAQmcmQm1Hmcmmin bAx min r cm2cm2mmbAxm 2 b2Qm1e1Qm1Hmcm 2 T 1/2 m b Q e Q H c b Q e Q H c 2 m1 1 m1 m m 2 m1 1 m1 m m beTcTHTQTQ beHc 21 mm m1m1 21 mm T 1/2 b e H c b e H c b e H c 21mm21mm21mm 21/2Find the optimized solutionmin b e H c AQR cm 21 mm2 RxQTbbAxRecall Page 84 of Lecture 3. This is a linear regression problem.The shape of Hm is (m + 1)m. It can be solved by QR decomposition: Hm mRmm : orthogonal matrixRm : upper (right) triangular matrixR c b T e mm2m1Since Rm is upper (right) triangular matrix. cm can be solved by backward substitution qH21,H22 H R c x 2qH22222qH31, H32, H33 H R c x 3qH33333Krylov subspace method Generalized minimal residual method (GMRES)H Ax bq 11 H R c x1qH11111 21NoIfbAx12 <?No32 IfbAx22 <?x*=x Yes IfbAx <? 33243SciPy : Generate a random sparse matriximport numpy as npimport matplotlib.pyplot as pltimport scipy.sparseM = 5; N = 10 # dimensions of matrixA = scipy.sparse.random(M, N, 0.5)# M*N matrix. Randomly pick 50% elements to be 0. The other 50% elements are random numbers in [0, 1)print(A) # note how the sparse matrix is stored print(A.toarray()) # convert A to format of NumPy array fig, axs = plt.subplots()axs.imshow(A.toarray(), cmap=’rainbow’)axs.set_xlabel(‘col index’)axs.set_ylabel(‘row index’)plt.tight_layout()plt.show() SciPy : Solve system of linear equations by GMRESimport numpy as npimport matplotlib.pyplot as pltimport scipy.sparseimport scipy.sparse.linalgN = 100A = scipy.sparse.random(N, N, 0.5) + 4*scipy.sparse.eye(N) b = np.random.rand(N, 1)fig, axs = plt.subplots()axs.imshow(A.toarray(), cmap=’rainbow’)plt.show()guess = np.zeros(shape=b.shape)x, exitCode = scipy.sparse.linalg.gmres(A, b, x0=guess, tol=1e-10, maxiter=10000)# x0 : the initial guess of the solution# tol : tolerance |b – Ax|/|b|# maxiter : maximum steps allowedprint(x)print(exitCode)# =0: successful exit; >0: number of iterations but tolerance has not been achieved; <0: illegal input or breakdownKrylov subspace method Conjugate gradient method (CG)Ax bOnly consider the case where A is symmetric and positivedefiniteIf the initial guess is x0, then the residual is r bAx00A(xx )r 00A x b xK (A,b)span 2M1 bRecallxK (A,b)xx K (A,r )spanr ,Ar ,A2r ,,AM1r M b,Ab,A b,,AM 0M00000xx span 02 M1r ,Ar ,A r ,,A r 0000So, exact solution x x0 is linear combination of r0, Ar0, A2r0, , AM1r0Conjugate gradient method (CG)Iterative improvementx xWe want to approach this exact solution by iterationThe exact solution is0 span00002 M1r ,Ar ,A r ,,A r xxr n1 n nnHow to determine rn? How to determine n?Conditions for the choice of residuals {ri}x x span 2 n1 K (A,r )r,Ar,Ar,,A rn0 000 0n0r bAx r A(x x )span 2 n K (A,r ) r ,Ar ,A r ,,A r n n 0 n 0 0 0 0 0 n1 0The 1st condition: r Kii i1Let k i + 1rK spanA2r 0Conjugate gradient method (CG) Conditions for the choice of residuals {ri}xix0Ki andxi1x0Ki1r x x (x x )(x x )Kii i1 i i1 0 i 0 i1k1 k r,Ar,,A r,Ark k1 0 0 0 0r2K3rk is outside Kk space. E.g. if Kk is a plane, rk must be a vector pointing out of the plane (the additional basis Akr0 leads to the out ofplane component)Ar0 r1 K2Idea of CG Most efficient search is rk all basis vectors of Kk:r0r(K).Sor (K) k k k1 k Conjugate gradient method (CG) Conditions for the choice of residuals {ri}Ki+1(Kk) (Ki+1) KkLet k i + 1 A(r)Ax Ax (bAx )bAx r r(K)(K )kkk1k k1 kk1kki1The 2nd condition:Both conditions must be satisfied whenA( r )(K ) kk i1Aconjugate: (r)T A( r )0rTAr 0 (ki) iikkik 1 2x2A 1T1xA1b 2 AConjugate gradient method (CG)Understanding CG assumption ({ri} are Aconjugate) from geometryDefine Anorm:Recall that A is symmetric and positivedefiniteSeek x s.t. the Anorm of the residual is minimum: xargmin1 xA1b2x xTAx A2xA1b AxA1b2xTAxxTbbTATAxbTATb ATA1T TTT11T T1T1 2 x A x x b b x b A b 2 x A x x b 2 b A b constantSo, solving system of linear equations is equivalent to finding minimumof the functionf (x) 1 xT Ax xT b 2How does the function f(x) look like? Consider 2D case (N = 2):Conjugate gradient method (CG)Understanding CG assumption ({ri} are Aconjugate) from geometry Seek x to minimize the function:f (x) 1 xT Ax xT b 2f(x ,x )1x x A A x x 1 2 1 2 11 12 1 1x b 2 1 2AAxb 21 22 2 2 1A x21(A A )xx1A x2bxbx 2 11 1 2 12 21 1 2 2 22 2 1 1 2 2 Example:A A 1 0, b 0 11 12 1 AA01b0 21 22 2 f(x ,x )1×2 1×2 122122 Example:A A 1 0, b 0 11 12 1 AA01b0 21 22 2 Example:A A 2.5 0 , b 0 11 12 1 A A 0 0.4 b 0 21 22 2 Example:A A 1 0.8, b 0 11 12 1 A A 0.8 1 b 0 21 22 2 In 2D, N = 2, and A = I: converge by at most 2 steps. New search direction is along gradient (rapid decrease) and orthogonal to old onex01r1xx1 0r0rTr0, rf(x) 101xx1r bAx 00So, direction of r0 is arbitrary. But, always we only need 2 stepsIf A = I and still use orthogonal vectors, rTr0, rf(x)xxk(ki)x01r1 2r2 x3 4r40r0x1kikThen, we need much more steps for convergenceThe orthogonal vectors are not good basis vectors for searching the minimumx23r3x4x If A = I and use Aorthogonal vectors, rTAr0 (ki)kiThen, again we only need 2 steps for convergence in 2D. Why?0r0 x0x11r1xIf A = I and use Aorthogonal vectors, rTAr0 (ki)kiThen, again we only need 2 steps for convergence in 2D. Why?f xT Ax f xT (FT AF1)xtransformationofcoordinate: xFx If pick F s.t. FTAF1 = I A = FTF, after transformation by F and under the new coordinate system (primed), the search directions should beorthogonalWhat does this require for {ri}?rTr0 kirTrrTFTFr rT Ar 0 kikikiAconjugateIn the coordinate system with A = I (primed), the search directions should be orthogonal, while in the coordinate system with A = I (unprimed), the search direction should be Aorthogonal Conjugate gradient method (CG)How to construct Aorthogonal basis vectors? GramSchmidt processGiven {r0, r1, r2, }Find Aorthogonal basis set {q0, q1, q2, }ur, q 000u0 u0u r (qTAr)q , q 11010 1A u1ur(qTAr)q(qTAr)q, q 220201212u2 A k 1ur (qTAr)q, quk ukkki0iki ku1 u2AA qT k1Aq k1 k qT Ar (qT Ar )(qT Aq )Conjugate gradient method (CG)How to construct Aorthogonal basis vectors? GramSchmidt processCheck if {q0, q1, q2, } constructed by GramSchmidt process are AorthogonalTkk qTAq uTAu 1kk2ukqTAqq0Au1 1 qTAr(qTAr)(qTAq)001u1 u1010100 AAAk1qAu 1 Tk k1 ki k k1 i i0 uk A uk A 1 qT Ar(qT Ar)(qT Aq )0uk k1 k k1 k k1 k1 Akki0i k i kConjugate gradient method (CG)r bAx 00u r, q u0 x x q 000u 10000Ar bAx bAx Aq r Aq1 1 000000ur(qTAr)q, q u1 xxq 110101u 21111A r r Aq2111ur(qTAr)q(qTAr)q, q u2 xxq 220201212u 3222 k 1How to determine {i}?ur (qTAr)q,quk ukxxq k1 k k kA2AConjugate gradient method (CG)How to determine {i}? Search k s.t. f(xk+1) is minimizedf(x )f(x q)1(x q)TA(x q)(x q)Tb k1 k kk 2k kk k kk k kk1qTAq 2 1xTAq 1qTAx qTb 1xTAx xTb kkkkkkkkkkkk2222 11f(xq)qTAq xTAq qTAxqTb0 kkk kkk2kk2kkkkqTAq qTb xTAq qTAx1 1 1 kkkkkkkk22T 1TT 1T ATA T Tq b q A x q A x q ( b A x ) q r k2kk2kkkkkk kqrTkk qk AqkT Conjugate gradient method (CG)r bAx 00Tu r, q u0000u 10000A 1000r r Aq ur(qTAr)q, q u1Txxq 110101u 2111r r Aq 21111AT2 220201212u 3222ur(qTAr)q(qTAr)q, q u2 2AT k 1Tur (qTAr)q,qkxxq k1 k k kkki0i k i ku ukTA 0qr 1qr kqr00q0Aq0 x x qT11 q1Aq1Tqr22 q2Aq2kk qk Aqkxxq N = 100SciPy : Solve system of linear equations by CGimport numpy as npimport matplotlib.pyplot as pltimport scipy.sparseimport scipy.sparse.linalg# note that A must be a symmetric, positive-definite matrixA = scipy.sparse.random(N, N, 0.5)A = A.toarray()A = A + np.transpose(A)A = scipy.sparse.csr_matrix(A) # convert numpy array to a compressed sparse row matrix (csr)print(A)A += 4*scipy.sparse.eye(N)b = np.random.rand(N, 1)guess = np.zeros(shape=b.shape)x, exitCode = scipy.sparse.linalg.cg(A, b, x0=guess, tol=1e-10, maxiter=10000)print(x)print(exitCode) System of linear equationsAx bHow to solve a system of linear equations? Direct method(i) Gaussian elimination (ii) LU decomposition (iii) QR decompositionIterative methodBasics 1: Normed spaceBasics 2: Banach fixed point theorem Stationary methods (old)(i) Jacobi method(ii) GaussSeidel method(iii) Relaxation methodKrylov subspace methods (modern)(i) Generalized minimal residual (GMRES) (ii) Conjugate gradient (CG)