1function [x,stats,N,NT] = spqr_pinv (A, varargin) 2%SPQR_PINV approx pseudoinverse solution to min(norm(B-A*X)) 3% 4% usage: [x,stats,N,NT] = spqr_pinv (A,B,opts) ; 5% 6% This function returns an approximate psuedoinverse solution to 7% min || B - A x || (1) 8% for rank deficient matrices A. The psuedoinverse solution is the minimum 9% norm solution to the least squares problem (1). Also, optionally, the 10% routine returns statistics including the numerical rank of the matrix A for 11% tolerance tol (i.e. the number of singular values > tol) and other statistics 12% (see below), as well as, if requested, an orthonormal basis for the numerical 13% null space to A and an orthonormal basis for the null space of A transpose. 14% The psuedoinverse solution is approximate since the algorithm allows small 15% perturbations in A (columns of A may be changed by no more than a user 16% defined value in opts.tol). 17% 18% Input: A -- an m by n matrix 19% B -- an m by p right hand side matrix 20% opts (optional) -- see below 21% Output: x -- this n by p matrix contains psuedoinverse solutions to (1). 22% If B is empty then x will also be empty. 23% stats (optional) -- statistics including: 24% 25% Examples: 26% A = sparse (gallery('kahan',100)) ; 27% B = randn (100,1) ; B = B / norm(B) ; 28% x = spqr_pinv(A,B) ; 29% x_pinv = pinv(full(A))*B ; 30% rel_error_in_x = norm (x - x_pinv) / norm (x_pinv) 31% % or 32% [x,stats,N,NT] = spqr_pinv (A,B) ; 33% norm_A_times_N = norm (full(spqr_null_mult(N,A,3))) 34% norm_N_transpose_times_A = norm (full(spqr_null_mult(NT,A,0))) 35% % or 36% opts = struct('tol',1.e-5) ; 37% [x,stats] = spqr_pinv (A,B,opts) ; 38% stats 39% 40% See also spqr_basic, spqr_null, spqr_pinv, spqr_cod. 41 42% Copyright 2012, Leslie Foster and Timothy A Davis. 43 44% Algorithm: a basic solution is calculated using spqr_basic. Following 45% this an orthogonal basis, stored in N, for the numerical null space 46% of A is calculated using spqr_null. The psuedoinverse solution is 47% then x - N*(N'*x) which is calculated using the routine spqr_null_mult. 48 49%------------------------------------------------------------------------------- 50% set tolerance and number of singular values to estimate 51%------------------------------------------------------------------------------- 52 53[B,opts,stats,start_tic,ok] = spqr_rank_get_inputs (A, 1, varargin {:}) ; 54 55if (~ok || nargout > 4) 56 error ('usage: [x,stats,N,NT] = spqr_pinv (A,B,opts)') ; 57end 58 59% get the options 60get_details = opts.get_details ; 61 62if (get_details == 1) 63 stats.opts_used = opts ; 64end 65 66% set the order of the stats fields 67% stats.flag, stats.rank, stats.rank_spqr, stats.rank_spqr (if get_details 68% >= 1), stats.tol, stats.tol_alt, stats.normest_A (if calculated), 69% stats.est_sval_upper_bounds, stats.est_sval_lower_bounds, and 70% stats.sval_numbers_for_bounds already initialized in spqr_rank_get_inputs 71if nargout >= 3 72 stats.est_norm_A_times_N = -1 ; 73end 74if nargout == 4 75 stats.est_norm_A_transpose_times_NT = -1 ; 76end 77% order for the additional stats fields needed when get_details is 1 will be 78% set using spqr_rank_order_fields at the end of the routine 79 80%------------------------------------------------------------------------------- 81% find basic solution 82%------------------------------------------------------------------------------- 83 84if get_details == 0 85 % opts2 used to include a few extra statistics in stats_ssi 86 opts2 = opts; 87 opts2.get_details = 2; 88else 89 opts2 = opts; 90end 91if nargout == 4 92 % save basis for null space of A', if there are four output 93 % parameters. This will require more memory. 94 [x,stats_spqr_basic,NT] = spqr_basic (A, B, opts2) ; 95else 96 % do not save the basis for the null space of A'. Saves memory. 97 [x,stats_spqr_basic] = spqr_basic (A, B, opts2) ; 98end 99 100stats.rank = stats_spqr_basic.rank ; 101 102% save the stats 103if (get_details == 1 || get_details == 2) 104 stats.rank_spqr = stats_spqr_basic.rank_spqr ; 105end 106 107%------------------------------------------------------------------------------- 108% check for early return 109%------------------------------------------------------------------------------- 110 111if stats_spqr_basic.flag == 4 112 % overflow in spqr_ssi, called by spqr_basic. 113 % spqr_basic has already issued a warning regarding overflow 114 [stats x N NT] = spqr_failure (4, stats, get_details, start_tic) ; 115 return 116end 117 118%------------------------------------------------------------------------------- 119% calculate basis for the null space of A 120%------------------------------------------------------------------------------- 121 122% spqr_null calls spqr_ssi. Set the block size in spqr_ssi based on the 123% difference in the numerical rank calculated by spqr_basic and by spqr. Note 124% that this overwrites the user-provided value of opts.ssi_min_block, if any. 125opts.ssi_min_block = max (2, stats_spqr_basic.rank_spqr - stats.rank + 1); 126 127opts.ssi_min_block = ... 128 min (opts.ssi_min_block, stats_spqr_basic.stats_ssi.ssi_max_block_used) ; 129 130[N,stats_spqr_null] = spqr_null (A, opts2) ; 131 132if (get_details == 1) 133 stats.stats_spqr_basic = stats_spqr_basic ; 134 stats.stats_spqr_null = stats_spqr_null ; 135end 136 137%------------------------------------------------------------------------------- 138% check for early return 139%------------------------------------------------------------------------------- 140 141% spqr_basic and spqr_null both return error flags. Since both must be correct 142% to calculate the psuedoinverse solution, choose the largest (worst) error. 143 144stats.flag = max (stats_spqr_basic.flag, stats_spqr_null.flag) ; 145 146if stats_spqr_basic.flag == 0 && stats_spqr_null.flag == 0 && ... 147 stats_spqr_basic.rank ~= stats_spqr_null.rank 148 % Rank from spqr_basic and spqr_null differ. This is so rare that we know 149 % of no matrices that trigger this condition. This block of code is thus 150 % untested. 151 error ('spqr_rank:inconsistent', 'inconsistent rank estimates') ; % untested 152 % an alternative, which would cause a return in the code below. 153 % stats.flag = 5 ; 154end 155 156if stats.flag >= 4 157 % early return: overflow in inverse power method in ssi, 158 % or inconsistent rank estimates by spqr_basic and spqr_null. 159 [stats x N NT] = spqr_failure (stats.flag, stats, get_details, start_tic) ; 160 return 161end 162 163%------------------------------------------------------------------------------- 164% calculate the psuedoinverse solution 165%------------------------------------------------------------------------------- 166 167x = x - spqr_null_mult (N, spqr_null_mult (N,x,0), 1) ; 168 169%------------------------------------------------------------------------------- 170% select from two estimates of numerical rank and two sets of bounds 171%------------------------------------------------------------------------------- 172 173% Strategy -- the choice of stats.flag above was the 174% max (stats_spqr_basic.flag, stats_spqr_null.flag) 175% We will choose the bounds and rank corresponding to this choice of flag. 176% When the flags are the same we will use spqr_basic results for 177% the bounds, except when the flags are both 1 we will choose the bounds 178% corresponding to the maximum of stats_spqr_basic.tol_alt and 179% stats_spqr_null.tol_alt (the conservative choice). 180 181if stats_spqr_basic.flag == 1 && stats_spqr_null.flag == 1 182 if stats_spqr_basic.tol_alt >= stats_spqr_null.tol_alt 183 st = stats_spqr_basic ; 184 else 185 st = stats_spqr_null ; 186 end 187elseif stats_spqr_basic.flag >= stats_spqr_null.flag 188 st = stats_spqr_basic ; 189else 190 st = stats_spqr_null ; 191end 192stats.rank = st.rank ; 193if isfield(st,'tol_alt') 194 stats.tol_alt = st.tol_alt ; 195end 196stats.sval_numbers_for_bounds = st.sval_numbers_for_bounds ; 197stats.est_sval_upper_bounds = st.est_sval_upper_bounds ; 198stats.est_sval_lower_bounds = st.est_sval_lower_bounds ; 199 200%------------------------------------------------------------------------------- 201% return statistics 202%------------------------------------------------------------------------------- 203 204if (get_details == 1 || nargout >= 3 ) 205 stats.est_norm_A_times_N = stats_spqr_null.est_norm_A_times_N ; 206end 207if (get_details == 1) 208 stats.est_err_bound_norm_A_times_N = ... 209 stats_spqr_null.est_err_bound_norm_A_times_N ; 210end 211if nargout == 4 212 stats.est_norm_A_transpose_times_NT = ... 213 stats_spqr_basic.est_norm_A_transpose_times_NT ; 214 if (get_details == 1) 215 stats.est_err_bound_norm_A_transpose_times_NT = ... 216 stats_spqr_basic.est_err_bound_norm_A_transpose_times_NT ; 217 end 218end 219 220if (get_details == 1) 221 stats.time_svd = stats_spqr_basic.time_svd + ... 222 stats_spqr_basic.stats_ssi.time_svd + ... 223 stats_spqr_null.time_svd + ... 224 stats_spqr_null.stats_ssi.time_svd + ... 225 stats_spqr_null.stats_ssp_N.time_svd ; 226 if nargout == 4 227 stats.time_svd = stats.time_svd+stats_spqr_basic.stats_ssp_NT.time_svd ; 228 end 229 stats.time_basis = stats_spqr_basic.time_basis + stats_spqr_null.time_basis; 230end 231 232if stats.tol_alt == -1 233 stats = rmfield(stats, 'tol_alt') ; 234end 235 236if get_details == 1 237 % order the fields of stats in a convenient order (the fields when 238 % get_details is 0 or 2 are already in a good order) 239 stats = spqr_rank_order_fields(stats); 240end 241 242if (get_details == 1) 243 stats.time = toc (start_tic) ; 244end 245 246