octave_packages/m/linear-algebra/condest.m

   1 ## Copyright (C) 2007-2012 Regents of the University of California
   2 ##
   3 ## This file is part of Octave.
   4 ##
   5 ## Octave is free software; you can redistribute it and/or modify it
   6 ## under the terms of the GNU General Public License as published by
   7 ## the Free Software Foundation; either version 3 of the License, or (at
   8 ## your option) any later version.
   9 ##
  10 ## Octave is distributed in the hope that it will be useful, but
  11 ## WITHOUT ANY WARRANTY; without even the implied warranty of
  12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13 ## General Public License for more details.
  14 ##
  15 ## You should have received a copy of the GNU General Public License
  16 ## along with Octave; see the file COPYING.  If not, see
  17 ## <http://www.gnu.org/licenses/>.
  18
  19 ## -*- texinfo -*-
  20 ## @deftypefn  {Function File} {} condest (@var{A})
  21 ## @deftypefnx {Function File} {} condest (@var{A}, @var{t})
  22 ## @deftypefnx {Function File} {[@var{est}, @var{v}] =} condest (@dots{})
  23 ## @deftypefnx {Function File} {[@var{est}, @var{v}] =} condest (@var{A}, @var{solve}, @var{solve_t}, @var{t})
  24 ## @deftypefnx {Function File} {[@var{est}, @var{v}] =} condest (@var{apply}, @var{apply_t}, @var{solve}, @var{solve_t}, @var{n}, @var{t})
  25 ##
  26 ## Estimate the 1-norm condition number of a matrix @var{A}
  27 ## using @var{t} test vectors using a randomized 1-norm estimator.
  28 ## If @var{t} exceeds 5, then only 5 test vectors are used.
  29 ##
  30 ## If the matrix is not explicit, e.g., when estimating the condition
  31 ## number of @var{A} given an LU@tie{}factorization, @code{condest} uses the
  32 ## following functions:
  33 ##
  34 ## @table @var
  35 ## @item apply
  36 ## @code{A*x} for a matrix @code{x} of size @var{n} by @var{t}.
  37 ##
  38 ## @item apply_t
  39 ## @code{A'*x} for a matrix @code{x} of size @var{n} by @var{t}.
  40 ##
  41 ## @item solve
  42 ## @code{A \ b} for a matrix @code{b} of size @var{n} by @var{t}.
  43 ##
  44 ## @item solve_t
  45 ## @code{A' \ b} for a matrix @code{b} of size @var{n} by @var{t}.
  46 ## @end table
  47 ##
  48 ## The implicit version requires an explicit dimension @var{n}.
  49 ##
  50 ## @code{condest} uses a randomized algorithm to approximate
  51 ## the 1-norms.
  52 ##
  53 ## @code{condest} returns the 1-norm condition estimate @var{est} and
  54 ## a vector @var{v} satisfying @code{norm (A*v, 1) == norm (A, 1) * norm
  55 ## (@var{v}, 1) / @var{est}}.  When @var{est} is large, @var{v} is an
  56 ## approximate null vector.
  57 ##
  58 ## References:
  59 ## @itemize
  60 ## @item
  61 ## N.J. Higham and F. Tisseur, @cite{A Block Algorithm
  62 ## for Matrix 1-Norm Estimation, with an Application to 1-Norm
  63 ## Pseudospectra}. SIMAX vol 21, no 4, pp 1185-1201.
  64 ## @url{http://dx.doi.org/10.1137/S0895479899356080}
  65 ##
  66 ## @item
  67 ## N.J. Higham and F. Tisseur, @cite{A Block Algorithm
  68 ## for Matrix 1-Norm Estimation, with an Application to 1-Norm
  69 ## Pseudospectra}. @url{http://citeseer.ist.psu.edu/223007.html}
  70 ## @end itemize
  71 ##
  72 ## @seealso{cond, norm, onenormest}
  73 ## @end deftypefn
  74
  75 ## Code originally licensed under
  76 ##
  77 ##  Copyright (c) 2007, Regents of the University of California
  78 ##  All rights reserved.
  79 ##
  80 ##  Redistribution and use in source and binary forms, with or without
  81 ##  modification, are permitted provided that the following conditions
  82 ##  are met:
  83 ##
  84 ##     * Redistributions of source code must retain the above copyright
  85 ##       notice, this list of conditions and the following disclaimer.
  86 ##
  87 ##     * Redistributions in binary form must reproduce the above
  88 ##       copyright notice, this list of conditions and the following
  89 ##       disclaimer in the documentation and/or other materials provided
  90 ##       with the distribution.
  91 ##
  92 ##     * Neither the name of the University of California, Berkeley nor
  93 ##       the names of its contributors may be used to endorse or promote
  94 ##       products derived from this software without specific prior
  95 ##       written permission.
  96 ##
  97 ##  THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS''
  98 ##  AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
  99 ##  TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
 100 ##  PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS AND
 101 ##  CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
 102 ##  SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
 103 ##  LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
 104 ##  USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
 105 ##  ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
 106 ##  OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
 107 ##  OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
 108 ##  SUCH DAMAGE.
 109
 110 ## Author: Jason Riedy <ejr@cs.berkeley.edu>
 111 ## Keywords: linear-algebra norm estimation
 112 ## Version: 0.2
 113
 114 function [est, v] = condest (varargin)
 115
 116   if (nargin < 1 || nargin > 6)
 117     print_usage ();
 118   endif
 119
 120   default_t = 5;
 121
 122   have_A = false;
 123   have_t = false;
 124   have_solve = false;
 125
 126   if (ismatrix (varargin{1}))
 127     A = varargin{1};
 128     if (! issquare (A))
 129       error ("condest: matrix must be square");
 130     endif
 131     n = rows (A);
 132     have_A = true;
 133
 134     if (nargin > 1)
 135       if (isscalar (varargin{2}))
 136         t = varargin{2};
 137         have_t = true;
 138       elseif (nargin > 2)
 139         solve = varargin{2};
 140         solve_t = varargin{3};
 141         have_solve = true;
 142         if (nargin > 3)
 143           t = varargin{4};
 144           have_t = true;
 145         endif
 146       else
 147         error ("condest: must supply both SOLVE and SOLVE_T");
 148       endif
 149     endif
 150   elseif (nargin > 4)
 151     apply = varargin{1};
 152     apply_t = varargin{2};
 153     solve = varargin{3};
 154     solve_t = varargin{4};
 155     have_solve = true;
 156     n = varargin{5};
 157     if (! isscalar (n))
 158       error ("condest: dimension argument of implicit form must be scalar");
 159     endif
 160     if (nargin > 5)
 161       t = varargin{6};
 162       have_t = true;
 163     endif
 164   else
 165     error ("condest: implicit form of condest requires at least 5 arguments");
 166   endif
 167
 168   if (! have_t)
 169     t = min (n, default_t);
 170   endif
 171
 172   if (! have_solve)
 173     if (issparse (A))
 174       [L, U, P, Pc] = lu (A);
 175       solve = @(x) Pc' * (U \ (L \ (P * x)));
 176       solve_t = @(x) P' * (L' \ (U' \ (Pc * x)));
 177     else
 178       [L, U, P] = lu (A);
 179       solve = @(x) U \ (L \ (P*x));
 180       solve_t = @(x) P' * (L' \ (U' \ x));
 181     endif
 182   endif
 183
 184   if (have_A)
 185     Anorm = norm (A, 1);
 186   else
 187     Anorm = onenormest (apply, apply_t, n, t);
 188   endif
 189
 190   [Ainv_norm, v, w] = onenormest (solve, solve_t, n, t);
 191
 192   est = Anorm * Ainv_norm;
 193   v = w / norm (w, 1);
 194
 195 endfunction
 196
 197 %!demo
 198 %!  N = 100;
 199 %!  A = randn (N) + eye (N);
 200 %!  condest (A)
 201 %!  [L,U,P] = lu (A);
 202 %!  condest (A, @(x) U\ (L\ (P*x)), @(x) P'*(L'\ (U'\x)))
 203 %!  condest (@(x) A*x, @(x) A'*x, @(x) U\ (L\ (P*x)), @(x) P'*(L'\ (U'\x)), N)
 204 %!  norm (inv (A), 1) * norm (A, 1)
 205
 206 ## Yes, these test bounds are really loose.  There's
 207 ## enough randomization to trigger odd cases with hilb().
 208
 209 %!test
 210 %!  N = 6;
 211 %!  A = hilb (N);
 212 %!  cA = condest (A);
 213 %!  cA_test = norm (inv (A), 1) * norm (A, 1);
 214 %!  assert (cA, cA_test, -2^-8);
 215
 216 %!test
 217 %!  N = 6;
 218 %!  A = hilb (N);
 219 %!  solve = @(x) A\x; solve_t = @(x) A'\x;
 220 %!  cA = condest (A, solve, solve_t);
 221 %!  cA_test = norm (inv (A), 1) * norm (A, 1);
 222 %!  assert (cA, cA_test, -2^-8);
 223
 224 %!test
 225 %!  N = 6;
 226 %!  A = hilb (N);
 227 %!  apply = @(x) A*x; apply_t = @(x) A'*x;
 228 %!  solve = @(x) A\x; solve_t = @(x) A'\x;
 229 %!  cA = condest (apply, apply_t, solve, solve_t, N);
 230 %!  cA_test = norm (inv (A), 1) * norm (A, 1);
 231 %!  assert (cA, cA_test, -2^-6);
 232
 233 %!test
 234 %!  N = 12;
 235 %!  A = hilb (N);
 236 %!  [rcondA, v] = condest (A);
 237 %!  x = A*v;
 238 %!  assert (norm(x, inf), 0, eps);