;+ ; 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 ; ; $Id: mcholdc.pro,v 1.6 2010/11/13 09:45:55 cmarkwar Exp $ ; ;- ; Copyright (C) 2001, 2002, 2010, Craig Markwardt ; This software is provided as is without any warranty whatsoever. ; Permission to use, copy, modify, and distribute modified or ; unmodified copies is granted, provided this copyright and disclaimer ; are included unchanged. ;- pro mcholdc, a, d, e, outfull=full, sparse=sparse, pivot=pivot, $ permute=pp, invpermute=ipp, tau=tau0, cholsol=cholsol if n_params() EQ 0 then begin message, 'USAGE: MCHOLDC, A, D, E [, /OUTFULL, /SPARSE, /PIVOT, '+$ '/CHOLSOL, TAU=, PERMUTE=, INVPERMUTE= ]', /info return endif ;; Test for proper dimensions sz = size(a) if sz(0) EQ 0 then begin d = a a = a(0)*0 + 1 e = a(0)*0 return endif if sz(0) NE 2 then $ message, 'ERROR: Matrix A must be two dimensional' if sz(1) NE sz(2) then $ message, 'ERROR: Matrix A must be square' if n_elements(tau0) GT 0 then tau = tau0(0) n = sz(1) zero = a(0)*0 one = zero + 1 eps = zero + 1e-6 ;; Gamma and Xi are the max diagonal and off-diagonal components diag = lindgen(n)*n + lindgen(n) gamma = max(abs(a(diag))) xi = zero for i = 0L, n-2 do begin xi = max( [xi, max(abs(a(i+1:*,i)))] ) endfor eps1 = max([gamma,xi]) * eps del = max([eps, eps1]) ;; Compute the bound on the diagonal elements bound = max([gamma, xi/sqrt(n^2-1), eps]) ;; Compute output arrays d = replicate(zero, n) e = d pp = lindgen(n) for j = 0, n-1 do begin ;; Pivoting - search for max element of abs(diag(A)) if keyword_set(pivot) then begin maxa = max(abs(a(diag(j:*))), jmax) if jmax GT 0 then begin jmax = jmax + j nn = lindgen(n) j1 = (nn*n+j ) < (nn+j*n) j2 = (nn*n+jmax) < (nn+jmax*n) temp = j1(j) & j1(j) = j1(jmax) & j1(jmax) = temp ;; Exchange the row/columns temp = a(j1) & a(j1) = a(j2) & a(j2) = temp ;; Change the permutation vector temp = pp(j) & pp(j) = pp(jmax) & pp(jmax) = temp endif endif ;; Compute update to the L matrix, in place if j GT 0 AND j LT n-1 then begin ;; The sparse path computes the same thing, but is faster if ;; very few components of the matrix are non-zero. if keyword_set(sparse) then begin ajk = reform(a(j,0:j-1)*d(0:j-1), /overwrite) wh = where(ajk NE 0, nk) if nk GT 0 then $ a(j+1:*,j) = a(j+1:*,j) - a(j+1:*,wh) # ajk(wh) endif else begin a(j+1:*,j) = a(j+1:*,j) - $ a(j+1:*,0:j-1) # reform(a(j,0:j-1)*d(0:j-1), /overwrite) endelse endif ;; Compute correction to diagonal element thj2 = max(a(j:*,j))^2 ;; Compute the unusual modified factorization of Xie and ;; Schlick, or else default to the standard Gill, Murray & Wright if n_elements(tau) GT 0 then begin ww = a(j,j) + tau(0) if ww GT del then d(j) = max([ww, thj2/bound]) $ else if abs(ww) LE del then d(j) = del $ else d(j) = min([ww,-thj2/bound]) endif else begin d(j) = max([del, abs(a(j,j)), thj2/bound]) endelse e(j) = d(j) - a(j,j) ;; Apply corrections if j LT n-1 then begin i = lindgen(n-1-j)+j+1 a(i,i) = a(i,i) - a(i,j)^2/d(j) a(j+1:*,j) = a(j+1:*,j) / d(j) endif endfor ;; Invert the permutation vector ipp = sort(pp) ;; Expand to a full matrix if requested if keyword_set(cholsol) then begin d = sqrt(d) for j = 0, n-2 do a(j+1:*,j) = a(j+1:*,j)*d(j) a = transpose(temporary(a)) endif else if keyword_set(full) then begin for j = 0, n-1 do a(0:j,j) = 0 a(diag) = 1 dd = a*0 ee = dd dd(diag) = d ee(diag) = e d = dd e = ee endif end