Page 54 - Numerical Methods for Chemical Engineering

Matrix factorization 43

From the (1, 2) element,

a 12 = L 11 × L 21 L 21 = a 12 /L 11 (1.214)
Similarly, we compute the first column of L for j = 3, 4,..., N,
(1.215)
L j1 = a 1 j /L 11
Next, we move to the second column, and from the (2, 2) element obtain
2 1/2
L 22 = a 22 − L (1.216)
a 22 = L 21 L 21 + L 22 L 22
21
and for j = 3, 4,..., N,
(1.217)
a 2 j = L 21 L j1 + L 22 L j2 L j2 = (a 2 j − L 21 L j1 )/L 22
Continuing this process for columns 3, 4, etc. yields the algorithm:
for i = 1, 2,..., N; iterate over each column of L
1/2

i−1
2
L ii ← a ii − L ik
k=1
for j = i + 1, i + 2,..., N; iterate over elements below the diagonal

i−1
1
L ji ← a ij − L ik L jk
L ii
k=1
end j = i + 1, i + 2,..., N
end i = 1, 2,..., N
2
As there are only two nested loops, the FLOPs required scale only as N .
In MATLAB, Cholesky decomposition is invoked by chol; however, this function returns
T
T
not a lower-triangular matrix L with A = LL but an upper-triangular matrix R = L with
T
A = R R. For the positive-definite matrix
 
211
A =  142  (1.218)
126
T
the Cholesky factorization A = R R is performed by the code,
A=[211;142;126];
R = chol(A),
R=
1.4142 0.7071 0.7071
0 1.8708 0.8018
0 0 2.2039
T
We see that A = R R by
∗
R R,
ans =
2.0000 1.0000 1.0000
1.0000 4.0000 2.0000
1.0000 2.0000 6.0000