Mercurial > octave
view scripts/linear-algebra/normest1.m @ 22330:53e246fd8124
doc: Spellcheck documentation ahead of 4.2 release.
* aspell-octave.en.pws: Add new words to Octave dictionary.
* func.txi, xzip.cc, publish.m, normest1.m, version.m, light.m:
Spellcheck documentation ahead of 4.2 release.
author | Rik <rik@octave.org> |
---|---|
date | Wed, 17 Aug 2016 08:20:26 -0700 |
parents | 9fc91bb2aec3 |
children | 7b190a2f11cb |
line wrap: on
line source
## Copyright (C) 2016 Marco Caliari ## ## This file is part of Octave. ## ## Octave is free software; you can redistribute it and/or modify it ## under the terms of the GNU General Public License as published by ## the Free Software Foundation; either version 3 of the License, or (at ## your option) any later version. ## ## Octave is distributed in the hope that it will be useful, but ## WITHOUT ANY WARRANTY; without even the implied warranty of ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU ## General Public License for more details. ## ## You should have received a copy of the GNU General Public License ## along with Octave; see the file COPYING. If not, see ## <http://www.gnu.org/licenses/>. ## -*- texinfo -*- ## @deftypefn {} {@var{c} =} normest1 (@var{a}) ## @deftypefnx {} {@var{c} =} normest1 (@var{a}, @var{t}) ## @deftypefnx {} {@var{c} =} normest1 (@var{a}, @var{t}, @var{x0}) ## @deftypefnx {} {@var{c} =} normest1 (@var{afun}, @var{t}, @var{x0}, @var{p1}, @var{p2}, @dots{}) ## @deftypefnx {} {[@var{c}, @var{v}] =} normest1 (@var{a}, @dots{}) ## @deftypefnx {} {[@var{c}, @var{v}, @var{w}] =} normest1 (@var{a}, @dots{}) ## @deftypefnx {} {[@var{c}, @var{v}, @var{w}, @var{it}] =} normest1 (@var{a}, @dots{}) ## Estimate the 1-norm of the matrix @var{a} using a block algorithm. ## ## For a medium size matrix @var{a}, @code{norm (@var{a}, 1)} should be ## used instead. For a large sparse matrix, when only an estimate of the norm ## is needed, @code{normest1 (@var{a})} might be faster. Moreover, it can be ## used for the estimate of the 1-norm of a linear ## operator @var{a} when matrix-vector products @code{@var{a} * @var{x}} and ## @code{@var{a}' * @var{x}} can be cheaply computed. In this case, ## instead of the matrix @var{a}, a function ## @code{@var{afun} (@var{flag}, @var{x})} can be used. It should return: ## ## @itemize @bullet ## @item ## the dimension @var{n} of @var{a}, if @var{flag} is @qcode{"dim"} ## ## @item ## true if @var{a} is a real operator, if @var{flag} is @qcode{"real"} ## ## @item ## the result @code{@var{a} * @var{x}}, if @var{flag} is @qcode{"notransp"} ## ## @item ## the result @code{@var{a}' * @var{x}}, if @var{flag} is @qcode{"transp"} ## @end itemize ## ## A typical case is @var{a} defined by @code{@var{b} ^ @var{m}}, ## in which the result @code{@var{a} * @var{x}} can be computed without ## even forming explicitly @code{@var{b} ^ @var{m}} by: ## ## @example ## @group ## @var{y} = @var{x}; ## for @var{i} = 1:@var{m} ## @var{y} = @var{b} * @var{y}; ## endfor ## @end group ## @end example ## ## The parameters @var{p1}, @var{p2}, @dots{} are arguments of ## @code{@var{afun} (@var{flag}, @var{x}, @var{p1}, @var{p2}, @dots{})}. ## ## The default value for @var{t} is 2. The algorithm requires ## matrix-matrix products with sizes @var{n} x @var{n} and ## @var{n} x @var{t}. ## ## The initial matrix @var{x0} should have columns of unit 1-norm. ## The default initial matrix @var{x0} has the first column ## @code{ones (@var{n}, 1) / @var{n}} ## and, if @var{t} > 1, the remaining columns with random elements ## @code{-1 / @var{n}}, @code{1 / @var{n}}, divided by @var{n}. ## Therefore, if consistent results are required, the @qcode{"state"} of the ## random generator should be fixed before invoking @code{normest1}. ## ## On output, @var{c} is the desired estimate, @var{v} and @var{w} ## vectors such that @code{@var{w} = @var{a} * @var{v}}, with ## @code{norm (@var{w}, 1)} = @code{@var{c} * norm (@var{v}, 1)}. ## @var{it} contains in @code{@var{it}(1)} the number of iterations ## (the maximum number is hardcoded to 5) and in @code{@var{it}(2)} ## the total number of products @code{@var{a} * @var{x}} or ## @code{@var{a}' * @var{x}} performed by the algorithm. ## ## Reference: @nospell{N. J. Higham and F. Tisseur}, ## @cite{A block algorithm for matrix 1-norm estimation, with and ## application to 1-norm @nospell{pseudospectra}}, ## @nospell{SIAM J. Matrix Anal. Appl.}, ## pp. 1185--1201, Vol 21, No. 4, 2000. ## ## @seealso{normest, rand} ## @end deftypefn ## Ideally, we would set t and X to their default values but Matlab ## compatibility would require we set the default even when they are ## empty. function [est, v, w, k] = normest1 (A, t = [], X = [], varargin) if (nargin < 1) print_usage (); endif if (isempty (t)) t = 2; endif ## FIXME: t < 0 should print trace information if (isnumeric (A) && issquare (A)) Aisnum = true; n = rows (A); if ((n <= 4) || (t == n)) ## compute directly [est, idx] = max (sum (abs (A), 1), [] ,2); v = zeros (n, 1); v(idx) = 1; w = A(:, idx); ## Matlab incompatible on purpose. Matlab returns k as a row vector ## for this special case, but a column vector in all other cases. ## This is obviously a bug in Matlab that we don't reproduce. k = [0; 1]; return else realm = isreal (A); endif elseif (is_function_handle (A)) Aisnum = false; n = A ("dim", [], varargin{:}); realm = A ("real", [], varargin{:}); Afun = @(x) A ("notransp", x, varargin{:}); A1fun = @(x) A ("transp", x, varargin{:}); else error ("normest1: A must be a square matrix or a function handle"); endif t = min (t, n); if (isempty (X)) X = [ones(n, 1), sign(2 * rand(n, t - 1) - 1)]; i = 2; imax = min (t, 2 ^ (n - 1)); ## There are at most 2^(n-1) unparallel columns, see later. while (i <= imax) if (any (abs (X(:,i)' * X(:,1:i-1)) == n)) ## column i is parallel to a colum 1:i-1. Change it. X(:,i) = sign (2 * rand (n, 1) - 1); else i++; endif endwhile X /= n; elseif (columns (X) < t) error ("normest1: X should have %d columns", t); endif itmax = 5; ind_hist = zeros (n, 1); est_old = 0; ind = zeros (n, 1); S = zeros (n, t); k = [0; 0]; conv = false; while ((! conv) && (k(1) < itmax)) k(1)++; if (Aisnum) Y = A * X; else Y = Afun (X); endif k(2)++; [est, j] = max (sum (abs (Y), 1), [], 2); if ((est > est_old) || (k(1) == 2)) ind_best = ind(j); w = Y(:, j); # there is an error in Algorithm 2.4 endif if ((est <= est_old) && (k(1) >= 2)) # (1) of Algorithm 2.4 est = est_old; break; # while endif est_old = est; Sold = S; S = sign (Y); S(S==0) = 1; possible_break = false; if (realm) ## test parallel (only real case) if (all (any (abs (Sold' * S) == n))) # (2) of Algorithm 2.4 ## all columns of S parallel to a column of Sold, exit possible_break = true; conv = true; else if (t > 1) ## at least two columns of S are not parallel i = 1; imax = min (t, 2 ^ (n - 1)); while (i <= imax) ## The maximum number of parallel columns of length n with entries ## {-1,1} is 2^(n-1). Therefore, if the number of columns of S is ## greater than 2^(n-1), for sure some of them are parallel to some ## columns of Sold. Don't even try to change them (i <= 2^(n-1)). ## Now, check if S(:,i) is parallel to any previous column of S p = (any (abs (S(:,i)' * S(:,1:i-1)) == n)); if (p || (any (abs (S(:,i)' * Sold) == n))) ## i-th column of S parallel to a previous or to a ## column of Sold: change it. S(:,i) = sign (2*rand (n, 1)-1); else i++; endif endwhile endif endif endif if (! possible_break) if (Aisnum) Z = A' * S; else Z = A1fun (S); # (3) of Algorithm 2.4 endif k(2)++; h = max (abs (Z), [], 2); ind = (1:n)'; if ((k(1) >= 2) && (max (h, [], 1) == h(ind_best))) # (4) of Alg. 2.4 break; # while endif [h, ind] = sort (h, "descend"); if (t > 1) if (all (ind_hist(ind(1:t)))) # (5) of Algorithm 2.4 break; # while endif ind = ind(! ind_hist(ind)); ## length(ind) could be less than t, especially if t is not << n. ## This is not considered in point (5) of Algorithm 2.4. tmax = min (numel (ind), t); ind = ind(1:tmax); else tmax = 1; endif X = zeros (n, tmax); X(sub2ind (size (X), ind(1:tmax), (1:tmax)')) = 1; ind_hist(ind(1:tmax)) = 1; endif endwhile v = zeros (n, 1); v(ind_best) = 1; endfunction %!function z = afun_A (flag, x, A, n) %! switch flag %! case {"dim"} %! z = n; %! case {"real"} %! z = isreal (A); %! case {"transp"} %! z = A' * x; %! case {"notransp"} %! z = A * x; %! endswitch %!endfunction %!function z = afun_A_P (flag, x, A, m) %! switch flag %! case "dim" %! z = length (A); %! case "real" %! z = isreal (A); %! case "transp" %! z = x; for i = 1:m, z = A' * z;, endfor %! case "notransp" %! z = x; for i = 1:m, z = A * z;, endfor %! endswitch %!endfunction %!test %! A = reshape ((1:16)-8, 4, 4); %! assert (normest1 (A), norm (A, 1), eps) ## test t=1 %!test %! A = rand (4); # for positive matrices always work %! assert (normest1 (A, 1), norm (A,1), 2 * eps) ## test t=3 %!test %! A = [-0.21825 0.16598 0.19388 0.75297 %! -1.47732 0.78443 -1.04254 0.42240 %! 1.39857 -0.34046 2.28617 0.68089 %! 0.31205 1.50529 -0.75804 -1.22476]; %! X = [1,1,-1;1,1,1;1,1,-1;1,-1,-1]/3; %! assert (normest1 (A, 3, X), norm (A, 1), 2 * eps) ## test afun %!test %! A = rand (10); %! n = length (A); %! afun = @(flag, x) afun_A (flag, x, A, n); %! assert (normest1 (afun), norm (A, 1), 2 * eps) ## test afun with parameters %!test %! A = rand (10); %! assert (normest1 (@afun_A_P, [], [], A, 3), norm (A ^ 3, 1), 1000 * eps) ## test output %!test %! A = reshape (1:16,4,4); %! [c, v, w, it] = normest1 (A); %! assert (norm (w, 1), c * norm (v, 1), eps) ## test output %!test %! A = rand (100); %! A(A <= 1/3) = -1; %! A(A > 1/3 & A <= 2/3) = 0; %! A(A > 2/3) = 1; %! [c, v, w, it] = normest1 (A, 10); %! assert (w, A * v, eps) %!test %! A = rand (5); %! c = normest1 (A, 6); %! assert (c, norm (A,1), eps) %!test %! A = rand (5); %! c = normest1 (A, 2, ones (5, 2) / 5); %! assert (c, norm (A,1), eps) %!test %! N = 10; %! A = ones (N); %! [nm1, v1, w1] = normest1 (A); %! [nminf, vinf, winf] = normest1 (A', 6); %! assert (nm1, N, -2*eps) %! assert (nminf, N, -2*eps) %! assert (norm (w1, 1), nm1 * norm (v1, 1), -2*eps) %! assert (norm (winf, 1), nminf * norm (vinf, 1), -2*eps) %!test %! N = 5; %! A = hilb (N); %! [nm1, v1, w1] = normest1 (A); %! [nminf, vinf, winf] = normest1 (A', 6); %! assert (nm1, norm (A, 1), -2*eps) %! assert (nminf, norm (A, inf), -2*eps) %! assert (norm (w1, 1), nm1 * norm (v1, 1), -2*eps) %! assert (norm (winf, 1), nminf * norm (vinf, 1), -2*eps) ## Only likely to be within a factor of 10. %!test %! old_state = rand ("state"); %! unwind_protect %! rand ("state", 42); # Initialize to guarantee reproducible results %! N = 100; %! A = rand (N); %! [nm1, v1, w1] = normest1 (A); %! [nminf, vinf, winf] = normest1 (A', 6); %! assert (nm1, norm (A, 1), -.1) %! assert (nminf, norm (A, inf), -.1) %! assert (norm (w1, 1), nm1 * norm (v1, 1), -2*eps) %! assert (norm (winf, 1), nminf * norm (vinf, 1), -2*eps) %! unwind_protect_cleanup %! rand ("state", old_state); %! end_unwind_protect ## Check IT is always a column vector. %!test %! [~, ~, ~, it] = normest1 (rand (3), 3); %! assert (iscolumn (it)) %! [~, ~, ~, it] = normest1 (rand (50), 20); %! assert (iscolumn (it))