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
Deal inputs to outputs. | |
Display message and abort function. | |
Extract lower triangular part. | |
Display warning message; disable or enable warning messages. |
Random vectors from a multivariate Normal distribution |
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
Comments or corrections to the WAFO group