function [R, indef, E] = cholmod(A) % CHOLMOD Modified Cholesky factorization % R = cholmod(A) returns the upper Cholesky factor of A (same as CHOL) % if A is (sufficiently) positive definite, and otherwise returns a % modified factor R with diagonal enries >= sqrt(delta) and % offdiagonal entries <= beta in absolute value, % where delta and beta are defined in terms of size of diagonal and % offdiagonal entries of A and the machine precision; see below. % The idea is to ensure that E = A - R'*R is reasonably small if A is % not too far from being indefinite. If A is sparse, so is R. % The output parameter indef is set to 0 if A is sufficiently positive % definite and to 1 if the factorization is modified. % % The point of modified Cholesky is to avoid computing eigenvalues, % but for dense matrices, MODCHOL takes longer than calling the built-in % EIG, because of the cost of interpreting the code, even though it % only has one loop and uses vector operations. % reference: Nocedal and Wright, Algorithm 3.4 and subsequent discussion % (not Algorithm 3.5, which is more complicated) % original algorithm is due to Gill and Murray, 1974 % written by M. Overton, overton@cs.nyu.edu, last modified Feb 2005 % convenient to work with A = LDL' where D is diagonal, L is unit % lower triangular, and so R = (LD^(1/2))' % if sum(sum(abs(A-A'))) > 0 error('A is not symmetric') end % set parameters governing bounds on L and D (eps is machine epsilon) n = length(A); diagA = diag(A); gamma = max(abs(diagA)); % max diagonal entry xi = max(max(abs(A - diag(diagA)))); % max offidagonal entry delta = eps*(max([gamma+xi, 1])); beta = sqrt(max([gamma, xi/n, eps])); indef = 0; % initialize d and L d = zeros(n,1); if issparse(A) L = speye(n); % sparse identity else L = eye(n); end % there are no inner for loops, everything implemented with % vector operations for a reasonable level of efficiency for j = 1:n K = 1:j-1; % column index: all columns to left of diagonal % d(K) doesn't work in case K is empty djtemp = A(j,j) - L(j,K)*(d(K,1).*L(j,K)'); % C(j,j) in book if j < n I = j+1:n; % row index: all rows below diagonal Ccol = A(I,j) - L(I,K)*(d(K,1).*L(j,K)'); % C(I,j) in book theta = max(abs(Ccol)); % guarantees d(j) not too small and L(I,j) not too big % in sufficiently positive definite case, d(j) = djtemp d(j) = max([abs(djtemp), (theta/beta)^2, delta]); L(I,j) = Ccol/d(j); else d(j) = max([abs(djtemp), delta]); end if d(j) > djtemp % A was not sufficiently positive definite indef = 1; end end % convert to usual output format: replace L by L*sqrt(D) and transpose for j=1:n L(:,j) = L(:,j)*sqrt(d(j)); % L = L*diag(sqrt(d)) bad in sparse case end; R = L'; if nargout == 3 E = A - R'*R; end