1######################################################################## 2## 3## Copyright (C) 2006-2021 The Octave Project Developers 4## 5## See the file COPYRIGHT.md in the top-level directory of this 6## distribution or <https://octave.org/copyright/>. 7## 8## This file is part of Octave. 9## 10## Octave is free software: you can redistribute it and/or modify it 11## under the terms of the GNU General Public License as published by 12## the Free Software Foundation, either version 3 of the License, or 13## (at your option) any later version. 14## 15## Octave is distributed in the hope that it will be useful, but 16## WITHOUT ANY WARRANTY; without even the implied warranty of 17## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 18## GNU General Public License for more details. 19## 20## You should have received a copy of the GNU General Public License 21## along with Octave; see the file COPYING. If not, see 22## <https://www.gnu.org/licenses/>. 23## 24######################################################################## 25 26## -*- texinfo -*- 27## @deftypefn {} {@var{s} =} svds (@var{A}) 28## @deftypefnx {} {@var{s} =} svds (@var{A}, @var{k}) 29## @deftypefnx {} {@var{s} =} svds (@var{A}, @var{k}, @var{sigma}) 30## @deftypefnx {} {@var{s} =} svds (@var{A}, @var{k}, @var{sigma}, @var{opts}) 31## @deftypefnx {} {[@var{u}, @var{s}, @var{v}] =} svds (@dots{}) 32## @deftypefnx {} {[@var{u}, @var{s}, @var{v}, @var{flag}] =} svds (@dots{}) 33## 34## Find a few singular values of the matrix @var{A}. 35## 36## The singular values are calculated using 37## 38## @example 39## @group 40## [@var{m}, @var{n}] = size (@var{A}); 41## @var{s} = eigs ([sparse(@var{m}, @var{m}), @var{A}; 42## @var{A}', sparse(@var{n}, @var{n})]) 43## @end group 44## @end example 45## 46## The eigenvalues returned by @code{eigs} correspond to the singular values 47## of @var{A}. The number of singular values to calculate is given by @var{k} 48## and defaults to 6. 49## 50## The argument @var{sigma} specifies which singular values to find. When 51## @var{sigma} is the string @qcode{'L'}, the default, the largest singular 52## values of @var{A} are found. Otherwise, @var{sigma} must be a real scalar 53## and the singular values closest to @var{sigma} are found. As a corollary, 54## @code{@var{sigma} = 0} finds the smallest singular values. Note that for 55## relatively small values of @var{sigma}, there is a chance that the 56## requested number of singular values will not be found. In that case 57## @var{sigma} should be increased. 58## 59## @var{opts} is a structure defining options that @code{svds} will pass 60## to @code{eigs}. The possible fields of this structure are documented in 61## @code{eigs}. By default, @code{svds} sets the following three fields: 62## 63## @table @code 64## @item tol 65## The required convergence tolerance for the singular values. The default 66## value is 1e-10. @code{eigs} is passed @code{@var{tol} / sqrt (2)}. 67## 68## @item maxit 69## The maximum number of iterations. The default is 300. 70## 71## @item disp 72## The level of diagnostic printout (0|1|2). If @code{disp} is 0 then 73## diagnostics are disabled. The default value is 0. 74## @end table 75## 76## If more than one output is requested then @code{svds} will return an 77## approximation of the singular value decomposition of @var{A} 78## 79## @example 80## @var{A}_approx = @var{u}*@var{s}*@var{v}' 81## @end example 82## 83## @noindent 84## where @var{A}_approx is a matrix of size @var{A} but only rank @var{k}. 85## 86## @var{flag} returns 0 if the algorithm has successfully converged, and 1 87## otherwise. The test for convergence is 88## 89## @example 90## @group 91## norm (@var{A}*@var{v} - @var{u}*@var{s}, 1) <= @var{tol} * norm (@var{A}, 1) 92## @end group 93## @end example 94## 95## @code{svds} is best for finding only a few singular values from a large 96## sparse matrix. Otherwise, @code{svd (full (@var{A}))} will likely be more 97## efficient. 98## @seealso{svd, eigs} 99## @end deftypefn 100 101function [u, s, v, flag] = svds (A, k, sigma, opts) 102 103 persistent root2 = sqrt (2); 104 105 if (nargin < 1 || nargin > 4) 106 print_usage (); 107 endif 108 109 if (ndims (A) > 2) 110 error ("svds: A must be a 2-D matrix"); 111 endif 112 113 if (nargin < 4) 114 opts.tol = 0; # use ARPACK default 115 opts.disp = 0; 116 opts.maxit = 300; 117 else 118 if (! isstruct (opts)) 119 error ("svds: OPTS must be a structure"); 120 endif 121 if (! isfield (opts, "tol")) 122 opts.tol = 0; # use ARPACK default 123 else 124 opts.tol = opts.tol / root2; 125 endif 126 if (isfield (opts, "v0")) 127 if (! isvector (opts.v0) || (length (opts.v0) != sum (size (A)))) 128 error ("svds: OPTS.v0 must be a vector with rows (A) + columns (A) entries"); 129 endif 130 endif 131 endif 132 133 if (nargin < 3 || strcmp (sigma, "L")) 134 if (isreal (A)) 135 sigma = "LA"; 136 else 137 sigma = "LR"; 138 endif 139 elseif (isscalar (sigma) && isnumeric (sigma) && isreal (sigma)) 140 if (sigma < 0) 141 error ("svds: SIGMA must be a positive real value"); 142 endif 143 else 144 error ("svds: SIGMA must be a positive real value or the string 'L'"); 145 endif 146 147 [m, n] = size (A); 148 max_a = max (abs (nonzeros (A))); 149 if (isempty (max_a)) 150 max_a = 0; 151 endif 152 ## Must initialize variable value, otherwise it may appear to interpreter 153 ## that code is trying to call flag() colormap function. 154 flag = 0; 155 156 if (max_a == 0) 157 s = zeros (k, 1); # special case of zero matrix 158 else 159 if (nargin < 2) 160 k = min ([6, m, n]); 161 else 162 k = min ([k, m, n]); 163 endif 164 165 ## Scale everything by the 1-norm to make things more stable. 166 b = A / max_a; 167 b_opts = opts; 168 ## Call to eigs is always a symmetric matrix by construction 169 b_opts.issym = true; 170 b_sigma = sigma; 171 if (! ischar (b_sigma)) 172 b_sigma /= max_a; 173 endif 174 175 if (b_sigma == 0) 176 ## Find the smallest eigenvalues 177 ## The eigenvalues returns by eigs for sigma=0 are symmetric about 0. 178 ## As we are only interested in the positive eigenvalues, we have to 179 ## double k and then throw out the k negative eigenvalues. 180 ## Separately, if sigma is nonzero, but smaller than the smallest 181 ## singular value, ARPACK may not return k eigenvalues. However, as 182 ## computation scales with k we'd like to avoid doubling k for all 183 ## scalar values of sigma. 184 b_k = 2 * k; 185 else 186 b_k = k; # Normal case, find just the k largest eigenvalues 187 endif 188 189 if (nargout > 1) 190 [V, s, flag] = eigs ([sparse(m,m), b; b', sparse(n,n)], 191 b_k, b_sigma, b_opts); 192 s = diag (s); 193 else 194 s = eigs ([sparse(m,m), b; b', sparse(n,n)], b_k, b_sigma, b_opts); 195 endif 196 197 if (ischar (sigma)) 198 norma = max (s); 199 else 200 norma = normest (A); 201 endif 202 ## We wish to exclude all eigenvalues that are less than zero as these 203 ## are artifacts of the way the matrix passed to eigs is formed. There 204 ## is also the possibility that the value of sigma chosen is exactly 205 ## a singular value, and in that case we're dead!! So have to rely on 206 ## the warning from eigs. We exclude the singular values which are 207 ## less than or equal to zero to within some tolerance scaled by the 208 ## norm since if we don't we might end up with too many singular 209 ## values. 210 if (b_sigma == 0) 211 if (sum (s>0) < k) 212 ## It may happen that the number of positive s is less than k. 213 ## In this case, take -s (if s in an eigenvalue, so is -s), 214 ## flipped upside-down. 215 s = flipud (-s); 216 endif 217 endif 218 tol = norma * opts.tol; 219 ind = find (s > tol); 220 if (length (ind) < k) 221 ## Too few eigenvalues returned. Add in any zero eigenvalues of B, 222 ## including the nominally negative ones. 223 zind = find (abs (s) <= tol); 224 p = min (length (zind), k - length (ind)); 225 ind = [ind; zind(1:p)]; 226 elseif (length (ind) > k) 227 ## Too many eigenvalues returned. Select according to criterion. 228 if (b_sigma == 0) 229 ind = ind(end+1-k:end); # smallest eigenvalues 230 else 231 ind = ind(1:k); # largest eigenvalues 232 endif 233 endif 234 s = s(ind); 235 236 if (length (s) < k) 237 warning ("svds: returning fewer singular values than requested"); 238 if (! ischar (sigma)) 239 warning ("svds: try increasing the value of sigma"); 240 endif 241 endif 242 243 s *= max_a; 244 endif 245 246 if (nargout < 2) 247 u = s; 248 else 249 if (max_a == 0) 250 u = eye (m, k); 251 s = diag (s); 252 v = eye (n, k); 253 else 254 u = root2 * V(1:m,ind); 255 s = diag (s); 256 v = root2 * V(m+1:end,ind); 257 endif 258 259 if (nargout > 3) 260 flag = (flag != 0); 261 endif 262 endif 263 264endfunction 265 266 267%!shared n, k, A, u, s, v, opts, rand_state, randn_state, tol 268%! n = 100; 269%! k = 7; 270%! A = sparse ([3:n,1:n,1:(n-2)],[1:(n-2),1:n,3:n],[ones(1,n-2),0.4*n*ones(1,n),ones(1,n-2)]); 271%! [u,s,v] = svd (full (A)); 272%! s = diag (s); 273%! [~, idx] = sort (abs (s)); 274%! s = s(idx); 275%! u = u(:, idx); 276%! v = v(:, idx); 277%! rand_state = rand ("state"); 278%! rand ("state", 42); 279%! opts.v0 = rand (2*n,1); # Initialize eigs ARPACK starting vector 280%! # to guarantee reproducible results 281 282%!testif HAVE_ARPACK 283%! [u2,s2,v2,flag] = svds (A,k); 284%! s2 = diag (s2); 285%! assert (flag, ! 1); 286%! tol = 15 * eps * norm (s2, 1); 287%! assert (s2, s(end:-1:end-k+1), tol); 288 289%!testif HAVE_ARPACK, HAVE_UMFPACK 290%! [u2,s2,v2,flag] = svds (A,k,0,opts); 291%! s2 = diag (s2); 292%! assert (flag, ! 1); 293%! tol = 15 * eps * norm (s2, 1); 294%! assert (s2, s(k:-1:1), tol); 295 296%!testif HAVE_ARPACK, HAVE_UMFPACK 297%! idx = floor (n/2); 298%! % Don't put sigma right on a singular value or there are convergence issues 299%! sigma = 0.99*s(idx) + 0.01*s(idx+1); 300%! [u2,s2,v2,flag] = svds (A,k,sigma,opts); 301%! s2 = diag (s2); 302%! assert (flag, ! 1); 303%! tol = 15 * eps * norm (s2, 1); 304%! assert (s2, s((idx+floor (k/2)):-1:(idx-floor (k/2))), tol); 305 306%!testif HAVE_ARPACK 307%! [u2,s2,v2,flag] = svds (zeros (10), k); 308%! assert (u2, eye (10, k)); 309%! assert (s2, zeros (k)); 310%! assert (v2, eye (10, 7)); 311%! 312%!testif HAVE_ARPACK 313%! s = svds (speye (10)); 314%! assert (s, ones (6, 1), 8*eps); 315 316%!testif HAVE_ARPACK <57185> 317%! z = complex (ones (10), ones (10)); 318%! s = svds (z); 319%! assert (isreal (s)); 320 321%!test 322%! ## Restore random number generator seed at end of tests 323%! rand ("state", rand_state); 324