1function [U, S, V, stats] = spqr_ssi (R, varargin) 2%SPQR_SSI block power method or subspace iteration applied to inv(R) 3% to estimate rank, smallest singular values, and left/right singular vectors. 4% 5% [U,S,V,stats] = spqr_ssi (R,opts); 6% 7% Uses the block power method or subspace iteration applied to the inverse of 8% the triangular matrix R (implicitly) to estimate the numerical rank, smallest 9% singular values, and left and right singular vectors R. The algorithm will 10% be efficient in the case the nullity of R is relatively low and inefficient 11% otherwise. 12% 13% Input: 14% 15% R -- an n by n nonsingular triangular matrix 16% opts (optional) -- see 'help spqr_rank_opts' for details. 17% 18% Output: 19% 20% U -- n by k matrix containing estimates of the left singular vectors of R 21% corresponding to the singular values in S. When stats.flag is 0 (and 22% when the optional parameter opts.nsvals_large has its default value 1), 23% r = n - k + 1 is the estimated rank of R. Also NT=U(:,2:k) is an 24% orthonormal basis for the numerical null space of R'. 25% S -- A k by k diagonal matrix whose diagonal entries are the estimated 26% smallest k singular values of R, with k as described above. For i=1:k, 27% S(i,i) is an an estimate of singular value (r + i - 1) of R. 28% Note that, unless stats.flag = 3 (see below), S(1,1)) > tol and for 29% i =2:k, S(i,i) <= tol. 30% V -- n by k matrix containing estimates of the right singular vectors of R 31% corresponding to the singular values in S. When stats.flag is 0, 32% N=V(:,2:k) is an orthonormal basis for the numerical null space of R. 33% stats -- statistics, type 'help spqr_rank_stats' for details. 34% 35% Note that U' * R = S * V' (in exact arithmetic) and that R * V is 36% approximately equal to U * S. 37% 38% output (for one or two output arguments): [s,stats] = spqr_ssi (...) 39% s -- diag(S), with S as above. 40% 41% Example: 42% R = sparse (gallery ('kahan',100)) ; 43% [U,S,V] = spqr_ssi (R) ; 44% norm_of_residual = norm( U' * R - S * V' ) % should be near zero 45% % or 46% [U,S,V,stats] = spqr_ssi (R) ; 47% N = V(:,2:end); % orthonormal basis for numerical null space of R 48% NT = U(:,2:end); % orthonormal basis for numerical null space of R' 49% norm_R_times_N = norm(R*N) 50% norm_R_transpose_times_NT = norm(R'*NT) 51% % or 52% opts = struct('tol',1.e-5,'nsvals_large',3) ; 53% [s,stats] = spqr_ssi (R,opts); 54% stats % information about several singular values 55% 56% See also spqr_basic, spqr_null, spqr_pinv, spqr_cod. 57 58% Copyright 2012, Leslie Foster and Timothy A Davis. 59 60% Outline of algorithm: 61% let b = initial block size 62% U = n by b random matrix with orthonormal columns 63% repeat 64% V1 = R \ U 65% [V,D1,X1]= svd(V1,0) = the compact svd of V1 66% U1 = R' \ V 67% [U,D2,X2] = svd(U1,0) = the compact svd of U1 68% sdot = diag(D2) 69% s = 1 ./ sdot = estimates of singular values of R 70% if s(end) <= tol 71% increase the number of columns of U and repeat loop 72% else 73% repeat loop until stopping criteria (see code) is met 74% end 75% end 76% k = smallest i with sdot(i) > tol 77% V = V*X2 (see code) 78% V = first k columns of V in reverse order 79% U = first k columns of U in reverse order 80% s = 1 ./ (first k entries in sdot in reverse order) 81% S = diag(s) 82 83% Note for i = 1, 2, ..., k, S(i,i) is an upper bound on singular 84% value (r + i - 1) of R in exact arithmetic but not necessarily 85% in finite precision arithmetic. 86 87% start block power method to estimate the smallest singular 88% values and the corresponding right singular vectors (in the 89% r by nblock matrix V) and left singular vectors (in the r by 90% nblock matrix U) of the n by n matrix R 91 92% disable nearly-singular matrix warnings, and save the previous state 93warning_state = warning ('off', 'MATLAB:nearlySingularMatrix') ; 94 95%------------------------------------------------------------------------------- 96% get input options 97%------------------------------------------------------------------------------- 98 99[ignore,opts,stats,start_tic,ok] = spqr_rank_get_inputs (R,0,varargin{:}) ; %#ok 100clear ignore 101 102if (~ok || nargout > 4) 103 error ('usage: [U,S,V,stats] = spqr_ssi (R,opts)') ; 104end 105 106[m,n] = size (R) ; 107if m ~= n 108 error('R must be square') 109end 110 111% see spqr_rank_get_inputs.m for defaults: 112tol = opts.ssi_tol ; % default opts.tol 113min_block = opts.ssi_min_block ; % default 3 114max_block = opts.ssi_max_block ; % default 10 115min_iters = opts.ssi_min_iters ; % default 3 116max_iters = opts.ssi_max_iters ; % default 100 117nblock_increment = opts.ssi_nblock_increment ; % default 5 118convergence_factor = opts.ssi_convergence_factor ; % default 0.1 119nsvals_large = opts.nsvals_large ; % default 1 120get_details = opts.get_details ; % default false 121repeatable = opts.repeatable ; % default false 122 123private_stream = spqr_repeatable (repeatable) ; 124 125%------------------------------------------------------------------------------- 126% adjust the options to reasonable ranges 127%------------------------------------------------------------------------------- 128 129% cannot compute more than n large singular values, where R is n-by-n 130nsvals_large = min (nsvals_large, n) ; 131 132% make the block size large enough so that there is a possiblity of 133% calculating all nsvals_large singular values 134max_block = max (max_block, nsvals_large) ; 135 136% max_block cannot be larger than n 137max_block = min (max_block, n) ; 138 139% min_block cannot be larger than n 140min_block = min (min_block, n) ; 141 142% start with nblock = min_block; 143nblock = min_block ; 144 145%------------------------------------------------------------------------------- 146% initializations 147%------------------------------------------------------------------------------- 148 149% stats.flag and stats.rank initialized in spqr_rank_get_inputs 150% set the order of the remaining stats fields 151stats.tol = tol ; 152stats.tol_alt = -1 ; % removed later if remains -1 153if get_details == 1 && isfield(stats,'normest_A') 154 normest_A = stats.normest_A ; 155 stats = rmfield( stats,'normest_A' ) ; % to place in proper order 156 stats.normest_R = normest_A ; 157end 158stats.est_svals_of_R = 0 ; 159stats.est_error_bounds = -1 ; 160stats.sval_numbers_for_bounds = -1 ; 161if get_details == 2 162 stats.ssi_max_block_used = -1 ; 163 stats.ssi_min_block_used = -1 ; 164end 165if (get_details == 1) 166 stats.norm_R_times_N = -1 ; 167 stats.norm_R_transpose_times_NT = -1 ; 168 stats.iters = -1 ; 169 stats.nsvals_large_found = 0 ; 170 stats.final_blocksize = -1 ; 171 stats.ssi_max_block_used = -1 ; 172 stats.ssi_min_block_used = -1 ; 173 stats.opts_used = opts ; 174 stats.time = 0 ; 175 time_initialize = stats.time_initialize ; 176 stats = rmfield(stats,'time_initialize') ; % to place in proper order 177 stats.time_initialize = time_initialize ; 178 stats.time_iters = 0 ; 179 stats.time_est_error_bounds = 0 ; 180 stats.time_svd = 0 ; 181end 182 183 184stats.rank = [ ] ; 185stats.tol = tol ; 186if (get_details == 1 || get_details == 2) 187 stats.ssi_max_block_used = max_block ; 188 stats.ssi_min_block_used = nblock ; 189end 190 191if (get_details == 1) 192 start_iters_tic = tic ; 193end 194 195if (~isempty (private_stream)) 196 U = randn (private_stream, n, nblock) ; 197else 198 U = randn (n, nblock) ; 199end 200[U,ignore] = qr (U,0) ; %#ok 201clear ignore 202% est_error_bound_calculated = 0 ; % set to 1 later if bound calculated 203flag_overflow = 0 ; % set to 1 later if overflow occurs 204 205%------------------------------------------------------------------------------- 206% iterations 207%------------------------------------------------------------------------------- 208 209for iters = 1:max_iters 210 211 U0= U; 212 V1 = R \ U ; 213 if min(min(isfinite(V1))) == 0 214 flag_overflow = 1; 215 break % *************>>>>> early exit from for loop 216 end 217 218 if (get_details == 1), time_svd = tic ; end 219 220 [V,D1,X1] = svd (V1,0) ; 221 222 if (get_details == 1) 223 stats.time_svd = stats.time_svd + toc( time_svd ) ; 224 end 225 226 % d1 = diag(D1)'; 227 U1 = ( V' / R )'; % uses less memory than U1 = R' \ V; 228 % Note: with the inverse power method overflow is a potential concern 229 % in extreme cases (SJid = 137 or UFid = 1239 is close) 230 if min (min (isfinite (U1))) == 0 231 % *************>>>>> early exit from for loop 232 % We know of no matrix that triggers this condition, so the next 233 % two lines are untested. 234 flag_overflow = 1; % untested 235 break % untested 236 end 237 238 if (get_details == 1), time_svd = tic ; end 239 240 [U,D2,X2] = svd (U1,0) ; 241 242 if (get_details == 1) 243 stats.time_svd = stats.time_svd + toc(time_svd) ; 244 end 245 246 d2 = diag (D2)' ; 247 k = find(d2*tol<1); 248 249 if isempty(k) && (nblock == n) 250 % success since the numerical rank of R is zero 251 break % *************>>>>> early exit from for loop 252 end 253 254 if ~isempty(k) 255 k = k(1); % k equals the calculated nullity + 1, k corresponds to 256 % the first singular value of R that is greater than tol 257 % reduce nsvals_large if necessary due to max_block limit: 258 nsvals_large_old = nsvals_large; 259 nsvals_large = min(nsvals_large,max_block - k + 1); 260 261 if nblock >= k + nsvals_large - 1 262 % estimate error bound for singular value n - k + 1 of R: 263 est_error_bound = norm( ... 264 U0*(X1*(D1 \ X2(:,k))) - U(:,k)/D2(k,k) ) / sqrt(2); 265 % est_error_bound_calculated = 1; 266 k2 = k + nsvals_large - 1; 267 % When nsvals_large is 1, k and k2 are the same. When 268 % nsvals_large > 1, k2 corresponds to the largest singular 269 % value of R that is returned. 270 if k2 ~= k 271 % estimate error bound for singular value n - k2 + 1 of R: 272 est_error_bound2 = norm( ... 273 U0*(X1*(D1 \ X2(:,k2))) - U(:,k2)/D2(k2,k2) ) / sqrt(2); 274 else 275 est_error_bound2 = est_error_bound; 276 end 277 end 278 279 % Note that 280 % [ 0 R ] [ U ] = [ U0 * X1 * ( D1 \ X2 ) ] 281 % [ R' 0 ] [ V*X2 ] [ V * ( X2 / D2 ) ] 282 % It follows that, by Demmel's Applied Numerical Linear Algebra, 283 % Theorem 5.5, some eigenvalue of 284 % B = [ 0 R ] 285 % [ R' 0 ] 286 % will be within a distance of est_error_bound of 1 / d2(k), where s(1) 287 % = 1 / d2(k) is our estimate of singular value n - k + 1 of R. 288 % Typically the eigenvalue of B is singular value number n - k + 1 of 289 % R, although this is not guaranteed. est_error_bound is our estimate 290 % of a bound on the error in using s(1) to approximate singular value 291 % number n - k + 1 of R. 292 293 if nblock >= k + nsvals_large - 1 && ... 294 est_error_bound <= convergence_factor*abs( 1/d2(k)-tol ) && ... 295 est_error_bound2 <= convergence_factor*abs( 1/d2(k2) ) && ... 296 iters >= min_iters 297 % Motivation for the tests: 298 % The first test in the if statement is an attempt to insure 299 % that nsvals_large singular values of R larger than tol are 300 % calculated. 301 % 302 % Goal of the power method is to increase nblock until we find 303 % sigma = (singular value n - k + 1 of R) is > tol. If 304 % this is true it is guaranteed that n - k + 1 is the 305 % correct numerical rank of R. If we let s(1) = 1 / d2(k) 306 % then s(1) > tol. However s(1) is only an estimate of sigma. 307 % However, in most cases 308 % | s(1) - sigma | <= est_error_bound (1) 309 % and to be conservative assume only that 310 % |s(1) - sigma|<= (1/convergence_factor)*est_error_bound (2) 311 % where convergence_factor<=1. By the second test in the if 312 % statement 313 % est_error_bound <= convergence_factor * | s(1) - tol | (3) 314 % Equations (2) and (3) imply that 315 % | s(1) - sigma | <= | s(1) - tol | 316 % This result and s(1) > tol imply that sigma > tol, as 317 % desired. Thus the second test in the if statement attempts 318 % to carry out enough iterations to insure that the calculated 319 % numerical rank is correct. 320 % 321 % The third test in the if statement checks on the accuracy of 322 % the estimate for singular values n - k2 + 1. Let sigma2 be 323 % singular value n - k2 + 1 of R. Usually it is true 324 % that 325 % | s( k2 ) - sigma2 | <= est_error_bound2. (4) 326 % Assuming (4) and the third test in the if statement it 327 % follows that the estimated relative 328 % error in s(k2), as measured by est_error_bound2 / s( k2) , 329 % is less that or equal to convergence_factor. Therefore 330 % the third test in the if statement attempts to insure that 331 % the largest singular value returned by ssi has a relative 332 % error <= convergence_factor. 333 % 334 % SUCCESS!!!, found singular value or R larger than tol 335 break % *************>>>>> early exit from for loop 336 end 337 nsvals_large = nsvals_large_old; % restore original value 338 end 339 340 if nblock == max_block && iters >= min_iters && isempty(k) 341 % reached max_block block size without encountering any 342 % singular values of R larger than tolerance 343 break % *************>>>>> early exit from for loop 344 end 345 346 if (1 <= iters && iters < max_iters) && ( isempty(k) || ... 347 ( ~isempty(k) && nblock < k(1) + nsvals_large - 1) ) 348 % increase block size 349 nblock_prev = nblock; 350 if isempty(k) 351 nblock = min(nblock + nblock_increment, max_block); 352 else 353 nblock = min( k(1) + nsvals_large - 1, max_block ); 354 end 355 if (nblock > nblock_prev) 356 if (~isempty (private_stream)) 357 Y = randn (private_stream, n, nblock-nblock_prev) ; 358 else 359 Y = randn (n, nblock-nblock_prev) ; 360 end 361 Y = Y - U*(U'*Y); 362 [Y,ignore]=qr(Y,0) ; %#ok 363 clear ignore 364 U = [U, Y]; %#ok 365 end 366 end 367end 368 369if (get_details == 1) 370 stats.final_blocksize = nblock ; % final block size 371 stats.iters = iters ; % number of iterations taken in ssi 372 stats.time_iters = toc (start_iters_tic) ; 373end 374 375%------------------------------------------------------------------------------- 376% check for early return 377%------------------------------------------------------------------------------- 378 379if flag_overflow == 1 380 warning ('spqr_rank:overflow', 'overflow in inverse power method') ; 381 warning (warning_state) ; 382 [stats U S V] = spqr_failure (4, stats, get_details, start_tic) ; 383 if (nargout == 2) 384 S = stats ; 385 end 386 return 387end 388 389%------------------------------------------------------------------------------- 390% determine estimated singular values of R 391%------------------------------------------------------------------------------- 392 393est_error_bounds = [ ] ; 394 395if ~isempty(k) 396 % Note: in this case nullity = k - 1 and rank = n - k + 1 397 s = 1.0 ./ d2(min(nblock,k+nsvals_large-1):-1:1); 398 est_error_bounds (nsvals_large) = est_error_bound; 399 numerical_rank = n - k + 1; 400else 401 k = nblock; 402 if nblock == n 403 numerical_rank = 0; 404 % success since numerical rank is 0 405 else 406 % in this case rank not successfully determined 407 % Note: In this case k is a lower bound on the nullity and 408 % n - k is an upper bound on the rank 409 numerical_rank = n - k; %upper bound on numerical rank 410 nsvals_large = 0 ; % calculated no singular values > tol 411 end 412 s = 1.0 ./ d2(nblock:-1:1); 413end 414 415stats.rank = numerical_rank ; 416 417nkeep = length (s) ; % number of cols kept in U and V 418S = diag (s) ; 419 420if (get_details == 1) 421 stats.nsvals_large_found = nsvals_large ; 422end 423 424 425%------------------------------------------------------------------------------- 426% adjust V so that R'*U = V*S (in exact arithmetic) 427%------------------------------------------------------------------------------- 428 429V = V*X2; 430% reverse order of U and V and keep only nkeep singular vectors 431V = V(:,nkeep:-1:1); 432U = U(:,nkeep:-1:1); 433 434if (get_details == 1) 435 t = tic ; 436end 437 438if nsvals_large > 0 439 % this recalculates est_error_bounds(nsvals_large) 440 U0 = R * V(:,1:nsvals_large) - ... 441 U(:,1:nsvals_large) * S(1:nsvals_large,1:nsvals_large) ; 442 U0 = [ U0; 443 R' * U(:,1:nsvals_large) - V(:,1:nsvals_large)* ... 444 S(1:nsvals_large,1:nsvals_large)]; 445 est_error_bounds (1:nsvals_large) = sqrt(sum(U0 .* conj(U0) )) / sqrt(2) ; 446end 447 448% this code calculates estimated error bounds for singular values 449% nsvals_large+1 to nkeep 450ibegin = nsvals_large+1; 451U0 = R * V(:,ibegin:nkeep) - U(:,ibegin:nkeep)* ... 452 S(ibegin:nkeep,ibegin:nkeep); 453U0 = [ U0; 454 R' * U(:,ibegin:nkeep) - V(:,ibegin:nkeep)* ... 455 S(ibegin:nkeep,ibegin:nkeep)]; 456est_error_bounds (ibegin:nkeep) = sqrt(sum(U0 .* conj(U0) )) /sqrt(2); 457% Note that 458% [ 0 R ] [ U ] - [ U ] * S = [ R * V - U * S ] 459% [ R' 0 ] [ V ] - [ V ] [ 0 ]. 460% It follows, by Demmel's Applied Numerical Linear Algebra, 461% Theorem 5.5, for i = 2, . . ., k, some eigenvalue of 462% B = [ 0 R ] 463% [ R' 0 ] 464% will be within a distance of the norm of the ith column of 465% [ R * V - U * S ; R'*U - V*S] / sqrt(2) from S(i,i). Typically this 466% eigenvalue will be singular value number n - k + i of R, 467% although it is not guaranteed that the singular value number 468% is correct. est_error_bounds(i) is the norm of the ith column 469% of [ R * V - U * S; R' * U - V * S ] / sqrt(2). 470 471 472if (get_details == 1) 473 % Note that stats.time_est_error_bounds includes the time for the error 474 % bound calculations done outside of the subspace iterations loop 475 stats.time_est_error_bounds = toc (t) ; 476end 477 478stats.sval_numbers_for_bounds = ... 479 numerical_rank - nsvals_large + 1 : ... 480 numerical_rank - nsvals_large + length (est_error_bounds) ; 481 482%------------------------------------------------------------------------------- 483% compute norm R*N and R'*NT 484%------------------------------------------------------------------------------- 485 486% compute norm R*N where N = V(:,nsvals_large+1:end) is the approximate 487% null space of R and 488% norm R'*NT where NT = U(:,nsvals_large+1:end) is the approximate 489% null space of R' 490 491if (get_details == 1), t = tic ; end 492 493% svals_R_times_N = svd(R*V(:,nsvals_large+1:end))'; 494norm_R_times_N = norm (R * V(:,nsvals_large+1:end)) ; 495% svals_R_transpose_times_NT = svd(R'*U(:,nsvals_large+1:end))'; 496norm_R_transpose_times_NT = norm (R' * U(:,nsvals_large+1:end)) ; 497 498if (get_details == 1) 499 stats.norm_R_times_N = norm_R_times_N; 500 stats.norm_R_transpose_times_NT = norm_R_transpose_times_NT; 501 stats.time_svd = stats.time_svd + toc (t) ; 502end 503 504% Note: norm_R_times_N is an upper bound on sing. val. rank1+1 of R 505% and norm_R_transpose_times_NT is also an upper bound on sing. val. 506% rank1+1 of R 507max_norm_RN_RTNT = max (norm_R_times_N,norm_R_transpose_times_NT); 508% choose max here to insure that both null spaces are good 509 510%------------------------------------------------------------------------------- 511% determine flag indicating the accuracy of the rank calculation 512%------------------------------------------------------------------------------- 513 514if numerical_rank == 0 515 % numerical rank is 0 in this case 516 stats.flag = 0; 517elseif (numerical_rank == n || max_norm_RN_RTNT <= tol) && ... 518 (nsvals_large > 0 && ... 519 s(nsvals_large) - est_error_bounds(nsvals_large) > tol ) 520 % in this case, assuming est_error_bounds(nsvals_large) is a true 521 % error bound, then the numerical rank is correct. Also 522 % N = V(:,nsvals_large+1:end) and NT = U(:,nsvals_large+1:end) 523 % are bases for the numerical null space or R and R', respectively 524 stats.flag = 0; 525elseif ( nsvals_large > 0 && numerical_rank == n ) || ... 526 ( nsvals_large > 0 && ... 527 s(nsvals_large) - est_error_bounds(nsvals_large) > ... 528 max_norm_RN_RTNT ) 529 % in this case, assuming est_error_bounds(nsvals_large) is a true 530 % error bound, then the numerical rank is correct with a modified 531 % tolerance. This is a rare case. 532 stats.flag = 1; 533 tol_alt = ( s(nsvals_large) - est_error_bounds(nsvals_large) ); 534 tol_alt = tol_alt - eps(tol_alt); % so that tol_alt satisfies the > 535 % portion of the inequality below 536 % tol_alt = max_norm_RN_RTNT; 537 % Note: satisfactory values of tol_alt are in the range 538 % s(nsvals_large) - est_error_bounds(nsvals_large) > tol_alt 539 % >= max_norm_RN_RTNT 540 stats.tol_alt = tol_alt; 541elseif nsvals_large > 0 && s(nsvals_large) > tol && ... 542 max_norm_RN_RTNT <= tol && ... 543 s(nsvals_large) - est_error_bounds(nsvals_large) ... 544 <= max_norm_RN_RTNT 545 % in this case, assuming est_error_bounds(nsvals_large) is a true 546 % error bound, the error bound is too large to determine the 547 % numerical rank. This case is very rare. 548 stats.flag = 2; 549else 550 % in this case all the calculated singular values are 551 % smaller than tol or either N is not a basis for the numerical 552 % null space of R with tolerance tol or NT is not such a basis for R'. 553 % stats.rank is an upper bound on the numerical rank. 554 stats.flag = 3; 555end 556 557%------------------------------------------------------------------------------- 558% return results 559%------------------------------------------------------------------------------- 560 561% restore the warnings 562warning (warning_state) ; 563 564stats.est_svals_of_R = diag (S)' ; 565stats.est_error_bounds = est_error_bounds; 566 567if (get_details == 1) 568 stats.time = toc (start_tic) ; 569end 570 571if stats.tol_alt == -1 572 stats = rmfield(stats, 'tol_alt') ; 573end 574 575if (get_details == 1) 576 stats.time = toc (start_tic) ; 577end 578 579if (nargout <= 2) 580 U = diag (S) ; 581 S = stats ; 582end 583 584