function R = cholesky_inf(A, tolerance) % % Peform the Cholesky-Infinity factorization of ``A``. This differs % from the ordinary Cholesky factorization in that it handles % matrices which are not positive-definite. % % If a non-positive diagonal entry appears, it is replaced with % ``inf`` and the rest of that row is set to zero. This will force a % zero in the corresponding row of the solution vector. % % This function is based on the description of MATLAB's `cholinc(A, % 'inf')` function which has been deprecated and removed. % % INPUT: % % - ``A`` -- The matrix to factor. % % - ``tolerance`` -- (default: eps) The zero tolerance used to test % the diagonal entries. If `A(i,i) <= tolerance`, % then we consider `A(i,i) to be zero. % % OUTPUT: % % - ``R`` -- The upper-triangular Cholesky factor of ``A``. ``R`` % will be sparse if ``A`` is. % if (nargin < 2) tolerance = eps; end [m,n] = size(A); if (issparse(A)) % If `A` is sparse, `R` should be too. R = sparse(m,n); else R = zeros(m,n); end for i = [ 1 : n ] col_above = R(1:i-1, i); R(i,i) = A(i,i) - col_above'*col_above; if (R(i,i) <= tolerance) % A is not positive-definite. Rather than throw an error here, % we set this pivot to 'inf' and the rest of the row to zero. R(i,i) = inf; R(i, i+1 : end) = 0; else % Normal case with the pivot greater than zero. R(i,i) = sqrt(R(i,i)); for j = [ i+1 : n ] row1 = R(1:i-1, i); row2 = R(1:i-1, j); R(i,j) = (A(i,j) - row1'*row2)/R(i,i); end end end end