Nash - Compact Numerical Methods for Computers (523163), страница 20
Текст из файла (страница 20)
Since (7.5) and (7.7)define the first column of L, we have a stepwise procedure for computing itsremaining columns, and furthermore this can be arranged so that L overwrites A84The Choleski decomposition85within the computer. It remains to be shown that the procedure is stable and thatfor i = m the right-hand side of (7.8) is positive, so no square roots of negativenumbers are required.Firstly, A is positive definite ifx T A x >0for all x 0.(7.9)An equivalent statement is that all the eigenvalues of A are positive. From (7.9) itfollows by setting x to any column of the unit matrix that no diagonal element ofA can be non-positive. Likewise, by taking only xi and xi non-zero(7.10)which requires that the quadratic equationz2Aii + 2zAij + Ajj = 0(7.11)has only complex roots.
This occurs if(7.12)Consider now the ith step of the Choleski decomposition. For the momentsuppose that only rows 1, 2, . . . , (i - 1) of L have been computed, giving asubmatrix Li-1 which is a decomposition of the submatrix Ai-1 of A; hence(7.13)Following Kowalik and Osborne (1968), we haveLi - l c = a(7.14)or(7.15)where Li-1 is assumed non-singular.
In fact, it is positive definite providing thepositive square root is chosen in the computation of each of its diagonal elementsvia (7.8). Consider now the choice in (7.9) of an x such that the first ( i – 1)elements of x are given byxi = -1, and xj = 0 for j > i. This choice, using(7.13) gives(7.16)which reduces toAii - cTc > 0.(7.17)But a comparison of this with (7.8) shows that it implies the square of eachdiagonal element of L is positive, so that all the elements of L are real providing Ais positive definite.
Furthermore, an analysis similar to that used in (7.10), (7.11)and (7.12) demands that(7.18)86Compact numerical methods for computers(Again, the diagonal elements must be chosen to be positive in the decomposition.) Equations (7.17) and (7.18) give bounds to the size of the subdiagonalelements of L, which suggests the algorithm is stable. A much more completeanalysis which confirms this conjecture is given by Wilkinson (1961) who showsthe matrix LLT as computed is always close in some norm to A.Once the Choleski decomposition has been performed, linear equationsAx=LLT x = b(7.19)can be solved by a combination of a forward- and a back-substitution, that isLv = b(7.20)followed byRx = LT x = v(7.21)Twhere we have used R to emphasise the fact that L is upper-triangular.
In acomputer program, b, v and x can all occupy the same storage vector, that is, voverwrites b, and x overwrites v. The solution of (7.20) is termed forwardsubstitution because the triangular structure defines the elements v j in the order1, 2, . . . , n, that is(7.22)v l = b l/ L 1 1andfor j = 2, 3, . . . , n.(7.23)Likewise, the solution elements xj of (7.21) are obtained in the backward order n,(n – 1), . . . , 1 from(7.24)xn= vn / Ln n(7.25)(7.26)7.2. EXTENSION OF THE CHOLESKI DECOMPOSITION TONON-NEGATIVE DEFINITE MATRICESWhen A is non-negative definite, the inequality (7.9) becomesx TAx > 0for all x0(7.27)and inequalities (7.10), (7.12), (7.17) and (7.18) must be amended similarly.There is no difficulty in performing the decomposition unless a zero diagonalelement appears in L before the decomposition is complete. For L mm = 0, equations (7.3) and (7.8) are satisfied by any values of Lim for i > m.
However, if wedesire L to be non-negative definite and wish to satisfy the amended form of(7.18), that is(7.28)The Choleski decomposition87we should set Lim = 0 for i > m. This is a relatively trivial modification to thedecomposition algorithm. Lmm is, of course, found as the square root of a quantitywhich arises by subtraction, and in practice this quantity may be slightly negativedue to rounding error. One policy, adopted here, is to proceed with the decomposition, setting all Lim = 0 for i > m even if this quantity is negative, thusassuming the matrix A is non-negative definite.
Since there may be very goodreasons for presupposing A to be non-negative definite, for instance if it has beenformed as a sum-of-squares and cross-products matrix in a least-squares regression calculation, this is not as dangerous as it appears. Furthermore the decisionto continue the decomposition in §7.1 when the computed square of Lmm ispositive, rather than greater than some tolerance for zero, has been made afterthe following considerations.(i) The decomposition is valid to the precision available in the arithmetic beingused. When zero diagonal elements arise in L they reflect linear dependencies inthe set of equations for which a solution is to be found, and any of the infinity ofsolutions which exist is taken to be acceptable.
However, in recognition of thepossibility that there may only be a near-linear dependence, it does not seem wiseto presume a small number is zero, since the Choleski decomposition, unlike thesingular-value decomposition (algorithm 1), does not allow the user to decide atthe time solutions are computed which small numbers are to be assumed zero.(ii) The size of the computed square of Lmm is dependent on the scale of thematrix A. Unless the tolerance for zero is proportional to a norm of A, itsapplication has an effect which is not consistent from one problem to another.If the computed square of Lmm is non-positive, the mth column of L istherefore set to zeroLim = 0for i = m, (m + 1), .
. . , n.(7.29)The forward- and back-substitutions to solve the linear equationsAx = b(2.2)are unfortunately not now possible since division by zero occurs. If, however, theequations are consistent, that is, if b belongs to the column space of A, at least onesolution x exists (see, for example, Finkbeiner 1966, p 98).Consider the forward-substitution to solveLv = b(7.30)for v. If v k is set to zero whenever Lkk = 0, then solutionsLT x = v(7.31)are solutions to (2.2) for arbitrary values of those xk for which Lkk = 0. This is, ofcourse, only possible if(7.32)which is another way of stating the requirement of consistency for the equations.88Compact numerical methods for computersFor a specific example, consider that Lmm = 0 as above. Thus, in the forwardsubstitution v m is always multiplied by zero and could have arbitrary value, exceptthat in the back-substitution the mth row of LT is null.
Denoting this by the vector(7.33)it is easily seen that(7.34)so that vm must be zero or the equation (7.34) is not satisfied. From (7.30) and(7.31) one hasLv = LLT x= b =Ax.(7.35)Since xm is arbitrary by virtue of (7.34), the value chosen will influence thevalues of all xi, i < m, so some standard choice here is useful when, for instance,an implementation of the algorithm is being tested. I have always chosen to setxm = 0 (as below in step 14 of algorithm 8).The importance of the above ideas is that the solution of linear least-squaresproblems by the normal equationsBT Bx = BT y(7.36)provides a set of consistent linear equations with a symmetric non-negativedefinite matrix A = BTB, that is(7.37)(B is presumed to be m by n).
The equations (7.36) are always consistent sincethe vector BT y belongs to the row space of B, which is identical to the columnspace of BT B.There remain two organisational details: (a) in any program designed to savestorage it is useful to keep only one triangle of a symmetric matrix and further tostore only the non-zero parts of triangular matrices, and (b) the forward- andback-substitutions can be organised to overwrite υ on b and then x on v.These ideas are incorporated into the following algorithms for the Choleskidecomposition of a non-negative definite symmetric matrix (see Healy (1968) for aFORTRAN implementation) and solution of consistent linear equations by thesubstitutions described above.Algorithm 7.
Choleski decomposition in compact storageprocedure 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 incompact row-order form.
a[i*(i-1)/2+j] = A[i,j]Copyright 1988 J.C.Nash!The Choleski decomposition89Algorithm 7. Choleski decomposition in compact storage (cont.)}vari,j,k,m,q: integer;s : real; {accumulator}beginsingmat := false; {singmat will be set true if matrix foundcomputationally singular}for j := 1 to n do {STEP 1}begin {STEP 2}q := j*(i+1) div 2; {index of the diagonal element of row j}if j>1 then {STEP 3}begin {prepare for the subtraction in Eqn. (7.8). This is not neededfor the first column of the matrix.}for i := j to n do {STEP 4}beginm := (i*(i-1) div 2)+j; s := a[m];for k := 1 to (j-1) do s := s-a[m-k]*a[q-k];a[m] := s;end; {loop on i}end; {of STEP 4}if a[q]<=0.0 then {STEP 5}begin {matrix singular}singmat := true;a[q] := 0.0; {since we shall assume matrix is non-negative definite)}end;s := sqrt(a[q]); {STEP 7}for i := j to n do {STEP 8}beginm := (i*(i-1) div 2)+j;if s=0.0 then a[m] := 0 {to zero column elements in singular case}else a[m] := a[m]/s; {to perform the scaling}end; (loop on i)end; {loop on j -- end-loop is STEP 9}end; {alg07.pass == Choleski decomposition choldcmp}This completes the decomposition.
The lower-triangular factor L is left in the vector a inrow-wise storage mode.Algorithm 8. Choleski back-substitutionprocedure cholback(n: integer; {order of problem}a: smatvec; {the decomposed matrix}var x: rvector); {the right hand side}{alg08.pass ==Choleski back substitution for the solution of consistent sets oflinear equations with symmetric coefficient matrices.Copyright 1988 J.C.Nash90Compact numerical methods for computersAlgorithm 8. Choleski back-substitution (cont.)}vari,j,q : integer;begin {Forward substitution phase -- STEP 1}if a[1]=0.0 then x[1]:=0.0 {to take care of singular case}else x[1]:=x[1]/a[1];{STEP 2}if n>1 then {do Steps 3 to 8; otherwise problem is trivial}begin{STEP 3}q:=1; {to initialize the index of matrix elements}for i:=2 to n do {STEP 4}begin {STEP 5}for j:=1 to (i-1) dobeginq:=q+1;x[i]:=x[i]-a[q]*x[j];end; {loop on j}q:=q+1; {STEP 6 -- to give index of diagonal element of row i}if a[q]=0.0 then x[i]:=0.0 {to handle singular case}else x[i]:=x[i]/a[q]; {STEP 7}end; {loop on i -- STEP 8}end; {non-trivial case.
This completes the forward substitution}{STEP 9 -- Back substitution phase}if a[n*(n+1) div 2]=0.0 then x[n]:=0.0 {for singular case}else x[n]:=x[n]/a[n*(n+1) div 2];if n>1 then {STEP 10} {test for trivial case; otherwise do steps 11 to 15}begin {STEP 11}for i:=n down to 2 dobegin {STEP 12}q:=i*(i-1) div 2; {to give base index for row i}for j:=1 to (i-1) do x[j]:=x[j]-x[i]*a[q+j]; {STEP 13}if a[q]=0.0 then x[i-1]:=0.0 {for singular case}else x[i-1]:=x[i-1]/a[q]; {STEP 14}end; {loop on i -- STEP 15}end; {non-trivial case -- STEP 16}end; {alg08.pas == Choleski back-substitution cholback}Note that this algorithm will solve consistent sets of equations whose coefficient matrices aresymmetric and non-negative definite. It will not detect the cases where the equations are notconsistent or the matrix of coefficients is indefinite.7.3.
SOME ORGANISATIONAL DETAILSThe algorithm given for the Choleski decomposition uses the most direct form forcomputing the array indices. In practice, however, this may not be the mostefficient way to program the algorithm. Indeed, by running a FORTRAN version ofalgorithm 7 against the subroutine of Healy (1968) and that of Kevin Price(Nash 1984a, pp 97-101) on an IBM 370/168 it became quite clear that the latter twocodes were markedly faster in execution (using the FORTRAN G compiler) than myown. On the Data General NOVA and Hewlett-Packard 9830, however, this91The Choleski decompositionfinding seemed to be reversed, probably because these latter machines are run asinterpreters, that is, they evaluate the code as they encounter it in a form ofsimultaneous translation into machine instructions whereas a compiler makes thetranslation first then discards the source code.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
