1function [U, S, V, stats] = spqr_ssp (A, varargin) 2%SPQR_SSP block power method or subspace iteration applied to A or A*N 3% 4% [U,S,V,stats] = spqr_ssp (A, N, k, opts) ; 5% 6% Uses the block power method or subspace iteration to estimate the largest k 7% singular values, and left and right singular vectors an m-by-n matrix A or, 8% if the n-by-p matrix N is included in the input, of A*N. 9% 10% Alternate input syntaxes (opts are set to defaults if not present): 11% 12% spqr_ssp (A) k = 1, N is [ ] 13% spqr_ssp (A,k) k is a scalar 14% spqr_ssp (A,N) N is a matrix or struct 15% spqr_ssp (A,N,k) k is a scalar 16% spqr_ssp (A,N,opts) opts is a struct 17% 18% Input: 19% A -- an m-by-n matrix 20% N (optional) -- an n-by-p matrix representing the null space of A, stored 21% either explicitly as a matrix or implicitly as a structure, and 22% returned by spqr_basic, spqr_null, spqr_cod, or spqr_pinv. 23% N is ignored if [ ]. 24% k (optional) -- the number of singular values to calculate. Also the block 25% size in the block power method. Also can be specified as opts.k. 26% Default: 1. 27% opts (optional) -- type 'help spqr_rank_opts' for details 28% 29% Output: 30% U -- m-by-k matrix containing estimates of the left singular vectors of A 31% (or A*N if N is included in the input) corresponding to the singular 32% values in S. 33% S -- k-by-k diagonal matrix whose diagonal entries are the estimated 34% singular values of A (or A*N). Also for i = 1:nsval, S(i,i) is a lower 35% bound on singular value i of A (or A*N). 36% V -- n-by-k matrix containing estimates of the right singular vectors of A 37% corresponding to the singular values in S. 38% stats -- type 'help spqr_rank_stats' 39% 40% Note that A*V = U*S (or A*N*V = U*S) and that U'*A (or U'*A*N) 41% is approximately equal to S*V'. 42% 43% Output (for one or two output arguments): [s,stats]= spqr_ssp (...) 44% s -- diag(S), with S as above 45% 46% Example: 47% A = sparse (gallery ('kahan',100)) ; 48% [U,S,V] = spqr_ssp (A) ; % singular triple for largest singular value 49% % or 50% [U,S,V,stats] = spqr_ssp (A,4) ; % same for largest 4 singular values 51% % or 52% [s,stats] = spqr_ssp (A,4) ; 53% stats % stats has information for largest 4 singular values 54% % or 55% N = spqr_null ( A ) ; % N has a basis for the null space of A 56% s = spqr_ssp (A,N) ; 57% s % s has estimate of largest singular value (or norm) of A*N 58% 59% See also spqr_basic, spqr_null, spqr_pinv, spqr_cod. 60 61% Copyright 2012, Leslie Foster and Timothy A Davis. 62 63% Outline of algorithm: 64% Let B = A or B = A*N, if the optional input N is included 65% where A is m-by-n and N is n-by-p. Note that the code is written so 66% that B is not explicitly formed. Let n = p if N is not input. 67% Let k = block size = min(k,m,p) 68% U = m-by-k random matrix with orthonormal columns 69% repeat 70% V1 = B' * U 71% [V,D1,X1]= svd(V1,0) = the compact svd of V1 72% U1 = B * V 73% [U,D2,X2] = svd(U1,0) = the compact svd of U1 74% s = diag(D2) 75% kth_est_error_bound = norm(B'*U(:,k) - V*(X2(:,k)*D2(k,k))) / sqrt(2); 76% if kth_est_error_bound <= convergence_factor*s(k) && iters >= min_iters 77% flag = 0 78% exit loop 79% elseif number of iteration > max number of iterations 80% flag = 1 81% exit loop 82% else 83% continue loop 84% end 85% end 86% V = V*X2; % adjust V so that A*V = U*S 87% calculate estimated error bounds for singular values 1:k-1 88% as discussed in the code 89% S = diag(s); 90 91% start block power method to estimate the largest singular values and the 92% corresponding right singular vectors (in the n-by-k matrix V) and left 93% singular vectors (in the m-by-k matrix U) of the m-by-n matrix A 94 95%------------------------------------------------------------------------------- 96% get input options 97%------------------------------------------------------------------------------- 98[N,opts,stats,start_tic,ok] = spqr_rank_get_inputs (A, 2, varargin {:}) ; 99 100if (~ok || nargout > 4) 101 error ('usage: [U,S,V,stats] = spqr_ssp (A,N,k,opts)') ; 102end 103 104k = max (opts.k, 0) ; 105min_iters = opts.ssp_min_iters ; 106max_iters = opts.ssp_max_iters ; 107convergence_factor = opts.ssp_convergence_factor ; 108get_details = opts.get_details ; 109repeatable = opts.repeatable ; 110 111%------------------------------------------------------------------------------- 112% initializations 113%------------------------------------------------------------------------------- 114 115% the order of all the fields of stats are set in spqr_rank_get_inputs 116 117[m,n] = size (A) ; 118 119if (isempty (N)) 120 % do not use N 121 use_N = 0 ; 122 p = n ; 123elseif (isstruct (N)) 124 % N is a struct 125 use_N = 2 ; 126 p = size (N.X, 2) ; 127else 128 % N is a matrix 129 use_N = 1 ; 130 p = size (N, 2) ; 131end 132 133k = min ([k m p]) ; % number of singular values to compute 134 135stats.flag = 1 ; 136stats.est_error_bounds = zeros (1, max (k,1)) ; 137stats.sval_numbers_for_bounds = 1:k ; 138 139if (k <= 0) 140 % quick return. This is not an error condition. 141 stats.flag = 0 ; 142 stats.est_svals = 0 ; 143 stats.est_error_bounds = 0 ; 144 stats.sval_numbers_for_bounds = 0 ; 145 U = zeros (m, 0) ; 146 S = 0 ; 147 V = zeros (n, 0) ; 148 if (nargout == 1) 149 % usage: s = spqr_ssp ( ... ) 150 S = 0 ; 151 elseif (nargout == 2) 152 % usage: [s,stats] = spqr_ssp ( ... ) 153 U = 0 ; 154 S = stats ; 155 end 156 return 157end 158 159private_stream = spqr_repeatable (repeatable) ; 160 161if (~isempty (private_stream)) 162 U = randn (private_stream, m, k) ; 163else 164 U = randn (m, k) ; 165end 166 167[U, ignore] = qr (U, 0) ; %#ok 168clear ignore 169 170%------------------------------------------------------------------------------- 171% block power iterations 172%------------------------------------------------------------------------------- 173 174if (get_details == 1) 175 t = tic ; 176end 177 178for iters = 1:max_iters 179 180 if use_N == 0 181 V1 = A' * U ; 182 elseif use_N == 1 183 V1 = N' * (A' * U) ; 184 else 185 V1 = spqr_null_mult (N, (A' * U), 0) ; 186 end 187 188 if issparse (V1) 189 V1 = full (V1) ; % needed if U is 1-by-1 190 end 191 192 if (get_details == 1), time_svd = tic ; end 193 % [V,D1,X1] = svd (V1, 0) ; 194 [V,ignore1,ignore2] = svd (V1, 0) ; %#ok 195 clear ignore1 ignore2 196 if (get_details == 1) 197 stats.time_svd = stats.time_svd + toc(time_svd) ; 198 end 199 200 % d1 = diag(D1)'; 201 202 if use_N == 0 203 U1 = A * V; 204 elseif use_N == 1 205 U1 = A * (N * V); 206 else 207 U1 = A * spqr_null_mult (N,V,1) ; 208 end 209 210 if issparse(U1) 211 U1 = full (U1) ; % needed if V is 1-by-1 212 end 213 214 if (get_details == 1), time_svd = tic ; end 215 [U,D2,X2] = svd (U1, 0) ; 216 if (get_details == 1) 217 stats.time_svd = stats.time_svd + toc(time_svd) ; 218 end 219 220 d2 = diag (D2)' ; 221 222 % estimate error bound on kth singular value 223 if use_N == 0 224 V1 = A' * U(:,k) ; 225 elseif use_N == 1 226 V1 = N' * (A' * U(:,k)) ; 227 else 228 V1 = spqr_null_mult(N,(A' * U(:,k)),0); 229 end 230 kth_est_error_bound = norm( V1 - V *( X2(:,k)*D2(k,k) ) ) / sqrt(2); 231 232 % Note that 233 % [ 0 B ] [ U ] = [ U * D2 ] 234 % [ B' 0 ] [ V*X2 ] [ V * ( X2 * D2 ) ] 235 % where B = A or B = A*N (if N is included in the input). It follows that, 236 % by Demmel's Applied Numerical Linear Algebra, Theorem 5.5, some 237 % eigenvalue of 238 % C = [ 0 B ] 239 % [ B' 0 ] 240 % will be within a distance of kth_est_error_bound of D2(k,k). Typically 241 % the eigenvalue of C is the kth singular value of B, although this is not 242 % guaranteed. kth_est_error_bound is our estimate of a bound on the error 243 % in using D2(k,k) to approximate kth singular value of B. 244 245 if kth_est_error_bound <= convergence_factor* d2(k) && iters >= min_iters 246 % The test indicates that the estimated relative error in 247 % the estimate, D2(k,k) for the kth singular value is smaller than 248 % convergence_factor, which is the goal. 249 stats.flag = 0 ; 250 break 251 end 252end 253 254if (get_details == 1) 255 stats.time_iters = toc (t) ; 256end 257 258s = d2 (1:k); 259stats.est_error_bounds (k) = kth_est_error_bound ; 260S = diag (s) ; 261stats.est_svals = diag (S)' ; 262 263 264% adjust V so that A*V = U*S 265V = V*X2 ; 266 267%------------------------------------------------------------------------------- 268% estimate error bounds for the 1st through (k-1)st singular values 269%------------------------------------------------------------------------------- 270 271if (get_details == 1) 272 t = tic ; 273end 274 275if k > 1 276 277 if use_N == 0 278 V1 = A' * U(:,1:k-1) ; 279 elseif use_N == 1 280 V1 = N'* (A' * U(:,1:k-1)) ; 281 else 282 V1 = spqr_null_mult(N,(A' * U(:,1:k-1)),0); 283 end 284 285 U0 = V1 - V(:,1:k-1)*S(1:k-1,1:k-1); 286 stats.est_error_bounds (1:k-1) = sqrt (sum (U0 .* conj(U0))) / sqrt (2) ; 287 288 % Note (with V adjusted as above) that 289 % [ 0 B ] [ U ] - [ U ] * S = [ 0 ] 290 % [ B' 0 ] [ V ] - [ V ] [ A'*U - V*S ] 291 % where B = A or B = A*N (if N is included in the input). It follows, by 292 % Demmel's Applied Numerical Linear Algebra, Theorem 5.5, for i = 1:k, some 293 % eigenvalue of 294 % C = [ 0 B ] 295 % [ B' 0 ] 296 % will be within a distance of the norm of the ith column of [ B' * U - V * 297 % S ] / sqrt(2) from S(i,i). Typically this eigenvalue will be the ith 298 % singular value number of B, although it is not guaranteed that the 299 % singular value number is correct. stats.est_error_bounds(i) is the norm 300 % of the ith column of [ B' * U - V * S ] / sqrt(2). 301 302end 303 304%------------------------------------------------------------------------------- 305% return results 306%------------------------------------------------------------------------------- 307 308if (get_details == 1) 309 stats.time_est_error_bounds = toc (t) ; 310 stats.iters = iters ; 311end 312 313if (get_details == 1) 314 stats.time = toc (start_tic) ; 315end 316 317if (nargout <= 2) 318 % usage: [s,stats] = spqr_ssp ( ... ) 319 U = diag (S) ; 320 S = stats ; 321end 322 323