sparseNull returns computes the sparse Null basis of a matrix N = sparseNull(S, tol) Computes a basis of the null space for a sparse matrix. For sparse matrixes this is much faster than using null. It does however have lower numerical accuracy. N is itself sparse and not orthonormal. So in this way it is like using N = null(S, 'r'), except of course much faster. Jan Schellenberger 10/20/2009 based on this: http://www.mathworks.com/matlabcentral/fileexchange/11120-null-space-of-a-sparse-matrix
0001 function N = sparseNull(S, tol) 0002 % sparseNull returns computes the sparse Null basis of a matrix 0003 % 0004 % N = sparseNull(S, tol) 0005 % 0006 % Computes a basis of the null space for a sparse matrix. For sparse 0007 % matrixes this is much faster than using null. It does however have lower 0008 % numerical accuracy. N is itself sparse and not orthonormal. So in this 0009 % way it is like using N = null(S, 'r'), except of course much faster. 0010 % 0011 % Jan Schellenberger 10/20/2009 0012 % based on this: 0013 % http://www.mathworks.com/matlabcentral/fileexchange/11120-null-space-of-a-sparse-matrix 0014 if nargin <2 0015 tol = 1e-9; 0016 end 0017 [SpLeft, SpRight] = spspaces(S,2, tol); 0018 N = SpRight{1}(:,SpRight{3}); 0019 N(abs(N) < tol) = 0; 0020 0021 0022 %%%%%%%%%%%%%% code from website. I did not write this myself. -Jan 0023 %%%% 0024 0025 function [SpLeft, SpRight] = spspaces(A,opt,tol) 0026 % PURPOSE: finds left and right null and range space of a sparse matrix A 0027 % 0028 % --------------------------------------------------- 0029 % USAGE: [SpLeft, SpRight] = spspaces(A,opt,tol) 0030 % 0031 % INPUT: 0032 % A a sparse matrix 0033 % opt spaces to calculate 0034 % = 1: left null and range space 0035 % = 2: right null and range space 0036 % = 3: both left and right spaces 0037 % tol uses the tolerance tol when calculating 0038 % null subspaces (optional) 0039 % 0040 % OUTPUT: 0041 % SpLeft 1x4 cell. SpLeft = {} if opt =2. 0042 % SpLeft{1} an invertible matrix Q 0043 % SpLeft{2} indices, I, of rows of the matrix Q that 0044 % span the left range of the matrix A 0045 % SpLeft{3} indices, J, of rows of the matrix Q that 0046 % span the left null space of the matrix A 0047 % Q(J,:)A = 0 0048 % SpLeft{4} inverse of the matrix Q 0049 % SpRight 1x4 cell. SpRight = {} if opt =1. 0050 % SpLeft{1} an invertible matrix Q 0051 % SpLeft{2} indices, I, of rows of the matrix Q that 0052 % span the right range of the matrix A 0053 % SpLeft{3} indices, J, of rows of the matrix Q that 0054 % span the right null space of the matrix A 0055 % AQ(:,J) = 0 0056 % SpLeft{4} inverse of the matrix Q 0057 % 0058 % COMMENTS: 0059 % uses luq routine, that finds matrices L, U, Q such that 0060 % 0061 % A = L | U 0 | Q 0062 % | 0 0 | 0063 % 0064 % where L, Q, U are invertible matrices, U is upper triangular. This 0065 % decomposition is calculated using lu decomposition. 0066 % 0067 % This routine is fast, but can deliver inaccurate null and range 0068 % spaces if zero and nonzero singular values of the matrix A are not 0069 % well separated. 0070 % 0071 % WARNING: 0072 % right null and rang space may be very inaccurate 0073 % 0074 % Copyright (c) Pawel Kowal (2006) 0075 % All rights reserved 0076 % LREM_SOLVE toolbox is available free for noncommercial academic use only. 0077 % pkowal3@sgh.waw.pl 0078 0079 if nargin<3 0080 tol = max(max(size(A)) * norm(A,1) * eps,100*eps); 0081 end 0082 0083 switch opt 0084 case 1 0085 calc_left = 1; 0086 calc_right = 0; 0087 case 2 0088 calc_left = 0; 0089 calc_right = 1; 0090 case 3 0091 calc_left = 1; 0092 calc_right = 1; 0093 end 0094 0095 [L,U,Q] = luq(A,0,tol); 0096 0097 if calc_left 0098 if ~isempty(L) 0099 LL = L^-1; 0100 else 0101 LL = L; 0102 end 0103 S = max(abs(U),[],2); 0104 I = find(S>tol); 0105 if ~isempty(S) 0106 J = find(S<=tol); 0107 else 0108 J = (1:size(S,1))'; 0109 end 0110 SpLeft = {LL,I,J,L}; 0111 else 0112 SpLeft = {}; 0113 end 0114 if calc_right 0115 if ~isempty(Q) 0116 QQ = Q^-1; 0117 else 0118 QQ = Q; 0119 end 0120 S = max(abs(U),[],1); 0121 I = find(S>tol); 0122 if ~isempty(S) 0123 J = find(S<=tol); 0124 else 0125 J = (1:size(S,2))'; 0126 end 0127 SpRight = {QQ,I,J,Q}; 0128 else 0129 SpRight = {}; 0130 end 0131 0132 0133 function [L,U,Q] = luq(A,do_pivot,tol) 0134 % PURPOSE: calculates the following decomposition 0135 % 0136 % A = L |Ubar 0 | Q 0137 % |0 0 | 0138 % 0139 % where Ubar is a square invertible matrix 0140 % and matrices L, Q are invertible. 0141 % 0142 % --------------------------------------------------- 0143 % USAGE: [L,U,Q] = luq(A,do_pivot,tol) 0144 % INPUT: 0145 % A a sparse matrix 0146 % do_pivot = 1 with column pivoting 0147 % = 0 without column pivoting 0148 % tol uses the tolerance tol in separating zero and 0149 % nonzero values 0150 % 0151 % OUTPUT: 0152 % L,U,Q matrices 0153 % 0154 % COMMENTS: 0155 % based on lu decomposition 0156 % 0157 % Copyright (c) Pawel Kowal (2006) 0158 % All rights reserved 0159 % LREM_SOLVE toolbox is available free for noncommercial academic use only. 0160 % pkowal3@sgh.waw.pl 0161 0162 [n,m] = size(A); 0163 0164 if ~issparse(A) 0165 A = sparse(A); 0166 end 0167 0168 %-------------------------------------------------------------------------- 0169 % SPECIAL CASES 0170 %-------------------------------------------------------------------------- 0171 if size(A,1)==0 0172 L = speye(n); 0173 U = A; 0174 Q = speye(m); 0175 return; 0176 end 0177 if size(A,2)==0 0178 L = speye(n); 0179 U = A; 0180 Q = speye(m); 0181 return; 0182 end 0183 0184 %-------------------------------------------------------------------------- 0185 % LU DECOMPOSITION 0186 %-------------------------------------------------------------------------- 0187 if do_pivot 0188 [L,U,P,Q] = lu(A); 0189 Q = Q'; 0190 else 0191 [L,U,P] = lu(A); 0192 Q = speye(m); 0193 end 0194 p = size(A,1)-size(L,2); 0195 LL = [sparse(n-p,p);speye(p)]; 0196 L = [P'*L P(n-p+1:n,:)']; 0197 U = [U;sparse(p,m)]; 0198 0199 %-------------------------------------------------------------------------- 0200 % FINDS ROWS WITH ZERO AND NONZERO ELEMENTS ON THE DIAGONAL 0201 %-------------------------------------------------------------------------- 0202 if size(U,1)==1 || size(U,2)==1 0203 S = U(1,1); 0204 else 0205 S = diag(U); 0206 end 0207 I = find(abs(S)>tol); 0208 Jl = (1:n)'; 0209 Jl(I) = []; 0210 Jq = (1:m)'; 0211 Jq(I) = []; 0212 0213 Ubar1 = U(I,I); 0214 Ubar2 = U(Jl,Jq); 0215 Qbar1 = Q(I,:); 0216 Lbar1 = L(:,I); 0217 0218 %-------------------------------------------------------------------------- 0219 % ELININATES NONZEZO ELEMENTS BELOW AND ON THE RIGHT OF THE 0220 % INVERTIBLE BLOCK OF THE MATRIX U 0221 % 0222 % UPDATES MATRICES L, Q 0223 %-------------------------------------------------------------------------- 0224 if ~isempty(I) 0225 Utmp = U(I,Jq); 0226 X = Ubar1'\U(Jl,I)'; 0227 Ubar2 = Ubar2-X'*Utmp; 0228 Lbar1 = Lbar1+L(:,Jl)*X'; 0229 0230 X = Ubar1\Utmp; 0231 Qbar1 = Qbar1+X*Q(Jq,:); 0232 Utmp = []; 0233 X = []; 0234 end 0235 0236 %-------------------------------------------------------------------------- 0237 % FINDS ROWS AND COLUMNS WITH ONLY ZERO ELEMENTS 0238 %-------------------------------------------------------------------------- 0239 I2 = find(max(abs(Ubar2),[],2)>tol); 0240 I5 = find(max(abs(Ubar2),[],1)>tol); 0241 0242 I3 = Jl(I2); 0243 I4 = Jq(I5); 0244 Jq(I5) = []; 0245 Jl(I2) = []; 0246 U = []; 0247 0248 %-------------------------------------------------------------------------- 0249 % FINDS A PART OF THE MATRIX U WHICH IS NOT IN THE REQIRED FORM 0250 %-------------------------------------------------------------------------- 0251 A = Ubar2(I2,I5); 0252 0253 %-------------------------------------------------------------------------- 0254 % PERFORMS LUQ DECOMPOSITION OF THE MATRIX A 0255 %-------------------------------------------------------------------------- 0256 [L1,U1,Q1] = luq(A,do_pivot,tol); 0257 0258 %-------------------------------------------------------------------------- 0259 % UPDATES MATRICES L, U, Q 0260 %-------------------------------------------------------------------------- 0261 Lbar2 = L(:,I3)*L1; 0262 Qbar2 = Q1*Q(I4,:); 0263 L = [Lbar1 Lbar2 L(:,Jl)]; 0264 Q = [Qbar1; Qbar2; Q(Jq,:)]; 0265 0266 n1 = length(I); 0267 n2 = length(I3); 0268 m2 = length(I4); 0269 U = [Ubar1 sparse(n1,m-n1);sparse(n2,n1) U1 sparse(n2,m-n1-m2);sparse(n-n1-n2,m)];