Mercurial > octave
view scripts/linear-algebra/normest1.m @ 32060:1203a2d81a42
Add BIST tests for Matlab compatibility for functions which accept negative dimensions.
* data.cc (Fones, Finf, Feps, Feye): Add BIST test for negative dimensions.
* rand.cc (Frand): Add BIST test for negative dimensions.
* ov-cell.cc (Fcell): Add BIST test for negative dimensions.
author | Rik <rik@octave.org> |
---|---|
date | Wed, 26 Apr 2023 16:26:32 -0700 |
parents | 597f3ee61a48 |
children | 2e484f9f1f18 |
line wrap: on
line source
######################################################################## ## ## Copyright (C) 2016-2023 The Octave Project Developers ## ## See the file COPYRIGHT.md in the top-level directory of this ## distribution or <https://octave.org/copyright/>. ## ## 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 ## <https://www.gnu.org/licenses/>. ## ######################################################################## ## -*- texinfo -*- ## @deftypefn {} {@var{nest} =} normest1 (@var{A}) ## @deftypefnx {} {@var{nest} =} normest1 (@var{A}, @var{t}) ## @deftypefnx {} {@var{nest} =} normest1 (@var{A}, @var{t}, @var{x0}) ## @deftypefnx {} {@var{nest} =} normest1 (@var{Afcn}, @var{t}, @var{x0}, @var{p1}, @var{p2}, @dots{}) ## @deftypefnx {} {[@var{nest}, @var{v}] =} normest1 (@var{A}, @dots{}) ## @deftypefnx {} {[@var{nest}, @var{v}, @var{w}] =} normest1 (@var{A}, @dots{}) ## @deftypefnx {} {[@var{nest}, @var{v}, @var{w}, @var{iter}] =} normest1 (@var{A}, @dots{}) ## Estimate the 1-norm of the matrix @var{A} using a block algorithm. ## ## @code{normest1} is best for large sparse matrices where only an estimate of ## the norm is required. For small to medium sized matrices, consider using ## @code{norm (@var{A}, 1)}. In addition, @code{normest1} 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{Afcn} (@var{flag}, @var{x})} is used; it must 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{Afcn} (@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}. ## ## On output, @var{nest} is the desired estimate, @var{v} and @var{w} ## are vectors such that @code{@var{w} = @var{A} * @var{v}}, with ## @code{norm (@var{w}, 1)} = @code{@var{c} * norm (@var{v}, 1)}. @var{iter} ## contains in @code{@var{iter}(1)} the number of iterations (the maximum is ## hardcoded to 5) and in @code{@var{iter}(2)} the total number of products ## @code{@var{A} * @var{x}} or @code{@var{A}' * @var{x}} performed by the ## algorithm. ## ## Algorithm Note: @code{normest1} uses random numbers during evaluation. ## Therefore, if consistent results are required, the @qcode{"state"} of the ## random generator should be fixed before invoking @code{normest1}. ## ## 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, norm, cond, condest} ## @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 [nest, v, w, iter] = normest1 (A, t = [], x0 = [], 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)) Aismat = true; Aisreal = isreal (A); n = rows (A); if (n <= 4 || t == n) ## small input, compute directly [nest, idx] = max (sum (abs (A), 1), [] , 2); v = zeros (n, 1); v(idx) = 1; w = A(:, idx); ## Matlab incompatible on purpose. Matlab returns iter 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. iter = [0; 1]; return; endif elseif (is_function_handle (A)) Aismat = false; Aisreal = A ("real", [], varargin{:}); n = A ("dim", [], varargin{:}); Afcn = @(x) A ("notransp", x, varargin{:}); A1fcn = @(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 (x0)) 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 column 1:i-1. Change it. X(:,i) = sign (2 * rand (n, 1) - 1); else i++; endif endwhile X /= n; else if (columns (x0) < t) error ("normest1: X0 must have %d columns", t); endif X = x0; endif itmax = 5; idx_hist = zeros (n, 1); nest_old = 0; idx = zeros (n, 1); S = zeros (n, t); iter = [0; 0]; converged = false; while (! converged && (iter(1) < itmax)) iter(1)++; if (Aismat) Y = A * X; else Y = Afcn (X); endif iter(2)++; [nest, j] = max (sum (abs (Y), 1), [], 2); if ((nest > nest_old) || (iter(1) == 2)) idx_best = idx(j); w = Y(:, j); # there is an error in Algorithm 2.4 endif if (nest <= nest_old && iter(1) >= 2) # (1) of Algorithm 2.4 nest = nest_old; break; # while endif nest_old = nest; Sold = S; S = sign (Y); S(S==0) = 1; possible_break = false; if (Aisreal) ## 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; converged = true; else if (t > 1) ## at least two columns of S are not parallel i = 1; ## The maximum number of unparallel columns of length n with ## entries {-1,1} is 2^(n-1). n of them are already in Sold. imax = min (t, 2 ^ (n - 1) - n); while (i <= imax) if (any (abs (S(:,i)' * S(:,1:i-1)) == n) || any (abs (S(:,i)' * Sold) == n)) ## i-th column of S is parallel to a previous column ## 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 (Aismat) Z = A' * S; else Z = A1fcn (S); # (3) of Algorithm 2.4 endif iter(2)++; h = max (abs (Z), [], 2); idx = (1:n)'; if (iter(1) >= 2 && (max (h, [], 1) == h(idx_best))) # (4) of Alg. 2.4 break; # while endif [h, idx] = sort (h, "descend"); if (t > 1) if (all (idx_hist(idx(1:t)))) # (5) of Algorithm 2.4 break; # while endif idx = idx(! idx_hist(idx)); ## length(idx) 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 (idx), t); idx = idx(1:tmax); else tmax = 1; endif X = zeros (n, tmax); X(sub2ind (size (X), idx(1:tmax), (1:tmax)')) = 1; idx_hist(idx(1:tmax)) = 1; endif endwhile v = zeros (n, 1); v(idx_best) = 1; endfunction %!function z = afcn_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 = afcn_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 Afcn %!test %! A = rand (10); %! n = length (A); %! Afcn = @(flag, x) afcn_A (flag, x, A, n); %! assert (normest1 (Afcn), norm (A, 1), 2 * eps); ## test Afcn with parameters %!test %! A = rand (10); %! assert (normest1 (@afcn_A_P, [], [], A, 3), norm (A ^ 3, 1), 1000 * eps); ## test output %!test %! A = reshape (1:16,4,4); %! [nest, v, w, iter] = normest1 (A); %! assert (norm (w, 1), nest * 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; %! [nest, v, w, iter] = normest1 (A, 10); %! assert (w, A * v, eps); %!test %! A = rand (5); %! nest = normest1 (A, 6); %! assert (nest, norm (A,1), eps); %!test %! A = rand (5); %! nest = normest1 (A, 2, ones (5, 2) / 5); %! assert (nest, 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 ITER is always a column vector. %!test %! [~, ~, ~, it] = normest1 (rand (3), 3); %! assert (iscolumn (it)); %! [~, ~, ~, it] = normest1 (rand (50), 20); %! assert (iscolumn (it)); ## Test input validation %!error <Invalid call> normest1 () %!error <A must be a square matrix or a function handle> normest1 ({1}) %!error <A must be a square matrix> normest1 ([1 2]) %!error <X0 must have 2 columns> normest1 (magic (5), :, [1])