%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Constructs an LU decomposition of a symmetric positive definite matrix
% Input:
% A - matrix
% Output:
% L - a matrix such that L * L' = A
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function L = cholesky(A)
% check that A is a square matrix
if size(A,1) ~= size(A,2)
error( 'Error: Matrix is not square' );
end
% check if A is symmetric
if norm(A-A')>0.00001
error( 'Error: Matrix is not symmetric' );
end
n = size(A,1); % size of A
A0=A;
L=eye(n);
for i=1:n
% check if A is is positively definite
if A0(i,i)<=0
error('Matrix A is not positive definite or singular');
end
bi=A0(i+1:n,i);
Li=eye(n);
Li(i,i)=sqrt(A0(i,i));
Li(i+1:n,i)=bi/Li(i,i);
A1=eye(n);
A1(i+1:n,i+1:n) = A0(i+1:n,i+1:n) - (bi*bi')/A0(i,i);
A0=A1;
L = L * Li;
end
end