sparseNull

PURPOSE ^

sparseNull returns computes the sparse Null basis of a matrix

SYNOPSIS ^

function N = sparseNull(S, tol)

DESCRIPTION ^

 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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

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

Generated on Thu 21-Jun-2012 15:39:23 by m2html © 2003