From a253483b81f5302103142955dda782089cb5624d Mon Sep 17 00:00:00 2001 From: Michael Orlitzky Date: Sat, 11 May 2013 15:46:12 -0400 Subject: [PATCH] Add the cholesky_inf() function. --- cholesky_inf.m | 62 ++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 62 insertions(+) create mode 100644 cholesky_inf.m diff --git a/cholesky_inf.m b/cholesky_inf.m new file mode 100644 index 0000000..4adbf34 --- /dev/null +++ b/cholesky_inf.m @@ -0,0 +1,62 @@ +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 -- 2.43.2