Code covered by the BSD License  

Highlights from
inv4spd.m

inv4spd.m

by

 

Inversion algorithm for symmetric, positive definite matrix M. Minv = inv4spd(M).

inv4spd(A)
function P = inv4spd(A)
% find the inverse of a symmetric positive definite matrix, 
% inv4spd.m 6/24/2009 developed by D. Hsu, revised on 8/05/2009
% A = L*U = L*L'  with U = L'
% P = inv(A) = inv(U)*inv(L)
% (1) find lower triangular factor of A --- L
n = length(A);
for k = 1:n
  A(k:n,k) = A(k:n,k)/sqrt(A(k,k));
  for j=k+1:n
      A(j:n,j) = A(j:n,j) - A(j:n,k)*A(j,k);
  end
end
L = tril(A);  % lower triangular factor of A
% (2) inverse of lower triangular factor --- inv(L)
for j=1:n
   % if L(j,j)==0  % or if abs(L(j,j))<1.e-12
   %     error('L could not be inverted')
   % end
    L(j,j) = 1/L(j,j);      
    i = 1:j-1;
    L(j,i) = -L(j,i)*L(i,i)*L(j,j) 
end
% take transpose of inv(L) --- inv(U) 
U = L';    %%
% (3) take product P = inv(U)*inv(L), P is Upper Triangular 
for j=n:-1:1
      P(1:j,j) = U(1:j,j:n)*U(j,j:n)';  %%  [j x (n-j)] x [(n-j) x 1]
end
% form symmetric matrix from upper triangular matrix P
for j=1:n-1
  P(j+1:n,j)=P(j,j+1:n)';
end

 

Contact us