octave_packages/m/sparse/svds.m

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