octave-lyh: scripts/sparse/svds.m annotate

annotate scripts/sparse/svds.m @ 8846:492801cd65eb

Fix small error in previous commit, which was two changes that should not have been combined; the first change was to improve breakpoints in debugging and was not mentioned in the commit message. Ugh.

author	John W. Eaton <jwe@octave.org>
date	Mon, 23 Feb 2009 16:37:40 -0500
parents	5a6db6bd1a02
children	eb63fbe60fab

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
8507 cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8417 diff changeset	83 persistent root2 = sqrt (2);
cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8417 diff changeset	84
8417 654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	85 if (nargin < 1 \|\| nargin > 4)
8507 cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8417 diff changeset	86 print_usage ();
8417 654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	87 endif
654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	88
654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	89 if (nargin < 4)
8507 cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8417 diff changeset	90 opts.tol = 1e-10 / root2;
8417 654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	91 opts.disp = 0;
654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	92 opts.maxit = 300;
654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	93 else
8507 cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8417 diff changeset	94 if (!isstruct (opts))
cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8417 diff changeset	95 error ("svds: opts must be a structure");
8417 654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	96 endif
8507 cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8417 diff changeset	97 if (!isfield (opts, "tol"))
cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8417 diff changeset	98 opts.tol = 1e-10 / root2;
8417 654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	99 endif
654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	100 endif
654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	101
8507 cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8417 diff changeset	102 if (nargin < 3 \|\| strcmp (sigma, "L"))
cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8417 diff changeset	103 if (isreal (a))
8417 654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	104 sigma = "LA";
654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	105 else
654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	106 sigma = "LR";
654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	107 endif
8507 cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8417 diff changeset	108 elseif (isscalar (sigma) && isreal (sigma))
cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8417 diff changeset	109 if (sigma < 0)
cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8417 diff changeset	110 error ("svds: sigma must be a positive real value");
8417 654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	111 endif
654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	112 else
8507 cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8417 diff changeset	113 error ("svds: sigma must be a positive real value or the string 'L'");
8417 654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	114 endif
654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	115
8507 cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8417 diff changeset	116 max_a = max (abs (a(:)));
cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8417 diff changeset	117 if (max_a == 0)
cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8417 diff changeset	118 u = eye (m, k);
cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8417 diff changeset	119 s = zeros (k, k);
cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8417 diff changeset	120 v = eye (n, k);
8417 654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	121 else
8507 cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8417 diff changeset	122 [m, n] = size (a);
8417 654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	123 if (nargin < 2)
8507 cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8417 diff changeset	124 k = min ([6, m, n]);
8417 654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	125 else
8507 cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8417 diff changeset	126 k = min ([k, m, n]);
8417 654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	127 endif
654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	128
654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	129 ## Scale everything by the 1-norm to make things more stable.
8507 cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8417 diff changeset	130 b = a / max_a;
cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8417 diff changeset	131 b_opts = opts;
cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8417 diff changeset	132 b_opts.tol = opts.tol / max_a;
cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8417 diff changeset	133 b_sigma = sigma;
cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8417 diff changeset	134 if (!ischar (b_sigma))
cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8417 diff changeset	135 b_sigma = b_sigma / max_a;
8417 654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	136 endif
654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	137
8507 cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8417 diff changeset	138 if (!ischar (b_sigma) && b_sigma == 0)
8417 654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	139 ## 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	140 ## 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	141 ## 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	142 ## 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	143 ## k for all scalar value of sigma...
8507 cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8417 diff changeset	144 [V, s, flag] = eigs ([sparse(m,m), b; b', sparse(n,n)],
cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8417 diff changeset	145 2 * k, b_sigma, b_opts);
8417 654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	146 else
8507 cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8417 diff changeset	147 [V, s, flag] = eigs ([sparse(m,m), b; b', sparse(n,n)],
cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8417 diff changeset	148 k, b_sigma, b_opts);
8417 654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	149 endif
8507 cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8417 diff changeset	150 s = diag (s);
8417 654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	151
8507 cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8417 diff changeset	152 if (ischar (sigma))
cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8417 diff changeset	153 norma = max (s);
8417 654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	154 else
8507 cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8417 diff changeset	155 norma = normest (a);
8417 654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	156 endif
8507 cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8417 diff changeset	157 V = root2 * V;
8417 654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	158 u = V(1:m,:);
654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	159 v = V(m+1:end,:);
654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	160
654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	161 ## 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	162 ## 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	163 ## 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	164 ## 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	165 ## 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	166 ## 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	167 ## 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	168 ## values. What is appropriate for the tolerance?
654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	169 tol = norma * opts.tol;
654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	170 ind = find(s > tol);
8507 cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8417 diff changeset	171 if (length (ind) < k)
8417 654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	172 ## 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	173 ## nominally negative.
8507 cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8417 diff changeset	174 zind = find (abs (s) <= tol);
cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8417 diff changeset	175 p = min (length (zind), k - length (ind));
cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8417 diff changeset	176 ind = [ind; zind(1:p)];
cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8417 diff changeset	177 elseif (length (ind) > k)
8417 654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	178 ind = ind(1:k);
654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	179 endif
654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	180 u = u(:,ind);
654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	181 s = s(ind);
654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	182 v = v(:,ind);
654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	183
8507 cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8417 diff changeset	184 if (length (s) < k)
cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8417 diff changeset	185 warning ("returning fewer singular values than requested");
cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8417 diff changeset	186 if (!ischar (sigma))
cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8417 diff changeset	187 warning ("try increasing the value of sigma");
8417 654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	188 endif
654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	189 endif
654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	190
8507 cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8417 diff changeset	191 s = s * max_a;
8417 654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	192 endif
654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	193
654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	194 if (nargout < 2)
654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	195 u = s;
654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	196 else
654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	197 s = diag(s);
654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	198 if (nargout > 3)
8507 cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8417 diff changeset	199 flag = norm (av - us, 1) > root2 * opts.tol * norm (a, 1);
8417 654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	200 endif
654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	201 endif
654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	202 endfunction
654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	203
654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	204 %!shared n, k, a, u, s, v, opts
654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	205 %! n = 100;
654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	206 %! k = 7;
654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	207 %! a = sparse([3:n,1:n,1:(n-2)],[1:(n-2),1:n,3:n],[ones(1,n-2),0.4nones(1,n),ones(1,n-2)]);
654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	208 %! %%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	209 %! [u,s,v] = svd(full(a));
654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	210 %! s = diag(s);
654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	211 %! [dum, idx] = sort(abs(s));
654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	212 %! s = s(idx);
654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	213 %! u = u(:,idx);
654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	214 %! v = v(:,idx);
8846 492801cd65eb Fix small error in previous commit, which was two changes that should not have been combined; the first change was to improve breakpoints in debugging and was not mentioned in the commit message. Ugh. John W. Eaton <jwe@octave.org> parents: 8845 diff changeset	215 %! randn('state',42)
8585 e6497be3f3d6 Skip tests if ARPACK is missing. Ben Abbott <bpabbott@mac.com> parents: 8507 diff changeset	216 %!testif HAVE_ARPACK
8417 654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	217 %! [u2,s2,v2,flag] = svds(a,k);
654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	218 %! s2 = diag(s2);
654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	219 %! assert(flag,!1);
654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	220 %! assert(s(end:-1:end-k+1), s2, 1e-10);
8585 e6497be3f3d6 Skip tests if ARPACK is missing. Ben Abbott <bpabbott@mac.com> parents: 8507 diff changeset	221 %!testif HAVE_ARPACK
8417 654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	222 %! [u2,s2,v2,flag] = svds(a,k,0);
654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	223 %! s2 = diag(s2);
654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	224 %! assert(flag,!1);
654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	225 %! assert(s(k:-1:1), s2, 1e-10);
8585 e6497be3f3d6 Skip tests if ARPACK is missing. Ben Abbott <bpabbott@mac.com> parents: 8507 diff changeset	226 %!testif HAVE_ARPACK
8417 654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	227 %! idx = floor(n/2);
654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	228 %! % 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	229 %! sigma = 0.99s(idx) + 0.01s(idx+1);
654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	230 %! [u2,s2,v2,flag] = svds(a,k,sigma);
654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	231 %! s2 = diag(s2);
654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	232 %! assert(flag,!1);
654bcfb937bf Add the eigs and svds functions David Bateman <dbateman@free.fr> parents: diff changeset	233 %! assert(s((idx+floor(k/2)):-1:(idx-floor(k/2))), s2, 1e-10);

Mercurial > hg > octave-lyh

annotate scripts/sparse/svds.m @ 8846:492801cd65eb