rcond, condest, and a block 1-norm estimator

Jason Riedy <[email protected]> Mon, 19 Nov 2007 17:56:25 -0800
Newsgroups gmane.comp.gnu.octave.sources
Organization CS Div, EECS Dept, Univ. of California, Berkeley
Message-ID <[email protected]>
--=-=-=

Appended are the rcond(), condest(), and block_onenorm_est()
routines I've been using. I don't think they exactly match the
MATLAB(tm) routines, but I also don't care. ;) The files also
are living at
  http://www.cs.berkeley.edu/~ejr/dev/linalg/
for now.

They have been moderately beaten around, but there still may be
bugs lurking. I'm not 100% sure I'm handling complex matrices
correctly, but I also haven't run into any noticable problems.

Improvements and bug fixes gleefully accepted.

Jason


--=-=-=
Content-Type: text/plain; charset=utf-8
Content-Disposition: inline; filename=rcond.m
Content-Description: rcond.m
Content-Transfer-Encoding: quoted-printable
X-MIME-Autoconverted: from 8bit to quoted-printable by mail.cae.wisc.edu id lALKneKa029180

## Copyright (C) 2007, Regents of the University of California
## -*- mode: octave; -*-
##
## This program is free software distributed under the "modified" or
## 3-clause BSD license appended to this file.

function [est, v, w] =3D rcond (varargin)
  ## -*- texinfo -*-
  ## @deftypefn {Function File} {[@var{est}, @var{v}] =3D } rcond (@var{A=
}, @var{t} =3D min (size (A, 1), 5))=20
  ## @deftypefnx {Function File} {[@var{est}, @var{v}] =3D } rcond (@var{=
A}, @var{SOLVE}, @var{SOVLE_T}, @var{t} =3D min (n, 5))
  ## @deftypefnx {Function File} {[@var{est}, @var{v}] =3D } rcond (@var{=
APPLY}, @var{APPLY_T}, @var{SOLVE}, @var{SOVLE_T}, @var{n}, @var{t} =3D m=
in (n, 5))
  ##
  ## Estimate the inverse of the 1-norm condition number of a matrix
  ## matrix @var{A} using @var{t} test vectors pusing the randomized
  ## 1-norm estimator in block_onenorm_est.  If the matrix is not
  ## explicit, e.g. when estimating the condition number of @var{A}
  ## given an LU factorization, rcond uses the following functions:
  ##
  ## @table @var
  ## @item APPLY
  ## @code{A*x} for a matrix @code{x} of size @var{n} by @var{t}.
  ## @item APPLY_T
  ## @code{A'*x} for a matrix @code{x} of size @var{n} by @var{t}.
  ## @item SOLVE
  ## @code{A \ b} for a matrix @code{b} of size @var{n} by @var{t}.
  ## @item SOLVE_T
  ## @code{A' \ b} for a matrix @code{b} of size @var{n} by @var{t}.
  ## @end table
  ##
  ## The implicit version requires an explicit dimension @var{n}.
  ##
  ## @code{rcond} uses a randomized algorithm to approximate
  ## the 1-norms.
  ##
  ## @code{rcond} returns the inverse of the 1-norm condition estimate
  ## @var{est} and a vector @var{v} satisfying @code{norm
  ## (@var{A}*@var{v}, 1) =3D=3D norm (@var{A}, 1) * norm (@var{v}, 1) *
  ## @var{est}}. When @var{est} is large, @var{v} is an approximate null
  ## vector.
  ##
  ## References:=20
  ## @itemize
  ## @item Nicholas J. Higham and Fran=C3=A7oise Tisseur, "A Block Algori=
thm
  ## for Matrix 1-Norm Estimation, with an Application to 1-Norm
  ## Pseudospectra." SIMAX vol 21, no 4, pp 1185-1201.
  ## http://dx.doi.org/10.1137/S0895479899356080
  ## @item Nicholas J. Higham and Fran=C3=A7oise Tisseur, "A Block Algori=
thm
  ## for Matrix 1-Norm Estimation, with an Application to 1-Norm
  ## Pseudospectra." http://citeseer.ist.psu.edu/223007.html
  ## @end itemize
  ##
  ## @seealso{block_onenorm_est, condest, norm, cond}
  ##
  ## @end deftypefn

  ## Author: Jason Riedy <[email protected]>
  ## Keywords: linear-algebra norm estimation
  ## Version: 0.2

%!demo
%!  N =3D 100;
%!  A =3D randn (N) + eye (N);
%!  rcond (A)
%!  [L,U,P] =3D lu (A);
%!  rcond (A, @(x) U\ (L\ (P*x)), @(x) P'*(L'\ (U'\x)))
%!  rcond (@(x) A*x, @(x) A'*x, @(x) U\ (L\ (P*x)), @(x) P'*(L'\ (U'\x)),=
 N)
%!  1 / (norm (inv (A), 1) * norm (A, 1))

  if size (varargin, 2) < 1 || size (varargin, 2) > 5,
    usage("rcond: Incorrect arguments.");
  endif

  DEFAULT_T =3D 5;
  ITMAX =3D 10;

  if ismatrix (varargin{1}),
    n =3D size (varargin{1}, 1);
    if n !=3D size (varargin{1}, 2),
      usage("Matrix must be square.");
    endif
    A =3D varargin{1};

    if size (varargin, 2) > 1,
      if isscalar (varargin{2}),
	t =3D varargin{2};
      else
	if size (varargin, 2) < 3,
	  usage("Must supply both SOLVE and SOLVE_T.");
	else
	  SOLVE =3D varargin{2};
	  SOLVE_T =3D varargin{3};
	  if size (varargin, 2) > 3,
	    t =3D varargin{4};
	  endif
	endif
      endif
    endif
  else
    if size (varargin, 2) < 5,
      usage("Implicit form of rcond requires at least 5 arguments.");
    endif
    APPLY =3D varargin{1};
    APPLY_T =3D varargin{2};
    SOLVE =3D varargin{3};
    SOLVE_T =3D varargin{4};
    n =3D varargin{5};
    if !isscalar (n),
      usage("Dimension argument of implicit form must be scalar.");
    endif
    if size (varargin, 2) > 5,
      t =3D varargin{6};
    endif
  endif

  if !exist ("t", "var"),
    t =3D min (n, DEFAULT_T);
  endif

  if !exist ("SOLVE", "var"),
    if issparse (A),
      colord =3D colamd(A);
      Pc =3D speye(n);
      Pc(:,colord) =3D Pc;
      [L,U,P] =3D lu(A(:,Pc));
      SOLVE =3D @(x) Pc' * (U\ (L\ (P*x)));
      SOLVE_T =3D @(x) P'*(L'\ (U'\ (Pc*x)));
    else
      [L,U,P] =3D lu(A);
      SOLVE =3D @(x) U\ (L\ (P*x));
      SOLVE_T =3D @(x) P' * (L'\ (U'\x));
    endif
  endif

  if exist ("A", "var"),
    Anorm =3D norm(A, 1);
  else
    Anorm =3D block_onenorm_est(APPLY, APPLY_T, n, t);
  endif

  [Ainv_norm, v, w] =3D block_onenorm_est(SOLVE, SOLVE_T, n, t);

  est =3D 1/Anorm;
  est /=3D Ainv_norm;
  v =3D w / norm (w, 1);

endfunction

## Yes, these test bounds are really loose.  There's
## enough randomization to trigger odd cases with hilb().

%!test
%!  N =3D 6;
%!  A =3D hilb (N);
%!  cA =3D 1 / rcond (A);
%!  cA_test =3D norm (inv (A), 1) * norm (A, 1);
%!  assert (cA, cA_test, 2**-12);

%!test
%!  N =3D 6;
%!  A =3D hilb (N);
%!  SOLVE =3D @(x) A\x; SOLVE_T =3D @(x) A'\x;
%!  cA =3D 1 / rcond (A, SOLVE, SOLVE_T);
%!  cA_test =3D norm (inv (A), 1) * norm (A, 1);
%!  assert (cA, cA_test, 2**-12);

%!test
%!  N =3D 6;
%!  A =3D hilb (N);
%!  APPLY =3D @(x) A*x; APPLY_T =3D @(x) A'*x;
%!  SOLVE =3D @(x) A\x; SOLVE_T =3D @(x) A'\x;
%!  cA =3D 1 / rcond (APPLY, APPLY_T, SOLVE, SOLVE_T, N);
%!  cA_test =3D norm (inv (A), 1) * norm (A, 1);
%!  assert (cA, cA_test, 2**-6);

%!test
%!  N =3D 12;
%!  A =3D hilb (N);
%!  [rcondA, v] =3D rcond (A);
%!  x =3D A*v;
%!  assert (norm(x, inf), 0, eps);

## Copyright (c) 2007, Regents of the University of California
## All rights reserved.
## Redistribution and use in source and binary forms, with or without
## modification, are permitted provided that the following conditions are=
 met:
##
##     * Redistributions of source code must retain the above copyright
##       notice, this list of conditions and the following disclaimer.
##     * Redistributions in binary form must reproduce the above copyrigh=
t
##       notice, this list of conditions and the following disclaimer in =
the
##       documentation and/or other materials provided with the distribut=
ion.
##     * Neither the name of the University of California, Berkeley nor t=
he
##       names of its contributors may be used to endorse or promote prod=
ucts
##       derived from this software without specific prior written permis=
sion.
##
## THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AN=
D ANY
## EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPL=
IED
## WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
## DISCLAIMED. IN NO EVENT SHALL THE REGENTS AND CONTRIBUTORS BE LIABLE F=
OR ANY
## DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAM=
AGES
## (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SER=
VICES;
## LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSE=
D AND
## ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR =
TORT
## (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE =
OF THIS
## SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

--=-=-=
Content-Type: text/plain; charset=utf-8
Content-Disposition: inline; filename=block_onenorm_est.m
Content-Description: block_onenorm_est.m
Content-Transfer-Encoding: quoted-printable
X-MIME-Autoconverted: from 8bit to quoted-printable by mail.cae.wisc.edu id lALKneKa029180

## Copyright (C) 2007, Regents of the University of California
## -*- mode: octave; -*-
##
## This program is free software distributed under the "modified" or
## 3-clause BSD license appended to this file.

function [est, v, w, iter] =3D block_onenorm_est (varargin)
  ## -*- texinfo -*-
  ## @deftypefn {Function File} {[@var{est}, @var{v}, @var{w}, @var{iter}=
] =3D } block_onenorm_est (@var{A}, @var{t} =3D min (size (A, 1), 5))=20
  ## @deftypefnx {Function File} {[@var{est}, @var{v}, @var{w}, @var{iter=
}] =3D } block_onenorm_est (@var{APPLY}, @var{APPLY_T}, @var{n}, @var{t} =
=3D min (n, 5))
  ##
  ## Apply Higham and Tisseur's randomized block 1-norm estimator to
  ## matrix @var{A} using @var{t} test vectors.  If the matrix is not
  ## explicit, e.g. when estimating the norm of @code{inv(@var{A})} given=
 an
  ## LU factorization, block_onenorm_est applies @var{A} and its conjugat=
e
  ## transpose through a pair of functions @var{APPLY} and @var{APPLY_T},
  ## respectively, to a dense matrix of size @var{n} by @var{t}. The
  ## implicit version requires an explicit dimension @var{n}.
  ##
  ## Returns the norm estimate @var{est}, two vectors @var{v} and
  ## @var{w} related by norm
  ## @code{(@var{w}, 1) =3D @var{est} * norm (@var{v}, 1)},
  ## and the number of iterations @var{iter}.  The number of
  ## iterations is limited to 10 and is at least 2.
  ##
  ## References:=20
  ## @itemize
  ## @item Nicholas J. Higham and Fran=C3=A7oise Tisseur, "A Block Algori=
thm
  ## for Matrix 1-Norm Estimation, with an Application to 1-Norm
  ## Pseudospectra." SIMAX vol 21, no 4, pp 1185-1201.
  ## http://dx.doi.org/10.1137/S0895479899356080
  ## @item Nicholas J. Higham and Fran=C3=A7oise Tisseur, "A Block Algori=
thm
  ## for Matrix 1-Norm Estimation, with an Application to 1-Norm
  ## Pseudospectra." http://citeseer.ist.psu.edu/223007.html
  ## @end itemize
  ##
  ## @seealso{condest, norm, cond}
  ##
  ## @end deftypefn

  ## Author: Jason Riedy <[email protected]>
  ## Keywords: linear-algebra norm estimation
  ## Version: 0.2

%!demo
%!  N =3D 100;
%!  A =3D randn(N) + eye(N);
%!  [L,U,P] =3D lu(A);
%!  nm1inv =3D block_onenorm_est(@(x) U\(L\(P*x)), @(x) P'*(L'\(U'\x)), \
%!                             N, 30)
%!  norm(inv(A), 1)

  if size (varargin, 2) < 1 || size (varargin, 2) > 4,
    print_usage();
  endif

  DEFAULT_T =3D 5;
  ITMAX =3D 10;

  if ismatrix (varargin{1}),
    n =3D size (varargin{1}, 1);
    if n !=3D size (varargin{1}, 2),
      error("Matrix must be square.");
    endif
    APPLY =3D @(x) varargin{1} * x;
    APPLY_T =3D @(x) varargin{1}' * x;
    if size (varargin) > 1,
      t =3D varargin{2};
    else
      t =3D min (n, DEFAULT_T);
    endif
  else
    if size (varargin, 2) < 3,
      print_usage();
    endif
    n =3D varargin{3};
    APPLY =3D varargin{1};
    APPLY_T =3D varargin{2};
    if size (varargin) > 3,
      t =3D varargin{4};
    else
      t =3D DEFAULT_T;
    endif
  endif

  ## Initial test vectors X.
  X =3D rand (n, t);
  X =3D X ./ (ones (n,1) * sum (abs (X), 1));

  been_there =3D zeros (n, 1); # Track if a vertex has been visited.
  est_old =3D 0; # To check if the estimate has increased.
  S =3D zeros (n, t); # Normalized vector of signs.  The normalization is=
=20

  for iter=3D1:ITMAX+1,
    Y =3D feval (APPLY, X);

    ## Find the initial estimate as the largest A*x.
    [est, ind_best] =3D max (sum (abs (Y), 1));
    if (est > est_old || iter =3D=3D 2),
      w =3D Y(:,ind_best);
    endif
    if (iter >=3D 2 && est < est_old),
      ## No improvement, so stop.
      est =3D est_old;
      break;
    endif

    est_old =3D est;
    S_old =3D S;
    if (iter > ITMAX),
      ## Gone too far.  Stop.
      break;
    endif

    S =3D sign (Y);

    ## Test if any of S are approximately parallel to previous S
    ## vectors or current S vectors.  If everything is parallel,
    ## stop. Otherwise, replace any parallel vectors with
    ## rand{-1,+1}.
    partest =3D any (abs (S_old' * S - n) < 4*eps*n);
    if all (partest),
      ## All the current vectors are parallel to old vectors.
      ## We've hit a cycle, so stop.
      break;
    endif
    if any (partest),
      ## Some vectors are parallel to old ones and are cycling,
      ## but not all of them.  Replace the parallel vectors with
      ## rand{-1,+1}.
      numpar =3D sum (partest);
      replacements =3D 2*(rand (n,numpar) < 0.5) - 1;
      S(:,partest) =3D replacements;
    endif
    ## Now test for parallel vectors within S.
    partest =3D any ( (S' * S - eye (t)) =3D=3D n );
    if any (partest),
      numpar =3D sum (partest);
      replacements =3D 2*(rand (n,numpar) < 0.5) - 1;
      S(:,partest) =3D replacements;
    endif
   =20
    Z =3D feval (APPLY_T, S);

    ## Now find the largest non-previously-visted index per
    ## vector.
    h =3D max (abs (Z),2);
    [mh, mhi] =3D max (h);
    if iter >=3D 2 && mhi =3D=3D ind_best,
      ## Hit a cycle, stop.
      break;
    endif
    [h, ind] =3D sort (h, 'descend');
    if t > 1,
      firstind =3D ind(1:t);
      if all (been_there(firstind)),
	## Visited all these before, so stop.
	break;
      endif
      ind=3Dind(!been_there(ind));
      if length (ind) < t,
	## There aren't enough new vectors, so we're practically
	## in a cycle. Stop.
	break;
      endif
    endif

    ## Visit the new indices.
    X =3D zeros (n, t);
    for zz =3D 1:t,
      X(ind(zz),zz) =3D 1;
    endfor
    been_there(ind(1:t)) =3D 1;

  endfor

  ## The estimate est and vector w are set in the loop above. The
  ## vector v selects the ind_best column of A.
  v =3D zeros (n, 1);
  v(ind_best) =3D 1;
endfunction

%!test
%!  N =3D 10;
%!  A =3D ones (N);
%!  [nm1, v1, w1] =3D block_onenorm_est (A);
%!  [nminf, vinf, winf] =3D block_onenorm_est (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 =3D 10;
%!  A =3D ones (N);
%!  [nm1, v1, w1] =3D block_onenorm_est (@(x) A*x, @(x) A'*x, N, 3);
%!  [nminf, vinf, winf] =3D block_onenorm_est (@(x) A'*x, @(x) A*x, N, 3)=
;
%!  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 =3D 5;
%!  A =3D hilb (N);
%!  [nm1, v1, w1] =3D block_onenorm_est (A);
%!  [nminf, vinf, winf] =3D block_onenorm_est (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
%!  N =3D 100;
%!  A =3D rand (N);
%!  [nm1, v1, w1] =3D block_onenorm_est (A);
%!  [nminf, vinf, winf] =3D block_onenorm_est (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)

## Copyright (c) 2007, Regents of the University of California
## All rights reserved.
## Redistribution and use in source and binary forms, with or without
## modification, are permitted provided that the following conditions are=
 met:
##
##     * Redistributions of source code must retain the above copyright
##       notice, this list of conditions and the following disclaimer.
##     * Redistributions in binary form must reproduce the above copyrigh=
t
##       notice, this list of conditions and the following disclaimer in =
the
##       documentation and/or other materials provided with the distribut=
ion.
##     * Neither the name of the University of California, Berkeley nor t=
he
##       names of its contributors may be used to endorse or promote prod=
ucts
##       derived from this software without specific prior written permis=
sion.
##
## THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AN=
D ANY
## EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPL=
IED
## WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
## DISCLAIMED. IN NO EVENT SHALL THE REGENTS AND CONTRIBUTORS BE LIABLE F=
OR ANY
## DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAM=
AGES
## (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SER=
VICES;
## LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSE=
D AND
## ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR =
TORT
## (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE =
OF THIS
## SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

--=-=-=
Content-Type: text/plain; charset=utf-8
Content-Disposition: inline; filename=condest.m
Content-Description: condest.m
Content-Transfer-Encoding: quoted-printable
X-MIME-Autoconverted: from 8bit to quoted-printable by mail.cae.wisc.edu id lALKneKa029180

## Copyright (C) 2007, Regents of the University of California
## -*- mode: octave; -*-
##
## This program is free software distributed under the "modified" or
## 3-clause BSD license appended to this file.

function [est, v] =3D condest (varargin)
  ## -*- texinfo -*-
  ## @deftypefn {Function File} {[@var{est}, @var{v}] =3D } condest (@var=
{A}, @var{t} =3D min (size (A, 1), 5))=20
  ## @deftypefnx {Function File} {[@var{est}, @var{v}] =3D } condest (@va=
r{A}, @var{SOLVE}, @var{SOVLE_T}, @var{t} =3D min (n, 5))
  ## @deftypefnx {Function File} {[@var{est}, @var{v}] =3D } condest (@va=
r{APPLY}, @var{APPLY_T}, @var{SOLVE}, @var{SOVLE_T}, @var{n}, @var{t} =3D=
 min (n, 5))
  ##
  ## Estimate the 1-norm condition number of a matrix matrix @var{A}
  ## using @var{t} test vectors pusing the randomized 1-norm estimator in
  ## block_onenorm_est.  If the matrix is not explicit, e.g. when
  ## estimating the condition number of @var{A} given an LU
  ## factorization, condest uses the following functions:
  ##
  ## @table @var
  ## @item APPLY
  ## @code{A*x} for a matrix @code{x} of size @var{n} by @var{t}.
  ## @item APPLY_T
  ## @code{A'*x} for a matrix @code{x} of size @var{n} by @var{t}.
  ## @item SOLVE
  ## @code{A \ b} for a matrix @code{b} of size @var{n} by @var{t}.
  ## @item SOLVE_T
  ## @code{A' \ b} for a matrix @code{b} of size @var{n} by @var{t}.
  ## @end table
  ##
  ## The implicit version requires an explicit dimension @var{n}.
  ##
  ## @code{condest} uses a randomized algorithm to approximate
  ## the 1-norms.
  ##
  ## @code{condest} returns the 1-norm condition estimate @var{est} and
  ## a vector @var{v} satisfying @code{norm (@var{A}*@var{v}, 1) =3D=3D n=
orm
  ## (@var{A}, 1) * norm (@var{v}, 1) / @var{est}}. When @var{est} is
  ## large, @var{v} is an approximate null vector.
  ##
  ## References:=20
  ## @itemize
  ## @item Nicholas J. Higham and Fran=C3=A7oise Tisseur, "A Block Algori=
thm
  ## for Matrix 1-Norm Estimation, with an Application to 1-Norm
  ## Pseudospectra." SIMAX vol 21, no 4, pp 1185-1201.
  ## http://dx.doi.org/10.1137/S0895479899356080
  ## @item Nicholas J. Higham and Fran=C3=A7oise Tisseur, "A Block Algori=
thm
  ## for Matrix 1-Norm Estimation, with an Application to 1-Norm
  ## Pseudospectra." http://citeseer.ist.psu.edu/223007.html
  ## @end itemize
  ##
  ## @seealso{rcond, block_onenorm_est, norm, cond}
  ##
  ## @end deftypefn

  ## Author: Jason Riedy <[email protected]>
  ## Keywords: linear-algebra norm estimation
  ## Version: 0.2

%!demo
%!  N =3D 100;
%!  A =3D randn (N) + eye (N);
%!  condest (A)
%!  [L,U,P] =3D lu (A);
%!  condest (A, @(x) U\ (L\ (P*x)), @(x) P'*(L'\ (U'\x)))
%!  condest (@(x) A*x, @(x) A'*x, @(x) U\ (L\ (P*x)), @(x) P'*(L'\ (U'\x)=
), N)
%!  norm (inv (A), 1) * norm (A, 1)

  [est, v] =3D rcond(varargin{:});
  est =3D 1/est;
endfunction

## These tests are the same as in rcond().

%!test
%!  N =3D 6;
%!  A =3D hilb (N);
%!  cA =3D condest (A);
%!  rcA =3D 1 / rcond (A);
%!  assert (cA, rcA, 2**-12);

## Copyright (c) 2007, Regents of the University of California
## All rights reserved.
## Redistribution and use in source and binary forms, with or without
## modification, are permitted provided that the following conditions are=
 met:
##
##     * Redistributions of source code must retain the above copyright
##       notice, this list of conditions and the following disclaimer.
##     * Redistributions in binary form must reproduce the above copyrigh=
t
##       notice, this list of conditions and the following disclaimer in =
the
##       documentation and/or other materials provided with the distribut=
ion.
##     * Neither the name of the University of California, Berkeley nor t=
he
##       names of its contributors may be used to endorse or promote prod=
ucts
##       derived from this software without specific prior written permis=
sion.
##
## THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AN=
D ANY
## EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPL=
IED
## WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
## DISCLAIMED. IN NO EVENT SHALL THE REGENTS AND CONTRIBUTORS BE LIABLE F=
OR ANY
## DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAM=
AGES
## (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SER=
VICES;
## LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSE=
D AND
## ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR =
TORT
## (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE =
OF THIS
## SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

--=-=-=
Content-Type: text/plain; charset="us-ascii"
MIME-Version: 1.0
Content-Disposition: inline
Content-Transfer-Encoding: 7bit

_______________________________________________
Octave-sources mailing list
[email protected]
https://www.cae.wisc.edu/mailman/listinfo/octave-sources

--=-=-=--