From: Michael Orlitzky <michael@orlitzky.com>
Date: Sat, 11 May 2013 19:46:12 +0000 (-0400)
Subject: Add the cholesky_inf() function.
X-Git-Url: https://gitweb.michael.orlitzky.com/?a=commitdiff_plain;h=a253483b81f5302103142955dda782089cb5624d;p=octave.git

Add the cholesky_inf() function.
---

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