REGSOL returns the regularized Cholesky factorization of a given n-by-n symmetric and positive semi-defnite matrix A given a vector b, i.e. points out to a regularized solution of the system Ax=b.

regsol.m

function [L, J, NJ, D, x] = regsol(A,b)
%REGSOL returns the regularized Cholesky factorization of a given
%n-by-n symmetric and positive semi-defnite matrix A given a vector b, i.e.
%points out to a regularized solution of the system Ax=b.
%
%Syntax:
% [L, J, NJ, x, D] = regsol(A,b)
%
%Input arguments:
%A,b : A positive semidefinite matrix and a vector, of compatible
%dimensions. We need to determine a regularized solution for the linear
%system Ax=b.
%
%Output arguments:
%L : A regularized Choleasky factor of A
%x : A regularized solution for the linear system Ax=b, in the sense that
%    there exists a set of indices I, such that A(I,I) has the same rank as
%    A and is positive definite.
%J : A set of indices such that x(J)=b(J)
%NJ: The complement of J.
%
%Reference:
%Wu Li and John Swetis, "Regularized Newton Methods for Minimizaton of
%Convex Quadratic Splines with Singular Hessians," Kluwer Academic
%Publishers, 1998
%
%


%% Check Input
n=size(A,1);

assert(n==size(A,2),'A is not a square matrix');
assert(size(b,2)==1,'b is not a column-vector');
assert(n==size(b,1),'Dimensions of A and b are inconsistent');

epsTOL=1e-10;

%% Step 1
J=zeros(1,0);
jind=0;
L=zeros(n,n);
if A(1,1) <= epsTOL
    jind=1;
    J(1)=1;
    L(1,1)=1;
    L(2:n)=0;
else
    L(1,1)=sqrt(A(1,1));
    for i=2:n
        L(i,1) = A(i,1)/L(1,1);
    end
end


%% Step 2
for k=2:n-1
    g=A(k,k);
    for j=1:k-1
        g = g - L(k,j)^2;
    end    
    if g<=epsTOL %instead of g<=0 !!! works better...
        jind=jind+1;
        J(jind)=k; %J = [J k];
        L(k,k)=1;
        for i=1:k-1
            L(k,i)=0;
        end
        for i=k+1:n
            L(i,k)=0;
        end
    else
        L(k,k) = sqrt(g);
        for i=k+1:n
            h=A(i,k);
            for j=1:k-1
                h = h - L(k,j)*L(i,j);
            end
            L(i,k) = h/L(k,k);
        end
    end
end


%% Step 3
f = A(n,n);
for j=1:n-1
    f = f - L(n,j)^2;
end
if f<=epsTOL
    jind=jind+1;
    J(jind)=n;
    J = [J n];
    L(n,n)=1;
    for i=1:n-1
        L(n,i) = 0;
    end
else
    L(n,n) = sqrt(f);
end


%% Step 4
%b_(J) = b(J);
%b_(NJ) = b(NJ)-A(NJ,J)*b(J);

b_ = b;
NJ=setdiff(1:n,J);

% iterate over the indices in NJ:
for q = 1:length(NJ)
    sum = 0;
    for i = 1:length(J)
        sum = sum + A(NJ(q),J(i)) * b(J(i));
    end
    b_(NJ(q)) = b(NJ(q)) - sum;
end

D=zeros(n,1);
D(NJ)=1;
D=diag(D);
%% Solve (Regularized) Cholesky System
% L is a lower-triangular matrix.
% We have to solve: L*L'*x = b_
%   Lz = b
%  L'x = z
if nargout<=4
    return
end

z=zeros(n,1);
z(1)=b_(1)/L(1,1);
for k=2:n
    sum = b_(k);
    for j=1:k-1
        sum = sum - L(k,j)*z(j);
    end
    z(k)=sum/L(k,k);
end

x=zeros(n,1);
x(n)=z(n)/L(n,n);
for k=n-1:-1:1
    sum = z(k);
    for j=k+1:n
        sum = sum - L(j,k)*x(j);
    end
    x(k) = sum/L(k,k);
end



end  % end of function