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