cholesky_inf.m

   1 function R = cholesky_inf(A, tolerance)
   2   %
   3   % Peform the Cholesky-Infinity factorization of ``A``. This differs
   4   % from the ordinary Cholesky factorization in that it handles
   5   % matrices which are not positive-definite.
   6   %
   7   % If a non-positive diagonal entry appears, it is replaced with
   8   % ``inf`` and the rest of that row is set to zero. This will force a
   9   % zero in the corresponding row of the solution vector.
  10   %
  11   % This function is based on the description of MATLAB's `cholinc(A,
  12   % 'inf')` function which has been deprecated and removed.
  13   %
  14   % INPUT:
  15   %
  16   %   - ``A`` -- The matrix to factor.
  17   %
  18   %   - ``tolerance`` -- (default: eps) The zero tolerance used to test
  19   %                      the diagonal entries. If `A(i,i) <= tolerance`,
  20   %                      then we consider `A(i,i) to be zero.
  21   %
  22   % OUTPUT:
  23   %
  24   %   - ``R`` -- The upper-triangular Cholesky factor of ``A``. ``R``
  25   %              will be sparse if ``A`` is.
  26   %
  27   if (nargin < 2)
  28     tolerance = eps;
  29   end
  30
  31   [m,n] = size(A);
  32
  33   if (issparse(A))
  34      % If `A` is sparse, `R` should be too.
  35     R = sparse(m,n);
  36   else
  37     R = zeros(m,n);
  38   end
  39
  40   for i = [ 1 : n ]
  41     col_above = R(1:i-1, i);
  42     R(i,i) = A(i,i) - col_above'*col_above;
  43
  44     if (R(i,i) <= tolerance)
  45       % A is not positive-definite. Rather than throw an error here,
  46       % we set this pivot to 'inf' and the rest of the row to zero.
  47       R(i,i) = inf;
  48       R(i, i+1 : end) = 0;
  49     else
  50       % Normal case with the pivot greater than zero.
  51       R(i,i) = sqrt(R(i,i));
  52
  53       for j = [ i+1 : n ]
  54         row1 = R(1:i-1, i);
  55         row2 = R(1:i-1, j);
  56         R(i,j) = (A(i,j) - row1'*row2)/R(i,i);
  57       end
  58     end
  59
  60   end
  61
  62 end