[SOLVED] 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: 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

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

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
teelementmethod
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 1u(1) daT u(1)
(1) u(1) x(1)
a (2) u(2) b (3) u(3)
Determine c0 and c1:
ccxu 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

xx(1) u(1) naT (x)u(1)
(2) (1) a 1xxL
a (2) (2) L Lu u
a a (2) (2) dx LLuu
(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) 2TdxA (1) 2 u(1)
a a a a aT a a K A L E d d A L E

aaa kAE 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 1LL
aaaa a
Ub 2 u(2) u(3) Kbu , 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 UUa Ub 2 u(1)
u(1) 1 u(2) Kau(2)2 u(2)
x(3)
u(2)
ka ka 0u(1) 1
2 u(1) u(2) u(3) ka ka 0u(2)
0 0 0u(1) 1 (1) (2) (3) b b (2)

2 u u u 0 k k u

u(3) Kbu(3)
0 0 0u(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:
0Kuf0 Kuf
u k
uiKijuj ij
ka ka 0u(1)f(1) a a b b(2)(2)
u uiKikuk ukKkjuj
kkkkuf 0 kb kbu(2)f(2)
kij uK Ku2 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) fka
ka ka kb kb
x(3) kbu(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 0u(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,3N2 1,3N1 1,3N 1,x 1,x
KKK K Kuf 22 23 2,3N2 2,3N1 2,3N 1,y 1,y
KK K Kuf 33 3,3N2 3,3N1 3,3N 1,z 1,z

K3N2,3N1 K3N2,3N uN,x fN,x K Kuf
sym
K3N2,3N2 Axb

3N1,3N1
3N1,3N N,y N,y K u f
FEM is to solve the largescale 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 3N3N 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:
yy whenx(1) x(1),x(2) x(2),,x(n) x(n) 0000
yy whenx(1) x(1),x(2) x(2),,x(n) x(n) 1111

yy 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)
00 0 yy1 1 x(1) x(n) a(1)
Xa
11

y 1x(1) x(n) a(n)
m (m1)1 m m (m1)(n1) (n1)1

Example 2 : Linear regression
y 1x(1) x(n) a(0) 000
yy1 ?1 x(1) x(n) a(1) 11
Xa
y 1x(1) x(n) a(n)
m (m1)1 m m (m1)(n1) (n1)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 Xay 2 1(Xay)(Xay)1a(XTX)aaXTy1yy 22222
J(a) 0 (XT X)a XT y 0 (XT X)a XT y a
XT X a XT y Normal (n1)(n1) (n1)1 (n1)1 equations
Axb

Solve the system of linear equations:
x1 22 33 1
2
21 22 x3 0
4
41 52 73 2 Forward
Gaussian elimination
elimination
x1 22 33 1
0 22 53 2 3
0 32 53 2 Forward
2
elimination
x 2x 3x 1 12315
022 53 2 x20
Backward substitution
005×1 x2 23 35
x 1

Gaussian elimination
Forward elimination (reduce the number of unknowns)
AxAxAxb A
111 122 133 1 21 A
A x A x A x b A 31
211 222 233 2 11 AxAxAxb
A
31 1 32 2 33 3 3
AxAxAxb 11 1 12 2 13 3 1
21 12 21 13 21 1 0 AAAx AAAx bAb
22 A2 23 A3 2 A 1111 11
31 12 31 13 31 1 0 AAAx AAAx bAb
32 A2 33 A3 3 A 1111 11
11

Gaussian elimination
Forward elimination (reduce the number of unknowns)
AxAxAxb A
111 122 133 1 21 A
A x A x A x b A 31
211 222 233 2 11 AxAxAxb
A
31 1 32 2 33 3 3
AxAxAxb
111 122A133A1b
AA22 AA23 Ab2 0 A 21 12 x A 21 13 x b 21 1
22 A2 23 A3 2 A 11A11A 11b
AA32 AA33 Ab3 0 A 31 12 x A 31 13 x b 31 1
32 A2 33 A3 3 A 1111 11
AxAxAxb 11 1 12 2 13 3 1
AxAxb 22 2 23 3 2
AxAxb 32 2 33 3 3
11

Forward elimination (reduce the number of unknowns)
AxAxAxb
11 1
12 2 13 3 1
AxAxb
222 233 2 A
AxAxb 32 2 33 3 3
32
AxAxAxb 111 122 133 1
AxAxb 22 2 23 3 2
32 23 32 2 0 AAAx bAb
33A3 3A 22 22
A
22

Forward elimination (reduce the number of unknowns)
AxAxAxb
11 1
12 2 13 3 1
AxAxb
222 233 2 A
AxAxb 32 2 33 3 3
32
AxAxAxb 111 122 133 1
AxAxb
222 233A2b
AA33 Ab3 0 A 32 23 x b 32 2
33A3 3A 22 22
AxAxAxb 11 1 12 2 13 3 1
AxAxb 22 2 23 3 2
Axb 333 3
A
22

11 1
12 2 13 3 1
Backward substitution
AxAxAxb
AxAxb
22 2
23 3 2
b
AxAxAxb
11 1
12 2 13 3 1 AxAxb x2 233
Axbx 333 3 3
3 33
222 233 2 2 Axb
bAx A
111 122 133 1 1 AxAxb
A
333 3
AxAxAxb xbAxAx
22 2 23 3 2
Axb 333 3
A
22
1 12 2 13 3
11

Forward elimination: Step (1) AxAxAxb
11 1 12 2 13 3 1
AxAxAxb 21 1 22 2 23 3 2
AxAxAxb 31 1 32 2 33 3 3
AxAxAxb
111 122A133A1b
AA22 AA23 Ab2 0 A 21 12 x A 21 13 x b 21 1
22 A2 23 A3 2 A 11A11A 11b
AA32 AA33 Ab3 0 A 31 12 x A 31 13 x b 31 1
32 A2 33 A3 3 A 1111 11
AxAxAxb 11 1 12 2 13 3 1
AxAxb 22 2 23 3 2
AxAxb 32 2 33 3 3

Forward elimination: Step (1) jm mk
A A A A (m1;j,k2,3) jk jk
AxAxAxb 11 1 A12 2 13 3 1
mm
AxAxAxb 21 1 A b22 2 23 3 2
bb jm m (m1;j2,3)
Aj x j AA x A x b 31 1 32 2 33 3 3
mm
AxAxAxb
111 122A133A1b
AA22 AA23 Ab2 0 A 21 12 x A 21 13 x b 21 1
22 A2 23 A3 2 A 11A11A 11b
AA32 AA33 Ab3 0 A 31 12 x A 31 13 x b 31 1
32 A2 33 A3 3 A 1111 11
AxAxAxb 11 1 12 2 13 3 1
AxAxb 22 2 23 3 2
AxAxb 32 2 33 3 3

Forward elimination: Step (2) AxAxAxb
11 1 12 2 13 3 1
AxAxb 22 2 23 3 2
AxAxb 32 2 33 3 3
AxAxAxb 111 122 133 1
AxAxb
222 233A2b
AA33 Ab3 0 A 32 23 x b 32 2
33A3 3A 22 22
AxAxAxb 11 1 12 2 13 3 1
AxAxb 22 2 23 3 2
Axb 333 3

Forward elimination: Step (2)
AA
AxAxAxb
11jk1
A A jm mk (m2;j,k3)
jk 12 2 13 3 1
A
A xmm A x b
22 2 23 3 2
Ab
Ajm mx A x b
(m2;j3) AxAxAxb
b b jj
32 2 33 3 3
A mm
111 122 133 1
AxAxb
222 233A2b
AA33 Ab3 0 A 32 23 x b 32 2
33A3 3A 22 22
AxAxAxb 11 1 12 2 13 3 1
AxAxb 22 2 23 3 2
Axb 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,km1,,N)
jk jk
A mm
jm m bbA b
jj
mm
(jm1,,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,km1,,N)
jk jk
A mm
jm m bbA 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)
(jm1,,N) Python code

Backward substitution:
Step (3)
x
A
3 3
Step (2)
xbAx 2 233
2
22
Step (1)
xbAxAx 1 12 2 13 3
1
A
Step (m) (m = 2, 1)
x m
m
km1
mk k
b
A
3 bAx
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
km1
mk k
x N
b N
Step (m) (m = N 1, , 1)
Algorithm
A NN
N bAx
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
km1
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 bAx
A mm

Number of operations for input of size N = ?
A A 21 22
A aAaAaA a N 2N 2 211 222 2NN
Computational complexity
Example: multiplication of a matrix with a vector
Assume that multiplication is much more timeconsuming than addition
A A 11 12
A a AaAaA 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 106 sec per multiplication, we will get result after 106N2 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)
AAAAxb 21 22 23 2N22
11
A
A N1
A A A Axb 31 32 33 3N 3 3
11
11

A
A A A Axb N1 N2 N3 NN N N
A A A Axb 11 12 13 1n11
Step(2):(N1)2 1
0 A A Axb 22 23 2n 2 2
Step (N 1): 12 1 N2
0 A A Axb 32 33 3n 3 3

(k1) k1
Computational complexity
Forward elimination
(N 1) Step(1):(N+1)(N1)=N2 1
0 A A Axb
N2 N3 NNN N(N1)(N1)

N2,N2
N2,N1 N2,N N2 N2 N1,N1 N1,N N1 N1

A Axb Axb

NN N N 1
b Step(N): xN
Computational complexity
Backward substitution
A A A A xb
11
1,N2 1,N1 1N 1 1

A A Axb
N
x b N1,N N
Step (N 1): N1
N1 A A
Step(1):
1122 1NN x b A x A x
N
N k1
k
1
AAA
A
NN Ax
2 Step(N2): x N2 N2,N1 N1 N2,N N 3
N2
AAA
N1,N1 N1,N1 bAxAx
N2,N2
N2,N2
N2,N2

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
k1 k1
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 1010 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)
AxAxAxb A
111 122 133 1 21 A
A x A x A x b A 31 211 222 233 2 11 A
AxAxAxb 11 31 1 32 2 33 3 3
1 00AAA 11 12 13
LAAA 10A A A 1 21 11 21 22 23
A A 0 1 A A A 3111 313233
AAA 1112 13
0AAAAAAAA 22 2112 11 23 2113 11
0AAAAAAAA
32 3112 11 33 3113 11
AAA 11 12 13
0AA 22 23
0AA
32 33

AxAxAxb
11 1
12 2 13 3 1
1
0A A A 11 12 13
Step (2)
AxAxb
222 233 2 A
AxAxb
32
32 2
33 3 3
A
0
LLA0 1 00AA
21 22 23 0 A A 1 0 A A
32 22 32 33 AAA
1112 13 0AA
22 23 00AAAA
33321322 AAA
11 12 13 0AA
22 23
00A
33
22

1 0 01 0 0A 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 0a 1 0A
21 22 23 22 23

U U 0c1b01AAA 00U
A A A1 0 0U U U 11 12 13 11 12 13
A A Aa 1 00 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 1A A A A
N111 N1N2N3 AAAA
NN
11 12 13 1N 0AAA
2223 2N 0AAA
3233 3N
0AAA N2N3 NN

L2L1A
0 1
0 A A 1
0AAA 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
0A A 010 A A A N2 22 N2 N3 NN
AAAA 11 12 13 1N
0AAA 2223 2N
00AA 33 3N

00AA 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 00AA
33 3N
0 0 A A 1 0 0 A A N333 N3 NN
AAAA 11 12 13 1N
0AAA 2223 2N
00AA 33 3N

000A NN

LN1L3L2L1A 1
AAA A 11 12 1,N1 1N

0 1
0 A A A 22 2,N1 2N
00
1
00 A A
If A is a N N matrix Gaussian elimination: Step (N 1)

N,N1 N1,N1 N,N1 NN

N1,N1 N1,N 00AA 100A A
AAAA 11 12 1,N1 1N
0AA A 22 2,N1 2N
00U AA
N1,N1 N1,N
000A NN

1
UUUU 11 12 13 1N
LU decomposition LN1L3L2L1A U
A (LN1L3L2L1)1ULU
L1UUU 21 2223 2N
LLL1UUU 3132 333N

LLL1 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
LL L 1 U 31 32
U U 33 34
LLL1 414243
U 44

U14
L41U11

U A (j1,,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
LUALAU (i2,,N) i111 i1 i1 i1 11
L41U11

U A (j1,,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 (i2,,N) i1 i1 11
i1 12 i2 22
i2 i2 i2 i1 12
22
LU decomposition: Doolittle algorithm
LU U A U A LU (j2,,N) (2)21 1j 2j 2j 2j 2j 21 1j
LU LU A L (A LU )U (i3,,N)
L41U11

U A (j1,,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 (i2,,N) i1 i1 11
LU decomposition: Doolittle algorithm
U A L U (j2,,N) (2)2j 2j 211j
L (A LU )U (i3,,N) i2 i2 i112 22
LU LU U A U A LU LU (j3,,N) (3)31 1j 32 2j 3j 3j 3j 3j 31 1j 32 2j
LU LU LU A L (A LU LU )U (i4,,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
m1

im im
ikkm mm

k1
(im1,,N)
LA LUU
Step (1)
U A (j1,,N) 1j 1j
L A U (i2,,N) i1 i1 11
Step (m)
(m = 2, , N 1)
m1
U A L U (jm,,N)
Step (N) N1
UALU
NN NN
Nk kN
k1

Solving system of linear equations by LU decomposition
Axb LUxb

L
L
U
x=z
A x=b
U
x=b
z=b
Lzb Ux z

Solving system of linear equations by LU decomposition
Forward substitution
Axb LUxb

1 z b
zb 11
1 1 21 22

L 1 zb
z b L z 2 2 211
LL1zb
N1 zbLz
N1 N2 Backward substitution

NN UUUxz
NN
Nk k k1
1kk 11 11 12 1N11

2N 2 2
k3

2kk 22

UUxz N
22
x2z2 UxU
NNNN xzU
N N NN
Uxz
Lzb Ux z
N x1z1 UxU
k2

Solving system of linear equations by LU decomposition
Forward substitution
Axb LUxb

1 z b
zb 11
1 1 21 22

L 1 zb
z b L z 2 2 211
LL1zb N1 N2 NN
N1 zbLz
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]
Lzb Ux z
NN
Nk k k1

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
Axb LUxb

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
Lzb Ux z

LU decomposition: Doolittle algorithm
mj mj mkkj k 1
m1 LA LUU

im im
ikkm mm

k1
(im1,,N)
Step (1)
U A (j1,,N) 1j 1j
L A U (i2,,N) i1 i1 11
It is not guaranteed thatallUmm =0! Pivoting
Step (m)
(m = 2, , N 1)
m1
U A L U (jm,,N)
Large Umm reduces roundoff error
P i v o t i n g
Step (N)
N 1 UALU
NN NN
k1
Nk kN

Pivoting
AxAxAxb 11 1 12 2 13 3 1
Exchange 2 equations
AxAxAxb 21 1 22 2 23 3 2
AxAxAxb 31 1 32 2 33 3 3
AxAxAxb 21 1 22 2 23 3 2
AxAxAxb 11 1 12 2 13 3 1
AxAxAxb 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 22 33 1
2x 5x x 0 1223
41 52 73 2
Idea of pivoting: an example

x1 22 33 1 2x 5x x 0
1223
41 52 73 2
Idea of pivoting: an example

x1 22 33 1 2x 5x x 0
1223
41 52 73 2
41 52 73 2 2x 5x x 0
1223
x1 22 33 1
Idea of pivoting: an example
Pivoting

x1 22 33 1 2x 5x x 0
1223
41 52 73 2
Idea of pivoting: an example
Pivoting
41 52 73 2 2x5xx0 4 1
2 12234
x1 22 33 1

x1 22 33 1 2x 5x x 0
1223
41 52 73 2
Idea of pivoting: an example
Pivoting
2 x 2x 3x 1 0 3x 5x 1
41 52 73 2 2x 5x x 0
41 52
1223 23
1 2 3 Forward 42 43 2
elimination
73
0 0 5x 1

41 52 73 2 2x 5x x 0
41 52 73 0 0 5x
2 1
Idea of pivoting: an example
x1 22 33 1 41 52 73 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

33 1 2x 5x x 0
41 52 73 2 0 3x 5x 1
Idea of pivoting: an example
x1 22
1223 42432
Solution: x2, x0, x1 35215
4x 5x 7x 2 0 0 5x 1 12323
41 52 73 2 2x 5x x 0
41 52 73 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 1a a 11 12
a a a 13 31 32
a 33
1 0 0a a
a a a 23 11 12
a 13

0 1
21 22 0a31 a32
a23
P
a33 a21 a22 0 0 1b b
1 3 1 0 0b 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
AMN QMMRMN Q : unitary matrix: QQ = 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 bRx 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 AHA 1 0
0

0 0 RHHHA 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)1H1H1H1
N1 21 12 N1
But we still dont know (1) how to construct these {Hk} and (2) why Q is
an orthogonal matrix

Householder reflector
Hzz2e1, Hz20 z0
3 Householder reflector H
transform a vector z to z2e1, a vector of same length
e2
z zTz 2
z
z1 z z zz
z1 z 2
Hz
e1

3 z2 z2 z2 (length)
1232 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
22 Projztovz cosvz zvvz
e2

zv v vTz
v vv vTv
z22
Hz
e1
vv 22

Householder reflector
e2
21 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
Hzz2 T vz2 T zI2 T zHI2 T vvvvvvvv
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 = z2e1
e1
Householder reflector Alternatively, Hz z 2 e1
In this case, need to calculate z1 + z2; since z1 > 0, this is addition. In previous case, need to calculate z1 z2; since z1 > 0, this is subtraction. Avoid calculating difference of 2 floatingpoint numbers which are close (why?)
e2
Hzsgn(z1) z 2 e1 v z sgn(z1 ) z 2 e1
z
HI2vvT vTv

AAA 11 12 13
AAA 11 12 13
AA A A 21 22 23
HA0AA 1 22 23
AAA 11 12 13
vT v
A11 A12 A13
QR decomposition
AAA 0AA
31 32 33 a1
32 33 a2
vasgn(A)a e 111121
vasgn(A)a e 2 22221
H I2vvT 1 vT v
1 0
HHA0AAR HA0AA 21 2223
12223
00A
0AA 33 32 33

H2 0 I2vvT
Q (H H )1 H1H1 2112

QR decomposition
Q (H H )1 H1H1 2112
Check if Q is orthogonal. Property of Householder reflector I0I0
H vvT HT vvT H 0 I2 0 I2
vT v vT v
I 0I 0I0 HHT vvT vvT I
0 I2 0 I2 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
HHH T

k kk1
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
T11 1 RHN1H2H1A, Q H1 H2 HN1 HN1H2H1I
R A QT I
0 ,0 (k1,,N1)
R HR QT HQT
k kk1 Ik1
0 T
vasgn(A)a e,H
k kkk21 k0I2vv
R:R ,QT:QT N1 N1
Nk1
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 115 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 ,vasgn(A)a e 2
1 1 1
1 1 3 1 0 1
a1
HI2vvT I213 1 1 11 5 1 1
3 2222
1 4 vTv 4 121 1 1 5 1 2666
11115 2666
242 1 1 1 1 2222
33TT2666 RHR0 4 8,QHQ1 5 1 1
1 1002161 101 1 5 1 332666
04 20 1 1 1 5 332666
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 , va sgn(A )a e 2202
4 1 10 10216 2 2222133
33 44 0 4 20 3 03
33 a2
10 212 vvT3 3 3
H T,I2 1 14 2 20I2vv3 T315 15
3 vv2211 vTv 3 15 15
242 1111 2222
TT2222 R HR 0 2 8, Q HQ 1 1 1 1
2 2100162 2171 7 1 5 10 10 10 10
0 0 12 1 7 1 7 5 10101010

0 2 8 20016 3
16 1 36 3332112012
5 0 0 12
55
5 a3
QR decomposition
Step (3)
R ,vasgn(A)a e5 4 5
2 4 2
I2 0
H ,I2I5361255
vvT
5 36 4 3
3 4 555
3 0I2vvT2 vTv2 144125 5 2 vTv
242 1 1 1 1 2222
TT2222 R HR 0 2 8, Q HQ 1 1 1 1
3 320043 321 111 2222
000 1111 2222

A QR
11 22
1 10 2 8 22
Q(1)
R(1) Q(2)
Q(1)R(1)
11 2 2
1 1242 2 2
1 1 2 2
1 10 0 4 2 2
0
QR decomposition
Result
242 1111 2222
TT2222 RR 0 2 8, Q Q 1 1 1 1
3004 31 111 2222
000 1111 2222
1111000 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()
11111 1 1242 2 2 2 2
A1 3 3 1 1 1 10 2 8QR 2222
1 1 5 1 1 1 10 0 4 2 2 2 2
1371111000
2 2 2 2

Solve system of linear equations by QR decomposition Ax QRx bRx 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 MN 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 MN 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 MN 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 AQR T T
x(A A) A b(QR) (QR)1(QR) b RT QT QR1 RT QT b RT R1 RT QT b
R1RTRTQTbR1QTb
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 (j1,,N) 0 (1) 1j 1j
N (k1)(Nk1)
(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 (i2,,N) 0 i1 i1 11
k1 N
U A L U (j2,,N) (N1)
(k1)(Nk) k1
(2)2j 2j 211j
L (A L U ) U
(i3,,N) (N2)
i2 i2 i1 12 22 U A L U L U
N3N2N
(j3,,N) 2(N2)
L (A L U L U ) U (i4,,N) 2(N3)
3
2
6
L41U11
L41U13 + L42U23 + L43U33
L41U14 + L42U24 + L43U34 + U44

Cholesky decomposition Ax b
If A is symmetric and positivedefinite, what is LU for A? 1 0 0U U U
11 12 13 AL 1 00 U U
21 22 23 LL100U
31 32 33
1 00U 0 01UU UU
11 12 11 13 11 L 100U 00 1 UU
21 22 2322 LL100U001
313233
AAT
1 00U001LL
L D LT
T 11 21 31
AUU 1 00U 001L 12 11 22 32
UUUU100U001 13 11 23 22 33

11 AL 1 00 U 0L 1 0
1 L
0 0 U11 0 1 00 U
0 00 U
0 T1 0L
0 0T 1 0
22 L L 1
22
21

21

L L 1

Cholesky decomposition
A LDL
1 0 0U 0 01 0 0T
21 22 21 LL100ULL1
31 32 AispositivedefiniteUii 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
RR 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 RARR ,RARR ,,RARR
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
RARR ,RARR ,,RARR 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
RARRRR ,RARRRR ,,RARRRR 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 fork1,,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 k1
21 31 N1

R A R2
j1
11 11 11
32 32R21 31 42 42R21 41 N2 N2R21 N1 RARR ,RARR ,,RARR
R A R2
k1 22 22 22
ik ik
kj ji kk RARRRR ,RARRRR ,,RARRRR
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

j1
RA RRT R (ik1,,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
(N2) RARRRR ,RARRRR ,,RARRRR 2(N3)
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()
Axb RRTxb

# 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
Rzb 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 (j1,,N) (1) 1j 1j
LUALAU (i2,,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 (j1,,N) (1) 1j 1j
LUALAU (i2,,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 (j2,,N)
(2)21 1j 2j 2j 2j 2j 21 1j
LU LU A L (A LU )U (i3,,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 (j2,,N)
(2)21 1j 2j 2j 2j 2j 21 1j
LU LU A L (A LU )U (i3,,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 (j3,,N)
(3)31 1j 32 2j 3j 3j 3j 3j 31 1j 32 2j
LU LU LU A L (A LU LU )U (i4,,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 (j3,,N)
(3)31 1j 32 2j 3j 3j 3j 3j 31 1j 32 2j
LU LU LU A L (A LU LU )U (i4,,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