Add the cholesky_inf() function.
[octave.git] / 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