Mercurial > octave-antonio

annotate scripts/sparse/svds.m @ 8417:654bcfb937bf

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add the eigs and svds functions
author David Bateman <dbateman@free.fr>
date Tue, 23 Dec 2008 08:28:23 +0100
parents
children cadc73247d65
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
8417
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
1 ## Copyright (C) 2006 David Bateman
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
2 ##
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
3 ## This program is free software; you can redistribute it and/or modify
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
4 ## it under the terms of the GNU General Public License as published by
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
5 ## the Free Software Foundation; either version 2 of the License, or
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
6 ## (at your option) any later version.
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
7 ##
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
8 ## This program is distributed in the hope that it will be useful,
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
9 ## but WITHOUT ANY WARRANTY; without even the implied warranty of
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
10 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
11 ## GNU General Public License for more details.
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
12 ##
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
13 ## You should have received a copy of the GNU General Public License
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
14 ## along with this program; If not, see <http://www.gnu.org/licenses/>.
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
15
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
16 ## -*- texinfo -*-
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
17 ## @deftypefn {Function File} {@var{s} =} svds (@var{a})
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
18 ## @deftypefnx {Function File} {@var{s} =} svds (@var{a}, @var{k})
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
19 ## @deftypefnx {Function File} {@var{s} =} svds (@var{a}, @var{k}, @var{sigma})
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
20 ## @deftypefnx {Function File} {@var{s} =} svds (@var{a}, @var{k}, @var{sigma}, @var{opts})
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
21 ## @deftypefnx {Function File} {[@var{u}, @var{s}, @var{v}, @var{flag}] =} svds (@dots{})
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
22 ##
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
23 ## Find a few singular values of the matrix @var{a}. The singular values
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
24 ## are calculated using
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
25 ##
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
26 ## @example
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
27 ## @group
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
28 ## [@var{m}, @var{n}] = size(@var{a})
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
29 ## @var{s} = eigs([sparse(@var{m}, @var{m}), @var{a}; @dots{}
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
30 ## @var{a}', sparse(@var{n}, @var{n})])
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
31 ## @end group
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
32 ## @end example
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
33 ##
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
34 ## The eigenvalues returned by @code{eigs} correspond to the singular
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
35 ## values of @var{a}. The number of singular values to calculate is given
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
36 ## by @var{k}, whose default value is 6.
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
37 ##
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
38 ## The argument @var{sigma} can be used to specify which singular values
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
39 ## to find. @var{sigma} can be either the string 'L', the default, in
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
40 ## which case the largest singular values of @var{a} are found. Otherwise
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
41 ## @var{sigma} should be a real scalar, in which case the singular values
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
42 ## closest to @var{sigma} are found. Note that for relatively small values
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
43 ## of @var{sigma}, there is the chance that the requested number of singular
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
44 ## values are not returned. In that case @var{sigma} should be increased.
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
45 ##
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
46 ## If @var{opts} is given, then it is a structure that defines options
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
47 ## that @code{svds} will pass to @var{eigs}. The possible fields of this
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
48 ## structure are therefore determined by @code{eigs}. By default three
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
49 ## fields of this structure are set by @code{svds}.
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
50 ##
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
51 ## @table @code
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
52 ## @item tol
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
53 ## The required convergence tolerance for the singular values. @code{eigs}
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
54 ## is passed @var{tol} divided by @code{sqrt(2)}. The default value is
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
55 ## 1e-10.
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
56 ##
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
57 ## @item maxit
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
58 ## The maximum number of iterations. The defaut is 300.
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
59 ##
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
60 ## @item disp
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
61 ## The level of diagnostic printout. If @code{disp} is 0 then there is no
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
62 ## printout. The default value is 0.
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
63 ## @end table
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
64 ##
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
65 ## If more than one output argument is given, then @code{svds} also
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
66 ## calculates the left and right singular vectors of @var{a}. @var{flag}
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
67 ## is used to signal the convergence of @code{svds}. If @code{svds}
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
68 ## converges to the desired tolerance, then @var{flag} given by
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
69 ##
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
70 ## @example
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
71 ## @group
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
72 ## norm (@var{a} * @var{v} - @var{u} * @var{s}, 1) <= @dots{}
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
73 ## @var{tol} * norm (@var{a}, 1)
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
74 ## @end group
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
75 ## @end example
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
76 ##
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
77 ## will be zero.
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
78 ## @end deftypefn
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
79 ## @seealso{eigs}
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
80
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
81 function [u, s, v, flag] = svds (a, k, sigma, opts)
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
82
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
83 if (nargin < 1 || nargin > 4)
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
84 error ("Incorrect number of arguments");
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
85 endif
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
86
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
87 if (nargin < 4)
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
88 opts.tol = 1e-10 / sqrt(2);
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
89 opts.disp = 0;
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
90 opts.maxit = 300;
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
91 else
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
92 if (!isstruct(opts))
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
93 error("opts must be a structure");
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
94 endif
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
95 if (!isfield(opts,"tol"))
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
96 opts.tol = 1e-10 / sqrt(2);
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
97 endif
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
98 endif
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
99
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
100 if (nargin < 3 || strcmp(sigma,"L"))
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
101 if (isreal(a))
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
102 sigma = "LA";
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
103 else
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
104 sigma = "LR";
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
105 endif
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
106 elseif (isscalar(sigma) && isreal(sigma))
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
107 if ((sigma < 0))
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
108 error ("sigma must be a positive real value");
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
109 endif
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
110 else
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
111 error ("sigma must be a positive real value or the string 'L'");
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
112 endif
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
113
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
114 maxA = max(max(abs(a)));
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
115 if (maxA == 0)
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
116 u = eye(m, k);
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
117 s = zeros(k, k);
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
118 v = eye(n, k);
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
119 else
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
120 [m, n] = size(a);
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
121 if (nargin < 2)
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
122 k = min([6, m, n]);
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
123 else
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
124 k = min([k, m, n]);
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
125 endif
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
126
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
127 ## Scale everything by the 1-norm to make things more stable.
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
128 B = a / maxA;
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
129 Bopts = opts;
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
130 Bopts.tol = opts.tol / maxA;
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
131 Bsigma = sigma;
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
132 if (!ischar(Bsigma))
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
133 Bsigma = Bsigma / maxA;
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
134 endif
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
135
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
136 if (!ischar(Bsigma) && Bsigma == 0)
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
137 ## The eigenvalues returns by eigs are symmetric about 0. As we
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
138 ## are only interested in the positive eigenvalues, we have to
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
139 ## double k. If sigma is smaller than the smallest singular value
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
140 ## this can also be an issue. However, we'd like to avoid double
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
141 ## k for all scalar value of sigma...
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
142 [V, s, flag] = eigs ([sparse(m,m), B; B', sparse(n,n)],
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
143 2 * k, Bsigma, Bopts);
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
144 else
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
145 [V, s, flag] = eigs ([sparse(m,m), B; B', sparse(n,n)],
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
146 k, Bsigma, Bopts);
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
147 endif
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
148 s = diag(s);
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
149
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
150 if (ischar(sigma))
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
151 norma = max(s);
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
152 else
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
153 norma = normest(a);
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
154 endif
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
155 V = sqrt(2) * V;
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
156 u = V(1:m,:);
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
157 v = V(m+1:end,:);
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
158
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
159 ## We wish to exclude all eigenvalues that are less than zero as these
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
160 ## are artifacts of the way the matrix passed to eigs is formed. There
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
161 ## is also the possibility that the value of sigma chosen is exactly
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
162 ## a singular value, and in that case we're dead!! So have to rely on
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
163 ## the warning from eigs. We exclude the singular values which are
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
164 ## less than or equal to zero to within some tolerance scaled by the
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
165 ## norm since if we don't we might end up with too many singular
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
166 ## values. What is appropriate for the tolerance?
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
167 tol = norma * opts.tol;
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
168 ind = find(s > tol);
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
169 if (length(ind) < k)
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
170 ## Find the zero eigenvalues of B, Ignore the eigenvalues that are
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
171 ## nominally negative.
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
172 zind = find(abs(s) <= tol);
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
173 p = min(length(zind), k-length(ind));
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
174 ind = [ind;zind(1:p)];
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
175 elseif (length(ind) > k)
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
176 ind = ind(1:k);
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
177 endif
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
178 u = u(:,ind);
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
179 s = s(ind);
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
180 v = v(:,ind);
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
181
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
182 if (length(s) < k)
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
183 warning("returning fewer singular values than requested.");
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
184 if (!ischar(sigma))
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
185 warning("try increasing the value of sigma");
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
186 endif
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
187 endif
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
188
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
189 s = s * maxA;
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
190 endif
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
191
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
192 if (nargout < 2)
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
193 u = s;
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
194 else
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
195 s = diag(s);
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
196 if (nargout > 3)
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
197 flag = norm(a*v - u*s, 1) > sqrt(2) * opts.tol * norm(a, 1);
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
198 endif
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
199 endif
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
200 endfunction
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
201
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
202 %!shared n, k, a, u, s, v, opts
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
203 %! n = 100;
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
204 %! k = 7;
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
205 %! 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)]);
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
206 %! %%a = sparse([3:n,1:n,1:(n-2)],[1:(n-2),1:n,3:n],[ones(1,n-2),1:n,-ones(1,n-2)]);
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
207 %! [u,s,v] = svd(full(a));
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
208 %! s = diag(s);
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
209 %! [dum, idx] = sort(abs(s));
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
210 %! s = s(idx);
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
211 %! u = u(:,idx);
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
212 %! v = v(:,idx);
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
213 %! randn('state',42)
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
214 %!test
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
215 %! [u2,s2,v2,flag] = svds(a,k);
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
216 %! s2 = diag(s2);
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
217 %! assert(flag,!1);
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
218 %! assert(s(end:-1:end-k+1), s2, 1e-10);
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
219 %!test
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
220 %! [u2,s2,v2,flag] = svds(a,k,0);
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
221 %! s2 = diag(s2);
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
222 %! assert(flag,!1);
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
223 %! assert(s(k:-1:1), s2, 1e-10);
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
224 %!test
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
225 %! idx = floor(n/2);
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
226 %! % Don't put sigma right on a singular value or there are convergence
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
227 %! sigma = 0.99*s(idx) + 0.01*s(idx+1);
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
228 %! [u2,s2,v2,flag] = svds(a,k,sigma);
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
229 %! s2 = diag(s2);
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
230 %! assert(flag,!1);
654bcfb937bf Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
231 %! assert(s((idx+floor(k/2)):-1:(idx-floor(k/2))), s2, 1e-10);