Problem with Cholesky's decomposition of a positive semi-definite matrix

Hello everyone. I need to perform the Cholesky decomposition of a positive semi-definite matrix (M) as M=R’R. The usual chol function does not work for me, since it only works with positive definite matrices.

(Removed cholesky decomposition code.)

I also found the following code, which performs another decomposition over the matrix, but instead of providing the R matrix as in the previous paragraph, it gives two matrices such that M=LDL’. If someone could tell me how to adapt this function to return the matrix R instead of L and D I would be extremely thankful.

 function [L, DMC, P, D] = modchol_ldlt(A,delta)
%modchol_ldlt  Modified Cholesky algorithm based on LDL' factorization.
%   [L D,P,D0] = modchol_ldlt(A,delta) computes a modified
%   Cholesky factorization P*(A + E)*P' = L*D*L', where 
%   P is a permutation matrix, L is unit lower triangular,
%   and D is block diagonal and positive definite with 1-by-1 and 2-by-2 
%   diagonal blocks.  Thus A+E is symmetric positive definite, but E is
%   not explicitly computed.  Also returned is a block diagonal D0 such
%   that P*A*P' = L*D0*L'.  If A is sufficiently positive definite then 
%   E = 0 and D = D0.  
%   The algorithm sets the smallest eigenvalue of D to the tolerance
%   delta, which defaults to sqrt(eps)*norm(A,'fro').
%   The LDL' factorization is compute using a symmetric form of rook 
%   pivoting proposed by Ashcraft, Grimes and Lewis.

%   Reference:
%   S. H. Cheng and N. J. Higham. A modified Cholesky algorithm based
%   on a symmetric indefinite factorization. SIAM J. Matrix Anal. Appl.,
%   19(4):1097-1110, 1998. doi:10.1137/S0895479896302898,

%   Authors: Bobby Cheng and Nick Higham, 1996; revised 2015.

if ~ishermitian(A), error('Must supply symmetric matrix.'), end
if nargin < 2, delta = sqrt(eps)*norm(A,'fro'); end

n = max(size(A));

[L,D,p] = ldl(A,'vector'); 
DMC = eye(n);

% Modified Cholesky perturbations.
k = 1;
while k <= n

      if k == n || D(k,k+1) == 0 % 1-by-1 block

         if D(k,k) <= delta
            DMC(k,k) = delta;
         else
            DMC(k,k) = D(k,k);
         end
         k = k+1;
      
      else % 2-by-2 block

         E = D(k:k+1,k:k+1);
         [U,T] = eig(E);
         for ii = 1:2
             if T(ii,ii) <= delta
                T(ii,ii) = delta;
             end
         end
         temp = U*T*U';
         DMC(k:k+1,k:k+1) = (temp + temp')/2;  % Ensure symmetric.
         k = k + 2;

      end

end

if nargout >= 3, P = eye(n); P = P(p,:); end

The only idea that I have to do this by myself is to add a small value to the diagonal of the matrix M and then use chol. I don’t like this, since I don’t consider it very scientific and I have no idea on how the results are altered by this, so if someone can offer a different alternative to my problem which involves chol and not adding a differential value to the diagonal, I would be thankful too.

Thanks for reading.

[Matt J]

But which Cholesky decomposition should it give you? The Cholesky decomposition of a semi-definite matrix can be non-unique. How would you have the solution given to you resolve the ambiguities?

[Jaime de la Mota]

The code that I have posted implements what is described in this article, where a method is developed to add a perturbation to a non positive definite matrix making it possitive definite and then applying the LDL' transformation on said positive definite matrix. My question is as follows: Can I perform the R'R decomposition in said matrix instead of the LDL'?

[Matt J]

Can I perform the R'R decomposition in said matrix instead of the LDL'?

I understood your question. But you continue to talk about R as if it is unique, and my question was, how would you choose which "R" if the space of possible decompositions was non-unique, as it can be for a semi-definite matrix?

[Matt J]

For example, suppose the given matrix is

M =

     0     0
     0     1

All R of the form

R =

     0     a
     0     b

with a^2+b^2=1 give valid decompositions. Which of these infinite choices would be your choice in this case?

[John D'Errico]

Please don't post MATLAB supplied code in your questions.

[Bruno Luong]

But the modchol_ldlt code is not MATLAB supplied

Answer 1

[Bruno Luong]

M=LDL’. If someone could tell me how to adapt this function to return the matrix R instead of L and D I would be extremely thankful.

Wasn't simply

R = L*sqrtm(D)
% M = R*R' = L*D*L'

D is 1x1 and 2x2 block diagonal positive by construction,no?

EDIT: here is the code of this idea

% Generate A symmetric but almost positive
A = randn(10);
A = A*A';
emin = min(eig(A));
A = A-1.01*emin*eye(size(A))

% This will return flag > 0 and R incomplete Chol decomposition
[R,flag] = chol(A)

% Call Cheng and N. J. Higham approximation
[L, D, P] = modchol_ldlt(A); % https://github.com/higham/modified-cholesky

% Compute "sqrtD", the cholesky (square root) of D
n = size(A,1);
d0 = D(1:n+1:end);        % diagonal
d1 = D(2:n+1:end);        % off diagonal
sqrtD = diag(sqrt(d0));
i = find(d1);
j = i+(i-1)*n;            % linear index of subindexes (i,i)
j1 = j+1;                 % (i+1;i)
k = j1+n;                 % (i+1,i+1)
sqrtD(j1) = D(j1)./sqrtD(j);
sqrtD(k) = sqrt(D(k)-D(j1).^2./D(j)); % you might double protect with sqrt(max(...,0))

R = P'*L*sqrtD; % L*sqrtD is lower triangular

Apos = R*R'

% Check how close Apos to A
norm(Apos-A,'fro')/norm(A,'fro')

[Matt J]

In general, though, this cannot be a stable solution for semi-definite matrices because the space of decompositions can be non-unique and connected.

[Bruno Luong]

The purpose here is to find the spd of a matrix close to the original matrix.

There is a lot matrix decompositions out there that is not unique, EIG for starter.

I like Cheng & Higham approximation, it seems like a good and quick way of correcting SPD.

Answer 2

[Ahmed Mahfouz]

Hello Jamie,

You probably got your answer already or even moved on from this question, but I will leave this answer here for those who might encounter the same problem in the future.

The simple way to do this factorization is to use the following matlab function:

R = cholcov(M);

You need to make sure that M is indeed PSD, otherwise, the output of the aformentioned code will be an empty matrix.