1function [x,stats,N,NT]= spqr_cod(A,varargin) 2%SPQR_COD approximate pseudoinverse solution to min(norm(B-A*x) 3% for a rank deficient matrix A. 4% 5% [x,stats,N,NT] = spqr_cod (A,B,opts) 6% 7% This function returns an approximate pseudoinverse solution to 8% min || B - A x || (1) 9% for rank deficient matrices A. 10% The psuedoinverse solution is the min norm solution to the 11% least squares problem (1). 12% 13% Optionally returns statistics including the numerical rank of the matrix A 14% for tolerance tol (i.e. the number of singular values > tol), and 15% orthnormal bases for the numerical null spaces of A and A'. 16% The solution is approximate since the algorithm allows small perturbations in 17% A (columns of A may be changed by no more than opts.tol). 18% 19% This routine calculates an approximate complete orthogonal decomposition of 20% A. The routine can be more accurate -- and more expensive -- than spqr_pinv. 21% 22% Input: 23% A -- an m by n matrix 24% B -- an m by p right hand side matrix 25% opts (optional) -- type 'help spqr_rank_opts' for details. 26% 27% Output: 28% x -- this n by p matrix contains psuedoinverse solutions to (1). 29% stats -- statistics, type 'help spqr_rank_stats' for details. 30% N -- orthonormal basis for numerical null space of A. 31% NT -- orthonormal basis for numerical null space of A'. 32% 33% Example: 34% 35% A = sparse(gallery('kahan',100)); 36% B = randn(100,1); B = B / norm(B); 37% x = spqr_cod(A,B); 38% x_pinv = pinv(full(A))*B; 39% rel_error_in_x = norm(x - x_pinv) / norm(x_pinv) 40% % or 41% [x,stats,N,NT]=spqr_cod(A,B); 42% norm_A_times_N = norm(full(spqr_null_mult(N,A,3))) 43% norm_A_transpose_times_NT = norm(full(spqr_null_mult(NT,A,0))) 44% % or 45% opts = struct('tol',1.e-5) ; 46% [x,stats]=spqr_cod(A,B,opts); 47% stats 48% 49% See also spqr_basic, spqr_null, spqr_pinv, spqr_ssi, spqr_ssp 50 51% Copyright 2012, Leslie Foster and Timothy A Davis. 52 53% Algorithm: First spqr is used to construct a QR factorization 54% of A: 55% m by n matrix A: A*P1 = Q1*R where R' = [ R1' 0 ] + E1, 56% R1 is a k by n upper trapezoidal matrix 57% or A': 58% n by m matrix A': A'*P1 = Q1*R where R' = [ R1' 0 ] + E1, 59% R1 is a k by m upper trapezoidal matrix 60% and where E1 is a small error matrix. 61% 62% Next a QR factorization of R1' is calculated: 63% R1' * P2 = Q2 * R2 where R2' = [ T' 0 ] and T is k by k, upper 64% triangular. This determines an approximate complete orthogonal 65% decomposition of A: 66% [ P2 0 ] [ T' 0 ] 67% A = Q1 * [ 0 I ] * [ 0 0 ] * Q2' * P1' + (E1*P1'). (2) 68% or of A': 69% [ T 0 ] [ P2' 0 ] 70% A = P1 *Q2 * [ 0 0 ] *[ 0 I ] * Q1' + (P1*E1'). (2) 71% 72% This is then used to calculate an approximate pseudoinverse solution 73% to (1) and an orthogonal basis for the left and right null spaces of 74% A, when they are requested. Subspace iteration, using the routine 75% spqr_ssi, is applied to T to determine if the rank returned by spqr is 76% correct and also, often, to determine the correct numerical 77% rank. When the two ranks are different deflation (see SIAM SISC, 78% 11:519-530, 1990.) is used in the caculation of the psuedoinverse 79% solution. 80% 81% When opts.start_with_A_transpose is 1 (the default value is 0) then 82% initially spqr constructs a QR factorization of A'. By default, A 83% is factorized. 84 85%------------------------------------------------------------------------------- 86% get opts: tolerance and number of singular values to estimate 87%------------------------------------------------------------------------------- 88 89[B,opts,stats,start_tic,ok] = spqr_rank_get_inputs (A, 1, varargin {:}) ; 90 91if (~ok || nargout > 4) 92 error ('usage: [x,stats,N,NT] = spqr_cod (A,B,opts)') ; 93end 94 95% get the options 96tol = opts.tol ; 97nsvals_small = opts.nsvals_small ; 98nsvals_large = opts.nsvals_large ; 99get_details = opts.get_details ; 100start_with_A_transpose = opts.start_with_A_transpose ; 101 102% set the order of the stats fields 103% stats.flag, stats.rank, stats.rank_spqr, stats.rank_spqr (if get_details 104% >= 1), stats.tol, stats.tol_alt, stats.normest_A (if calculated), 105% stats.est_sval_upper_bounds, stats.est_sval_lower_bounds, and 106% stats.sval_numbers_for_bounds already initialized in spqr_rank_get_inputs 107if nargout >= 3 108 stats.est_norm_A_times_N = -1 ; 109end 110if nargout == 4 111 stats.est_norm_A_transpose_times_NT = -1 ; 112end 113if get_details == 2 114 stats.stats_ssi = -1 ; 115end 116% order for the additional stats fields for case where get_details is 1 will be 117% set using spqr_rank_order_fields at the end of the routine 118 119 120if (get_details == 1) 121 stats.opts_used = opts ; 122 stats.time_basis = 0 ; 123end 124 125[m,n] = size (A) ; 126 127%------------------------------------------------------------------------------- 128% first QR factorization of A or A', and initial estimate of num. rank, via spqr 129%------------------------------------------------------------------------------- 130 131if (start_with_A_transpose) 132 % Compute Q1*R = A(p1,:)', do not compute C, and keep Q1. 133 [Q1,R,C,p1,info_spqr1] = ... 134 spqr_wrapper (A', [ ], tol, 'keep Q', get_details) ; 135 B = B(p1,:); % or: p=1:m; p(p1) = 1:length(p1); B=B(p,:); 136else 137 % Compute Q1*R = A(:,p1), C=Q1'*B. Keep Q only if needed for NT 138 if nargout <= 3 139 [Q1,R,C,p1,info_spqr1] = ... 140 spqr_wrapper (A, B, tol, 'discard Q', get_details) ; 141 else 142 [Q1,R,C,p1,info_spqr1] = ... 143 spqr_wrapper (A, B, tol, 'keep Q', get_details) ; 144 end 145end 146 147% the next line is the same as: rank_spqr = size (R,1) ; 148rank_spqr = info_spqr1.rank_A_estimate; 149norm_E_fro = info_spqr1.norm_E_fro ; 150 151% save the stats 152if (get_details == 1 || get_details == 2) 153 stats.rank_spqr = rank_spqr ; 154end 155 156if (get_details == 1) 157 stats.info_spqr1 = info_spqr1 ; 158end 159 160%------------------------------------------------------------------------------- 161% second QR factorization of R', with zero tolerance 162%------------------------------------------------------------------------------- 163 164if (start_with_A_transpose) 165 % Compute Q2*R2 = R(p2,:)', C=Q2'*B, overwrite R with R2. Keep Q for NT 166 if nargout <= 3 167 [Q2,R,C,p2,info_spqr2] = ... 168 spqr_wrapper (R', B, 0, 'discard Q', get_details) ; 169 else 170 [Q2,R,C,p2,info_spqr2] = ... 171 spqr_wrapper (R', B, 0, 'keep Q', get_details) ; 172 end 173else 174 % Compute Q2*R2 = R(p2,:)', do not compute C, keep Q2, overwrite R with R2. 175 [Q2,R,ignore,p2,info_spqr2] = spqr_wrapper (R', [ ], 0, 'keep Q', ... 176 get_details) ; %#ok 177 clear ignore 178end 179 180if (get_details == 1) 181 stats.info_spqr2 = info_spqr2 ; 182end 183 184%------------------------------------------------------------------------------- 185% check if the numerical rank is consistent between the two QR factorizations 186%------------------------------------------------------------------------------- 187 188if rank_spqr ~= size (R,1) 189 % Approximate rank from two sparse QR factorizations are inconsistent. 190 % This should "never" happen. We know of no matrix that triggers this 191 % condition, and so the following line of code is untested. 192 error ('spqr_rank:inconsistent', 'inconsistent rank estimates') ; % untested 193 % rather than returning an error, we could do the following instead, 194 % but the code would be still untestable: 195 % warning ('spqr_rank:inconsistent', 'inconsistent rank estimates') ; 196 % [stats x N NT] = spqr_failure (5, stats, get_details, start_tic) ; 197 % return 198end 199 200% R is now square and has dimension rank_spqr 201 202%------------------------------------------------------------------------------- 203% use spqr_ssi to check and adjust numerical rank from spqr 204%------------------------------------------------------------------------------- 205 206R11 = R ; 207if get_details == 0 208 % opts2 used to include a few extra statistics in stats_ssi 209 opts2 = opts; 210 opts2.get_details = 2; 211else 212 opts2 = opts; 213end 214[U,S,V,stats_ssi] = spqr_ssi (R11, opts2) ; 215 216if (get_details == 1 || get_details == 2) 217 stats.stats_ssi = stats_ssi ; 218end 219if (get_details == 1) 220 stats.time_svd = stats_ssi.time_svd ; 221end 222 223 224%------------------------------------------------------------------------------- 225% check for early return 226%------------------------------------------------------------------------------- 227 228if stats_ssi.flag == 4 229 % overflow occurred during the inverse power method in spqr_ssi 230 [stats x N NT] = spqr_failure (4, stats, get_details, start_tic) ; 231 return 232end 233 234%------------------------------------------------------------------------------- 235% Estimate lower bounds on the singular values of A 236%------------------------------------------------------------------------------- 237 238% In equation (2) the Frobenius norm of E1*P1' is equal to norm_E_fro = 239% info_spqr1.norm_E_fro and therefore || E1*P1' || <= norm_E_fro. 240% By the pertubation theorem for singular values, for i = 1, 2, ..., 241% rank_spqr, singular value i of A differs at most by norm_E_fro 242% from singular value i of R. The routine spqr_ssi returns estimates of the 243% singular values of R in S and stats_ssi.est_error_bounds contains 244% estimates of error bounds on the entries in S. Therefore, 245% for i = 1:k, where S is k by k, estimated lower bounds on singular 246% values number (rank_spqr - k + i) of A are in est_sval_lower_bounds: 247% 248est_sval_lower_bounds = ... 249 max (diag(S)' - stats_ssi.est_error_bounds - norm_E_fro, 0) ; 250 251% lower bounds on the remaining singular values of A are zero 252est_sval_lower_bounds (length(S)+1:length(S)+min(m,n)-rank_spqr) = 0 ; 253 254numerical_rank = stats_ssi.rank ; 255 256% limit nsvals_small and nsvals_large due to number of singular values 257% available and calculated by spqr_ssi 258nsvals_small = min ([nsvals_small, min(m,n) - numerical_rank]) ; 259 260nsvals_large = min (nsvals_large, rank_spqr) ; 261nsvals_large = min ([nsvals_large, numerical_rank, ... 262 numerical_rank - rank_spqr + stats_ssi.ssi_max_block_used]) ; 263 264% return nsvals_large + nsvals_small of the estimates 265est_sval_lower_bounds = est_sval_lower_bounds (1:nsvals_large+nsvals_small) ; 266 267%------------------------------------------------------------------------------- 268% Estimate upper bounds on the singular values of A 269%------------------------------------------------------------------------------- 270 271% Again, by the pertubation theorem for singular values, for i = 1, 2,..., 272% rank_spqr, singular value i of A differs at most by norm_E_fro 273% from singular value i of R. The routine spqr_ssi returns estimates of the 274% singular values of R in S and stats_ssi.est_error_bounds contains 275% estimates of error bounds on the entries in S. Therefore, 276% for i = 1:k, where S is k by k, estimated lower bounds on singular 277% values number (rank_spqr - k + i) of A are in est_sval_upper_bounds: 278est_sval_upper_bounds = diag(S)' + stats_ssi.est_error_bounds + norm_E_fro; 279 280% upper bounds on the remaining singular values of A are norm_E_fro 281est_sval_upper_bounds(length(S)+1:length(S)+min(m,n)-rank_spqr) = norm_E_fro; 282 283% return nsvals_large + nsvals_small components of the estimates 284est_sval_upper_bounds = est_sval_upper_bounds(1:nsvals_large+nsvals_small); 285 286%------------------------------------------------------------------------------- 287% calculate orthonormal basis for null space of A 288%------------------------------------------------------------------------------- 289 290% always construct null space basis for A since this can produce better 291% estimated upper bounds on the singular values and should not require 292% significantly more work or memory 293 294call_from = 2; 295[N, stats, stats_ssp_N, est_sval_upper_bounds] = spqr_rank_form_basis(... 296 call_from, A, U, V ,Q1, rank_spqr, numerical_rank, stats, opts, ... 297 est_sval_upper_bounds, nsvals_small, nsvals_large, p1, Q2, p2) ; 298 299%------------------------------------------------------------------------------- 300% if requested, form null space basis of A' 301%------------------------------------------------------------------------------- 302 303if nargout == 4 304 305 call_from = 3; 306 [NT, stats, stats_ssp_NT, est_sval_upper_bounds] = spqr_rank_form_basis(... 307 call_from, A, U, V ,Q1, rank_spqr, numerical_rank, stats, opts, ... 308 est_sval_upper_bounds, nsvals_small, nsvals_large, p1, Q2, p2) ; 309 310end 311 312%------------------------------------------------------------------------------- 313% find psuedoinverse solution to (1): 314%------------------------------------------------------------------------------- 315 316call_from = 2 ; 317x = spqr_rank_deflation(call_from, R, U, V, C, m, n, rank_spqr, ... 318 numerical_rank, nsvals_large, opts, p1, p2, N, Q1, Q2) ; 319 320%------------------------------------------------------------------------------- 321% determine flag which indicates accuracy of the estimated numerical rank 322% and return statistics 323%------------------------------------------------------------------------------- 324 325call_from = 2; 326if nargout < 4 327 stats_ssp_NT = [ ]; 328end 329stats = spqr_rank_assign_stats(... 330 call_from, est_sval_upper_bounds, est_sval_lower_bounds, tol, ... 331 numerical_rank, nsvals_small, nsvals_large, stats, ... 332 stats_ssi, opts, nargout, stats_ssp_N, stats_ssp_NT, start_tic) ; 333 334