Page 108 - Compact Numerical Methods For Computers

The symmetric positive definite matrix again 97
8.2. THE GAUSS-JORDAN ALGORITHM FOR THE INVERSE
OF A SYMMETRIC POSITIVE DEFINITE MATRIX
Bauer and Reinsch (in Wilkinson and Reinsch 1971, p 45) present a very compact
algorithm for inverting a positive definite symmetric matrix in situ, that is,
overwriting itself. The principal advantages of this algorithm are as follows.
(i) No pivoting is required. This is a consequence of positive definiteness and
symmetry. Peters and Wilkinson (1975) state that this is ‘well known’, but I
believe the full analysis is as yet unpublished.
(ii) Only a triangular portion of the matrix need be stored due to symmetry,
though a working vector of length n, where n is the order of the matrix, is needed.

The algorithm is simply the substitution procedure outlined above. The
modifications which are possible due to symmetry and positive definiteness,
however, cause the computational steps to look completely different.
Consider an intermediate situation in which the first k of the elements x and b
have been exchanged in solving

Ax = b (8.6)

by the Gauss-Jordan algorithm. At this stage the matrix of coefficients will have
the form
W X
Y Z (8.7)

with W, k by k; X, k by (n – k); Y, (n – k) by k; and Z, (n – k) by (n – k). Then W
is the inverse of the corresponding block of A since the equations (8.6) are now
given by their equivalents

(8.8)

Thus by setting x = 0, for j = (k + 1), (k + 2), . . . , n, in both (8.6) and (8.8) the
j
required association of W and the leading k by k block of A -1 is established.
Likewise, by setting b = 0, for j = 1, . . . , k, in (8.6) and (8.8), Z is the inverse of
j -l
the corresponding block of A . (Note that W and Z are symmetric because A and
A -1 are symmetric.)
From these results, W and Z are both positive definite for all k since A is
positive definite. This means that the diagonal elements needed as pivots in the
Gauss-Jordan steps are always positive. (This follows directly from the definition
T
of positive definiteness for A, that is, that x Ax > 0 for all x 0. )
In order to avoid retaining elements above or below the diagonal in a program,
we need one more result. This is that

T
Y = –X (8.9)
in matrix (8.7). This can be proved by induction using the Gauss-Jordan