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))

## CROSS-REFERENCE INFORMATION

This function calls:
 deal Deal inputs to outputs. error Display message and abort function. tril Extract lower triangular part. warning Display warning message; disable or enable warning messages.
This function is called by:
 wmnormrnd Random vectors from a multivariate Normal distribution

## 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 %
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