Code covered by the BSD License  

Highlights from
UD Factorization & Kalman Filtering

UD Factorization & Kalman Filtering

by

 

15 Aug 2011 (Updated )

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

[uti,r]=utinv(ut,tolp,toln)
% UTINV : inverse of upper triangular ut
%
%         function [uti,r]=utinv(ut,tolp,toln)
%
%         ut  : U'DU factorization
%         tolp: tolerance for positiveness
%               default 1e-12, see also toln
%         toln: error generation if diagonal part < toln,  default: tolp
%               if toln < 0 all the diagonal parts < max(0,tolp) are set to zero
%
%         uti : inv(ut)
%         r   : rank(di)
%
%         Used to compute : inv(p)=inv(u)'*inv(u)
%         [uti]=utinv(ut) computes inv(u)
%         Then: inv(p)=utt2sym(utinv(sym2ud(p)))
%
%         See also utt2sym, ut2sym, sym2ut
%
%         References: Factorization methods for discrete sequential estimation
%                     1977, Gerald J. Bierman
%
% L.G. van Willigenburg, W.L. de Koning, Update August 2011

  function [uti,r]=utinv(ut,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;

  [n,m]=size(ut);
  if n~=m; error(' ut must be square'); end;
  if n==0; error('  Compatible but empty inputs'); end;

  r=n;
  for i=1:n
    if ut(i,i)<toln; error('  toln violated');
    elseif ut(i,i)<=tolp; uti(i,i)=0; r=r-1;
    else uti(i,i)=1/ut(i,i); end
  end;

  for j=2:n
    jm1=j-1;
    for k=1:jm1
      sum=0;
      for i=k:jm1
        sum=sum-uti(k,i)*ut(i,j);
      end
      uti(k,j)=sum*uti(j,j);
    end
  end

Contact us