Home > wafo > misc > genchol.m

genchol

PURPOSE ^

Generalized Cholesky factorization

SYNOPSIS ^

[L,P,r] = genchol(A,tol)

DESCRIPTION ^

 GENCHOL Generalized Cholesky factorization 
   
  CALL: [L,P,r] = genchol(A,tol); 
   
   L  = lower triangular matrix, Cholesky factor 
   P  = permutation vector 
   r  = matrix rank, i.e., the number of eigenvalues larger than tol. 
        This is an estimate of the number of linearly 
        independent rows or columns of a matrix A. 
   A  = symmetric and semi-positive definite matrix 
  tol = tolerance used in factorization, i.e., norm(A(P,P)- L*L.') <= tol   
  
  GENCHOL computes the generalized Cholesky decomposition, 
            
    A(P,P) = L*L.'  
   
  where L is lower triangular and P is a permutation vector. 
  
  Example 
   H = hilb(10); 
   tol   = 1e-6; 
   [L,P] = genchol(H,tol); 
   spy(L*L.'-H(P,P)) 
  
  See also  chol

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

001 function [L,P,r] = genchol(A,tol) 
002 %GENCHOL Generalized Cholesky factorization 
003 %  
004 % CALL: [L,P,r] = genchol(A,tol); 
005 %  
006 %  L  = lower triangular matrix, Cholesky factor 
007 %  P  = permutation vector 
008 %  r  = matrix rank, i.e., the number of eigenvalues larger than tol. 
009 %       This is an estimate of the number of linearly 
010 %       independent rows or columns of a matrix A. 
011 %  A  = symmetric and semi-positive definite matrix 
012 % tol = tolerance used in factorization, i.e., norm(A(P,P)- L*L.') <= tol   
013 % 
014 % GENCHOL computes the generalized Cholesky decomposition, 
015 %           
016 %   A(P,P) = L*L.'  
017 %  
018 % where L is lower triangular and P is a permutation vector. 
019 % 
020 % Example 
021 %  H = hilb(10); 
022 %  tol   = 1e-6; 
023 %  [L,P] = genchol(H,tol); 
024 %  spy(L*L.'-H(P,P)) 
025 % 
026 % See also  chol   
027    
028 %error(nargoutchk(2,2,,nargout)) 
029  
030 [m,n] = size(A);  
031 if (m ~= n),   error('input matrix must be square'); end 
032 if nargin<2|isempty(tol),   tol = 1e-12; end 
033 %tol = tol/(100*n*n);%(4*sqrt(n)); 
034  
035  
036  
037 L = tril(A);  
038 %D = diag(A);  
039 %Dmax = max(D); 
040 %tol = Dmax*tol; 
041 %localTolerance = max(eps*Dmax,tol*1e-3); 
042 localTolerance = 0;%tol/(20*n^2); 
043  
044 P = 1:n; % permuation vector 
045 %x = P; 
046  
047 k = 1; 
048 nullity = 0; 
049 while (k<=n) 
050    % Find next pivot 
051    D = diag(L,-nullity); 
052    k0 = k-nullity; 
053     n0 = n-nullity; 
054    [big,imax] = max(D(k0:n0)); 
055     
056    if (big>tol), %D(imax)>tol) 
057       imax = imax+k-1; 
058       if imax~=k 
059          % Swap the new pivot at imax with k 
060          P([imax,k]) = P([k,imax]);% [P(imax),P(k)] = deal(P(k),P(imax)); 
061          L = rowcolChange(L,k,imax,n,nullity);      
062       end 
063       L(k,k0) = sqrt(L(k,k0)); % STD(Xk|X1,X2,...,Xk-1) 
064       k1 = k; 
065       for i = k+1:n, 
066          %disp(i) 
067          %tmp = 0; 
068          %for j = 1:k-1,  
069          %Cov(Xi,Xj|X1,X2,...Xk-2)*Cov(Xj,Xk|X1,X2,...Xk-2)/Cov() 
070          %tmp = tmp + L(i,j) * L(k,j);  
071          %end 
072          %L(i,k) = L(i,k)-tmp; % Cov(Xi,Xk|X1,X2,...Xk-1) 
073          %L(i,k) = L(i,k) / L(k,k);             
074           
075          j =  1:k0-1; 
076          % Cov(Xi,Xk|X1,X2,...Xk-1)/STD(Xk|X1,X2,...Xk-1) 
077          L(i,k0) = (L(i,k0)-sum(L(i,j).*L(k1,j)))/L(k1,k0); 
078          % Var(Xi|X1,X2,...Xk) 
079          i0 = i-nullity; 
080          L(i,i0) = L(i,i0) - L(i,k0)*L(i,k0);         
081           
082          if (L(i,i0)<=localTolerance) % & norm(L(nnDet,k+1:nnDet),'inf')<=tol)) 
083             %if (D(i)<=-sqrt(tol)) % make sure we are not too restrictive 
084             %   warning('Matrix is not positive semi-definite!') 
085             %end 
086       
087             if (k+1<i) 
088                % Swap the singular pivot at i with k+1 
089                P([i,k+1]) = P([k+1,i]); 
090                L = rowcolChange(L,k+1,i,n,nullity);   
091             end 
092             % shift 
093             if nullity>0 
094                n0 = n-nullity; 
095                L(k+1:n,k0+1:n0) = L(k+1:n,k0+2:n0+1);  
096             else 
097                L(k+1:n,k+1:n-1) = L(k+1:n,k+2:n); 
098                L(n,n) = 0; 
099             end 
100             nullity = nullity+1; 
101             k = k+1; 
102          end 
103      
104       end       
105    else 
106       if (big<=-sqrt(tol)) 
107          warning('Matrix is not positive semi-definite!') 
108       end 
109       k0 = k-nullity; 
110       L(k:end,k0:end) = 0; 
111       nullity = n-k0+1; 
112       break; 
113    end 
114    k = k+1; 
115 end  %while k 
116 r = n-nullity; 
117 %flop 
118 %nDet 
119 return 
120  
121  
122  
123 function L = rowcolChange(L,i,j,n,nullity) 
124 %rowcolChange exchange column and rows of i and j, 
125 %but only the lower triangular part 
126  
127 if (i<j) 
128    k = i; 
129    imax = j; 
130 else 
131    k = j; 
132    imax = i; 
133 end     
134     
135    %for  
136    k0 = k-nullity; 
137    if 1<k0 
138        iz=1:k0-1; 
139       [L(k,iz),L(imax,iz)] = deal(L(imax,iz),L(k,iz)); 
140    end     
141     
142    imax0 = imax-nullity; 
143    [L(k,k0),L(imax,imax0)] = deal(L(imax,imax0),L(k,k0)); 
144    %for 
145    iz=k+1:imax-1; 
146    iz0= k0+1:imax0-1; 
147    if any(iz) 
148       [L(iz,k0),L(imax,iz0)] = deal(L(imax,iz0).',L(iz,k0).'); 
149    end 
150 %   for  
151 iz=imax+1:n; 
152 if any(iz) 
153       [L(iz,k0),L(iz,imax0)] = deal(L(iz,imax0),L(iz,k0)); 
154   end

Mathematical Statistics
Centre for Mathematical Sciences
Lund University with Lund Institute of Technology

Comments or corrections to the WAFO group


Generated on Thu 06-Oct-2005 02:21:16 for WAFO by m2html © 2003