Code covered by the BSD License

# UD Factorization & Kalman Filtering

### Gerard Van Willigenburg (view profile)

15 Aug 2011 (Updated )

UD and LD factorization of nonnegative matrices and associated Kalman filter implementations.

[l,r,inz]=sym2lt(p,tolp,toln)
```% SYM2LT: Lower Triangular Cholesky factorization P=LL', P nonnegative symmetric
%
%         [l,r,inz]=sym2lt(p,tolp,toln);
%
%         p     : square symmetric matrix >= 0 (see toln)
%         tolp  : tolerance for positiveness
%         toln  : error generation if to be square rooted part < toln
%                 default: tolp
%                 if toln < 0 all the to be square rooted
%                 parts < tolp are set to zero
%         l     : factorization matrix lower triangular, p=ll'
%         r     : rank of u
%         inz   : row with indices non zero columns u
%
%
%         References: Factorization methods for discrete sequential estimation
%                     1977, Gerald J. Bierman
%
% L.G. van Willigenburg, W.L. de Koning, Update August 2011

function [l,r,inz]=sym2lt(p,tolp,toln)

if nargin > 3; error('  one, two or three input arguments required'); end;
if nargin==1; tolp=1e-12; toln=tolp;
elseif nargin==2; toln=tolp; end; tolp=max(0,tolp);
if (toln>tolp); error('  toln > tolp'); end;
%toln,tolp
[n,m]=size(p);
if n~=m; error('  p must be square'); end;
if n==0; error('  Compatible but empty inputs'); end;
l=zeros(n,n); r=n; inz=[];

for j=1:n-1
l(j,j)=p(j,j);
if l(j,j)<toln; error('  toln or symmetry violated');
elseif l(j,j)<=tolp; l(j,j)=0; a=0; r=r-1; zf=1;
else l(j,j)=sqrt(l(j,j)); a=1./l(j,j); inz=[j,inz]; zf=0;
end
if ~zf
for k=n:-1:j+1;
l(k,j)=a*p(k,j); b=l(k,j);
for i=k:n
p(i,k)=p(i,k)-b*l(i,j);
end
end
else
for k=n:-1:j+1;
l(k,j)=0;
end
end
end
l(n,n)=p(n,n);
if l(n,n)<toln; error('  toln or symmetry violated');
elseif l(n,n)<=tolp; l(n,n)=0; r=r-1;
else l(n,n)=sqrt(l(n,n));
end```