## 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