1function [nfailures SJid_failures] = demo_spqr_rank (ids,args) 2%DEMO_SPQR_RANK lengthy demo for spqr_rank functions (requires SJget) 3% Usage: demo_spqr_rank(ids,args) 4% 5% This is a demonstration program for the routines spqr_basic, spqr_null, 6% spqr_pinv, spqr_cod discussed in the paper "Algorithm xxx: Reliable 7% Calculation of Numerical Rank, Null Space Bases, Pseudoinverse Solutions and 8% Basic Solutions using SuiteSparseQR" by Leslie Foster and Timothy Davis, 9% submitted ACM Transactions on Mathematical Software, 2011. Plots similar to 10% those in Figures 2 - 5 in the paper are reproduced. 11% 12% If the first argument ids is a negative scalar, than the smallest (-ids) 13% matrices from the SJ Collection are used for the tests. Otherwise, ids is 14% a list of matrix IDs to use for the tests. 15% 16% If demo_spqr_rank has a second input parameter the second parameter 17% controls options in the demonstration program. For example if the second 18% parameter is one then the routine produces plots similar to only those in 19% Figure 2 of the paper and if the second parameter is zero then no plots 20% are produced. These last two cases run more quickly than the default case. 21% The second parameter can also be a structure (for details see the comments 22% in the body of the code prior to the first executable statement). 23% 24% Examples: 25% demo_spqr_rank % test code for 100 matrices, create 4 figures 26% or 27% demo_spqr_rank(-200); % test code for 200 matrices, create 4 figures 28% or 29% demo_spqr_rank(-300,1); % test code for 300 matrices, create one figure 30% 31% See also spqr_basic, spqr_null, spqr_pinv, spqr_cod, SJget. 32 33% Copyright 2012, Leslie Foster and Timothy A Davis. 34 35% Potential run times: 36% demo_spqr_rank can require 20 seconds 37% demo_spqr_rank(-300,1) can require a couple of minutes 38% demo_spqr_rank(-300,2) can require 15 minutes 39% demo_spqr_rank(-640,1) can require 5 hours 40 41% The function creates a file save_samples_demo_spqr_rank containing many 42% of the local variable used in demo_spqr_rank. demo_spqr_rank(0) uses 43% the file to redraw plots of the last run of demo_spqr_rank. Also 44% 45% load save_samples_demo_spqr_rank 46% 47% can be used to load and examine these local variables in the main MATLAB 48% workspace. 49 50% The second parameter can be a structure with any of the following fields 51% demo_opts.figures (default 2) 52% 0 -- produce no plots 53% 1 -- produce plots similar to those in figure 2 54% 2 -- produce plots similar to those in figures 2 through 5 55% demo_opts.null_spaces (default is 1) 56% 0 -- return no null space bases (when allowed) 57% 1 -- return one null space basis (of A or of A') 58% 2 -- return null spaces bases of A and A' (when allowed) 59% demo_opts.doprint (default is 0) 60% -1 -- do not print out anything 61% 0 -- only print out a summary of the success of the programs 62% 1 -- in demo_spqr_rank print out a summary of the 63% purpose of the demo_spqr_rank, a discussion of the plots 64% produced and information about the success of the programs 65% demo_opts.start_with_A_transpose (default is 0) 66% 0 -- in spqr_cod start with a QR factorization of A 67% 1 -- in spqr_cod start with a QR factorization of A transpose 68% demo_opts.implicit_null_space_basis (default is 1) 69% 1 -- return the null space basis in implicit form as Householder 70% transformations 71% 0 -- return the null space basis as columns of an explicit matrix 72% demo_opts.repeatable (default is 1) 73% 1 -- reproduce exactly the same random numbers 74% for each run so that results are repeatable 75% 0 -- the random number stream will be different on repeated runs 76% demo_opts.nsvals (default is 1) 77% the number of small and large singular values to estimate 78% demo_opts.get_details (default is 0, but can be modified by demo_opts.figures) 79% 0 -- return basic statistics in spqr_basic, spqr_null, etc. 80% 1 -- return detailed statistics in spqr_basic, spqr_null, etc. 81% 2 -- return basic statistics and a few additional statistics needed 82% for communication between routines 83% demo_opts.tol_norm_type (default is 0) 84% 0 -- let tol = max( size(A) ) * eps( norm(A) ) where norm(A) is determined 85% using the precomputed singular values in the SJSingular Data Base 86% or (for figures = 2) by MATLAB's dense matrix SVD 87% 1 -- compute the default tolerance tol = max( size(A) ) * eps( norm (A,1) ) 88% 2 -- compute the default tolerance tol = max( size(A) ) * 89% eps( normest (A,0.01) ) 90 91if (nargin < 1) 92 ids = -100; 93end 94 95demo_opts = struct('figures', 2, ... 96 'null_spaces', 1, ... 97 'doprint', 1, ... 98 'start_with_A_transpose', 0, ... 99 'implicit_null_space_basis', 1, ... 100 'repeatable', 1 , ... 101 'nsvals', 1 , ... 102 'get_details', 0, ... 103 'tol_norm_type', 0) ; 104 105if (nargin == 2) 106 %override default values using second argument 107 if isreal(args) 108 demo_opts.figures = args ; 109 else 110 % assumes that args is a structure 111 names_args = fieldnames(args) ; 112 for i = 1 : length(names_args) 113 demo_opts.(names_args{i}) = args.(names_args{i}) ; 114 end 115 end 116end 117 118% values of get_details and null_spaces restricted by demo_opts.figures: 119if demo_opts.figures == 1 120 demo_opts.get_details = max(1, demo_opts.get_details) ; 121elseif demo_opts.figures == 2 122 demo_opts.get_details = 1 ; 123 demo_opts.null_spaces = max(1,demo_opts.null_spaces) ; 124end 125 126% SPQR_BASIC, SPQR_NULL, SPQR_PINV and SPQR_COD return null space bases stored 127% in an implicit form by default. To have the routines return null space bases 128% as explicit matrices, use args(4) = 0. To have the routines return a null 129% space bases in the form that requires less memory, use args(4) = 2. 130 131opts.get_details = demo_opts.get_details ; 132opts.repeatable = demo_opts.repeatable ; 133opts.implicit_null_space_basis = demo_opts.implicit_null_space_basis ; 134 135if (demo_opts.doprint > 0) 136 disp(' ') 137 disp(['This program demonstrates use of the routines SPQR_BASIC, '... 138 'SPQR_NULL, SPQR_PINV,']) 139 disp(['and SPQR_COD discussed in the paper "Algorithm xxx: ',... 140 'Reliable Calculation of']) 141 disp(['Numerical Rank, Null Space Bases, Pseudoinverse Solutions and ',...) 142 ' Basic Solutions']) 143 disp(['using SuiteSparseQR" by Leslie Foster and Timothy Davis, ',... 144 'submitted ACM']) 145 disp(['Transactions on Mathematical Software, 2011. Plots similar ',... 146 'to those in Figures']) 147 disp(['2 - 5 or, optionally, just Figure 2 in the paper are ',... 148 'reproduced, except the']) 149 disp(['sample set is restricted to small matrices so that the demo ',... 150 'runs quickly. The']) 151 disp(['matrices come from the San Jose State University Singular ',... 152 'Matrix Database.']) 153 disp(' ') 154 disp('The routines are designed to work with rank deficient matrices.') 155 disp('The primary use of each routine is:') 156 disp(' SPQR_BASIC -- determine a basic solution to min ||b - A x||') 157 disp([' SPQR_NULL -- determine an orthonormal basis for the ',... 158 'numerical nullspace']) 159 disp(' of A') 160 disp([' SPQR_PINV -- determine a pseudoinverse or mininimum ',... 161 'norm solutio to']) 162 disp(' min || b - A x||') 163 disp([' SPQR_COD -- determine a pseudoinverse or mininimum norm ',... 164 'solution to']) 165 disp([' min || b - A x|| using a complete orthogonal ',... 166 'decomposition.']) 167 disp('The demonstration program creates plots that illustrate the accuracy') 168 disp('of rank determination, the accuracy of the null space bases, the') 169 disp(['accuracy of the basic solutions and the accuracy of the ',... 170 'pseudoinverse']) 171 disp('solutions. The above routines are compared with the calculations') 172 disp('using MATLAB''s SVD, MATLAB''s dense matrix QR factorization and') 173 disp(['with SPQR_SOLVE, part of SuiteSparseQR. In the demonstration ',... 174 'program']) 175 disp('the tolerance defining the numerical rank is min(m,n)*eps(||A||)') 176 disp('where the matrix A is m by n.') 177 %disp(' ') 178 % if (demo_opts.dopause) 179 % disp('Press enter to begin demonstration') 180 % pause 181 % end 182 disp(' ') 183end 184 185mfilename ('fullpath') ; 186install_SJget ; 187 188% intitialze 189index = SJget; 190 191if (isscalar (ids) && ids < 0) 192 % test with matrices 1 to (-ids) 193 dim = max (index.nrows, index.ncols) ; 194 % R2009b introduced '~' to denote unused output arguments, but we avoid that 195 % feature so that this code can run on R2008a (and perhaps earlier). 196 [ignore,indexs] = sort (dim) ; %#ok 197 clear ignore 198 indexs = indexs (1:(-ids)) ; 199 200elseif (isscalar (ids) && (ids == 0)) 201 202 % the file save_samples_demo_spqr_rank is created 203 % when demo_spqr_rank is run 204 if exist('save_samples_demo_spqr_rank.mat','file') 205 load save_samples_demo_spqr_rank 206 else 207 error (['prior to running demo_spqr_rank(0) run ' ... 208 'demo_spqr_rank(ids) with ids < 0 or ids a list of IDs']) ; 209 end 210 % demo_spqr_rank(0) can be used to redraw plots of the last 211 % run of demo_spqr_rank. Also the command 212 % load save_samples_demo_spqr_rank 213 % can be used to load many of the local variables used in 214 % demo_spqr_rank into the main MATLAB workspace. 215 216else 217 % list of matrix ID's has been passed in directly 218 indexs = ids ; 219end 220 221cnt = 0; 222time_start = clock; 223 224%------------------------------------------------------------------------------- 225% allocate space for vectors containing statistics calculated 226%------------------------------------------------------------------------------- 227 228nothing = -ones (1,length (indexs)) ; 229rank_svd_v = nothing ; 230gap_v = nothing ; 231m_v = nothing ; 232n_v = nothing ; 233flag_spqr_basic_v = nothing ; 234rank_svd_basic_v = nothing ; 235rank_spqr_basic_v = nothing ; 236rank_spqr_from_spqr_basic_v = nothing ; 237flag_spqr_null_v = nothing ; 238rank_svd_null_v = nothing ; 239rank_spqr_null_v = nothing ; 240rank_spqr_from_spqr_null_v = nothing ; 241flag_spqr_pinv_v = nothing ; 242rank_svd_pinv_v = nothing ; 243rank_spqr_pinv_v = nothing ; 244rank_spqr_from_spqr_pinv_v = nothing ; 245flag_spqr_cod_v = nothing ; 246rank_svd_cod_v = nothing ; 247rank_spqr_cod_v = nothing ; 248rank_spqr_from_spqr_cod_v = nothing ; 249tol_v = nothing ; 250norm_A_v = nothing ; 251norm_A_N_svd_v = nothing ; 252norm_A_NT_svd_v = nothing ; 253norm_A_NT_spqr_basic_v = nothing ; 254norm_A_N_spqr_null_v = nothing ; 255norm_A_N_spqr_pinv_v = nothing ; 256norm_A_N_spqr_cod_v = nothing ; 257norm_x_pinv_v = nothing ; 258norm_x_QR_dense_v = nothing ; 259norm_r_QR_dense_v = nothing ; 260norm_x_spqr_basic_v = nothing ; 261norm_r_spqr_basic_v = nothing ; 262norm_x_spqr_solve_v = nothing ; 263norm_r_spqr_solve_v = nothing ; 264norm_x_spqr_pinv_minus_x_pinv_v = nothing ; 265norm_x_spqr_cod_minus_x_pinv_v = nothing ; 266cond1_pinv_v = nothing ; 267cond1_cod_v = nothing ; 268norm_w_cod_v = nothing ; 269norm_w_pinv_v = nothing ; 270 271% ignore warnings from inverse power method in ssi 272user_warning_state = warning ('off', 'spqr_rank:overflow') ; 273 274% begin the calculations 275 276if (demo_opts.doprint > 0) 277 disp('Begin calculations') 278 disp(' ') 279end 280 281for i = indexs 282 283 cnt = cnt + 1; 284 285 if (demo_opts.doprint > 0) 286 fprintf ('.') ; 287 if (mod (cnt,50) == 0) 288 fprintf ('\n') ; 289 end 290 end 291 292 %--------------------------------------------------------------------------- 293 % generate the problem 294 %--------------------------------------------------------------------------- 295 296 % select matrix from SJSU singular matrices 297 Problem = SJget(i,index) ; 298 A=Problem.A ; 299 [m,n]=size(A); 300 private_stream = spqr_repeatable (opts.repeatable) ; 301 if (~isempty (private_stream)) 302 b1 = randn (private_stream, m, 1) ; 303 x2 = randn (private_stream, n, 1) ; 304 else 305 b1 = randn (m,1) ; 306 x2 = randn (n,1) ; 307 end 308 b2 = A*x2; % consistent right hand side 309 b = [ b1, b2 ]; 310 311 %--------------------------------------------------------------------------- 312 % for demo_opts.figures == 2 find numerical rank, null space basis, 313 % pseudoinverse soln using MATLAB svd 314 % for demo_opts.figures <= 1 find numerical rank using precomputed sing. 315 % values in SJsingular database 316 %--------------------------------------------------------------------------- 317 318 if demo_opts.figures <= 1 319 % use the precomputed singular values from the SJsingular database 320 s = Problem.svals; 321 normA = max(s); 322 tol = max(m,n) * eps( normA ); 323 rank_svd = SJrank(Problem,tol); 324 else 325 [U,S,V] = svd(full(A)); 326 if m > 1, s = diag(S); 327 elseif m == 1 328 s = S(1); 329 else 330 s = 0; 331 end 332 normA= max(s); 333 tol = max(m,n) * eps(normA); 334 rank_svd = sum(s > tol); 335 end 336 337 if (rank_svd == 0) 338 gap = NaN; 339 else 340 if rank_svd < min(m,n) 341 gap = s(rank_svd) / abs(s(rank_svd + 1)); 342 else 343 gap = Inf; 344 end 345 end 346 347 if demo_opts.figures == 2 348 if (rank_svd == 0) 349 pseudoinverse_svd = zeros(size(A')); 350 N_svd = eye(n,n); 351 NT_svd = eye(m,m); 352 else 353 S_pinv = diag(ones(rank_svd,1)./s(1:rank_svd)); 354 pseudoinverse_svd = V(:,1:rank_svd)*S_pinv*U(:,1:rank_svd)'; 355 N_svd = V(:,rank_svd+1:end); 356 NT_svd = U(:,rank_svd+1:end); 357 end 358 x_pinv = pseudoinverse_svd * b(:,1); % the plots only use x_pinv 359 % for b(:,1), the random rhs 360 norm_x_pinv = norm(x_pinv); 361 norm_A_N_svd = norm(A*N_svd); 362 norm_A_transpose_NT_svd = norm(A'*NT_svd); 363 364 % find basic solution using MATLAB's dense QR routine 365 warning_state = warning ('off','MATLAB:rankDeficientMatrix') ; 366 if m ~= n 367 x = full(A) \ b; 368 else 369 x = full([A,zeros(m,1)]) \ b; 370 x = x(1:n,:); 371 end 372 warning (warning_state) ; % restore the warnings 373 374 norm_x_QR_dense = norm(x(:,1)) ; %the plots only use ||x|| 375 % for the random rhs b(:,1) 376 r = b(:,2)- A*x(:,2); % the plots only use r for 377 % the consistent rhs b(:,2) 378 norm_r_QR_dense = norm(r) / norm(b(:,2)); 379 end 380 381 %--------------------------------------------------------------------------- 382 % run spqr_basic, spqr_null, spqr_pinv and spqr_cod: 383 %--------------------------------------------------------------------------- 384 385 if demo_opts.tol_norm_type == 0 386 opts.tol = tol; 387 else 388 opts.tol_norm_type = demo_opts.tol_norm_type; 389 % calculate tol inside spqr_basic, spqr_null, spqr_pinv and spqr_cod 390 end 391 392 % SPQR_COD uses a complete orthogonal decomposition (COD) of A. By 393 % default the COD is constructed by first factoring A. To select the 394 % option for SPQR_COD which initially factors A', which can have 395 % different numerical properties, use args(1) = 1. 396 397 opts.start_with_A_transpose = demo_opts.start_with_A_transpose ; 398 opts.nsvals_small = demo_opts.nsvals; 399 opts.nsvals_large = demo_opts.nsvals; 400 401 if demo_opts.null_spaces == 2 402 [x_spqr_basic, stats_basic, NT_spqr_basic] = spqr_basic (A,b,opts) ; 403 [N_spqr_null, stats_null] = spqr_null (A,opts) ; 404 [x_spqr_pinv, stats_pinv, N_spqr_pinv, NT] = spqr_pinv (A,b,opts) ; %#ok 405 [x_spqr_cod, stats_cod, N_spqr_cod, NT] = spqr_cod (A,b,opts) ; %#ok 406 elseif demo_opts.null_spaces == 1 407 [x_spqr_basic, stats_basic, NT_spqr_basic] = spqr_basic (A,b,opts) ; 408 [N_spqr_null, stats_null] = spqr_null (A,opts) ; 409 [x_spqr_pinv, stats_pinv, N_spqr_pinv] = spqr_pinv (A,b,opts) ; 410 [x_spqr_cod, stats_cod, N_spqr_cod] = spqr_cod (A,b,opts) ; 411 else 412 [x_spqr_basic, stats_basic] = spqr_basic (A,b,opts) ; 413 [N_spqr_null, stats_null] = spqr_null (A,opts) ; 414 [x_spqr_pinv, stats_pinv] = spqr_pinv (A,b,opts) ; 415 [x_spqr_cod, stats_cod] = spqr_cod (A,b,opts) ; 416 end 417 418 %--------------------------------------------------------------------------- 419 % calculate and save results for figures displaying ranks: 420 %--------------------------------------------------------------------------- 421 422 rank_svd_v(cnt) = rank_svd; 423 gap_v(cnt) = gap; 424 m_v(cnt) = m; 425 n_v(cnt) = n; 426 427 % spqr_basic results: 428 flag_spqr_basic_v(cnt) = stats_basic.flag; 429 if stats_basic.flag == 1 430 % use tol_alt returned by spqr_basic to determine true num. rank 431 if demo_opts.figures == 2 432 rank_svd_basic_v(cnt)= sum(s > stats_basic.tol_alt); 433 else 434 rank_svd_basic_v(cnt)= SJrank(Problem,stats_basic.tol_alt); 435 end 436 else 437 % rank_svd_basic_v(cnt)= rank_svd_v(cnt); 438 % use tol returned by spqr_basic to determine true num. rank 439 if demo_opts.figures == 2 440 rank_svd_basic_v(cnt)= sum(s > stats_basic.tol); 441 else 442 rank_svd_basic_v(cnt)= SJrank(Problem,stats_basic.tol); 443 end 444 end 445 if stats_basic.flag <= 3 446 rank_spqr_basic_v(cnt) = stats_basic.rank ; 447 end 448 if demo_opts.figures >= 1 449 rank_spqr_from_spqr_basic_v(cnt) = stats_basic.rank_spqr ; 450 else 451 rank_spqr_from_spqr_null_v(cnt) = -1; % not calculated 452 end 453 454 % spqr_null results: 455 flag_spqr_null_v(cnt) = stats_null.flag; 456 if stats_null.flag == 1 457 % use tol_alt returned by spqr_null to determine true num. rank 458 if demo_opts.figures == 2 459 rank_svd_null_v(cnt) = sum(s > stats_null.tol_alt); 460 else 461 rank_svd_null_v(cnt) = SJrank(Problem,stats_null.tol_alt); 462 end 463 else 464 % rank_svd_null_v(cnt)= rank_svd_v(cnt); 465 % use tol returned by spqr_null to determine true num. rank 466 if demo_opts.figures == 2 467 rank_svd_null_v(cnt)= sum(s > stats_null.tol); 468 else 469 rank_svd_null_v(cnt)= SJrank(Problem,stats_null.tol); 470 end 471 end 472 if stats_null.flag <= 3 473 rank_spqr_null_v(cnt) = stats_null.rank ; 474 end 475 if demo_opts.figures >= 1 476 rank_spqr_from_spqr_null_v(cnt) = stats_null.rank_spqr ; 477 else 478 rank_spqr_from_spqr_null_v(cnt) = -1 ; % not calculated 479 end 480 481 % spqr_pinv results: 482 flag_spqr_pinv_v(cnt) = stats_pinv.flag; 483 if stats_pinv.flag == 1 484 % use tol_alt returned by spqr_pinv to determine true num. rank 485 if demo_opts.figures == 2 486 rank_svd_pinv_v(cnt) = sum(s > stats_pinv.tol_alt); 487 else 488 rank_svd_pinv_v(cnt) = SJrank(Problem,stats_pinv.tol_alt); 489 end 490 else 491 % rank_svd_pinv_v(cnt)= rank_svd_v(cnt); 492 % use tol returned by spqr_pinv to determine true num. rank 493 if demo_opts.figures == 2 494 rank_svd_pinv_v(cnt)= sum(s > stats_pinv.tol); 495 else 496 rank_svd_pinv_v(cnt)= SJrank(Problem,stats_pinv.tol); 497 end 498 499 end 500 if stats_pinv.flag <= 3 501 rank_spqr_pinv_v(cnt) = stats_pinv.rank ; 502 end 503 if demo_opts.figures >= 1 504 rank_spqr_from_spqr_pinv_v(cnt) = stats_pinv.rank_spqr ; 505 else 506 rank_spqr_from_spqr_pinv_v(cnt) = -1 ; % not calculated 507 end 508 509 % spqr_cod results: 510 flag_spqr_cod_v(cnt) = stats_cod.flag; 511 if stats_cod.flag == 1 512 % use tol_alt returned by spqr_cod to determine true num. rank 513 if demo_opts.figures == 2 514 rank_svd_cod_v(cnt)= sum(s > stats_cod.tol_alt); 515 else 516 rank_svd_cod_v(cnt) = SJrank(Problem,stats_cod.tol_alt); 517 end 518 else 519 % rank_svd_cod_v(cnt)= rank_svd_v(cnt); 520 % use tol returned by spqr_cod to determine true num. rank 521 if demo_opts.figures == 2 522 rank_svd_cod_v(cnt)= sum(s > stats_cod.tol); 523 else 524 rank_svd_cod_v(cnt)= SJrank(Problem,stats_cod.tol); 525 end 526 527 end 528 if stats_cod.flag <= 3 529 rank_spqr_cod_v(cnt) = stats_cod.rank ; 530 end 531 if demo_opts.figures >= 1 532 rank_spqr_from_spqr_cod_v(cnt) = stats_cod.rank_spqr ; 533 else 534 rank_spqr_from_spqr_cod_v(cnt) = -1 ; % not calculated 535 end 536 537 tol_v(cnt) = tol; 538 norm_A_v(cnt) = normA; 539 540 if demo_opts.figures == 2 541 542 %---------------------------------------------------------------------- 543 % calculate and save results for figures displaying null space accuracy 544 %---------------------------------------------------------------------- 545 546 norm_A_N_svd_v(cnt) = norm_A_N_svd; 547 norm_A_NT_svd_v(cnt) = norm_A_transpose_NT_svd; 548 549 if stats_basic.flag <= 3 550 % use spqr_null_mult to form NT_spqr_basic' * A 551 A_transpose_times_NT = spqr_null_mult(NT_spqr_basic,A,0); 552 norm_A_NT_spqr_basic_v(cnt) = norm(full( A_transpose_times_NT )); 553 % For large matrices, rather than forming A_transpose_times_NT, 554 % it is more efficient to use the estimate of 555 % ||A' * NT_spqr_basic|| calculated by SPQR_BASIC (using spqr_ssp) 556 % when SPQR_BASIC returns a null space bases: 557 % norm_A_NT_spqr_basic_v(cnt) = ... 558 % stats_basic.est_norm_A_transpose_times_NT; 559 end 560 561 if stats_null.flag <= 3 562 % use spqr_null_mult to form A * N_spqr_null 563 A_times_N = spqr_null_mult(N_spqr_null,A,3); 564 norm_A_N_spqr_null_v(cnt) = norm(full( A_times_N )); 565 % For large matrices, rather than forming A_times_N, 566 % it is more efficient to use the estimate of 567 % ||A * N|| calculated by SPQR_NULL (using spqr_ssp): 568 % norm_A_N_spqr_null_v(cnt) = stats_null.est_norm_A_times_N; 569 end 570 571 if stats_pinv.flag <= 3 572 % use spqr_null_mult to form A * N_spqr_pinv 573 A_times_N = spqr_null_mult(N_spqr_pinv,A,3); 574 norm_A_N_spqr_pinv_v(cnt) = norm(full( A_times_N )); 575 % For large matrices, rather than forming A_times_N, 576 % it is more efficient to use the estimate of 577 % ||A * N|| calculated by SPQR_PINV (using spqr_ssp): 578 % norm_A_N_spqr_pinv_v(cnt) = stats_pinv.est_norm_A_times_N; 579 end 580 581 if stats_cod.flag <= 3 582 % use spqr_null_mult to form A * N_spqr_cod 583 A_times_N = spqr_null_mult(N_spqr_cod,A,3); 584 norm_A_N_spqr_cod_v(cnt) = norm(full( A_times_N )); 585 % For large matrices, rather than forming A_times_N, 586 % it is more efficient to use the estimate of 587 % ||A * N|| calculated by SPQR_COD (using spqr_ssp): 588 % norm_A_N_spqr_cod_v(cnt) = stats_cod.est_norm_A_times_N; 589 end 590 591 %---------------------------------------------------------------------- 592 % calculate and save results for basic solutions plots 593 %---------------------------------------------------------------------- 594 595 % norm_x_pinv, norm_x_QR_dense and norm_r_QR_dense already computed 596 norm_x_pinv_v(cnt) = norm_x_pinv; 597 norm_x_QR_dense_v(cnt) = norm_x_QR_dense; 598 norm_r_QR_dense_v(cnt) = norm_r_QR_dense; 599 600 norm_x_spqr_basic_v(cnt) = norm(x_spqr_basic(:,1)) ; 601 norm_r_spqr_basic_v(cnt) = norm(b(:,2) - A * x_spqr_basic(:,2)) / ... 602 norm(b(:,2)); 603 opts_spqr_solve.tol = stats_basic.tol; 604 warning_state = warning ('off','MATLAB:rankDeficientMatrix') ; 605 x_spqr_solve = spqr_solve (sparse(A), b, opts_spqr_solve) ; 606 warning (warning_state) ; 607 norm_x_spqr_solve_v(cnt) = norm(x_spqr_solve(:,1)); 608 norm_r_spqr_solve_v(cnt) = norm(b(:,2) - A * x_spqr_solve(:,2)) / ... 609 norm(b(:,2)); 610 611 %---------------------------------------------------------------------- 612 % calculate and save results for pinv solutions plots 613 %---------------------------------------------------------------------- 614 615 if stats_pinv.flag <= 3 616 norm_x_spqr_pinv_minus_x_pinv_v(cnt) = ... 617 norm(x_spqr_pinv(:,1)-x_pinv) / norm(x_pinv); 618 end 619 620 if stats_cod.flag <= 3 621 norm_x_spqr_cod_minus_x_pinv_v(cnt) = ... 622 norm(x_spqr_cod(:,1)-x_pinv) / norm(x_pinv); 623 end 624 625 cond1_cod_v(cnt) = s(1) / s( rank_svd_cod_v(cnt) ); 626 cond1_pinv_v(cnt) = s(1) / s( rank_svd_pinv_v(cnt) ); 627 norm_w_cod_v(cnt) = stats_cod.info_spqr1.norm_E_fro ; 628 norm_w_pinv_v(cnt) = max( ... 629 stats_pinv.stats_spqr_basic.info_spqr1.norm_E_fro, ... 630 stats_pinv.stats_spqr_null.info_spqr1.norm_E_fro) ; 631 end 632 633 %--------------------------------------------------------------------------- 634 % save results to a *.mat file for future reference 635 %--------------------------------------------------------------------------- 636 637 if ( mod(cnt,10) == 0 || cnt == length(indexs) && cnt > 1) 638 % save the information needed to draw the plots; 639 % if later the code is run with ids = 0 640 % save_samples_demo_spqr_rank will be loaded 641 time_required = etime(clock,time_start); 642 save save_samples_demo_spqr_rank ... 643 rank_svd_basic_v rank_spqr_basic_v flag_spqr_basic_v ... 644 rank_svd_null_v rank_spqr_null_v flag_spqr_null_v ... 645 rank_svd_pinv_v rank_spqr_pinv_v flag_spqr_pinv_v ... 646 rank_svd_cod_v rank_spqr_cod_v flag_spqr_cod_v ... 647 rank_spqr_from_spqr_basic_v gap_v tol_v ... 648 rank_spqr_from_spqr_null_v ... 649 rank_spqr_from_spqr_pinv_v ... 650 rank_spqr_from_spqr_cod_v ... 651 norm_A_NT_svd_v norm_A_NT_spqr_basic_v ... 652 norm_A_N_svd_v norm_A_N_spqr_null_v ... 653 norm_A_N_spqr_pinv_v ... 654 norm_A_N_spqr_cod_v ... 655 norm_x_pinv_v norm_x_QR_dense_v norm_r_QR_dense_v ... 656 norm_x_spqr_basic_v norm_r_spqr_basic_v ... 657 norm_x_spqr_solve_v norm_r_spqr_solve_v ... 658 norm_x_spqr_pinv_minus_x_pinv_v ... 659 norm_x_spqr_cod_minus_x_pinv_v cond1_pinv_v cond1_cod_v ... 660 norm_w_cod_v norm_w_pinv_v norm_A_v ... 661 cnt time_required indexs m_v n_v rank_svd_v ... 662 demo_opts 663 end 664 665end 666 667% restore user's warning state 668warning (user_warning_state) ; 669 670% if (demo_opts.dopause) 671% disp('Press enter to see the first plot') 672% pause 673% end 674 675if demo_opts.figures >= 1 676 %--------------------------------------------------------------------------- 677 % plot information for calculated ranks 678 %--------------------------------------------------------------------------- 679 680 figure(1) 681 subplot(2,2,1) 682 plot_ranks(rank_svd_basic_v,rank_spqr_basic_v,... 683 rank_spqr_from_spqr_basic_v,flag_spqr_basic_v,gap_v,'SPQR\_BASIC') 684 subplot(2,2,2) 685 plot_ranks(rank_svd_null_v,rank_spqr_null_v,rank_spqr_from_spqr_null_v,... 686 flag_spqr_null_v,gap_v,'SPQR\_NULL') 687 subplot(2,2,3) 688 plot_ranks(rank_svd_pinv_v,rank_spqr_pinv_v,rank_spqr_from_spqr_pinv_v,... 689 flag_spqr_pinv_v,gap_v,'SPQR\_PINV') 690 subplot(2,2,4) 691 plot_ranks(rank_svd_cod_v,rank_spqr_cod_v,rank_spqr_from_spqr_cod_v,... 692 flag_spqr_cod_v,gap_v,'SPQR\_COD') 693 694 if (demo_opts.doprint > 0) 695 disp(' ') 696 if demo_opts.figures == 1 697 disp('In the figure for each of SPQR_BASIC, SPQR_NULL, ') 698 else 699 disp('In the first figure for each of SPQR_BASIC, SPQR_NULL, ') 700 end 701 disp('SPQR_PINV, SPQR_COD and for SPQR the plots summarize the percent') 702 disp('of matrices where the calculated numerical rank is correct and') 703 disp('the percent of the matrices where the warning flag indicates') 704 disp('that the calculated numerical rank is correct with a warning') 705 disp('flag either 0 or 1 versus the singular value gap, the ratio of') 706 disp('singular number r over singular value number r+1, where r is the') 707 disp('calculated numerical rank.') 708 disp(' ') 709 disp('Note that the percent of the matrices where the new routines ') 710 disp('calculate the correct numerical rank approaches 100 percent') 711 disp('as the singular value gap increases.') 712 disp(' ') 713 disp('The plot is best seen as a full screen plot.') 714 disp(' ') 715 % if (demo_opts.dopause) 716 % disp('Press enter to view the second figure') 717 % pause 718 % end 719 end 720end 721 722drawnow 723 724if demo_opts.figures == 2 725 726 %--------------------------------------------------------------------------- 727 % plot information for null spaces 728 %--------------------------------------------------------------------------- 729 730 figure(2) 731 subplot(2,2,1) 732 plot_null_spaces(norm_A_NT_svd_v,tol_v, ... 733 norm_A_NT_spqr_basic_v,flag_spqr_basic_v,'SPQR\_BASIC') 734 subplot(2,2,2) 735 plot_null_spaces(norm_A_N_svd_v,tol_v, ... 736 norm_A_N_spqr_null_v,flag_spqr_null_v,'SPQR\_NULL') 737 subplot(2,2,3) 738 plot_null_spaces(norm_A_N_svd_v,tol_v, ... 739 norm_A_N_spqr_pinv_v,flag_spqr_pinv_v,'SPQR\_PINV') 740 subplot(2,2,4) 741 plot_null_spaces(norm_A_N_svd_v,tol_v, ... 742 norm_A_N_spqr_cod_v,flag_spqr_cod_v,'SPQR\_COD') 743 744 if (demo_opts.doprint > 0) 745 disp(' ') 746 disp('In the second figure ||AN||, where N is a calculated') 747 disp('orthonormal basis for the numerical null space, or, in the case') 748 disp('of SPQR_BASIC, ||(transpose of A) N||, normalized by the ') 749 disp('tolerance defining the numerical rank, is plotted for null') 750 disp('space bases calculated by MATLAB''s SVD and by SPQR_BASIC,') 751 disp('SPQR_NULL, SPQR_PINV and SPQR_COD.') 752 disp(' ') 753 disp('Note that the null space bases calculated by the new routines') 754 disp('are generally as good as the bases calculated by MATLAB''s SVD.') 755 disp('The tolerance used for normalization in the plots is ') 756 disp('O(relative machine precision times ||A||).') 757 disp(' ') 758 disp('The plot is best seen as a full screen plot.') 759 disp(' ') 760 % if (demo_opts.dopause) 761 % disp('Press enter to view the third figure') 762 % pause 763 % end 764 end 765 766 drawnow 767 768 %--------------------------------------------------------------------------- 769 % plot information for basic solutions 770 %--------------------------------------------------------------------------- 771 772 figure(3) 773 plot_basic(norm_x_pinv_v,norm_x_QR_dense_v, norm_r_QR_dense_v, ... 774 norm_x_spqr_basic_v, norm_r_spqr_basic_v, ... 775 norm_x_spqr_solve_v, norm_r_spqr_solve_v,flag_spqr_basic_v) 776 777 if (demo_opts.doprint) 778 disp(' ') 779 disp('In the third figure the left plot pictures ||x|| / ||x_PINV|| ') 780 disp('where x is a basic solution to min || b - A x ||') 781 disp('calculated by MATLAB''s dense matrix QR algorithm,') 782 disp('by SPQR_SOLVE or by SPQR_BASIC. SPQR_SOLVE is part of') 783 disp('SuiteSparseQR and can be used to construct basic solutions to') 784 disp('min ||b - A x||. x_PINV is computed using MATLAB''s PINV. In the') 785 disp('left hand plot the vectors b in min || b - A x || are random') 786 disp('vectors. The right plot pictures ||r|| = || b - A x || for x') 787 disp('vectors calculated using MATLAB''s dense matrix QR algorithm,') 788 disp('by SPQR_SOLVE or by SPQR_BASIC. In the right hand plot the') 789 disp('vectors b in min || b - A x || are of the form b = Ax where x') 790 disp('is a random vector.') 791 disp(' ') 792 disp('In the left hand plot note that for most, but not all, of the') 793 disp('matrices the norm of the basic solution calculated by SPQR_BASIC') 794 disp('or by SPQR_SOLVE is the same order of magnitude as the norm of') 795 disp('the pseudoinverse solution. Also note that SPQR_SOLVE calculates') 796 disp('a large norm solution more frequently than does SPQR_BASIC.') 797 disp(' ') 798 disp('In the right hand plot note that for most, but not all, matrices') 799 disp('the residuals for solutions calculated by SPQR_BASIC or by ') 800 disp('SPQR_SOLVE are similar in size to the residual for solutions') 801 disp('calcluated by MATLAB''s dense QR factorization. ') 802 disp(' ') 803 disp('The plot is best seen as a full screen plot.') 804 disp(' ') 805 % if (demo_opts.dopause) 806 % disp('Press enter to view the fourth figure.') 807 % pause 808 % end 809 end 810 811 %--------------------------------------------------------------------------- 812 % plot information for pseudoinverse solutions 813 %--------------------------------------------------------------------------- 814 815 figure(4) 816 plot_pinv(norm_x_spqr_pinv_minus_x_pinv_v, ... 817 norm_x_spqr_cod_minus_x_pinv_v, cond1_pinv_v, cond1_cod_v, ... 818 norm_w_pinv_v, norm_w_cod_v, flag_spqr_pinv_v, ... % tol_v, ... 819 flag_spqr_cod_v, norm_A_v) 820 821 if (demo_opts.doprint > 0) 822 disp(' ') 823 disp('In the fourth figure the left graph plots || x - x_PINV ||') 824 disp(' / ||x_PINV|| for x produced by SPQR_PINV for the matrices') 825 disp('where SPQR_PINV returns a flag of 0. x_PINV is') 826 disp('calculated using MATLAB''s PINV routine. Also part of a') 827 disp('perturbation theory result from "Matrix Perturbation Theory" by') 828 disp('Stewart and Sun, page 157, is plotted. The right hand graph') 829 disp('is the same plot for x produced by SPQR_COD for the matrices') 830 disp('where SPQR_COD returns a flag of 0. For the plots the ') 831 disp('vectors b in min||b-Ax|| are random vectors') 832 disp(' ') 833 disp('Note that the accuracies of the pseudoinverse solutions') 834 disp('calculated by SPQR_COD and, in most cases, by SPQR_PINV are as') 835 disp('good as or nearly as good as predicted by the perturbation') 836 disp('theory.') 837 disp(' ') 838 disp('The plot is best seen as a full screen plot.') 839 end 840 841 drawnow 842 843end 844 845if (demo_opts.doprint >= 0) 846 fprintf ('\n') ; 847end 848failures = 0 ; 849 850%------------------------------------------------------------------------------- 851% check that numerical rank calculations are accurate 852%------------------------------------------------------------------------------- 853 854ifail_spqr_basic = find( rank_svd_basic_v ~= rank_spqr_basic_v & ... 855 (flag_spqr_basic_v == 0 | flag_spqr_basic_v == 1) ); 856ifail_spqr_null = find( rank_svd_null_v ~= rank_spqr_null_v & ... 857 (flag_spqr_null_v == 0 | flag_spqr_null_v == 1) ); 858ifail_spqr_pinv = find( rank_svd_pinv_v ~= rank_spqr_pinv_v & ... 859 (flag_spqr_pinv_v == 0 | flag_spqr_pinv_v == 1) ); 860ifail_spqr_cod = find( rank_svd_cod_v ~= rank_spqr_cod_v & ... 861 (flag_spqr_cod_v == 0 | flag_spqr_cod_v == 1) ); 862 863nfail_spqr_basic = length(ifail_spqr_basic); 864nfail_spqr_null = length(ifail_spqr_null); 865nfail_spqr_pinv = length(ifail_spqr_pinv); 866nfail_spqr_cod = length(ifail_spqr_cod); 867 868SJid_fail_spqr_basic = sort( indexs( ifail_spqr_basic ) ); 869SJid_fail_spqr_null = sort( indexs( ifail_spqr_null ) ); 870SJid_fail_spqr_pinv = sort( indexs( ifail_spqr_pinv ) ); 871SJid_fail_spqr_cod = sort( indexs( ifail_spqr_cod ) ); 872 873SJid_fail = union(SJid_fail_spqr_basic, SJid_fail_spqr_null) ; 874SJid_fail = union(SJid_fail, SJid_fail_spqr_pinv) ; 875SJid_fail = union(SJid_fail, SJid_fail_spqr_cod) ; 876 877if (demo_opts.doprint > 0) 878 disp(' ') 879 disp('Check that the routines reliably calculate the numerical rank.') 880end 881 882if (demo_opts.doprint >= 0) 883 fprintf ('SPQR_BASIC: ') ; 884end 885iflagis0_or_1 = find(flag_spqr_basic_v == 0 | flag_spqr_basic_v == 1); 886nflagis0_or_1 = length(iflagis0_or_1); 887 888failures = failures + nfail_spqr_basic ; 889 890if (demo_opts.doprint > 0) 891 disp([' ',int2str(nflagis0_or_1 - nfail_spqr_basic),... 892 ' matrices have the correct numerical rank from the set of']) 893 if ( nfail_spqr_basic == 0 ) 894 disp([' ',int2str(nflagis0_or_1),... 895 ' matrices with a warning flag of 0 or 1.']) 896 elseif ( nfail_spqr_basic == 1 ) 897 disp([' ',int2str(nflagis0_or_1),' matrices with a warning ',... 898 'flag of 0 or 1. Failure is for matrix from']) 899 disp([' ','the SJSU Singular Matrix Database with SJid = ',... 900 int2str(SJid_fail_spqr_basic),'.']) 901 elseif ( nfail_spqr_basic >= 2 ) 902 disp([' ',int2str(nflagis0_or_1),' matrices with a warning ',... 903 'flag of 0 or 1. Failure is for matrices from']) 904 disp([' ','the SJSU Singular Matrix Database with SJid = ',... 905 int2str(SJid_fail_spqr_basic),'.']) 906 end 907elseif (demo_opts.doprint >= 0) 908 if ( nfail_spqr_basic == 0 ) 909 fprintf ('%3d matrices, failed: %d\n', nflagis0_or_1, ... 910 nfail_spqr_basic) ; 911 else 912 fprintf ('%3d matrices, failed: %d with SJid =', ... 913 nflagis0_or_1, nfail_spqr_basic) ; 914 fprintf(' %d', SJid_fail_spqr_basic ) ; 915 fprintf('\n') ; 916 end 917end 918 919if (demo_opts.doprint >= 0) 920 fprintf ('SPQR_NULL: ') ; 921end 922iflagis0_or_1 = find(flag_spqr_null_v == 0 | flag_spqr_null_v == 1); 923nflagis0_or_1 = length(iflagis0_or_1); 924 925failures = failures + nfail_spqr_null ; 926 927if (demo_opts.doprint > 0) 928 disp([' ',int2str(nflagis0_or_1 - nfail_spqr_null),... 929 ' matrices have the correct numerical rank from the set of']) 930 if ( nfail_spqr_null == 0 ) 931 disp([' ',int2str(nflagis0_or_1),... 932 ' matrices with a warning flag of 0 or 1.']) 933 elseif ( nfail_spqr_null == 1 ) 934 disp([' ',int2str(nflagis0_or_1),' matrices with a warning ',... 935 'flag of 0 or 1. Failure is for matrix from']) 936 disp([' ','the SJSU Singular Matrix Database with SJid = ',... 937 int2str(SJid_fail_spqr_null),'.']) 938 elseif ( nfail_spqr_null >= 2 ) 939 disp([' ',int2str(nflagis0_or_1),' matrices with a warning ',... 940 'flag of 0 or 1. Failure is for matrices from']) 941 disp([' ','the SJSU Singular Matrix Database with SJid = ',... 942 int2str(SJid_fail_spqr_null),'.']) 943 end 944elseif (demo_opts.doprint >= 0) 945 if ( nfail_spqr_null == 0 ) 946 fprintf ('%3d matrices, failed: %d\n', nflagis0_or_1, ... 947 nfail_spqr_null) ; 948 else 949 fprintf ('%3d matrices, failed: %d with SJid =', ... 950 nflagis0_or_1, nfail_spqr_null) ; 951 fprintf(' %d', SJid_fail_spqr_null ) ; 952 fprintf('\n') ; 953 end 954end 955 956if (demo_opts.doprint >= 0) 957 fprintf ('SPQR_PINV: ') ; 958end 959iflagis0_or_1 = find(flag_spqr_pinv_v == 0 | flag_spqr_pinv_v == 1); 960nflagis0_or_1 = length(iflagis0_or_1); 961 962failures = failures + nfail_spqr_pinv ; 963 964if (demo_opts.doprint > 0) 965 disp([' ',int2str(nflagis0_or_1 - nfail_spqr_pinv),... 966 ' matrices have the correct numerical rank from the set of']) 967 if ( nfail_spqr_pinv == 0 ) 968 disp([' ',int2str(nflagis0_or_1),... 969 ' matrices with a warning flag of 0 or 1.']) 970 elseif ( nfail_spqr_pinv == 1 ) 971 disp([' ',int2str(nflagis0_or_1),' matrices with a warning ',... 972 'flag of 0 or 1. Failure is for matrix from']) 973 disp([' ','the SJSU Singular Matrix Database with SJid = ',... 974 int2str(SJid_fail_spqr_pinv),'.']) 975 elseif ( nfail_spqr_pinv >= 2 ) 976 disp([' ',int2str(nflagis0_or_1),' matrices with a warning ',... 977 'flag of 0 or 1. Failure is for matrices from']) 978 disp([' ','the SJSU Singular Matrix Database with SJid = ',... 979 int2str(SJid_fail_spqr_pinv),'.']) 980 end 981elseif (demo_opts.doprint >= 0) 982 if ( nfail_spqr_pinv == 0 ) 983 fprintf ('%3d matrices, failed: %d\n', nflagis0_or_1, ... 984 nfail_spqr_pinv) ; 985 else 986 fprintf ('%3d matrices, failed: %d with SJid =', ... 987 nflagis0_or_1, nfail_spqr_pinv) ; 988 fprintf(' %d', SJid_fail_spqr_pinv ) ; 989 fprintf('\n') ; 990 end 991end 992 993if (demo_opts.doprint >= 0) 994 fprintf ('SPQR_COD: ') ; 995end 996iflagis0_or_1 = find(flag_spqr_cod_v == 0 | flag_spqr_cod_v == 1); 997nflagis0_or_1 = length(iflagis0_or_1); 998 999failures = failures + nfail_spqr_cod ; 1000 1001if (demo_opts.doprint > 0) 1002 disp([' ',int2str(nflagis0_or_1 - nfail_spqr_cod),... 1003 ' matrices have the correct numerical rank from the set of']) 1004 if ( nfail_spqr_cod == 0 ) 1005 disp([' ',int2str(nflagis0_or_1),... 1006 ' matrices with a warning flag of 0 or 1.']) 1007 elseif ( nfail_spqr_cod == 1 ) 1008 disp([' ',int2str(nflagis0_or_1),' matrices with a warning ',... 1009 'flag of 0 or 1. Failure is for matrix from']) 1010 disp([' ','the SJSU Singular Matrix Database with SJid = ',... 1011 int2str(SJid_fail_spqr_cod),'.']) 1012 elseif ( nfail_spqr_cod >= 2 ) 1013 disp([' ',int2str(nflagis0_or_1),' matrices with a warning ',... 1014 'flag of 0 or 1. Failure is for matrices from']) 1015 disp([' ','the SJSU Singular Matrix Database with SJid = ',... 1016 int2str(SJid_fail_spqr_cod),'.']) 1017 end 1018elseif (demo_opts.doprint >= 0) 1019 if ( nfail_spqr_cod == 0 ) 1020 fprintf ('%3d matrices, failed: %d\n', nflagis0_or_1, ... 1021 nfail_spqr_cod) ; 1022 else 1023 fprintf ('%3d matrices, failed: %d with SJid =', ... 1024 nflagis0_or_1, nfail_spqr_cod) ; 1025 fprintf(' %d', SJid_fail_spqr_cod ) ; 1026 fprintf('\n') ; 1027 end 1028end 1029 1030%------------------------------------------------------------------------------- 1031% return results 1032%------------------------------------------------------------------------------- 1033 1034if (nargout > 0) 1035 nfailures = failures ; 1036end 1037 1038if nargout > 1 1039 SJid_failures = SJid_fail; 1040end 1041 1042if (failures > 0) 1043 fprintf ('demo_spqr_rank: %d failures\n', failures) ; 1044end 1045 1046%------------------------------------------------------------------------------- 1047% subfunctions 1048%------------------------------------------------------------------------------- 1049 1050 1051%**************************************************************** 1052%**** plot_ranks 1053%**************************************************************** 1054 1055 1056function plot_ranks(rank_svd_v,rank_spqr_cod_v,rank_spqr_v, ... 1057 flag_spqr_cod_v, gap_v, method) 1058% plot the percent of matrices with a gap larger than a specified gap 1059% that have the correct numerical rank, that have a 1060% warning flag = 0 or 1' and that have a correct rank for spqr 1061% versus gap in singular values 1062% the input is computed by demo_reliable_spqr 1063 1064% ncor = sum (rank_svd_v == rank_spqr_cod_v); 1065 1066gap_tol = 10 .^ ((0:32)/2); 1067z = zeros (1, length (gap_tol)) ; 1068per_cor = z ; 1069per_cor_spqr = z ; 1070per_flag_0_1 = z ; 1071 1072i = 0; 1073for gap = gap_tol 1074 i = i+1; 1075 igt = find( gap_v >= gap & rank_svd_v > 0 ); 1076 % for igt use rank_svd_v > 0 to exclude rare cases where SJrank returns -1 1077 ncorrect = sum( rank_svd_v(igt) == rank_spqr_cod_v(igt) ); 1078 per_cor(i) = 100*ncorrect / length(igt); 1079 ncorrect = sum( rank_svd_v(igt) == rank_spqr_v(igt) ); 1080 per_cor_spqr(i) = 100*ncorrect / length(igt); 1081 n_flag_0_1 = sum( flag_spqr_cod_v(igt) <= 1 ); 1082 per_flag_0_1(i) = 100 * n_flag_0_1 / length(igt); 1083end 1084 1085semilogx(gap_tol,per_cor_spqr,'ks--',gap_tol,per_flag_0_1, ... 1086 'bo--',gap_tol,per_cor,'rx--'); 1087axisv=axis; 1088axisv3 = axisv(3); 1089%axisv3=60; 1090axisv2=[axisv(1) gap_tol(end) axisv3 axisv(4)+.01]; 1091axis(axisv2); 1092fs = 12; 1093ylabel('Percent','fontsize',fs) 1094xlabel('Gap in the singular value spectrum bigger than','fontsize',fs) 1095title(['Percent numerical rank correct and warning flag = 0 or 1',char(10),... 1096 'versus gap in singular values'],'fontsize',fs) 1097legend('% SPQR rank correct',['% flag = 0 or 1 in ',method],... 1098 ['% ',method,' rank correct'],'location','se') 1099grid 1100set(gca,'fontsize',fs) 1101 1102 1103%**************************************************************** 1104%**** plot_ basic 1105%**************************************************************** 1106 1107function plot_basic(norm_x_pinv_v,norm_x_QR_dense_v, norm_r_QR_dense_v, ... 1108 norm_x_spqr_basic_v, norm_r_spqr_basic_v, ... 1109 norm_x_spqr_solve_v, norm_r_spqr_solve_v,flag_v) 1110% plot the quality of the basic solutions 1111% the input is computed by demo_spqr_rank 1112 1113iflagis0 = find( flag_v == 0 ); 1114nflagis0 = length(iflagis0); 1115 1116norm_x_pinv0 = norm_x_pinv_v(iflagis0); 1117norm_x_spqr_basic0 = norm_x_spqr_basic_v(iflagis0); 1118norm_x_QR_dense0 = norm_x_QR_dense_v(iflagis0); 1119norm_x_spqr_solve0 = norm_x_spqr_solve_v(iflagis0); 1120norm_ratio_spqr_basic = norm_x_spqr_basic0 ./ norm_x_pinv0; 1121norm_ratio_QR_dense = norm_x_QR_dense0 ./ norm_x_pinv0; 1122norm_ratio_spqr_solve = norm_x_spqr_solve0 ./ norm_x_pinv0; 1123[norm_ratio_spqr_basic,isort]=sort(norm_ratio_spqr_basic); 1124norm_ratio_QR_dense = norm_ratio_QR_dense(isort); 1125norm_ratio_spqr_solve = norm_ratio_spqr_solve(isort); 1126 1127 1128subplot(1,2,1) 1129 1130semilogy(1:nflagis0,norm_ratio_QR_dense, 'bo',... 1131 1:nflagis0,norm_ratio_spqr_solve, 'ks',... 1132 1:nflagis0,norm_ratio_spqr_basic, 'rx') 1133fs = 12; 1134ylabel(' || x || / ||x_{ PINV} ||','fontsize',fs) 1135xlabel('matrix: ordered by ||x_{SPQR\_BASIC} || / ||x_{PINV} ||','fontsize',fs) 1136title([' Comparison of the norms of basic ',... 1137 'solutions divided by ||x_{PINV}||',char(10),'for ', ... 1138 int2str(nflagis0) , ' matrices with flag = 0', ... 1139 ' in SPQR\_BASIC'],'fontsize',fs) 1140legend('dense QR','SPQR\_SOLVE',... 1141 ' SPQR\_BASIC ','location','best') 1142grid 1143set(gca,'fontsize',fs) 1144 1145subplot(1,2,2) 1146norm_r_spqr_basic0 = norm_r_spqr_basic_v(iflagis0); 1147norm_r_QR_dense0 = norm_r_QR_dense_v(iflagis0); 1148norm_r_spqr_solve0 = norm_r_spqr_solve_v(iflagis0); 1149[norm_r_spqr_basic,isort]=sort(norm_r_spqr_basic0); 1150norm_r_QR_dense = norm_r_QR_dense0(isort); 1151norm_r_spqr_solve = norm_r_spqr_solve0(isort); 1152 1153semilogy(1:nflagis0,norm_r_QR_dense,'bo',... 1154 1:nflagis0,norm_r_spqr_solve,'ks',... 1155 1:nflagis0,norm_r_spqr_basic,'rx') 1156 1157fs = 12; 1158ylabel(' || r || / || b ||','fontsize',fs) 1159xlabel('matrix: ordered by ||r_{SPQR\_BASIC} || / ||b||','fontsize',fs) 1160title([' Comparison of the norms of residuals, ',... 1161 'r = b-A*x, divided by ||b||',char(10),' for ', ... 1162 int2str(length(iflagis0)) , ' matrices with flag = 0', ... 1163 ' in SPQR\_BASIC'],'fontsize',fs) 1164legend('dense QR','SPQR\_SOLVE',... 1165 ' SPQR\_BASIC ','location','best') 1166set(gca,'fontsize',fs) 1167grid 1168 1169%**************************************************************** 1170%**** plot_null_spaces 1171%**************************************************************** 1172 1173function plot_null_spaces(norm_A_N_svd_v,tol_v, ... 1174 norm_A_N_v,flag_v,method) 1175% plot the quality of the null spaces 1176% the input is computed by demo_spqr_rank 1177 1178iflagis0 = find( flag_v == 0 ); 1179nflagis0 = length(iflagis0); 1180n_method_better = sum( norm_A_N_v(iflagis0) <= norm_A_N_svd_v(iflagis0) ); 1181percent_method_better = 100*n_method_better / length(iflagis0); 1182quality_method = norm_A_N_v(iflagis0) ./ tol_v(iflagis0); 1183quality_svd = norm_A_N_svd_v(iflagis0) ./ tol_v(iflagis0); 1184[ignore, isort]=sort(quality_method); %#ok 1185clear ignore 1186fs = 12; 1187semilogy(1:nflagis0,quality_svd(isort),'bo',1:nflagis0, ... 1188 quality_method(isort),'rx'); 1189%axis1 = axis; 1190%axis1(3)=1.e-8; 1191%axis(axis1) 1192if strcmp(method,'SPQR\_BASIC') 1193 ylabel('|| A^T N || / tolerance','fontsize',fs) 1194 xlabel(['matrix: ordered by ||A^TN|| / tolerance for ',method],... 1195 'fontsize',fs) 1196 title(['Null space quality when flag in ',method,' is 0.',char(10),... 1197 method,' has ||A^T N|| smaller in ',... 1198 int2str(percent_method_better),'% of cases.'],'fontsize',fs) 1199else 1200 ylabel('|| A N || / tolerance','fontsize',fs) 1201 xlabel(['matrix: ordered by ||AN|| / tolerance for ',method],... 1202 'fontsize',fs) 1203 title(['Null space quality when flag in ',method,' is 0.',char(10),... 1204 method,' has ||AN|| smaller in ',... 1205 int2str(percent_method_better),'% of cases.'],'fontsize',fs) 1206end 1207legend('SVD null space',[method,' null space'],'location','SE') 1208grid 1209set(gca,'fontsize',fs) 1210 1211%**************************************************************** 1212%**** plot_pinv 1213%**************************************************************** 1214 1215function plot_pinv(norm_x_spqr_pinv_minus_x_pinv_v, ... 1216 norm_x_spqr_cod_minus_x_pinv_v, cond1_pinv_v, cond1_cod_v, ... 1217 norm_w_pinv_v, norm_w_cod_v, flag_spqr_pinv_v, ... % tol_v, ... 1218 flag_spqr_cod_v, norm_A_v) 1219 1220% plot the quality of the pseudoinverse solutions 1221% the input is computed by demo_spqr_rank 1222 1223% draw plot for results from spqr_pinv 1224subplot(1,2,1) 1225iflagis0 = find(flag_spqr_pinv_v == 0); 1226nflagis0 = length(iflagis0); 1227%perturbation_theory_v = cond1_v .* (max(tol_v,norm_w_pinv_v) ./ norm_A_v ); 1228perturbation_theory_v = cond1_pinv_v .* max(10*eps,norm_w_pinv_v ./ norm_A_v ); 1229 1230perturbation_theory0 = perturbation_theory_v(iflagis0); 1231[perturbation_theory,isort] = sort(perturbation_theory0); 1232norm_x_spqr_pinv_minus_x_pinv0 = norm_x_spqr_pinv_minus_x_pinv_v(iflagis0); 1233norm_x_spqr_pinv_minus_x_pinv = norm_x_spqr_pinv_minus_x_pinv0(isort); 1234fs = 12; 1235loglog(perturbation_theory, norm_x_spqr_pinv_minus_x_pinv,'o',... 1236 [perturbation_theory(1),max(1,perturbation_theory(end))], ... 1237 [perturbation_theory(1),max(1,perturbation_theory(end))],... 1238 'r--','linewidth',2 ) 1239ylabel(' || x_{SPQR\_PINV } - x_{PINV} || / ||x_{PINV} ||','fontsize',fs) 1240xlabel(' ( \sigma_1(A) / \sigma_r(A) ) max(10 \epsilon, ||w|| / ||A|| ) ',... 1241 'fontsize',fs) 1242title([' Comparison of the pseudoinverse solutions ',... 1243 'returned by SPQR\_PINV ',char(10),' and MATLAB''s PINV for ', ... 1244 int2str(nflagis0) , ' matrices with flag = 0', ... 1245 ' in SPQR\_PINV'],'fontsize',fs) 1246legend('|| x_{SPQR\_PINV } - x_{PINV} || / ||x_{PINV} ||',... 1247 '( \sigma_1(A)/\sigma_r(A) ) max(10\epsilon, ||w|| / ||A|| ) ',... 1248 'location','best') 1249set(gca,'fontsize',fs) 1250grid 1251shg 1252 1253% draw plot for results from spqr_cod 1254subplot(1,2,2) 1255iflagis0 = find(flag_spqr_cod_v == 0); 1256nflagis0 = length(iflagis0); 1257%perturbation_theory_v = cond1_v .* (max(tol_v,norm_w_cod_v) ./ norm_A_v); 1258perturbation_theory_v = cond1_cod_v .* max(10*eps,norm_w_cod_v ./ norm_A_v ); 1259perturbation_theory0 = perturbation_theory_v(iflagis0); 1260[perturbation_theory,isort] = sort(perturbation_theory0); 1261norm_x_spqr_cod_minus_x_pinv0 = norm_x_spqr_cod_minus_x_pinv_v(iflagis0); 1262norm_x_spqr_cod_minus_x_pinv = norm_x_spqr_cod_minus_x_pinv0(isort); 1263fs = 12; 1264loglog(perturbation_theory, norm_x_spqr_cod_minus_x_pinv,'o',... 1265 [perturbation_theory(1),max(1,perturbation_theory(end))], ... 1266 [perturbation_theory(1),max(1,perturbation_theory(end))],... 1267 'r--','linewidth',2 ) 1268ylabel('|| x_{SPQR\_COD } - x_{PINV} || / ||x_{PINV} ||','fontsize',fs) 1269xlabel('( \sigma_1(A)/\sigma_r(A) ) max(10\epsilon, ||w|| / ||A||) ',... 1270 'fontsize',fs) 1271title([' Comparison of the pseudoinverse solutions ',... 1272 'returned by SPQR\_COD ',char(10),' and MATLAB''s PINV for ', ... 1273 int2str(nflagis0) , ' matrices with flag = 0', ... 1274 ' in SPQR\_COD'],'fontsize',fs) 1275legend('|| x_{SPQR\_COD } - x_{PINV} || / ||x_{PINV} ||',... 1276 '( \sigma_1(A) / \sigma_r(A) ) max(10 \epsilon, ||w|| / ||A||) ',... 1277 'location','best') 1278set(gca,'fontsize',fs) 1279grid 1280shg 1281