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