Page 99 - Compact Numerical Methods For Computers

88 Compact numerical methods for computers
For a specific example, consider that L = 0 as above. Thus, in the forward-
mm
substitution v is always multiplied by zero and could have arbitrary value, except
m
T
that in the back-substitution the mth row of L is null. Denoting this by the vector
(7.33)
it is easily seen that
(7.34)
so that v m must be zero or the equation (7.34) is not satisfied. From (7.30) and
(7.31) one has
T
Lv =LL x= b=Ax. (7.35)
Since x m is arbitrary by virtue of (7.34), the value chosen will influence the
values of all x , i < m, so some standard choice here is useful when, for instance,
i
an implementation of the algorithm is being tested. I have always chosen to set
x = 0 (as below in step 14 of algorithm 8).
m
The importance of the above ideas is that the solution of linear least-squares
problems by the normal equations
T T
B Bx = B y (7.36)
provides a set of consistent linear equations with a symmetric non-negative
T
definite matrix A = B B, that is
(7.37)

(B is presumed to be m by n). The equations (7.36) are always consistent since
T
the vector B y belongs to the row space of B, which is identical to the column
T
space of B B.
There remain two organisational details: (a) in any program designed to save
storage it is useful to keep only one triangle of a symmetric matrix and further to
store only the non-zero parts of triangular matrices, and (b) the forward- and
back-substitutions can be organised to overwrite u u on b and then x on v.
These ideas are incorporated into the following algorithms for the Choleski
decomposition of a non-negative definite symmetric matrix (see Healy (1968) for a
FORTRAN implementation) and solution of consistent linear equations by the
substitutions described above.

Algorithm 7. Choleski decomposition in compact storage
procedure choldcmp(n: integer; {order of problem}
var a: smatvec; {matrix to decompose}
var singmat: boolean); {singularity flag}

(alg07.pas ==
Choleski decomposition of symmetric positive definite matrix stored in
compact row-order form. a[i*(i-1)/2+j] = A[i,j]
Copyright 1988 J.C.Nash