CS计算机代考程序代写 python algorithm Computational Methods

Computational Methods
for Physicists and Materials Engineers
3
Systems of Linear Equations I Problem and Direct 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 

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

Displacement
Example 1 : Finite Element Methods (FEM)
Problem: A material of arbitrary shape. Force is applied on the surface and the points in the body. Displacement is applied to the surface. How will the material deform?
Surface force, t
Element
yy Body force,
b
Node
xx
Idea: discretize the material to the elements connected by nodes. Solve the displacement of the nodes

Stress of a wheel rim
Deformation of a metamaterial T. Frenzel, et al. Science 358 (2017)
https://www.comsol.com/multiphysics/fini
te‐element‐method
Human knee joint
Naghibi Beidokhti et al. (2017)
A composite under 0.12% strain
S. Das, et al., Construct. Build. Mater. 127 (2016)
Example 1 : FEM

Example 1 : FEM
AEa Ebf La Lb
How is each point along this bar displaced under the force f?
Discretize this rod into two elements, “a” and “b”. Two elements are terminated by three nodes, (1), (2) and (3). Solve the displacement u(x)
(1) u(1) a (2) u(2) b (3) u(3)
x First, focus on the element “a”. If we know u(1) and u(2), the
x(1) x(2) x(3)
displacement of any point on element “a” can be obtained by u(x)c0 c1x

Strain: Stress:
ε(x)du(x)1 1u(1) daT u(1)      
(1) u(1) x(1)
a (2) u(2) b (3) u(3)
Determine c0 and c1:
ccxu 0 (2) (1) a
01
c c x(2) u(2) 01
u(2) u(1) u(2) u(1) c 
u(x) x(2) x  a
x(2) x(3) u(x)c0 c1x
x
c  x(2)u(1)  x(1)u(2)  x(2)u(1)  x(1)u(2)
(1) (1) x x L

xx(1) u(1) naT (x)u(1) 
(2) (1) a 1xxL
a (2) (2) L Lu u
 a a (2) (2) dx LLuu
σ(x)Eaε(x)EadaT u(1)  u(2) 


2
u  2  a 1 1
u   ka
(1) u(1) a x(1)
(2) u(2) x(2)
b
(3) u(3)
Elastic energy in element “a”: x(2) 1
x(2) 1 Ua A (1) 2σTεdxA (1) 2 u(1)
a a a a aT a a K  A L E d d  A L E
 
aaa kAE L
 xx
 u(1) 1
(1) (2) a a a aT  (2)  (1) (2) a  (2) 
1
uuALEdd uuK
L Similarly, elastic energy in element “b”:
1  u(2) kb
kb AEb
 1 
u(2) daEadaT u dx
ka            k k 
L 1LL
aaaa a
Ub 2 u(2) u(3) Kbu , Kb k k withkb  L (3) b b b
x(3) u(1) 
x
(2)
 u(1)

(1) u(1) x(1)
a (2) u(2) x(2)
b (3) u(3)
1 UUa Ub 2 u(1)
 u(1) 1 u(2) Kau(2)2 u(2)
x(3)
 u(2)
 ka ka 0u(1)  1  
 2 u(1) u(2) u(3) ka ka 0u(2)  
0 0 0u(1) 1 (1) (2) (3)  b b  (2) 

2 u u u 0 k k u 

u(3) Kbu(3)  
0 0 0u(3)  
u(3)   
0 kb kb
 ka ka 0 u(1) 
1  1
2 u(1) u(2) u(3) ka ka kb kb u(2) 2uTKu
 0 kb kb u(3)  
x

(1) u(1) x(1)
a (2) u(2) x(2)
b
(3) u(3)
Free energy:
Equilibrium is reached when free energy is minimized:
π0Kuf0 Kuf
u k
uiKijuj ij
ka ka 0u(1)f(1) a a b b(2)(2)
 u uiKikuk ukKkjuj 
kkkkuf 0 kb kbu(2)f(2)
kij uK Ku2 Ku
  Boundary condition:


π  U  W  1 uT Ku  uT f 2
u(1) 0, f(2) 0,
f(3)  f
 u
i ik
i j j
x(3)
x
kj j
kj j

(1) u(1) a x(1)
(2) u(2) x(2)
b (3) u(3)
ka u(1) 0, f(2) 0, f(3) fka
ka ka kb kb
x(3) kbu(2)f(2)
ka kb 
kb u(2) 0
     kakb
 kb
(1) ka
 k a  k b u ( 3 )  f (2) kb (3)
kb u(3) This solution suggests spring model
f
x(1)
x(2)
x(3)
x

    kb u(2) f(2)
 0 
    kau(2)  f
x 0u(1) f(1)

Ku  f
u : displacements on all the nodes (to be determined)
f : applied forces on all the nodes (known) K : matrix of stiffness (known)
Example 1 : FEM
K K K  K K K  u  f   11 12 13 1,3N2 1,3N1 1,3N   1,x   1,x 
KKK K Kuf 22 23 2,3N2 2,3N1 2,3N 1,y 1,y
 KK K Kuf  33 3,3N2 3,3N1 3,3N   1,z   1,z 
 

K3N2,3N1 K3N2,3N  uN,x   fN,x  K Kuf
sym 
K3N2,3N2 Axb
 
3N1,3N1
3N1,3N   N,y   N,y  K  u  f 
FEM is to solve the large‐scale system of linear equations. Consider a cube. If the number of nodes along one edge is n, the total number of nodes is n3 and the number of equations is 3n3
3N,3N 3N3N  N,z 3N  N,z 3N

Example 2 : Linear regression
Problem: We want to find a function y = y(x(1), ∙∙∙, x(n))
E.g. y = lifetime of a material before failure x(1) = service temperature
x(2) = humidity of the environment
x(3), x(4) = hydrostatic and deviatoric stresses
x(5) = square of average grain size g2
x(6), ∙∙∙, x(p) = compositions of all species of impurities x(p+1), ∙∙∙ = (x(i))α, log(x(i)), sin(x(i)), exp(x(i)), ∙∙∙ (i = 1, ∙∙∙, p)
The simplest guess is the linear relationship
y a(0) a(1)x(1) a(2)x(2) a(n)x(n)
(Note: x(i) could be some term like xα , log(x), ∙∙∙ s.t. y is NOT linear with x) Now, the problem is to determine {a(i)}

Example 2 : Linear regression
y a(0) a(1)x(1) a(2)x(2) a(n)x(n) Now, the problem is how to determine {a(i)}
(ǂ)
Collect a set of data by measurement:
yy whenx(1) x(1),x(2) x(2),,x(n) x(n) 0000
yy whenx(1) x(1),x(2) x(2),,x(n) x(n) 1111

yy whenx(1) x(1),x(2) x(2),,x(n) x(n)
mmmm Determine {a(i)} s.t. all the data are consistent with Eq. (ǂ)
If Eq. (ǂ) is perfectly satisfied,
y 1x(1) x(n) a(0)
00 0 yy1  1 x(1)  x(n) a(1)
Xa
11
   
y 1x(1) x(n) a(n)
 m (m1)1  m m (m1)(n1)  (n1)1

Example 2 : Linear regression
y 1x(1) x(n) a(0) 000
yy1  ?1 x(1)  x(n) a(1) 11
Xa
 y 1x(1) x(n) a(n)
 m (m1)1  m m (m1)(n1)  (n1)1 However, usually m > n (the number of equations is larger than the number of unknowns); so, there is no solution for {a(i)} (i.e., a)
We need to find {a(i)} s.t. RHS is very close to LHS
Math language: minimize the cost function J w.r.t. a
J(a)1 Xay 2 1(Xay)(Xay)1a(XTX)aaXTy1yy 22222
J(a)  0  (XT X)a  XT y  0  (XT X)a  XT y a
XT X a  XT y Normal  (n1)(n1) (n1)1  (n1)1 equations
Axb

Solve the system of linear equations:
x1  2×2  3×3  1
2
2×1  2×2  x3  0
4
4×1  5×2  7×3  2 Forward
Gaussian elimination
elimination
x1  2×2  3×3  1
0  2×2  5×3  2 3
0  3×2  5×3  2 Forward
2
elimination
x  2x  3x  1 12315
02×2 5×3 2 x20
Backward substitution
005x1 x2 23 35
x 1

Gaussian elimination
Forward elimination (reduce the number of unknowns)
AxAxAxb A
111 122 133 1 21 A
A x  A x  A x  b A  31
211 222 233 2 11 AxAxAxb
A
31 1 32 2 33 3 3
AxAxAxb 11 1 12 2 13 3 1
21 12 21 13 21 1 0 AAAx AAAx bAb
22 A2 23 A3 2 A 1111 11
31 12 31 13 31 1 0 AAAx AAAx bAb
32 A2 33 A3 3 A 1111 11
11

Gaussian elimination
Forward elimination (reduce the number of unknowns)
AxAxAxb A
111 122 133 1 21 A
A x  A x  A x  b A  31
211 222 233 2 11 AxAxAxb
A
31 1 32 2 33 3 3
AxAxAxb
111 122A133A1b
 AA22  AA23 Ab2 0  A  21 12 x  A  21 13 x  b  21 1
22 A2 23 A3 2 A 11A11A 11b
 AA32  AA33 Ab3 0  A  31 12 x  A  31 13 x  b  31 1
32 A2 33 A3 3 A 1111 11
AxAxAxb 11 1 12 2 13 3 1
AxAxb 22 2 23 3 2
AxAxb 32 2 33 3 3
11

Forward elimination (reduce the number of unknowns)
AxAxAxb
11 1
12 2 13 3 1
AxAxb
222 233 2 A
AxAxb 32 2 33 3 3
32
AxAxAxb 111 122 133 1
AxAxb 22 2 23 3 2
32 23 32 2 0 AAAx bAb
33A3 3A 22 22
A
22

Forward elimination (reduce the number of unknowns)
AxAxAxb
11 1
12 2 13 3 1
AxAxb
222 233 2 A
AxAxb 32 2 33 3 3
32
AxAxAxb 111 122 133 1
AxAxb
222 233A2b
 AA33 Ab3 0  A  32 23 x  b  32 2
33A3 3A 22 22
AxAxAxb 11 1 12 2 13 3 1
AxAxb 22 2 23 3 2
Axb 333 3
A
22

11 1
12 2 13 3 1
Backward substitution
AxAxAxb
AxAxb
22 2
23 3 2
b
AxAxAxb
11 1
12 2 13 3 1 AxAxb x2 233
Axbx 333 3 3
3 33
222 233 2 2 Axb
bAx A
111 122 133 1 1 AxAxb
A
333 3
AxAxAxb xbAxAx
22 2 23 3 2
Axb 333 3
A
22
1 12 2 13 3
11

Forward elimination: Step (1) AxAxAxb
11 1 12 2 13 3 1
AxAxAxb 21 1 22 2 23 3 2
AxAxAxb 31 1 32 2 33 3 3
AxAxAxb
111 122A133A1b
 AA22  AA23 Ab2 0  A  21 12 x  A  21 13 x  b  21 1
22 A2 23 A3 2 A 11A11A 11b
 AA32  AA33 Ab3 0  A  31 12 x  A  31 13 x  b  31 1
32 A2 33 A3 3 A 1111 11
AxAxAxb 11 1 12 2 13 3 1
AxAxb 22 2 23 3 2
AxAxb 32 2 33 3 3

Forward elimination: Step (1) jm mk
A A A A (m1;j,k2,3) jk jk
AxAxAxb 11 1 A12 2 13 3 1
mm
AxAxAxb 21 1 A b22 2 23 3 2
bb jm m (m1;j2,3)
Aj x j  AA x  A x  b 31 1 32 2 33 3 3
mm
AxAxAxb
111 122A133A1b
 AA22  AA23 Ab2 0  A  21 12 x  A  21 13 x  b  21 1
22 A2 23 A3 2 A 11A11A 11b
 AA32  AA33 Ab3 0  A  31 12 x  A  31 13 x  b  31 1
32 A2 33 A3 3 A 1111 11
AxAxAxb 11 1 12 2 13 3 1
AxAxb 22 2 23 3 2
AxAxb 32 2 33 3 3

Forward elimination: Step (2) AxAxAxb
11 1 12 2 13 3 1
AxAxb 22 2 23 3 2
AxAxb 32 2 33 3 3
AxAxAxb 111 122 133 1
AxAxb
222 233A2b
 AA33 Ab3 0  A  32 23 x  b  32 2
33A3 3A 22 22
AxAxAxb 11 1 12 2 13 3 1
AxAxb 22 2 23 3 2
Axb 333 3

Forward elimination: Step (2)
AA
AxAxAxb
11jk1
A A  jm mk (m2;j,k3)
jk 12 2 13 3 1
A
A xmm  A x  b
22 2 23 3 2
Ab
Ajm mx  A x  b
(m2;j3) AxAxAxb
b b  jj
32 2 33 3 3
A mm
111 122 133 1
AxAxb
222 233A2b
 AA33 Ab3 0  A  32 23 x  b  32 2
33A3 3A 22 22
AxAxAxb 11 1 12 2 13 3 1
AxAxb 22 2 23 3 2
Axb 333 3

Forward elimination for N equations and N unknowns Step (m) (m = 1, ∙∙∙, N ‒ 1)
for m = 1, ···, N ‒ 1
for j = m + 1, ···, N
jm mk
A A A A (j,km1,,N)
jk jk
A mm
jm m bbA b
jj
mm
(jm1,,N) Algorithm
A
b[j] = b[j] – A[j][m] * b[m] / A[m][m]
for k = m + 1, ···, N
A[j][k] = A[j][k] – A[j][m] * A[m][k] / A[m][m]

Forward elimination for N equations and N unknowns Step (m) (m = 1, ∙∙∙, N ‒ 1)
def forward_elim(A, b):
N = b.size
A = A.astype(‘f’)
b = b.astype(‘f’)
for m in range(N-1):
# 0, …, N-2
# m+1, …, N-1
jm mk
A A A A (j,km1,,N)
jk jk
A mm
jm m bbA b
jj
mm
A
for j in range(m+1, N):
b[j] -= A[j,m] * b[m] / A[m,m]
for k in range(m+1, N):
# m+1, …, N-1
A[j,k] -= A[j,m] * A[m,k] / A[m,m]
return (A, b)
(jm1,,N) Python code

Backward substitution:
Step (3)
x
A
3 3
Step (2)
xbAx 2 233
2
22
Step (1)
xbAxAx 1 12 2 13 3
1
A
Step (m) (m = 2, 1)
x m
m
km1
mk k
b
A
3 bAx
33
11
A mm

Backward substitution for N equations and N unknowns Step (N)
x[N] = b[N] / A[N][N]
for m = N ‒ 1, ···, 1
x[m] = b[m]
for k = m + 1, ···, N
x[m] = x[m] – A[m][k] * x[k]
x[m] = x[m] / A[m][m]
x m
m
km1
mk k
x N
b N
Step (m) (m = N ‒ 1, ∙∙∙, 1)
Algorithm
A NN
N bAx
A mm

Backward substitution for N equations and N unknowns Step (N)
def backward_subs(A, b):
x N
b N
Step (m) (m = N ‒ 1, ∙∙∙, 1)
x m
m
km1
mk k
Python code
N = b.size; A = A.astype(‘f’); b = b.astype(‘f’) x = np.empty(shape=b.shape, dtype=’f’)
x[N-1] = b[N-1] / A[N-1, N-1]
for m in range(N-2, -1, -1):
return x
x[m] = b[m]
for k in range(m+1, N):
# N-2, N-3, …, 1, 0
# m+1, m+2, …, N-2, N-1
x[m] -= A[m,k] * x[k]
x[m] = x[m] / A[m,m]
A NN
N bAx
A mm

∙∙∙
Number of operations for input of size N = ?
A A 21 22
 A aAaAaA a N 2N 2 211 222 2NN
Computational complexity
Example: multiplication of a matrix with a vector
Assume that multiplication is much more time‐consuming than addition
A A  11 12
 A a AaAaA a N 1N  1   11 1 12 2 1N N 
N A A  A a  A a A a A a  N
     N1 N2 NN  N   N1 1 N2 2 NN N 
Computer needs to do N2 multiplications; complexity is O(N2) If 10‒6 sec per multiplication, we will get result after 10‒6N2 sec
For algorithm of Gaussian elimination,
Number of multiplication operations for Gaussian elimination on a matrix of size N N = ?

∙∙∙
A  11
A A
 A x  b  1N  1   1 
(N+1) (N+1) 21 A
12 13
A
A 31
(N + 1)
AAAAxb 21 22 23 2N22
11
A
A N1
A A A Axb 31 32 33 3N 3 3
11 
11
   
A
A A A Axb  N1 N2 N3 NN  N   N 
A A A Axb 11 12 13 1n11
Step(2):(N‒1)2 ‒1 ⁞
0 A A Axb 22 23 2n 2 2
Step (N ‒ 1): 12 ‒ 1 N2
0 A A Axb  32 33 3n  3   3 
 
(k1) k1
Computational complexity
Forward elimination
(N ‒ 1) Step(1):(N+1)(N‒1)=N2 ‒1
0 A A Axb
 N2 N3 NNN N(N‒1)(N‒1)


N2,N2
N2,N1 N2,N  N2   N2  N1,N1 N1,N  N1   N1 

A Axb Axb
 
NN  N   N  1
 b Step(N): xN
Computational complexity
Backward substitution
A  A A A xb
 11
1,N2 1,N1 1N  1   1     

A A Axb
N
x  b  N1,N N
Step (N ‒ 1): N1
N1 A A
Step(1):
1122 1NN x  b A x A x
N
N k1
k
1
AAA
A
NN Ax
2 Step(N‒2): x  N2  N2,N1 N1  N2,N N 3
N2
AAA
N1,N1 N1,N1 bAxAx
N2,N2
N2,N2
N2,N2
⁞
1111 11

Computational complexity
Number of multiplication operations for Gaussian elimination on a matrix of size N N =
N N N3 N (k2 1)k N2 
k1 k1
k, N = symbols(‘k N’, real=True, positive=True)
from sympy import *
pprint(summation(k**2 – 1 + k, (k, 1, N))) When N is large, only the highest order term matters
N3 O(N2) 3
Floating point operations per second (FLOPS)
For PC, about 1010 FLOPS per core about 10‒10 sec per multiplication If N = 104, time ≈ 30 sec
If N = 105, time ≈ 10 hr
If N = 106, time ≈ 1 yr
33

Matrix representation of Gaussian elimination
Step (1)
AxAxAxb A
111 122 133 1 21 A
A x  A x  A x  b A  31 211 222 233 2 11 A
AxAxAxb 11 31 1 32 2 33 3 3
1 00AAA   11 12 13
LAAA 10A A A 1 21 11 21 22 23
A A 0 1 A A A 3111 313233
AAA 1112 13
0AAAAAAAA  22 2112 11 23 2113 11
0AAAAAAAA
 32 3112 11 33 3113 11
AAA  11 12 13 
0AA  22 23
0AA
 32 33

AxAxAxb
11 1
12 2 13 3 1
1 
0A A A  11 12 13
Step (2)
AxAxb
222 233 2 A
AxAxb
32
32 2
33 3 3
A
0
LLA0 1 00AA
21   22 23 0 A A 1 0 A A
 32 22  32 33 AAA
1112 13 0AA
22 23 00AAAA
 33321322 AAA
 11 12 13  0AA
 22 23
00A
 33 
22

1 0 01 0 0A    11
A A  U 12 13  11
U U  12 13
LU decomposition L2L1 A  U
A A  0
  31 32 33  33
0 1 0a 1 0A
  21 22 23 22 23

U U  0c1b01AAA 00U
A A A1 0 0U U U  11 12 13   11 12 13
A A Aa 1 00 U U  21 22 23   22 23
A A A b c 1 0 0 U
 31 32 33 
 33
L: lower triangular matrix U: upper triangular matrix
A  LU

If A is a N  N matrix Gaussian elimination: Step (1)
L1A
 1 A A A  A 
 111213 1N A A 1 A A A  A 
21 11 21 22 23 2N A A 0 1 A A A  A 
3111 31 32 33 3N  
A A 0 0  1A A A  A 
N111 N1N2N3 AAAA
NN 
 11 12 13 1N  0AAA
2223 2N 0AAA
3233 3N 
0AAA N2N3 NN

L2L1A
0 1
0 A A 1
0AAA 2223 2N
If A is a N  N matrix Gaussian elimination: Step (2)
1 A A A  A   111213 1N
0 A A  A   32 33 3N
 32 22  
0A A 010 A A A  N2 22  N2 N3 NN
AAAA  11 12 13 1N 
0AAA 2223 2N
00AA  33 3N

00AA  N3 NN

If A is a N  N matrix Gaussian elimination: Step (3)
L3L2L1A
1 A A A  A 
 111213 1N 0 1 0 A A  A 
2223 2N 001 00AA
  33 3N  
0 0 A A  1 0 0 A  A   N333  N3 NN
AAAA  11 12 13 1N 
0AAA 2223 2N
00AA  33 3N

000A  NN 

LN1L3L2L1A 1
AAA A  11 12 1,N1 1N

0 1
0 A  A A 22 2,N1 2N
00
 1
00   A A
If A is a N  N matrix Gaussian elimination: Step (N ‒ 1)

 N,N1 N1,N1  N,N1 NN 

N1,N1 N1,N 00AA 100A A
AAAA  11 12 1,N1 1N 
0AA A 22 2,N1 2N
00U AA
 N1,N1 N1,N 
000A  NN 

1 
 UUUU  11 12 13 1N
LU decomposition LN1L3L2L1A U
A (LN1L3L2L1)1ULU
L1UUU 21 2223 2N
LLL1UUU 3132 333N
  
LLL1  N1 N2 N3  
U NN
A=
U L
Theorem: For a nonsingular matrix A, Gaussian elimination yields decomposition: A = LU
How can we obtain L and U?

Assume that A = LU
LU = U11
U12 L21U12 + U22
U13 L21U13 + U23
L21U11
L21U14 + U24 ∙∙∙∙ L31U14 + L32U24 + U34 ∙∙∙∙ L41U14 + L42U24 + L43U34 + U44 ∙∙∙∙
L31U11
L31U12 + L32U22 L41U12 + L42U22
L31U13 + L32U23 + U33 L41U13 + L42U23 + L43U33
LU decomposition: Doolittle algorithm
1  UUUU   11 12 13 14 
L1 UUU 21 22 23 24
LL L 1  U 31 32  
U U  33 34 
LLL1 414243 
U 44
    

U14 ∙∙∙∙
L41U11 ⁞⁞⁞⁞

U A (j1,,N) (1) 1j 1j
LU = U11
U12 L21U12 + U22
U13 L21U13 + U23
U14 ∙∙∙∙ L21U14 + U24 ∙∙∙∙ L31U14 + L32U24 + U34 ∙∙∙∙ L41U14 + L42U24 + L43U34 + U44 ∙∙∙∙
L21U11
L31U11
L31U12 + L32U22 L41U12 + L42U22
L31U13 + L32U23 + U33 L41U13 + L42U23 + L43U33
LU decomposition: Doolittle algorithm
LUALAU (i2,,N) i111 i1 i1 i1 11
L41U11 ⁞⁞⁞⁞

U A (j1,,N) (1) 1j 1j
LU = U11
U12 L21U12 + U22
U13 L21U13 + U23
U14 ∙∙∙∙ L21U14 + U24 ∙∙∙∙ L31U14 + L32U24 + U34 ∙∙∙∙ L41U14 + L42U24 + L43U34 + U44 ∙∙∙∙
L21U11
L31U11
L31U12 + L32U22 L41U12 + L42U22
L31U13 + L32U23 + U33 L41U13 + L42U23 + L43U33
L A U (i2,,N) i1 i1 11
i1 12 i2 22
i2 i2 i2 i1 12
22
LU decomposition: Doolittle algorithm
LU U A U A LU (j2,,N) (2)21 1j 2j 2j 2j 2j 21 1j
LU LU A L (A LU )U (i3,,N)
L41U11 ⁞⁞⁞⁞

U A (j1,,N) (1) 1j 1j
U11
U12 L21U12 + U22
U13 L21U13 + U23
U14
L21U11
L31U11
L31U12 + L32U22 L41U12 + L42U22
L31U13 + L32U23 + U33 L41U13 + L42U23 + L43U33
L A U (i2,,N) i1 i1 11
LU decomposition: Doolittle algorithm
U A L U (j2,,N) (2)2j 2j 211j
L (A LU )U (i3,,N) i2 i2 i112 22
LU LU U A U A LU LU (j3,,N) (3)31 1j 32 2j 3j 3j 3j 3j 31 1j 32 2j
LU LU LU A L (A LU LU )U (i4,,N)
i1 13 i2 23 LU =
i3 33
i3 i3
i3
i1 13
i2 23 33
L41U11 ⁞⁞⁞⁞
∙∙∙∙ L21U14 + U24 ∙∙∙∙ L31U14 + L32U24 + U34 ∙∙∙∙ L41U14 + L42U24 + L43U34 + U44 ∙∙∙∙



mj mj
 k  1
mk kj

LU decomposition: Doolittle algorithm
m1 

im im
ikkm mm

 k1 
(im1,,N)
LA LUU
Step (1)
U A (j1,,N) 1j 1j
L A U (i2,,N) i1 i1 11
Step (m)
(m = 2, ∙∙∙, N ‒ 1)
m1
U A  L U (jm,,N)
Step (N) N1
UALU
NN NN
Nk kN
 k1

Solving system of linear equations by LU decomposition
Axb  LUxb 

L
L
U
x=z
A x=b
U
x=b
z=b
Lzb Ux  z

Solving system of linear equations by LU decomposition
Forward substitution
Axb  LUxb 

1 z b
zb 11
 1   1  21 22

L 1 zb 
z b L z 2 2 211
LL1zb
N1 zbLz
N1 N2 Backward substitution

NN UUUxz
NN
Nk k k1
1kk 11 11 12 1N11 

2N  2  2 
 k3 

 
2kk 22

   
UUxz N
22
x2z2 UxU
 NNNN xzU
N N NN
Uxz
Lzb Ux  z
N x1z1 UxU
 k2 

Solving system of linear equations by LU decomposition
Forward substitution
Axb  LUxb 

1 z b
zb 11
 1   1  21 22

L 1 zb 
z b L z 2 2 211
LL1zb N1 N2 NN
N1 zbLz
def forward_subs(A, b):
N = b.size; A = A.astype(‘f’); b = b.astype(‘f’) x = np.empty(shape=b.shape, dtype=’f’)
x[0] = b[0] / A[0,0]
for m in range(1, N):
return x
x[m] = b[m]
for k in range(m):
# 1, …, N-1
# 0, …, m-1
x[m] -= A[m,k] * x[k]
x[m] /= A[m,m]
Lzb Ux  z
NN
Nk k k1


Why LU decomposition?
Ax  b
Sometimes we want to solve x for same A but different b. E.g. A describes materials property, x is response and b is external force; how does the material respond to different external forces?
Without use of LU, employ Gaussian elimination every time b is changed
Axb  LUxb 

By LU decomposition, only need to do the decomposition once. Then, L and U are fixed. For different b, we just need to solve 2 equations simply by substitution
Lzb Ux  z

LU decomposition: Doolittle algorithm
mj mj mkkj  k  1
m1  LA LUU

im im
ikkm mm

 k1 
(im1,,N)
Step (1)
U A (j1,,N) 1j 1j
L A U (i2,,N) i1 i1 11
It is not guaranteed thatallUmm ≠0! Pivoting
Step (m)
(m = 2, ∙∙∙, N ‒ 1)
 m1
U A L U (jm,,N)
Large Umm reduces roundoff error
P i v o t i n g
Step (N)
N  1 UALU
NN NN
k1
Nk kN

Pivoting
AxAxAxb 11 1 12 2 13 3 1
Exchange 2 equations
AxAxAxb 21 1 22 2 23 3 2
AxAxAxb 31 1 32 2 33 3 3
AxAxAxb 21 1 22 2 23 3 2
AxAxAxb 11 1 12 2 13 3 1
AxAxAxb 31 1 32 2 33 3 3
Reordering 2 equations does not change the solution Pivoting means such reordering

Solve the system of equations x1  2×2  3×3  1
2x  5x  x  0 1223
4×1  5×2  7×3  2
Idea of pivoting: an example

x1  2×2  3×3  1 2x  5x  x  0
1223
4×1  5×2  7×3  2
Idea of pivoting: an example

x1  2×2  3×3  1 2x  5x  x  0
1223
4×1  5×2  7×3  2
4×1  5×2  7×3  2 2x  5x  x  0
1223
x1  2×2  3×3  1
Idea of pivoting: an example
Pivoting

x1  2×2  3×3  1 2x  5x  x  0
1223
4×1  5×2  7×3  2
Idea of pivoting: an example
Pivoting
4×1  5×2  7×3  2 2x5xx0 4 1
2 12234
x1  2×2  3×3  1

x1  2×2  3×3  1 2x  5x  x  0
1223
4×1  5×2  7×3  2
Idea of pivoting: an example
Pivoting
 2 x  2x  3x  1 0  3x  5x  1
4×1  5×2  7×3  2 2x  5x  x  0
4×1  5×2
1223 23
1 2 3 Forward 42 43 2
elimination
 7×3
0  0  5x  1

4×1 5×2 7×3 2 2x  5x  x  0
4×1 5×2 7×3 0  0  5x
2  1
Idea of pivoting: an example
x1  2×2  3×3  1 4×1  5×2  7×3  2 2x  5x  x  0 0  3x  5x  1
1223 42432
4x  5x  7x  2 0  0  5x  1 12323
Pivoting
Pivoting
1223 23 x  2x  3x  1 0  3x  5x
 1 2
1 2 3 Forward 42 43 elimination

 3×3  1 2x  5x  x  0
4×1  5×2  7×3  2 0  3x  5x  1
Idea of pivoting: an example
x1  2×2
1223 42432
Solution: x2, x0, x1 35215
4x  5x  7x  2 0  0  5x  1 12323
4×1 5×2 7×3 2 2x  5x  x  0
4×1 5×2 7×3 0  0  5x
2  1
Pivoting
Pivoting
1223 23 x  2x  3x  1 0  3x  5x
 1 2
1 2 3 Forward 42 43 elimination
Backward substitution

Pivoting can be represented by permutation matrix P
0 0 1a a   11 12
a  a a 13   31 32
a  33 
1 0 0a a
a a a 23   11 12
a  13 

0 1
 21 22 0a31 a32
a23 
P
a33  a21 a22 0 0 1b b
  1   3  1 0 0b b 
 
 2   1  010b b
 3   2 
Permutating and then permutating back = no operation: PTP = I
P
Repeatedly exchanging multiple rows = permutating once: P1P2∙∙∙ = P LU decomposition with pivoting
Pivoting  LU decomp. of A Lz  PT b A x    P A x  b      P L U x  b   U x  z


A = np.array([[2, 5, 8, 7],
[5, 2, 2, 8],
[7, 5, 6, 6],
[5, 4, 4, 8]])
b = np.array([[1],
print(P); print(L); print(U)
# obtain the solution x
[0],
[1],
[0]])
Pivoted LU decomposition by SciPy
# pivoted LU decomposition: A = P L U
P, L, U = scipy.linalg.lu(A)
scipy.linalg.lu()
Pivoting  LU decomp. of A Lz  PT b A x    P A x  b      P L U x  b   U x  z
z = forward_subs(L, np.matmul(np.transpose(P), b)) print(backward_subs(U, z))


A = np.array([[2, 5, 8, 7],
[5, 2, 2, 8],
[7, 5, 6, 6],
Solve system of linear equations by SciPy
[5, 4, 4, 8]])
b = np.array([[1],
[0],
[1],
[0]])
scipy.linalg.solve()
SciPy provides function to solve system of linear equations by LU decomposition
print(scipy.linalg.solve(A, b))
# compare with the solution from other methods print(np.matmul(np.linalg.inv(A), b)) print(backward_subs(*forward_elim(A, b)))
# * : unpack the tuple

QR decomposition Any real square matrix A can be decomposed as
A  QR Q : orthogonal matrix: QTQ = QQT = I
R : upper triangular matrix (also called “right” triangular matrix)
A=Q
Any complex rectangular M  N (M > N) matrix A can be decomposed as
AMN QMMRMN Q : unitary matrix: Q†Q = QQ† = I
R : upper triangular matrix (also called “right” triangular matrix) R
A=Q
R

Assume that QR decomposition can be achieved
QR decomposition
Ax  b
QRx bRx QTb
Since R is upper triangular, x can be obtained by backward substitution Assume that QR decomposition can be achieved. True?
A x=b
QRx=b
R
x=QT b

QR decomposition
       
    0    AHA     1 0   
    0    
       
0    0    RHHHA  HHA  3 2 1 0 0   2 1 0 0  
0 0 0  0 0   
HN 1 H2H1 A  R  A  (HN 1 H2H1 )1 R  QR Q(H HH)1H1H1H1
N1 21 12 N1
But we still don’t know (1) how to construct these {Hk} and (2) why Q is
an orthogonal matrix

Householder reflector
Hzz2e1, Hz20 z0
3 Householder reflector H
transform a vector z to ‖z‖2e1, a vector of same length
e2
z  zTz 2
z
z1  z z zz
z1 z  2
Hz
e1

3 z2  z2  z2 (length)
1232 z
123

Householder reflector
Hz  z 2 e1
H reflects z about a mirror
Difference between z and Hz is v
e2
z
Hz
e1

Householder reflector
Hz  z 2 e1
Projectztov  v  v
2ˆ2ˆˆˆ Projztovz cosθvz zvvz 
e2

zv v vTz
v vv vTv
z22 θ
Hz
e1
vv 22

Householder reflector
e2
21 2 3 12 22 23 v  v v v v v v 
Hz  z 2 e1
Difference between z and Hz
vTz vvT  vvT vvT
Hzz2 T vz2 T zI2 T zHI2 T vvvvvvvv
v1  v1v1 v1v2 vvT v v v v v v v v
v1v3  v v 
3 13 23 33
z
2 vT z v vT v
Hz
e1

Hz = –‖z‖2e1
e1
Householder reflector Alternatively, Hz   z 2 e1
In this case, need to calculate z1 + ‖z‖2; since z1 > 0, this is addition. In previous case, need to calculate z1 ‒ ‖z‖2; since z1 > 0, this is subtraction. Avoid calculating difference of 2 floating‐point numbers which are close (why?)
e2
Hzsgn(z1) z 2 e1 v  z  sgn(z1 ) z 2 e1
z
HI2vvT vTv

AAA  11 12 13 
AAA  11 12 13 
AA A A  21 22 23 
HA0AA 1  22 23
AAA 11 12 13
 vT v 
A11 A12 A13 
QR decomposition
AAA 0AA
 31 32 33  a1
 32 33 a2
vasgn(A)a e 111121
vasgn(A)a e 2 22221
H I2vvT 1 vT v
1 0 
HHA0AAR HA0AA 21 2223
12223 
00A
0AA 33  32 33

H2 0 I2vvT  
Q  (H H )1  H1H1 2112

QR decomposition
Q  (H H )1  H1H1 2112
Check if Q is orthogonal. Property of Householder reflector I0I0
H vvT HT  vvT H 0 I2  0 I2 
  vT v   vT v 
I 0I 0I0 HHT  vvT  vvT  I
0 I2 0 I2  0 I
  
vT v  vT v  1 T

T  1 1   1 1  1 1
Q Q H1 H2 H1 H2 H2H1H1 H2 I So, Q is indeed an orthogonal matrix
HHH T

k kk1
At each step, we can construct Hk
QR decomposition
Iterative construction of R and QT (finally we need QT not Q) T
After iteration,
Solve the equations Rx = QTb by backward substitution
Algorithm
T11 1 RHN1H2H1A, Q  H1 H2 HN1 HN1H2H1I
R A QT I
0 ,0 (k1,,N1)
R HR QT HQT
k kk1 Ik1
0  T
vasgn(A)a e,H
k kkk21 k0I2vv
R:R ,QT:QT N1 N1
 Nk1 
 vT v 

Example Do QR decomposition on the matrix
QR decomposition
Step (0) 1 1 1
1 0 0 0
1 1 1
1 3 3 A  1 1 5
1 3 7 
0 1 0 0 R A , QT I  0 115 0 0010
1 3 3
1 3 7 0 0 0 1 

1 3 3
0 1 1 5
1 11
1 2 1
1 0 1 1 0 1
1 3 7  
 1 1 1 1
QR decomposition
Step (1)
R ,vasgn(A)a e 2 
1 1 1
1 1 3 1 0 1
a1
HI2vvT I213 1 1 11 5 1 1
3 2222
1 4 vTv 4 121 1 1 5 1  2666
11115 2666
242 1 1 1 1 2222
33TT2666 RHR0 4 8,QHQ1 5 1 1
1 1002161 101 1 5 1 332666
04 20 1 1 1 5 332666
2666

R0 = Matrix([[-1, -1, 1],
[ 1, 3, 3],
[-1, -1, 5],
[ 1, 3, 7]])
QT0 = eye(R0.shape[0])
QR decomposition
We can do the calculation by SymPy
# .shape returns shape of matrix (number of rows, number of cols) # eye(N): identity matrix of size NxN
a1 = R0[:, 0] # access a block of matrix
v = a1
v[0,0] = v[0,0] + np.sign(a1[0,0]) * a1.norm()
# .norm() returns Frobenius norm of matrix or L2 norm of vector print(“v = “); pprint(v)
H1 = eye(v.shape[0]) – 2 * v * v.T / (v.dot(v)) # .dot() dot product; returns a scalar print(“H1 = “); pprint(H1)
R1 = H1 * R0
print(“R1 = “); pprint(R1)
QT1 = H1 * QT0
print(“QT1 = “); pprint(QT1)

QR decomposition
Step (2)
048 3 3
2 4 2
R  3 3 , va sgn(A )a e 2202 
4 1 10 10216 2 2222133
33 44 0 4 20 3 03
33 a2
10 212 vvT3 3 3
H T,I2 1 14 2 20I2vv3 T315 15
3 vv2211  vTv 3 15 15
242 1111 2222
TT2222 R HR 0 2 8, Q HQ 1 1 1 1
2 2100162 2171 7 1  5 10 10 10 10
0 0 12 1 7 1 7 5 10101010

0 2 8 20016 3
16 1 36 3332112012
 5 0 0 12
55
5 a3
QR decomposition
Step (3)
R ,vasgn(A)a e5 4 5
2 4 2
I2 0 
H ,I2I5361255
vvT
5 36 4 3
3 4 555
3 0I2vvT2 vTv2 144125 5  2 vTv
242 1 1 1 1 2222
TT2222 R HR 0 2 8, Q HQ 1 1 1 1
3 320043 321 111 2222
000 1111 2222

A  QR 
11 22
1 10 2 8 22
 Q(1)
R(1)  Q(2) 
 Q(1)R(1)
11 2 2
1 1242 2 2 
1 1 2 2
1 10 0 4 2 2 
0 
QR decomposition
Result
242 1111 2222
TT2222 RR 0 2 8, Q Q 1 1 1 1
3004 31 111 2222
000 1111 2222
1111000 2 2 2 2 
Note that QR decomposition can be performed on rectangular M  N matrices A with M > N

import scipy.linalg
import numpy as np
A = np.array([[-1, -1, 1],
[ 1, 3, 3],
[-1, -1, 5],
# QR decomposition: A = Q R
Q, R = scipy.linalg.qr(A)
QR decomposition by SciPy
[ 1, 3, 7]])
print(Q); print(np.matmul(np.transpose(Q), Q)) print(R)
print(np.matmul(Q, R) – A)
scipy.linalg.qr()
11111 1 1242  2 2 2 2 
A1 3 3 1 1 1 10 2 8QR 2222
1 1 5 1 1 1 10 0 4  2 2 2 2 
1371111000
 2 2 2 2 

Solve system of linear equations by QR decomposition Ax QRx bRx QTb
import numpy as np
import scipy.linalg
from Gaussian_elimination import forward_subs, backward_subs
A = np.array([[2, 5, 8, 7],
[5, 2, 2, 8],
[7, 5, 6, 6],
[5, 4, 4, 8]])
b = np.array([[1],
# solve by LU decomposition
[0],
[1],
[0]])
P, L, U = scipy.linalg.lu(A)
z = forward_subs(L, np.matmul(np.transpose(P), b)) print(backward_subs(U, z))
# solve by QR decomposition
Q, R = scipy.linalg.qr(A)
print(backward_subs(R, np.matmul(np.transpose(Q), b)))

QR vs. LU
If A is of shape MN and M > N, the computational complexity (count of
multiplications) of QR is 2MN2/3 For square A, the complexity is:
LU: N3 O(N2), QR: 2N3 O(N2) 33
Both QR and LU are used to solve system of linear equations. QR requires twice the cost of LU. Why do we need QR?
We note that QR can be applied to rectangular A (LU only for square A)
If A is of shape MN and M > N, we have more equations than unknowns. The equality is not valid. The problem is to find x such that ‖Ax ‒ b‖ reaches minimum (linear regression). QR can be used to do linear regression. We will learn least squares fitting in Lecture 5

QR vs. LU
If A is of shape MN and M > N, we have more equations than unknowns. The equality is not valid. The problem is to find x such that ‖Ax ‒ b‖ reaches minimum (linear regression). QR can be used to do linear regression. Start from the normal equations:
AT Ax  AT b
T 1 T AQR  T  T
x(A A) A b(QR) (QR)1(QR) b  RT QT QR1 RT QT b  RT R1 RT QT b
 R1RTRTQTbR1QTb
 Rx  QT b
Recall that R is an upper (right) triangular matrix. We can solve this system of equations simply be backward substitution

System of linear equations
Why solving system of linear equations is important?
Finite element method, linear regression, etc.
How to solve a system of linear equations?
(i) Gaussian elimination (ii) LU decomposition
Direct method
Ax  b
Matrix form of Gaussian elimination
Pivoted LU decomposition (iii) QR decomposition

U A (j1,,N) 0 (1) 1j 1j
N (k1)(Nk1)
(3)3j 3j 31 1j 32 2j
i3 LU =
i3 i1 13
i2 23
33
U11 L21U11 L31U11
U12 L21U12 + U22
U13 L21U13 + U23
U14 L21U14 + U24
∙∙∙∙
∙∙∙∙
∙∙∙∙
∙∙∙∙
L31U12 + L32U22 L41U12 + L42U22
L31U13 + L32U23 + U33
L31U14 + L32U24 + U34
LU decomposition: Doolittle algorithm
L A U (i2,,N) 0 i1 i1 11
k1 N
U A L U (j2,,N) (N1)
(k1)(Nk) k1
(2)2j 2j 211j
L (A L U ) U
(i3,,N) (N2)
i2 i2 i1 12 22 U A L U L U
N3N2N
(j3,,N) 2(N2)
L (A L U L U ) U (i4,,N) 2(N3)
3
2
6
L41U11 ⁞⁞⁞⁞
L41U13 + L42U23 + L43U33
L41U14 + L42U24 + L43U34 + U44

Cholesky decomposition Ax  b
If A is symmetric and positive‐definite, what is LU for A? 1 0 0U U U
   11 12 13  AL 1 00 U U
 21   22 23  LL100U
 31 32   33 
1 00U 0 01UU UU
 11  12 11 13 11 L 100U 00 1 UU
21 22 2322 LL100U001
313233 
AAT 
1 00U001LL
L D LT
T  11  21 31
AUU 1 00U 001L 12 11  22  32
UUUU100U001 13 11 23 22  33 

11  AL 1 00 U 0L 1 0
1 L
0 0 U11 0 1 00 U
0 00 U
0 T1 0L
0 0T 1 0
 22 L L 1
 22
21
 
21
 
 
L L 1
 
Cholesky decomposition
A  LDL
1 0 0U 0 01 0 0T
21  22 21  LL100ULL1
 31 32  Aispositive‐definiteUii 0
33  31 32  0  U11
 31 32  0 0
R RT
U33  0 Cholesky decomposition: A = RRT
0
32
A=
RT R
T
U33   31

Assume that A = RRT R11
 R11 R21 R31 R41  RRR
RRT = R112
R11R21 R212 + R222
R11R31 R21R31 + R22R32
R11R41 ∙∙∙∙ R21R41 + R22R42 ∙∙∙∙
R11R21
R11R31
R21R31 + R22R32 R21R41 + R22R42
R312 + R322 + R332 R31R41 + R32R42 + R33R43
R31R41 + R32R42 + R33R43 ∙∙∙∙ R412 + R422 + R432 + R442 ∙∙∙∙
RR 21 22
223242 RT RR
RR R R  31 32 33
 3343   R44 
Cholesky decomposition
R41 R42 R43 R44
 

R11R41 ⁞⁞⁞⁞

R  A (1) 11 11
Cholesky decomposition
RRT = R112
R11R21 R212 + R222
R11R31 R21R31 + R22R32
R11R41 ∙∙∙∙ R21R41 + R22R42 ∙∙∙∙
R11R21
R11R31
R21R31 + R22R32 R21R41 + R22R42
R312 + R322 + R332 R31R41 + R32R42 + R33R43
R31R41 + R32R42 + R33R43 ∙∙∙∙ R412 + R422 + R432 + R442 ∙∙∙∙
R A , R A , , R A 21R31R N1R

21 31 N1
11 11 11
R11R41 ⁞⁞⁞⁞

R  A (1) 11 11
Cholesky decomposition
RRT = R112
R11R21 R212 + R222
R11R31 R21R31 + R22R32
R11R41 ∙∙∙∙ R21R41 + R22R42 ∙∙∙∙
R11R21
R11R31
R21R31 + R22R32 R21R41 + R22R42
R312 + R322 + R332 R31R41 + R32R42 + R33R43
R31R41 + R32R42 + R33R43 ∙∙∙∙ R412 + R422 + R432 + R442 ∙∙∙∙
R A , R A , , R A 21R31R N1R
21 31 N1

R  A R2
(2) 22 22 21
11 11 11
32R21 31 42 42R21 41 N2 N2R21 N1 RARR ,RARR ,,RARR
32
22 22 22
R11R41 ⁞⁞⁞⁞

R  A (1) 11 11
Cholesky decomposition
R A , R A , , R A 21R31R N1R
21 31 N1

R  A R2
11 11 11
(2) 22 22 21
32 32R21 31 42 42R21 41 N2 N2R21 N1
RARR ,RARR ,,RARR 22 22 22
R  A R2 R2 (3) 33 33 31 32
RRT = R112
R11R21 R212 + R222
R11R31 R21R31 + R22R32
R11R41 ∙∙∙∙ R21R41 + R22R42 ∙∙∙∙
R11R21
R11R31
R21R31 + R22R32 R21R41 + R22R42
R312 + R322 + R332 R31R41 + R32R42 + R33R43
R31R41 + R32R42 + R33R43 ∙∙∙∙ R412 + R422 + R432 + R442 ∙∙∙∙
RARRRR ,RARRRR ,,RARRRR 43 31 41 32 42 53 31 51 32 52 N3 31 N1 32 N2
43 R 53 R N3 R

33 33 33
R11R41 ⁞⁞⁞⁞

R  A
(1) 11 11 fork1,,N
(2) 22 22 21
kk kk
kj
RRT = R112
R11R21 R212 + R222
R11R31 R21R31 + R22R32
R11R41 ∙∙∙∙ R21R41 + R22R42 ∙∙∙∙
R11R21
R11R31
R21R31 + R22R32 R21R41 + R22R42
R312 + R322 + R332 R31R41 + R32R42 + R33R43
R31R41 + R32R42 + R33R43 ∙∙∙∙ R412 + R422 + R432 + R442 ∙∙∙∙
Cholesky decomposition
R A , R A , , R A
21 R 31 R N1 R k1
21 31 N1

R  A R2
 j1
11 11 11
32 32R21 31 42 42R21 41 N2 N2R21 N1 RARR ,RARR ,,RARR
R  A  R2
 k1  22 22 22
ik ik 
kj ji kk RARRRR ,RARRRR ,,RARRRR
R  A R2 R2 (3) 33 33 31 32
43 31 41 32 42 53 31 51 32 52 N3 31 N1 32 N2 43 R 53 R N3 R

33 33 33
R11R41 ⁞⁞⁞⁞

j1 
RA RRT R (ik1,,N)

R  A 0 (1) 11 11
Cholesky decomposition
RRT = R112
R11R21 R212 + R222
R11R31 R21R31 + R22R32
R11R41 ∙∙∙∙ R21R41 + R22R42 ∙∙∙∙
R11R21
R11R31
R21R31 + R22R32 R21R41 + R22R42
R312 + R322 + R332 R31R41 + R32R42 + R33R43
R31R41 + R32R42 + R33R43 ∙∙∙∙ R412 + R422 + R432 + R442 ∙∙∙∙
R A , R A , , R A 0
Number of multiplication operations=?
21R31R N1R 11 11 11
21 31 N1
R  A R2 1 (2) 22 22 21
R A R R , R A R R , , R A R R 32 32R21 31 42 42R21 41 N2 N2R21 N1
(N2) RARRRR ,RARRRR ,,RARRRR 2(N3)
22 22 22
R  A R2 R2 2 (3) 33 33 31 32
43 31 41 32 42 53 31 51 32 52 N3 31 N1 32 N2 43 R 53 R N3 R

33 33 33
R11R41 ⁞⁞⁞⁞

A = np.array([[12, 5, 1, 7],
[5, 12, 2, 8],
[1, 2, 16, 6],
[7, 8, 6, 18]])
b = np.array([[1],
print(R)
# obtain the solution x
[0],
[1],
[0]])
Cholesky decomposition by SciPy
scipy.linalg.cholesky() scipy.linalg.solve()
Axb  RRTxb 

# Cholesky decomposition: A = R R.T
R = scipy.linalg.cholesky(A, lower=True)
z = forward_subs(R, b)) print(backward_subs(np.transpose(R), z))
# directly solve the equations print(scipy.linalg.solve(A, b, assume_a=”pos”))
# assume_a=”pos” denotes that A is positive-definite
Rzb RT x  z

Exchange
LU =
L UA21 21 11
L U +U A22 21 12 22
L U +U 21 13
A23 23
L U +U A24 ∙∙∙∙ 21 14 24
L UA31 31 11
L U +L UA32 31 12 32 22
+U A33 31 13 32 23 33
L U +L U 31 14 32 24
+U A34 ∙∙∙∙ 34
max |A11| Pivoting
U A (j1,,N) (1) 1j 1j
LUALAU (i2,,N) i111 i1 i1 i1 11
UA11UA12 UA13 UA14∙∙∙∙ 11 12 13 14
L U +L U
LUA41 LU +LUA42 LU +LU +LUA43 LU +LU +LU +UA44 ∙∙∙∙
41 11 41 12 42 22 41 13 42 23 43 33 41 14 42 24 43 34 44
⁞⁞⁞⁞
The row with largest |Ap1|

LU =
max |A11| Pivoting
U A (j1,,N) (1) 1j 1j
LUALAU (i2,,N) i111 i1 i1 i1 11
UA11UA12 UA13 UA14∙∙∙∙ 11 12 13 14
L UA21 21 11
L U +U A22 21 12 22
L U +U 21 13
A23 23
L U +U A24 ∙∙∙∙ 21 14 24
L UA31 31 11
L U +L UA32 31 12 32 22
+U A33 31 13 32 23 33
L U +L U 31 14 32 24
+U A34 ∙∙∙∙ 34
L U +L U
LUA41 LU +LUA42 LU +LU +LUA43 LU +LU +LU +UA44 ∙∙∙∙
41 11 41 12 42 22 41 13 42 23 43 33 41 14 42 24 43 34 44 ⁞⁞⁞⁞

Exchange
LU =
max |A22|
LU U A U A LU (j2,,N)
(2)21 1j 2j 2j 2j 2j 21 1j
LU LU A L (A LU )U (i3,,N)
i1 12 i2 22 i2 i2 i2 i1 12 22
L UA21 21 11
L U +U A22 21 12 22
L U +U 21 13
A23 23
L U +U A24 ∙∙∙∙ 21 14 24
L UA31 31 11
L U +L UA32 31 12 32 22
+U A33 31 13 32 23 33
L U +L U 31 14 32 24
+U A34 ∙∙∙∙ 34
Pivoting
UA11UA12 UA13 UA14∙∙∙∙ 11 12 13 14
L U +L U
LUA41 LU +LUA42 LU +LU +LUA43 LU +LU +LU +UA44 ∙∙∙∙
41 11 41 12 42 22 41 13 42 23 43 33 41 14 42 24 43 34 44
⁞⁞⁞⁞
The row with largest |Ap2|

max |A22|
LU U A U A LU (j2,,N)
(2)21 1j 2j 2j 2j 2j 21 1j
LU LU A L (A LU )U (i3,,N)
LU =
i1 12 i2 22 i2 i2 i2 i1 12 22
Pivoting
UA11UA12 UA13 UA14∙∙∙∙ 11 12 13 14
L UA21 21 11
L U +U A22 21 12 22
L U +U 21 13
A23 23
L U +U A24 ∙∙∙∙ 21 14 24
L UA31 31 11
L U +L UA32 31 12 32 22
+U A33 31 13 32 23 33
L U +L U 31 14 32 24
+U A34 ∙∙∙∙ 34
L U +L U
LUA41 LU +LUA42 LU +LU +LUA43 LU +LU +LU +UA44 ∙∙∙∙
41 11 41 12 42 22 41 13 42 23 43 33 41 14 42 24 43 34 44 ⁞⁞⁞⁞

L UA21 21 11
L U +U A22 21 12 22
L U +U 21 13
A23 23
L U +U A24 ∙∙∙∙ 21 14 24
L UA31 31 11
L U +L UA32 31 12 32 22
+U A33 31 13 32 23 33
L U +L U 31 14 32 24
+U A34 ∙∙∙∙ 34
Pivoting
max |A33|
LU LU U A U A LU LU (j3,,N)
(3)31 1j 32 2j 3j 3j 3j 3j 31 1j 32 2j
LU LU LU A L (A LU LU )U (i4,,N)
i1 13 i2 23 i3 33 i3 i3 i3 i1 13 i2 23 33 LU =
UA11UA12 UA13 UA14∙∙∙∙ 11 12 13 14
L U +L U
LUA41 LU +LUA42 LU +LU +LUA43 LU +LU +LU +UA44 ∙∙∙∙
41 11 41 12 42 22 41 13 42 23 43 33 41 14 42 24 43 34 44
⁞⁞⁞⁞
The row with largest |Ap3|

Pivoting
max |A33|
LU LU U A U A LU LU (j3,,N)
(3)31 1j 32 2j 3j 3j 3j 3j 31 1j 32 2j
LU LU LU A L (A LU LU )U (i4,,N)
i1 13 i2 23 i3 33 i3 i3 i3 i1 13 i2 23 33 LU =
UA11UA12 UA13 UA14∙∙∙∙ 11 12 13 14
L UA21 21 11
L U +U A22 21 12 22
L U +U 21 13
A23 23
L U +U A24 ∙∙∙∙ 21 14 24
L UA31 31 11
L U +L UA32 31 12 32 22
+U A33 31 13 32 23 33
L U +L U 31 14 32 24
+U A34 ∙∙∙∙ 34
L U +L U
LUA41 LU +LUA42 LU +LU +LUA43 LU +LU +LU +UA44 ∙∙∙∙
41 11 41 12 42 22 41 13 42 23 43 33 41 14 42 24 43 34 44 ⁞⁞⁞⁞