rcond, condest, and a block 1-norm estimator

[Jason Riedy-2]

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


## 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] = rcond (varargin)
  ## -*- texinfo -*-
  ## @deftypefn {Function File} {[@var{est}, @var{v}] = } rcond (@var{A}, @var{t} = min (size (A, 1), 5))
  ## @deftypefnx {Function File} {[@var{est}, @var{v}] = } rcond (@var{A}, @var{SOLVE}, @var{SOVLE_T}, @var{t} = min (n, 5))
  ## @deftypefnx {Function File} {[@var{est}, @var{v}] = } rcond (@var{APPLY}, @var{APPLY_T}, @var{SOLVE}, @var{SOVLE_T}, @var{n}, @var{t} = min (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) == norm (@var{A}, 1) * norm (@var{v}, 1) *
  ## @var{est}}. When @var{est} is large, @var{v} is an approximate null
  ## vector.
  ##
  ## References:
  ## @itemize
  ## @item Nicholas J. Higham and Françoise Tisseur, "A Block Algorithm
  ## 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çoise Tisseur, "A Block Algorithm
  ## 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 <ejr@...>
  ## Keywords: linear-algebra norm estimation
  ## Version: 0.2

%!demo
%!  N = 100;
%!  A = randn (N) + eye (N);
%!  rcond (A)
%!  [L,U,P] = 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 = 5;
  ITMAX = 10;

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

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

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

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

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

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

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

endfunction

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

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

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

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

%!test
%!  N = 12;
%!  A = hilb (N);
%!  [rcondA, v] = rcond (A);
%!  x = 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 copyright
##       notice, this list of conditions and the following disclaimer in the
##       documentation and/or other materials provided with the distribution.
##     * Neither the name of the University of California, Berkeley nor the
##       names of its contributors may be used to endorse or promote products
##       derived from this software without specific prior written permission.
##
## THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND ANY
## EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
## WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
## DISCLAIMED. IN NO EVENT SHALL THE REGENTS AND CONTRIBUTORS BE LIABLE FOR ANY
## DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
## (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
## LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED 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.

## 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] = block_onenorm_est (varargin)
  ## -*- texinfo -*-
  ## @deftypefn {Function File} {[@var{est}, @var{v}, @var{w}, @var{iter}] = } block_onenorm_est (@var{A}, @var{t} = min (size (A, 1), 5))
  ## @deftypefnx {Function File} {[@var{est}, @var{v}, @var{w}, @var{iter}] = } block_onenorm_est (@var{APPLY}, @var{APPLY_T}, @var{n}, @var{t} = 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 conjugate
  ## 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) = @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:
  ## @itemize
  ## @item Nicholas J. Higham and Françoise Tisseur, "A Block Algorithm
  ## 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çoise Tisseur, "A Block Algorithm
  ## 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 <ejr@...>
  ## Keywords: linear-algebra norm estimation
  ## Version: 0.2

%!demo
%!  N = 100;
%!  A = randn(N) + eye(N);
%!  [L,U,P] = lu(A);
%!  nm1inv = 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 = 5;
  ITMAX = 10;

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

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

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

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

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

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

    S = 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 = 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 = sum (partest);
      replacements = 2*(rand (n,numpar) < 0.5) - 1;
      S(:,partest) = replacements;
    endif
    ## Now test for parallel vectors within S.
    partest = any ( (S' * S - eye (t)) == n );
    if any (partest),
      numpar = sum (partest);
      replacements = 2*(rand (n,numpar) < 0.5) - 1;
      S(:,partest) = replacements;
    endif
   
    Z = feval (APPLY_T, S);

    ## Now find the largest non-previously-visted index per
    ## vector.
    h = max (abs (Z),2);
    [mh, mhi] = max (h);
    if iter >= 2 && mhi == ind_best,
      ## Hit a cycle, stop.
      break;
    endif
    [h, ind] = sort (h, 'descend');
    if t > 1,
      firstind = ind(1:t);
      if all (been_there(firstind)),
        ## Visited all these before, so stop.
        break;
      endif
      ind=ind(!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 = zeros (n, t);
    for zz = 1:t,
      X(ind(zz),zz) = 1;
    endfor
    been_there(ind(1:t)) = 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 = zeros (n, 1);
  v(ind_best) = 1;
endfunction

%!test
%!  N = 10;
%!  A = ones (N);
%!  [nm1, v1, w1] = block_onenorm_est (A);
%!  [nminf, vinf, winf] = 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 = 10;
%!  A = ones (N);
%!  [nm1, v1, w1] = block_onenorm_est (@(x) A*x, @(x) A'*x, N, 3);
%!  [nminf, vinf, winf] = 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 = 5;
%!  A = hilb (N);
%!  [nm1, v1, w1] = block_onenorm_est (A);
%!  [nminf, vinf, winf] = 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 = 100;
%!  A = rand (N);
%!  [nm1, v1, w1] = block_onenorm_est (A);
%!  [nminf, vinf, winf] = 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 copyright
##       notice, this list of conditions and the following disclaimer in the
##       documentation and/or other materials provided with the distribution.
##     * Neither the name of the University of California, Berkeley nor the
##       names of its contributors may be used to endorse or promote products
##       derived from this software without specific prior written permission.
##
## THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND ANY
## EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
## WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
## DISCLAIMED. IN NO EVENT SHALL THE REGENTS AND CONTRIBUTORS BE LIABLE FOR ANY
## DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
## (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
## LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED 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.

## 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] = condest (varargin)
  ## -*- texinfo -*-
  ## @deftypefn {Function File} {[@var{est}, @var{v}] = } condest (@var{A}, @var{t} = min (size (A, 1), 5))
  ## @deftypefnx {Function File} {[@var{est}, @var{v}] = } condest (@var{A}, @var{SOLVE}, @var{SOVLE_T}, @var{t} = min (n, 5))
  ## @deftypefnx {Function File} {[@var{est}, @var{v}] = } condest (@var{APPLY}, @var{APPLY_T}, @var{SOLVE}, @var{SOVLE_T}, @var{n}, @var{t} = 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) == norm
  ## (@var{A}, 1) * norm (@var{v}, 1) / @var{est}}. When @var{est} is
  ## large, @var{v} is an approximate null vector.
  ##
  ## References:
  ## @itemize
  ## @item Nicholas J. Higham and Françoise Tisseur, "A Block Algorithm
  ## 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çoise Tisseur, "A Block Algorithm
  ## 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 <ejr@...>
  ## Keywords: linear-algebra norm estimation
  ## Version: 0.2

%!demo
%!  N = 100;
%!  A = randn (N) + eye (N);
%!  condest (A)
%!  [L,U,P] = 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] = rcond(varargin{:});
  est = 1/est;
endfunction

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

%!test
%!  N = 6;
%!  A = hilb (N);
%!  cA = condest (A);
%!  rcA = 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 copyright
##       notice, this list of conditions and the following disclaimer in the
##       documentation and/or other materials provided with the distribution.
##     * Neither the name of the University of California, Berkeley nor the
##       names of its contributors may be used to endorse or promote products
##       derived from this software without specific prior written permission.
##
## THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND ANY
## EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
## WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
## DISCLAIMED. IN NO EVENT SHALL THE REGENTS AND CONTRIBUTORS BE LIABLE FOR ANY
## DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
## (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
## LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED 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.

_______________________________________________
Octave-sources mailing list
Octave-sources@...
https://www.cae.wisc.edu/mailman/listinfo/octave-sources

[dbateman3]

I believe the rcond function should take rcond from the LAPACK
factorization, at least I believe is what matlab does. Matlab calls the
normest1 function, that is basically your block_onenorm_est, and then
multiples the estimated inverse norm by norm(A,1) to get the condition
number estimate.

condest is one of the annoying missing sparse functions in Octave and I
appreciate that your code treats sparse matrices as a special case. I'd
hesitate to include it though with the rcond function with its current
name...

D.




Jason Riedy wrote:

> 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
>
>
>
> ------------------------------------------------------------------------
>
> ## 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] = rcond (varargin)
>   ## -*- texinfo -*-
>   ## @deftypefn {Function File} {[@var{est}, @var{v}] = } rcond (@var{A}, @var{t} = min (size (A, 1), 5))
>   ## @deftypefnx {Function File} {[@var{est}, @var{v}] = } rcond (@var{A}, @var{SOLVE}, @var{SOVLE_T}, @var{t} = min (n, 5))
>   ## @deftypefnx {Function File} {[@var{est}, @var{v}] = } rcond (@var{APPLY}, @var{APPLY_T}, @var{SOLVE}, @var{SOVLE_T}, @var{n}, @var{t} = min (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) == norm (@var{A}, 1) * norm (@var{v}, 1) *
>   ## @var{est}}. When @var{est} is large, @var{v} is an approximate null
>   ## vector.
>   ##
>   ## References:
>   ## @itemize
>   ## @item Nicholas J. Higham and Françoise Tisseur, "A Block Algorithm
>   ## 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çoise Tisseur, "A Block Algorithm
>   ## 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 <ejr@...>
>   ## Keywords: linear-algebra norm estimation
>   ## Version: 0.2
>
> %!demo
> %!  N = 100;
> %!  A = randn (N) + eye (N);
> %!  rcond (A)
> %!  [L,U,P] = 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 = 5;
>   ITMAX = 10;
>
>   if ismatrix (varargin{1}),
>     n = size (varargin{1}, 1);
>     if n != size (varargin{1}, 2),
>       usage("Matrix must be square.");
>     endif
>     A = varargin{1};
>
>     if size (varargin, 2) > 1,
>       if isscalar (varargin{2}),
> t = varargin{2};
>       else
> if size (varargin, 2) < 3,
>  usage("Must supply both SOLVE and SOLVE_T.");
> else
>  SOLVE = varargin{2};
>  SOLVE_T = varargin{3};
>  if size (varargin, 2) > 3,
>    t = varargin{4};
>  endif
> endif
>       endif
>     endif
>   else
>     if size (varargin, 2) < 5,
>       usage("Implicit form of rcond requires at least 5 arguments.");
>     endif
>     APPLY = varargin{1};
>     APPLY_T = varargin{2};
>     SOLVE = varargin{3};
>     SOLVE_T = varargin{4};
>     n = varargin{5};
>     if !isscalar (n),
>       usage("Dimension argument of implicit form must be scalar.");
>     endif
>     if size (varargin, 2) > 5,
>       t = varargin{6};
>     endif
>   endif
>
>   if !exist ("t", "var"),
>     t = min (n, DEFAULT_T);
>   endif
>
>   if !exist ("SOLVE", "var"),
>     if issparse (A),
>       colord = colamd(A);
>       Pc = speye(n);
>       Pc(:,colord) = Pc;
>       [L,U,P] = lu(A(:,Pc));
>       SOLVE = @(x) Pc' * (U\ (L\ (P*x)));
>       SOLVE_T = @(x) P'*(L'\ (U'\ (Pc*x)));
>     else
>       [L,U,P] = lu(A);
>       SOLVE = @(x) U\ (L\ (P*x));
>       SOLVE_T = @(x) P' * (L'\ (U'\x));
>     endif
>   endif
>
>   if exist ("A", "var"),
>     Anorm = norm(A, 1);
>   else
>     Anorm = block_onenorm_est(APPLY, APPLY_T, n, t);
>   endif
>
>   [Ainv_norm, v, w] = block_onenorm_est(SOLVE, SOLVE_T, n, t);
>
>   est = 1/Anorm;
>   est /= Ainv_norm;
>   v = w / norm (w, 1);
>
> endfunction
>
> ## Yes, these test bounds are really loose.  There's
> ## enough randomization to trigger odd cases with hilb().
>
> %!test
> %!  N = 6;
> %!  A = hilb (N);
> %!  cA = 1 / rcond (A);
> %!  cA_test = norm (inv (A), 1) * norm (A, 1);
> %!  assert (cA, cA_test, 2**-12);
>
> %!test
> %!  N = 6;
> %!  A = hilb (N);
> %!  SOLVE = @(x) A\x; SOLVE_T = @(x) A'\x;
> %!  cA = 1 / rcond (A, SOLVE, SOLVE_T);
> %!  cA_test = norm (inv (A), 1) * norm (A, 1);
> %!  assert (cA, cA_test, 2**-12);
>
> %!test
> %!  N = 6;
> %!  A = hilb (N);
> %!  APPLY = @(x) A*x; APPLY_T = @(x) A'*x;
> %!  SOLVE = @(x) A\x; SOLVE_T = @(x) A'\x;
> %!  cA = 1 / rcond (APPLY, APPLY_T, SOLVE, SOLVE_T, N);
> %!  cA_test = norm (inv (A), 1) * norm (A, 1);
> %!  assert (cA, cA_test, 2**-6);
>
> %!test
> %!  N = 12;
> %!  A = hilb (N);
> %!  [rcondA, v] = rcond (A);
> %!  x = 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 copyright
> ##       notice, this list of conditions and the following disclaimer in the
> ##       documentation and/or other materials provided with the distribution.
> ##     * Neither the name of the University of California, Berkeley nor the
> ##       names of its contributors may be used to endorse or promote products
> ##       derived from this software without specific prior written permission.
> ##
> ## THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND ANY
> ## EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
> ## WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
> ## DISCLAIMED. IN NO EVENT SHALL THE REGENTS AND CONTRIBUTORS BE LIABLE FOR ANY
> ## DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
> ## (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
> ## LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED 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.
>
>
> ------------------------------------------------------------------------
>
> ## 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] = block_onenorm_est (varargin)
>   ## -*- texinfo -*-
>   ## @deftypefn {Function File} {[@var{est}, @var{v}, @var{w}, @var{iter}] = } block_onenorm_est (@var{A}, @var{t} = min (size (A, 1), 5))
>   ## @deftypefnx {Function File} {[@var{est}, @var{v}, @var{w}, @var{iter}] = } block_onenorm_est (@var{APPLY}, @var{APPLY_T}, @var{n}, @var{t} = 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 conjugate
>   ## 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) = @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:
>   ## @itemize
>   ## @item Nicholas J. Higham and Françoise Tisseur, "A Block Algorithm
>   ## 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çoise Tisseur, "A Block Algorithm
>   ## 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 <ejr@...>
>   ## Keywords: linear-algebra norm estimation
>   ## Version: 0.2
>
> %!demo
> %!  N = 100;
> %!  A = randn(N) + eye(N);
> %!  [L,U,P] = lu(A);
> %!  nm1inv = 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 = 5;
>   ITMAX = 10;
>
>   if ismatrix (varargin{1}),
>     n = size (varargin{1}, 1);
>     if n != size (varargin{1}, 2),
>       error("Matrix must be square.");
>     endif
>     APPLY = @(x) varargin{1} * x;
>     APPLY_T = @(x) varargin{1}' * x;
>     if size (varargin) > 1,
>       t = varargin{2};
>     else
>       t = min (n, DEFAULT_T);
>     endif
>   else
>     if size (varargin, 2) < 3,
>       print_usage();
>     endif
>     n = varargin{3};
>     APPLY = varargin{1};
>     APPLY_T = varargin{2};
>     if size (varargin) > 3,
>       t = varargin{4};
>     else
>       t = DEFAULT_T;
>     endif
>   endif
>
>   ## Initial test vectors X.
>   X = rand (n, t);
>   X = X ./ (ones (n,1) * sum (abs (X), 1));
>
>   been_there = zeros (n, 1); # Track if a vertex has been visited.
>   est_old = 0; # To check if the estimate has increased.
>   S = zeros (n, t); # Normalized vector of signs.  The normalization is
>
>   for iter=1:ITMAX+1,
>     Y = feval (APPLY, X);
>
>     ## Find the initial estimate as the largest A*x.
>     [est, ind_best] = max (sum (abs (Y), 1));
>     if (est > est_old || iter == 2),
>       w = Y(:,ind_best);
>     endif
>     if (iter >= 2 && est < est_old),
>       ## No improvement, so stop.
>       est = est_old;
>       break;
>     endif
>
>     est_old = est;
>     S_old = S;
>     if (iter > ITMAX),
>       ## Gone too far.  Stop.
>       break;
>     endif
>
>     S = 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 = 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 = sum (partest);
>       replacements = 2*(rand (n,numpar) < 0.5) - 1;
>       S(:,partest) = replacements;
>     endif
>     ## Now test for parallel vectors within S.
>     partest = any ( (S' * S - eye (t)) == n );
>     if any (partest),
>       numpar = sum (partest);
>       replacements = 2*(rand (n,numpar) < 0.5) - 1;
>       S(:,partest) = replacements;
>     endif
>    
>     Z = feval (APPLY_T, S);
>
>     ## Now find the largest non-previously-visted index per
>     ## vector.
>     h = max (abs (Z),2);
>     [mh, mhi] = max (h);
>     if iter >= 2 && mhi == ind_best,
>       ## Hit a cycle, stop.
>       break;
>     endif
>     [h, ind] = sort (h, 'descend');
>     if t > 1,
>       firstind = ind(1:t);
>       if all (been_there(firstind)),
> ## Visited all these before, so stop.
> break;
>       endif
>       ind=ind(!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 = zeros (n, t);
>     for zz = 1:t,
>       X(ind(zz),zz) = 1;
>     endfor
>     been_there(ind(1:t)) = 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 = zeros (n, 1);
>   v(ind_best) = 1;
> endfunction
>
> %!test
> %!  N = 10;
> %!  A = ones (N);
> %!  [nm1, v1, w1] = block_onenorm_est (A);
> %!  [nminf, vinf, winf] = 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 = 10;
> %!  A = ones (N);
> %!  [nm1, v1, w1] = block_onenorm_est (@(x) A*x, @(x) A'*x, N, 3);
> %!  [nminf, vinf, winf] = 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 = 5;
> %!  A = hilb (N);
> %!  [nm1, v1, w1] = block_onenorm_est (A);
> %!  [nminf, vinf, winf] = 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 = 100;
> %!  A = rand (N);
> %!  [nm1, v1, w1] = block_onenorm_est (A);
> %!  [nminf, vinf, winf] = 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 copyright
> ##       notice, this list of conditions and the following disclaimer in the
> ##       documentation and/or other materials provided with the distribution.
> ##     * Neither the name of the University of California, Berkeley nor the
> ##       names of its contributors may be used to endorse or promote products
> ##       derived from this software without specific prior written permission.
> ##
> ## THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND ANY
> ## EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
> ## WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
> ## DISCLAIMED. IN NO EVENT SHALL THE REGENTS AND CONTRIBUTORS BE LIABLE FOR ANY
> ## DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
> ## (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
> ## LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED 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.
>
>
> ------------------------------------------------------------------------
>
> ## 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] = condest (varargin)
>   ## -*- texinfo -*-
>   ## @deftypefn {Function File} {[@var{est}, @var{v}] = } condest (@var{A}, @var{t} = min (size (A, 1), 5))
>   ## @deftypefnx {Function File} {[@var{est}, @var{v}] = } condest (@var{A}, @var{SOLVE}, @var{SOVLE_T}, @var{t} = min (n, 5))
>   ## @deftypefnx {Function File} {[@var{est}, @var{v}] = } condest (@var{APPLY}, @var{APPLY_T}, @var{SOLVE}, @var{SOVLE_T}, @var{n}, @var{t} = 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) == norm
>   ## (@var{A}, 1) * norm (@var{v}, 1) / @var{est}}. When @var{est} is
>   ## large, @var{v} is an approximate null vector.
>   ##
>   ## References:
>   ## @itemize
>   ## @item Nicholas J. Higham and Françoise Tisseur, "A Block Algorithm
>   ## 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çoise Tisseur, "A Block Algorithm
>   ## 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 <ejr@...>
>   ## Keywords: linear-algebra norm estimation
>   ## Version: 0.2
>
> %!demo
> %!  N = 100;
> %!  A = randn (N) + eye (N);
> %!  condest (A)
> %!  [L,U,P] = 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] = rcond(varargin{:});
>   est = 1/est;
> endfunction
>
> ## These tests are the same as in rcond().
>
> %!test
> %!  N = 6;
> %!  A = hilb (N);
> %!  cA = condest (A);
> %!  rcA = 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 copyright
> ##       notice, this list of conditions and the following disclaimer in the
> ##       documentation and/or other materials provided with the distribution.
> ##     * Neither the name of the University of California, Berkeley nor the
> ##       names of its contributors may be used to endorse or promote products
> ##       derived from this software without specific prior written permission.
> ##
> ## THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND ANY
> ## EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
> ## WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
> ## DISCLAIMED. IN NO EVENT SHALL THE REGENTS AND CONTRIBUTORS BE LIABLE FOR ANY
> ## DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
> ## (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
> ## LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED 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.
>
>
> ------------------------------------------------------------------------
>
> _______________________________________________
> Octave-sources mailing list
> Octave-sources@...
> https://www.cae.wisc.edu/mailman/listinfo/octave-sources

_______________________________________________
Octave-sources mailing list
Octave-sources@...
https://www.cae.wisc.edu/mailman/listinfo/octave-sources

[Jason Riedy-2]

And David Bateman writes:
> I believe the rcond function should take rcond from the LAPACK
> factorization, at least I believe is what matlab does. [...]

For various reasons, I'm not referring to MATLAB(tm) itself or
any docs other than what is on-line and trivially findable.  I
have no idea if MATLAB's rcond applies to sparse matrices.  I
know LAPACK's does not.

That said, I have been losing the argument to use the block
estimator in LAPACK...  We may weaken the routine slightly
instead.  For "small" problems, computing RCOND is very, very
expensive right now.  The current norm estimator uses a very safe
but very slow triangular solve.  We likely will replace that with
the fast solve, followed by a check for overflow and re-running
with the slow solve if necessary.  That will catch and re-scale
overflows, but it won't catch possibly avoidable underflows
during estimation.

So it's up to y'all what rcond() should be.  The routine I
included computes the inverse of the 1-norm condition number of A
(iirc), and LAPACK's RCOND returned by xyySVX is the inverse of
the inf-norm condition number of the matrix that actually is
factored.  That matrix may have been equilibrated, so determining
how LAPACK's RCOND relates to your original matrix requires
examining a bunch of parameters MATLAB's rcond() doesn't return.
Without experimenting, I can't tell what MATLAB's rcond() really
does.

The only use for LAPACK's RCOND directly is to see if the
factorization algorithm came close to solving with a nearly a
singular matrix when forming a Schur complement.  It tells you
very little about the actual problem you are trying to solve or
even about any errors when using the factorization; RCOND may be
measured with an entirely different norm (scaled by
equilibration) than your problem.

And if you pick up our extra-precise refinement drivers (once we
actually ship the damned things, see LAWN 165), you will have
more and more relevant condition numbers from which to chose.

> Matlab calls the normest1 function, that is basically your
> block_onenorm_est, and then multiples the estimated inverse
> norm by norm(A,1) to get the condition number estimate.

>From MathWorks' pages, it seems normest has a very different
interface.  Hence the cumbersome block_onenorm_est name.

But your description is not the RCOND returned from LAPACK's
linear system drivers.  ;)

> condest is one of the annoying missing sparse functions in
> Octave...

Amen.  Limits Octave's applicability for many colleagues.

> I'd hesitate to include it though with the rcond function with
> its current name...

I have little ego about these things.  Use, change, and rename as
you see fit.  Just thought I should send along my good-enough
versions as a starting point.

I'm not in a position to get into the copyright assignment
argument here right now, so feel free to make sufficient,
copyrightable changes that are assigned and licensed
appropriately.  ;) These are just the routines I'm using right
now.

Jason
_______________________________________________
Octave-sources mailing list
Octave-sources@...
https://www.cae.wisc.edu/mailman/listinfo/octave-sources

[David Bateman]

Jason Riedy wrote:

> And David Bateman writes:
>  
>> I believe the rcond function should take rcond from the LAPACK
>> factorization, at least I believe is what matlab does. [...]
>>    
>
> For various reasons, I'm not referring to MATLAB(tm) itself or
> any docs other than what is on-line and trivially findable.  I
> have no idea if MATLAB's rcond applies to sparse matrices.  I
> know LAPACK's does not.
>  
>> help rcond

 RCOND  LAPACK reciprocal condition estimator.
    RCOND(X) is an estimate for the reciprocal of the
    condition of X in the 1-norm obtained by the LAPACK
    condition estimator. If X is well conditioned, RCOND(X)
    is near 1.0. If X is badly conditioned, RCOND(X) is
    near EPS.
 
    See also COND, NORM, CONDEST, NORMEST.

>> rcond(speye(1024))
??? Error using ==> rcond
Use rcond(full(S)) or condest(S).

So the matlab rcond function works only for full matrices and uses the
lapack rcond estimation from the factorization.. There needs to be a
little reorganization to get this to work optimally in Octave. Till then
a minimal rcond implementation for Octave is

function  rc = rcond (x)
  [dummy, rc] = inv (x);
endfunction

As Octave also calculates the condition number for inverse of sparse
matrices if requested, this works also for sparse matrices, though
forming the inverse of a sparse matrix is always a bad idea..

> That said, I have been losing the argument to use the block
> estimator in LAPACK...  We may weaken the routine slightly
> instead.  For "small" problems, computing RCOND is very, very
> expensive right now.  The current norm estimator uses a very safe
> but very slow triangular solve.  We likely will replace that with
> the fast solve, followed by a check for overflow and re-running
> with the slow solve if necessary.  That will catch and re-scale
> overflows, but it won't catch possibly avoidable underflows
> during estimation.
>
> So it's up to y'all what rcond() should be.  The routine I
> included computes the inverse of the 1-norm condition number of A
> (iirc), and LAPACK's RCOND returned by xyySVX is the inverse of
> the inf-norm condition number of the matrix that actually is
> factored.  That matrix may have been equilibrated, so determining
> how LAPACK's RCOND relates to your original matrix requires
> examining a bunch of parameters MATLAB's rcond() doesn't return.
> Without experimenting, I can't tell what MATLAB's rcond() really
> does.
>  

I believe its simpler than that.. I think matlab uses the xGETRF
functions of lapack that return an estimate of rcond, and matlab just
returns that estimate..
> The only use for LAPACK's RCOND directly is to see if the
> factorization algorithm came close to solving with a nearly a
> singular matrix when forming a Schur complement.  It tells you
> very little about the actual problem you are trying to solve or
> even about any errors when using the factorization; RCOND may be
> measured with an entirely different norm (scaled by
> equilibration) than your problem.
>  
xGETRF scales the inverse norm estimation by the 1-norm and so I don't
think this is the case. You're certainly right for the xGESVX LAPACK
routines though. The result is that if the matrix is near singular then
the rcond extimate will be incorrect. However, as rcond is only supposed
to be used to estimate if a matrix is singular, then this error probably
isn't an issue.

> And if you pick up our extra-precise refinement drivers (once we
> actually ship the damned things, see LAWN 165), you will have
> more and more relevant condition numbers from which to chose.
>  
Choice is good for a developer, though how to hide it from the user
behind simple driver routines is always an issue.

>  
>> Matlab calls the normest1 function, that is basically your
>> block_onenorm_est, and then multiples the estimated inverse
>> norm by norm(A,1) to get the condition number estimate.
>>    
>
> >From MathWorks' pages, it seems normest has a very different
> interface.  Hence the cumbersome block_onenorm_est name.
>  
normest is a 2-norm estimate using the power method to find the largest
singular value. Octave has that routine. No your block_onenorm_est is
like the normest1 routine from Matlab. The interface to normest1 has a
function handle to the solver passed as the first argument, whereas your
code has a solver and a solver for the conjugate transpose. Frankly I
don't think normest1 in itself is a useful function and so if its not
implemented who cares. We just need a matlab compatible interface to
condest..

> But your description is not the RCOND returned from LAPACK's
> linear system drivers.  ;)
>  
No, I was trying in my non-mathematicians manner to describe your code :-)

>  
>> condest is one of the annoying missing sparse functions in
>> Octave...
>>    
>
> Amen.  Limits Octave's applicability for many colleagues.
>  

Well, I definitely want condest and your code gives it. I see if I can
convert it to Octave style and sent it back to you for your comments..

>  
>> I'd hesitate to include it though with the rcond function with
>> its current name...
>>    
>
> I have little ego about these things.  Use, change, and rename as
> you see fit.  Just thought I should send along my good-enough
> versions as a starting point.
>
> I'm not in a position to get into the copyright assignment
> argument here right now, so feel free to make sufficient,
> copyrightable changes that are assigned and licensed
> appropriately.  ;) These are just the routines I'm using right
> now.
>  

If I used the BSD license you've chosen it allows relicensing under a
GPL, so I would suggest doing that and keeping you as the author..

Cheers, and thanks for the code.
D.




> Jason
>
>  


--
David Bateman                                David.Bateman@...
Motorola Labs - Paris                        [phone] (Ph)
Parc Les Algorithmes, Commune de St Aubin    [phone] (Mob)
91193 Gif-Sur-Yvette FRANCE                  [phone] (Fax)

The information contained in this communication has been classified as:

[x] General Business Information
[ ] Motorola Internal Use Only
[ ] Motorola Confidential Proprietary

_______________________________________________
Octave-sources mailing list
Octave-sources@...
https://www.cae.wisc.edu/mailman/listinfo/octave-sources