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