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

61 次查看（过去 30 天）
Jaime de la Mota
Jaime de la Mota 2019-8-26
回答： Ahmed Mahfouz 2023-10-30
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.
  6 个评论
显示 4更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

回答（2 个）

Bruno Luong
Bruno Luong 2019-8-27
编辑：Bruno Luong 2019-8-28
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')
  2 个评论
Matt J
Matt J 2019-8-27
编辑：Matt J 2019-8-27
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
Bruno Luong 2019-8-28
编辑：Bruno Luong 2019-8-28
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.

请先登录，再进行评论。

Ahmed Mahfouz
Ahmed Mahfouz 2023-10-30
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.

类别

在 Help Center 和 File Exchange 中查找有关 Linear Algebra 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by