[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: AxAxAxb
11 1 12 2 1N N 1 AxAxAxb
21 1 22 2 2N N 2 
AxAxAxb N1 1 N2 2 NN N N
AAAxb  11 12 1N   1   1 
21 22 2N 2 2 AAAxb
 A A Ax b
 N1 N2 NN   N   N  A xb

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 don’t 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 large‐N, sparse matrix E.g. FEM Ku = f
For aircraft flap actuator

Guess a solution x0 (probably it’s bad)
Residual
Correct the guess by residual
bAx0
x1 x0 C(bAx0)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(bAx1)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,yXandallα 
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
 p1/p
xx
p
i
 i1 

p i1  l norm is called maximum norm
i1,,N
Vector norm: “size” of a vector
2
i
 i1 
 N p 1/p x x 
p
 i1 N
i

x1 x i1
i
 21/2
x  x   xTx
N
l2 norm is called Euclidean norm because it is just the length of a vector
N
 p1/p
xlim x  maxx ii

2 x 3
i1
 3 2 1/2
Vector norm: “size” of a vector
Equivalence of the norms
Theorem: on a finite‐dimensional vector space, all norms are equivalent. When one norm goes to 0, all the other norms go to 0
 x1  xx
3
x1xi x1 x2 x3
x x   x2 x2 x2 (lengthofx)
2
i 123  i1 
x maxxmaxx,x,x  i1,2,3 i 1 2 3
Onlywhenx1,x2,x3 0atthesametime,‖x‖1 0.Andatthesame time ‖x‖2 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 2k1,,N i1ikN i1,,N k1A   maxN 1/2A2A2ik i,k1  MNik i1 k1AxAx xsupAx xA x x x1 xA ikE.g. in 3DAAA  11 12 13 N Recall x 1  xAAAA  21 22 23 Matrix norm: “size” of a matrix AAA i1i 31 32 33 A max A max A A A maxc ,c ,c  3  ik 1k 2k 3k 1 2 31 k1,2,3 i1 k1,2,3  1 1 13  ik i1 i2 i3 1 2 3A max A max A A A maxr ,r ,r i1,2,3 k1 i1,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,k1  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) xxInitial guess: Update: Update:∙∙∙f(x) x0f(x)‖x1 – x0‖ is small enough? NOUpdate: Solutionis∙∙∙xn+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 methodAxb(AI)xbx or (AI)1x(AI)1bx 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 xn1  xn  0 nExact solution:f (x)  0.1x  10  x 0.1x*10x*x* 10 11.111Iterative 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.0111Iterative 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 1‐equation/1‐variable to N‐equation/N‐variable f (x)  x  f(x)  xf(x,x,,x)  x 1 1 2 n 1f(x) is linear xx1  xf(x,x,,x)  x 212n2fN(x1,x2,,xN)  xN121 x1x22  x2f(x,x,,x)x 112 N1f(x) is nonlinear (next lecture)f(x,x,,x)x 212N20.5cosx1 0.5sinx2 0.5sinx1 0.5cosx2 x1  x2 f(x,x,,x) xN 1 2 NNWhether {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 nn0 1 2 3 4 5 6 7 8 9 10 11 12∙∙∙ 10000 ∙∙∙E.g. 11    lim  2n   n2n  0 x  21/n , limx  lim21/n 1n nnn      1nsin1 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 xy We may think of f as a functionfor all x, y  D 1‐variable example:f(x)0.1x10f (x)  f (y)  0.1 x  ySo, this f(x) is contractionf(x)10x10f (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 itself1‐variable example:f  100   0.1 100  10  Fixed point is x*  100f(x)0.1x10 999How to find the fixed pointStart from an arbitrary element x0. Define a sequence {xn} by recursion: xn1 :f(xn)We can find the fixed point bylimxn x* n 1‐variable example:f (x)  0.1x  3 (contraction)Exact solution for fixed point:0.1x*3x*x* 3 0.9xn 0.1xn1 3xn1 : f(xn)0.1xn 3 3)30.12 x (0.10 0.11)30.1(0.1x 0.12(0.1xn23)(0.10 0.11)30.13 xn2 n3n3(0.10 0.11 0.12)3 10.1nBanach fixed point theorem0.1 x0 (0.1 0.1 0.1 )30.1 x0  10.1 3n 0 1 n1 n1rn  Sum of a geometric series: ar0 ar1 ar2 arn1 a 1r limxn  3   n 0.9 Exact solution for fixed point:10 x * 3  x*  x*   3 9xn 10xn1 3xn1 : f(xn)10xn 3 3)3102 x (100 101)310(10x 102(10xn23)(100 101)3103 xn2 n3n3(100 101 102)3 110nBanach fixed point theorem1‐variable example:f (x)  10x  3 (not a contraction)10n x0 (100 101 10n1)310n x0  110 3limxn  n xn 0.1xn1 0.30.1(0.1x 0.12(0.1xn2 0.3)(0.10 0.11)0.30.13 xn2 n3n3(0.10 0.11 0.12)0.3 10.1nBanach fixed point theorem1‐variable 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)xn1 : f(xn)0.1xn 0.3 0.3)0.30.12 x (0.10 0.11)0.30.1 x0 (0.1 0.1 0.1 )0.30.1 x0  10.1 0.3n 0 1 n1 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.xn1 :f(xn) (n0,1,) f is contractionx  x  f ( x )  f ( x )       q x  x n1 n n n1 n n12 n nq f(xn1)f(xn2) q xn1 xn2 q x1 x0 0limxn1xn 0 nxn  nlimxn1xn lim n1 n0n nBut, obviously, {xn} does not converge with nStronger convergence condition:Banach fixed point theorem: prooflim xmxn0 n ,mIf {xn} satisfies this condition, it is called Cauchy sequence xn  nmα2n lim x x lim α2n n limα1 n n,m m n n nSo, {xn} is NOT Cauchy sequence. It does not converge xn1 : f(xn ) f is contractionProve that {xn} is Cauchy sequence, i.e.,lim xmxn 0 n ,m Banach fixed point theorem: proofxn1 : 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 11n nnx x qn x x n1n 10xm xn  xm xm1 xm1 xm2 xn2 xn1 xn1 xnxx x x x x x xtriangle inequalityn1m m1 m1 m2 n2 n1 m m1 n2 n1 n ,mn 1nlimxn lim1  e (irrationalnumber)n n nBut 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: proofxn1 :f(xn) For this sequence, there exists x* s.t.lim f(xn)f(x*) qlim xn x* 0 n nor limf(xn)f(x*) nSo, x* is the fixed pointxn U and limxn U nIf 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* nx*limxn1 limf(xn)f(x*) n n Priori (earlier) and posteriori (later) errorsLet f be a contraction with constant q < 1. The sequence constructed byxn1 :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* 1q x1 x0Posteriori error ‖xn ‒ xn‒1‖ bounds the exact error at step (n) viaProof is simplePickx asx n‒1 0xx*q xx n 1q 1 0xn xm lim xn xm(qn qn1 qm2 qm1) x1 x0 lim(qn qn1 qm2 qm1) x1 x0mm nxnx* q xnxn1 1q min n for convergenceIf the desired accuracy is ε s.t. ‖xn ‒ x*‖ ≤ ε, how many iterative steps are needed? (i.e., n = ?)qn ln(1q)1 ε xnx*1qx1x0 ε  n lnq lnqlnxx0.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)Bxzf(x)f(y)  B(xy)  B xy‖B‖ plays the role of constant q. f is contraction when ‖B‖ < 1 Hence, we can construct sequence by recursionxn1 :f(xn)Bxn zWith arbitrary initial guess x0, it converges to a unique fixed point x* s.t.Priori error estimateBn xnx*1B x1x0x*f(x*)Bx*z or (IB)x*zPosteriori error estimate B xnx*1B xnxn1A 11A22DA N1,N1  ANN 0AAA0 121,N11N 0AAA0 21Jacobi iteration Ax  b ADAL AU2,N1 2N AL AUAA00A N1,1 N1,2  N1,N A AA 0  0 N1N2 N,N1  A  11  A1A1 A 2222D1  D  A  A1   N1,N1  A N1,N1 1 NN 0 121,N11N 0AAA0 21Jacobi iteration Ax  b ADAL AUInverse of astraightforward  11 NN 0AAA2,N1 2N AL AUAA00A N1,1 N1,2  N1,N A AA 0  0 N1N2 N,N1 Adiagonal matrix isAssume that all diagonal elements are nonzero xD1(AL AU)xD1b(DAL AU)xb Dx(AL AU)xbRecall 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:xn1 :D1(AL AU)xn D1bxn1 :Bxn z x*Bx*zxn1 :xn D1(bAxn) residual xn1 :D1(AL AU)xn D1bik ix   N A x  b ( i  1 ,  , N )n1,iki ii iin,k k1A 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 element‐wise operationsxn1 :D1(AL AU)xn D1bik ix   N A x  b ( i  1 ,  , N )n1,iki ii iisum = 0for k in range(N):n,k k1A 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 xxn1 :D1(AL AU)xn D1bik ix   N A x  b ( i  1 ,  , N )n1,iki ii iin,k k1A 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 xD1(AL AU)xD1bRecap:A   max  ABzN i1,,N k1B  D1(AL AU) 1ikD1(AA) maxA A1 ork  1 , , N i k i1L Ui ki iN i1,,N k1 kiN A1max AJacobi iteration1 N A2A2D (ALAU) max A A 1 or i,k1 ik ii ik1 k1,,N i1 1/2ND1(AA) AA2 1 ik iiLU 2 i,k1 ik General trend: if {|Aii|} are large, convergence condition is satisfiedN 1/2 ik  Gauss‐Seidel iteration Axb, ADAL AU(DAL AU)xb (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.Gauss‐Seidel iteration:xn1 :Bxn zx*Bx*zxn1 :(DAL)1AUxn (DAL)1bxn1 :(DAL)1AUxn (DAL)1bIn practice, at Step (n), we solve the system of linear equations:(DAL)xn1 AUxn bSince D + AL is lower triangular matrix, it can be solved simply byforward substitution (recall LU decomposition)A  11xy xn1,1y1A11 n1,1   1   A A 21 22x y x yAx A n1,2  2 n1,2 2 21 n1,1 22   N1  AAAxyxyAxA N1 N2k1 NN  n1,N   N  n1,N  Ni1 N Nxn1,iyi Ax A, y A x b Ax bikn1,k ii iU,ik n,k i k1  k1 ki1ikn,k iik ik ix   i  1 A x  N A x  b ( i  1 ,  , N )n1,in1,kn,kk1 A ki1 A A ii ii iiNk n1,k NN  xikxikx i n,k(i1,,N)n1,in1,kii ii iii1A NA bk1 Aki1 A AStep (n) & (n+1) data can be stored in same array ‐‐ Method of successive displacementsalready updatednot updated yetGauss‐Seidel 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 breakxikxikx i n,k(i1,,N)n1,in1,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):i1A NA bk1 Aki1 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 element‐wise operationsGauss‐Seidel 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]element‐wiseif np.linalg.norm(x – xprev) < eps: break print(f”Number of iteractions is {count}”) return x xikxikx i n,k(i1,,N)n1,in1,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)i1A NA bk1 Aki1 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):Gauss‐Seidel 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:xn1 xn ωD1(bAxn)Small q or ‖B‖  fast convergencexn1 (IωD1A)xn ωD1bRelaxation method xn1 xn D1(bAxn)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 ‒D‒1(AL + AU), i.e., the B of the original Jacobi with ω = 1

Relaxation method (SOR method)
Start from Gauss‐Seidel method:
xn1 (DAL)1AUxn (DAL)1b
xn1 xn D bALxn1 (DAU)xn 1 
Gauss‐Seidel method with relaxation or Successive overrelaxation (SOR) method:
xn1 xn ωD bALxn1 (DAU)xn 1 
(DωAL)xn1 [(1ω)DωAU]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 (v∙Av > 0), SOR converges for all 0 < ω < 2 Stationary iterative methods vs. Krylov subspace methodsJacobi iteration Gauss‐Seidel iteration SORmethodMxn1 Nxn b M = DN = ‒(AL + AU)N = ‒AUN=(ω‒1 ‒1)D‒AUStationary iterative methodsM = D + AL M=ω‒1D+ALJacobi iteration, Gauss‐Seidel 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, e3e3e1,1  e2,1  e3,1  ee , ee , ee 1  1,2  e 2  2,2  e 3  3,2  e x 1,3eTe δ 1, 2,3 ik 3,3i k ik 0, ik e2Any vector in 3D space is a linear combination of basis vectorsxx1e1 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 xx1e1 x2e2x e2is a vector in the subspacee1Krylov subspace methode1,1  e2,1  e3,1  ee , ee , ee e31  1,2  e 2  2,2  e 3  3,2  e  3,3 span{e1, e2, e3, e4} is impossible, 1,3ei contains 3 components. So,e4e2 2,3Krylov subspace methodbecause e4 can be represented bye4 ae1 be2 ce3 e1The number of basis vectors ≤ dimension of basis vector K (A,v)  spanM v,Av,,AM1  vKrylov subspace methodThe m‐th Krylov subspace:K (A,v)span 2 m1 m v,Av,A v,,A v Basis vectors are v, Av, A2v, ∙∙∙, Am‒1vIf v is N1 and A is NN, m ≤ N (number of basis vectors ≤ dimension) K (A,v)span v K2(A,v)spanv,Av1K (A,v)span 2  3 v,Av,A vK1 K2 KM KM1 KNIt’s possible that all Km with m ≥ M are same K (A,v)  K (A,v) M M1MAK (A,v)  span Av,A2v,,AMvSo, KM is invariant w.r.t. AKrylov subspace method Ax  bK (A,x)span 2 M1  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 M1 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 c1bc2Abc3A2bcMAM1b We can construct the Krylov subspace:K (A,b)span 2M1  bKrylov subspace method Ax  bM b,Ab,A b,,A Then, construct the solution:x c1bc2Abc3A2bcMAM1bTwo 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? NO∙∙∙Krylov subspace method Generalized minimal residual method (GMRES)(i) Construction of Krylov subspaceAx  bVectors {b, Ab, A2b, ∙∙∙, Am–1b} are in general NOT orthogonal Gram‐Schmidt 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 m‐th Krylov subspace xm =c1q1 +c2q2 +∙∙∙+cmqmFind {c1, ∙∙∙, cm} s.t. ‖b – Axm‖2 is minimized Gram‐Schmidt processGiven {v1, v2, v3}They are linearly independent Anyvectorin3Dspacex=av1 +bv2 +cv3 Find orthonormal basis set {q1, q2, q3}v2v3v1 Gram‐Schmidt processGiven {v1, v2, v3}uv,q 111u1 u1 2q1v1 = u1Gram‐Schmidt processGiven {v1, v2, v3} uv(vTq)qvq(qTv), qu2 u2 22221121122q1Proj v2 to q1v2q2u2 Proj v3 to q2 q2Gram‐Schmidt processGiven {v1, v2, v3}u v (vTq )q (vTq )q v q (qTv )q (qTv ), q u3 u3 23 3311 3223113 2233q1Proj v3 to q1q3u3v3q2Gram‐Schmidt processGiven {v1, v2, v3}Find orthonormal basis set {q1, q2, q3}Anyvectorin3Dspacex=av1 +bv2 +cv3 Equivalently, x = a’q1 + b’q2 + c’q3q3q1v1v2v3 Gram‐Schmidt processGiven {v1, v2, v3, ∙∙∙}Find orthonormal basis set {q1, q2, q3, ∙∙∙}uv, q 111u1 u1 2u v q (qTv ), q  22112 2u2 u2 2uvq(qTv)q(qTv), q 331132233u3 u3 2ku v  q(qTv ),q  k1uk1 uk1 2k1 k1i1i i k1 Modified Gram‐Schmidt 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  32331132233u3 u3 2kv Aq, u v  q(qTv ),q  k1uk1 uk1 2k1 k k1 k1i1i i k1 Modified Gram‐Schmidt process (Arnoldi)Given {b, Ab, A2b, ∙∙∙}Find orthonormal basis set {q1, q2, q3, ∙∙∙}q1  qb b22Aq  q (qT Aq ) 1111Aq  q (qT Aq ) 11112Aq q(qTAq)q(qTAq)2112222 Aq q(qTAq)q(qTAq)q3 21122222Aq  kkq (qT Aq ) iiki1 qk1kkiikAq q (qT Aq )i12 A q 1  q 1 H 1 1  q 2 H 2 1Aq qH qH qH q Ti A q k ,  kjjkModified Gram‐Schmidt process (Arnoldi)Given {b, Ab, A2b, ∙∙∙}Find orthonormal basis set {q1, q2, q3, ∙∙∙}k1kkC2 112 222 332H   j1 2m1 Aqm qHi im i1k k1i1 i1i ik i1kk i1 i1Aq qk q (qT Aq ) Aq  iikkq (qT Aq ) iikAq Aqq C q(qTAq) qHq (qT Aq ) iikkk1 2iiki  kik Aqkq(qTAq),ik1 A A| |  11 1N   q1 qm A A| |Modified Gram‐Schmidt process (Arnoldi) N1 NN   H11 H12H1m  HHHH 21 22 2,m1 2m  | ||0 H32 H3,m1 H3mqqq    1 m m1 | ||00HH   m1,m1 m1,m  H1,m1 0 0  Hm,m1 Hmm 000Hm1,m AQm = Qm+1HmHessenberg matrixq1 b b2H qTAq 11 1 1Modified Gram‐Schmidt process (Arnoldi) How to obtain Hessenberg matrix Hm? q Ti A q k , ki  k Aqk q(qTAq) , ik1H ikjjk j1 2q:Aqq(qTAq)H q1 1 1 1 21 2 q1  q2 b b2H qTAq 11 1 1Modified Gram‐Schmidt process (Arnoldi) How to obtain Hessenberg matrix Hm? q Ti A q k , ki  k Aqk q(qTAq) , ik1H ikjjk j1 2q:Aqq(qTAq)H q1 1 1 1 21 2H qTAq, H qTAq 12 1 2 22 2 2q qTT2 q:Aqq(qAq)q(qAq)H q2 1 1 2 2 2 2 32 2 q1  q2 b b2H qTAq 11 1 1q3 q qTTT2q:Aqq(qAq)q(qAq)H q2 1 1 2 2 2 2 32 2Modified Gram‐Schmidt process (Arnoldi) How to obtain Hessenberg matrix Hm? q Ti A q k , ki  k Aqk q(qTAq) , ik1H ikjjk j1 2q:Aqq(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:Aqq(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 q2∙∙∙2q:Aqq(qAq)q(qAq)q(qAq)H q3 1 1 3 2 2 3 3 3 3 43 2Modified Gram‐Schmidt process (Arnoldi) How to obtain Hessenberg matrix Hm? q Ti A q k , ki  k Aqk q(qTAq) , ik1H ikjjk j1 2q:Aqq(qTAq)H q1 1 1 1 21 2H qTAq, H qTAq 12 1 2 22 2 2q qTT2q:Aqq(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  xKcqcqcqq qQcm m 11 22 mm1 mmm | |c Minimize the residual:bb2q1b2Qm1e1, AxmAQmcmQm1Hmcmmin bAx min r cm2cm2mmbAxm 2 b2Qm1e1Qm1Hmcm 2 T 1/2  m  b Q e Q H c b Q e Q H c   2 m1 1 m1 m m 2 m1 1 m1 m m  beTcTHTQTQ beHc 21 mm m1m1 21 mm T 1/2 b e H c b e H c   b e H c 21mm21mm21mm 21/2Find the optimized solutionmin b e H c AQR cm 21 mm2 RxQTbbAxRecall 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 ΩmRmΩm : 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 21NoIf‖b‒Ax1‖2 <ε?No32 If‖b‒Ax2‖2 <ε?x*=x Yes If‖b‒Ax‖ <ε? 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 positive‐definiteIf the initial guess is x0, then the residual is r bAx00A(xx )r 00A’ x’ b’ xK (A,b)span 2M1  bRecallxK (A,b)xx K (A,r )spanr ,Ar ,A2r ,,AM1r M b,Ab,A b,,AM 0M00000xx span 02 M1r ,Ar ,A r ,,A r 0000So, exact solution x ‒ x0 is linear combination of r0, Ar0, A2r0, ∙∙∙, AM‒1r0Conjugate gradient method (CG)Iterative improvementx  xWe want to approach this exact solution by iterationThe exact solution is0 span00002 M1r ,Ar ,A r ,,A r xxαr n1 n nnHow to determine rn? How to determine αn?Conditions for the choice of residuals {ri}x x span 2 n1 K (A,r )r,Ar,Ar,,A rn0 000 0n0r bAx 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 n1 0The 1st condition:α r Kii i1Let k ≥ i + 1rK spanA2r 0Conjugate gradient method (CG) Conditions for the choice of residuals {ri}xix0Ki andxi1x0Ki1αr x x (x x )(x x )Kii i1 i i1 0 i 0 i1k1 k r,Ar,,A r,Ark k1 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‐ of‐plane component)Ar0 r1 K2Idea of CG Most efficient search is rk  all basis vectors of Kk:r0r(K).Sor (K) k k k1 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 (bAx )bAx r r(K)(K )kkk1k k1 kk1kki1The 2nd condition:Both conditions must be satisfied whenA(α r )(K ) kk i1A‐conjugate: (αr)T A(α r )0rTAr 0 (ki) iikkik 1 2x2A 1T1xA1b 2 AConjugate gradient method (CG)Understanding CG assumption ({ri} are A‐conjugate) from geometryDefine A‐norm:Recall that A is symmetric and positive‐definiteSeek x s.t. the A‐norm of the residual is minimum: xargmin1 xA1b2x  xTAx A2xA1b AxA1b2xTAxxTbbTATAxbTATb ATA1T TTT11T T1T1   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 A‐conjugate) from geometry Seek x to minimize the function:f (x)  1 xT Ax  xT b 2f(x ,x )1x x A A x x 1 2 1 2  11 12  1  1x b  2  1 2AAxb 21 22  2   2 1A x21(A A )xx1A x2bxbx 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 onex0α1r1xx1 α0r0rTr0, rf(x) 101xx1r bAx 00So, direction of r0 is arbitrary. But, always we only need 2 stepsIf A ≠ I and still use orthogonal vectors, rTr0, rf(x)xxk(ki)x0α1r1 α2r2 x3 α4r4α0r0x1kikThen, we need much more steps for convergenceThe orthogonal vectors are not good basis vectors for searching the minimumx2α3r3x4x If A ≠ I and use A‐orthogonal vectors, rTAr0 (ki)kiThen, again we only need 2 steps for convergence in 2D. Why?α0r0 x0x1α1r1xIf A ≠ I and use A‐orthogonal vectors, rTAr0 (ki)kiThen, again we only need 2 steps for convergence in 2D. Why?f  xT Ax f  xT (FT AF1)xtransformationofcoordinate: xFx If pick F s.t. F–TAF–1 = I A = FTF, after transformation by F and under the new coordinate system (primed), the search directions should beorthogonalWhat does this require for {ri}?rTr0 kirTrrTFTFr rT Ar 0 kikikiA‐conjugateIn 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 A‐orthogonal Conjugate gradient method (CG)How to construct A‐orthogonal basis vectors? Gram‐Schmidt processGiven {r0, r1, r2, ∙∙∙}Find A‐orthogonal basis set {q0, q1, q2, ∙∙∙}ur, q 000u0 u0u r (qTAr)q , q  11010 1A u1ur(qTAr)q(qTAr)q, q 220201212u2 A k  1ur (qTAr)q, quk ukkki0iki ku1 u2AA qT k1Aq  k1 k  qT Ar  (qT Ar )(qT Aq )Conjugate gradient method (CG)How to construct A‐orthogonal basis vectors? Gram‐Schmidt processCheck if {q0, q1, q2, ∙∙∙} constructed by Gram‐Schmidt process are A‐orthogonalTkk qTAq uTAu 1kk2ukqTAqq0Au1 1 qTAr(qTAr)(qTAq)001u1 u1010100 AAAk1qAu 1  Tk k1 ki k k1 i i0 uk A uk A 1 qT Ar(qT Ar)(qT Aq )0uk k1 k k1 k k1 k1 Akki0i k i kConjugate gradient method (CG)r bAx 00u r, q  u0 x x αq 000u 10000Ar bAx bAx α Aq r α Aq1 1 000000ur(qTAr)q, q u1 xxαq 110101u 21111A r r α Aq2111ur(qTAr)q(qTAr)q, q u2 xxαq 220201212u 3222∙∙∙ ∙∙∙ k  1How to determine {αi}?ur (qTAr)q,quk ukxxαq k1 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 k1 k kk 2k kk k kk k kk1qTAq α2 1xTAq 1qTAx qTbα 1xTAx xTb kkkkkkkkkkkk2222 11f(xαq)qTAq α xTAq qTAxqTb0 αkkk kkk2kk2kkkkαqTAq qTb xTAq qTAx1 1 1  kkkkkkkk22T 1TT 1T ATA T Tq b  q A x  q A x    q ( b  A x )  q r k2kk2kkkkkkα kqrTkk qk AqkT Conjugate gradient method (CG)r bAx 00Tu r, q  u0000u 10000A 1000r r α Aq ur(qTAr)q, q u1Txxαq 110101u 2111r r α Aq 21111AαT2 220201212u 3222ur(qTAr)q(qTAr)q, q u2 2AT∙∙∙ ∙∙∙ k  1Tur (qTAr)q,qkxxαq k1 k k kkki0i k i ku ukTAα 0qrα 1qrα kqr00q0Aq0 x x αqT11 q1Aq1Tqr22 q2Aq2kk qk Aqkxxαq 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) Gauss‐Seidel method(iii) Relaxation methodKrylov subspace methods (modern)(i) Generalized minimal residual (GMRES) (ii) Conjugate gradient (CG)