MCHOLDC Name
MCHOLDC
Author
Craig B. Markwardt, NASA/GSFC Code 662, Greenbelt, MD 20770
craigm@lheamail.gsfc.nasa.gov Purpose
Modified Cholesky Factorization of a Symmetric Matrix
Major Topics
Linear Systems
Calling Sequence
MCHOLDC, A, D, E [, /OUTFULL, /SPARSE, /PIVOT, TAU=TAU, $
PERMUTE=PERMUTE, INVPERMUTE=INVPERMUTE ]
Description
Given a symmetric matrix A, the MCHOLDC procedure computes the
factorization:
A + E = TRANSPOSE(U) ## D ## U
where A is the original matrix (optionally permuted if the PIVOT
keyword is set), U is an upper triangular matrix, D is a diagonal
matrix, and E is a diagonal error matrix.
The standard Cholesky factorization is only defined for a positive
definite symmetric matrix. If the input matrix is positive
definite then the error term E will be zero upon output. The user
may in fact test the positive-definiteness of their matrix by
factoring it and testing that all terms in E are zero.
If A is *not* positive definite, then the standard Cholesky
factorization is undefined. In that case we adopt the "modified"
factorization strategy of Gill, Murray and Wright (p. 108), which
involves adding a diagonal error term to A in order to enforce
positive-definiteness. The approach is optimal in the sense that
it attempts to minimize E, and thus disturbing A as little as
possible. For optimization problems, this approximate
factorization can be used to find a direction of descent even when
the curvature is not positive definite.
The upper triangle of A is modified in place. By default, the
lower triangle is left unchanged, and the matrices D and E are
actually returned as vectors containing only the diagonal terms.
However, if the keyword OUTFULL is set then full matrices are
returned. This is useful when matrix multiplication will be
performed at the next step.
The modified Cholesky factorization is most stable when pivoting is
performed. If the keyword PIVOT is set, then pivoting is performed
to place the diagonal terms with the largest amplitude in the next
row. The permutation vectors returned in PERMUTE and INVPERMUTE
can be used to apply and reverse the pivoting.
[ i.e., (U(PP,*))(*,PP) applies the permutation and
(U(IPP,*))(*,IPP) reverses it, where PP and IPP are the
permutation and inverse permutation vectors. ]
If the matrix to be factored is very sparse, then setting the
SPARSE keyword may improve the speed of the computations. SPARSE
is more costly on a dense matrix, but only grows as N^2, where as
the standard computation grows as N^3, where N is the rank of the
matrix.
If the CHOLSOL keyword is set, then the output is slightly
modified. The returned matrix A that is returned is structured so
that it is compatible with the CHOLSOL built-in IDL routine. This
involves converting A to being upper to lower triangular, and
multiplying by SQRT(D). Users must be sure to check that all
elements of E are zero before using CHOLSOL.
Parameters
A - upon input, a symmetric NxN matrix to be factored.
Upon output, the upper triangle of the matrix is modified to
contain the factorization.
D - upon output, the diagonal matrix D.
E - upon output, the diagonal error matrix E.
Keyword Parameters
OUTFULL - if set, then A, D and E will be modified to be full IDL
matrices than can be matrix-multiplied. By default,
only the upper triangle of A is modified, and D and E
are returned as vectors.
PIVOT - if set, then diagonal elements of A are pivoted into place
and operated on, in decrease order of their amplitude.
The permutation vectors are returned in the PERMUTE and
INVPERMUTE keywords.
PERMUTE - upon return, the permutation vector which converts a
vector into permuted form.
INVPERMUTE - upon return, the inverse permutation vector which
converts a vector from permuted form back into
standard form.
SPARSE - if set, then operations optimized for sparse matrices are
employed. For large but very sparse matrices, this can
save a significant amount of computation time.
CHOLSOL - if set, then A and D are returned, suitable for input to
the built-in IDL routine CHOLSOL. CHOLSOL is mutually
exclusive with the FULL keyword.
TAU - if set, then use the Tau factor as described in the
"unconventional" modified Cholesky factorization, as
described by Xie & Schlick.
Default: the unconventional technique is not used.
Example
Example 1
---------
a = randomn(seed, 5,5) ;; Generate a random matrix
a = 0.5*(transpose(a)+a) ;; Symmetrize it
a1 = a ;; Make a copy
mcholdc, a1, d, e, /full ;; Factorize it
print, max(abs(e)) ;; This matrix is not positive definite
diff = transpose(a1) ## d ## a1 - e - a
;; Test the factorization by inverting
;; it and subtracting A
print, max(abs(diff)) ;; Differences are small
Example 2
---------
Solving a problem with MCHOLDC and CHOLSOL
a = [[6E,15,55],[15E,55,225],[55E,225,979]]
b = [9.5E,50,237]
mcholdc, a, d, e, /cholsol ;; Factorize matrix, compatible w/ CHOLSOL
if total(abs(e)) NE 0 then $
message, 'ERROR: Matrix A is not positive definite'
x = cholsol(a, d, b) ;; Solve with CHOLSOL
print, x
-0.500001 -0.999999 0.500000
which is within 1e-6 of the true solution.
References
Gill, P. E., Murray, W., & Wright, M. H. 1981
*Practical Optimization*, Academic Press
Schlick, T. & Fogelson, A., "TNPACK - A Truncated Newton
Minimization Package for Large- Scale Problems: I. Algorithm and
Usage," 1992, ACM TOMS, v. 18, p. 46-70. (Alg. 702)
Xie, D. & Schlick, T., "Remark on Algorithm 702 - The Updated
Truncated Newton Minimization Package," 1999, ACM TOMS, v. 25,
p. 108-122.
Modification History
Written, CM, Apr 2001
Added CHOLSOL keyword, CM, 15 Feb 2002
Fix bug when computing final correction factor (thanks to Aaron
Adcock), CM, 13 Nov 2010