Re: rcond, condest, and a block 1-norm estimator

[David Bateman]

I'm moving this to the maintainers list as this is now a patch against
Octave. What I've done in this patch relative to your code is

1) Relicense under the GPL. I kept the original license as a comment,
but perhaps that should be removed as well.

2) Change the style to be closer to the rest of Octave.
  - Limit the uppercase variables
  - Help text immediately following the copyright
  - write "if (a > 1)" and not "if a>1,"

3) Eliminate the rcond function as 1/condest(...) gives the same..

4) Change the name of block_onenorm_est to onenormest as something
shorter. I hesitated to make this an internal function, but thought that
a 1-norm estimator is useful in its own right. Also I didn't want to use
the name normest1 as the interface is not the same as the matlab
normest1 function.

5) Replace the block of code

      colord = colamd (A);
      Pc = speye (n);
      Pc(:, colord) = Pc;
      [L, U, P] = lu (A(:, Pc));

in condest.m with

      [L, U, P, Pc] = splu (A);

letting UMFPACK itself decide on the sparsity preserving permutation.

6) Add condest to the sparse documentation. Maybe onenormest should be
added as well.

Are you happy with the proposed changes? Can you suggest any further
changes?

D.

2007-11-22  Jason Riedy  <address@hidden>

        * linear-algebra/condest.m, linear-algebra/__onenormest__.m: New
        functions.
        * linear-algebra/Makefile.in (SOURCES): Add them to the sources.

2007-11-22  David Bateman  <address@hidden>

        * interpreter/sparse.txi: Document condest.

2007-11-22  David Bateman  <address@hidden>

        * PROJECTS: condest now implemented.

*** ./PROJECTS.orig23   2007-11-22 22:26:40.054952795 +0100
--- ./PROJECTS  2007-11-22 22:26:43.553792182 +0100
***************
*** 109,115 ****
        - colmmd      Superseded by colamd
        - treelayout
        - cholinc
-       - condest
        - bicg        Can this be taken from octave-forge?
        - bicgstab
        - cgs
--- 109,114 ----
*** ./scripts/linear-algebra/condest.m.orig23   2007-11-22 22:25:17.516741686 
+0100
--- ./scripts/linear-algebra/condest.m  2007-11-22 22:53:13.177820976 +0100
***************
*** 0 ****
--- 1,218 ----
+ ## Copyright (C) 2007, Regents of the University of California
+ ##
+ ## This file is part of Octave.
+ ##
+ ## Octave is free software; you can redistribute it and/or modify it
+ ## under the terms of the GNU General Public License as published by
+ ## the Free Software Foundation; either version 3 of the License, or (at
+ ## your option) any later version.
+ ##
+ ## Octave is distributed in the hope that it will be useful, but
+ ## WITHOUT ANY WARRANTY; without even the implied warranty of
+ ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ ## General Public License for more details.
+ ##
+ ## You should have received a copy of the GNU General Public License
+ ## along with Octave; see the file COPYING.  If not, see
+ ## <http://www.gnu.org/licenses/>.
+ 
+ ## -*- texinfo -*-
+ ## @deftypefn {Function File} address@hidden, @var{v}] =} condest (@var{A}, 
@var{t}) 
+ ## @deftypefnx {Function File} address@hidden, @var{v}] =} condest (@var{A}, 
@var{solve}, @var{solve_t}, @var{t})
+ ## @deftypefnx {Function File} address@hidden, @var{v}] =} condest 
(@var{apply}, @var{apply_t}, @var{solve}, @var{solve_t}, @var{n}, @var{t})
+ ##
+ ## Estimate the 1-norm condition number of a matrix matrix @var{A}
+ ## using @var{t} test vectors using a randomized 1-norm estimator.
+ ## If @var{t} exceeds 5, then only 5 test vectors are used.
+ ##
+ ## If the matrix is not explicit, e.g. when  estimating the condition 
+ ## number of @var{A} given an LU factorization, @code{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 (@address@hidden, 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.
+ ## @url{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." @url{http://citeseer.ist.psu.edu/223007.html}
+ ## @end itemize
+ ##
+ ## @seealso{norm, cond, onenormest}
+ ##
+ ## @end deftypefn
+ 
+ ## Code originally licensed under
+ ##
+ ##  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.
+ ##
+ ## Relicensed to GPL for inclusion in Octave.
+ 
+ ## Author: Jason Riedy <address@hidden>
+ ## Keywords: linear-algebra norm estimation
+ ## Version: 0.2
+ 
+ function [est, v] = condest (varargin)
+   if size (varargin, 2) < 1 || size (varargin, 2) > 5,
+     usage("condest: Incorrect arguments.");
+   endif
+ 
+   default_t = 5;
+ 
+   if (ismatrix (varargin{1}))
+     n = size (varargin{1}, 1);
+     if (n != size (varargin{1}, 2))
+       error ("condest: 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)
+         error ("condest: 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)
+       error ("condest: implicit form of condest 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))
+       error ("condest: 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))
+       [L, U, P, Pc] = splu (A);
+       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 = onenormest (apply, apply_t, n, t);
+   endif
+ 
+   [Ainv_norm, v, w] = onenormest (solve, solve_t, n, t);
+ 
+   est = Anorm * Ainv_norm;
+   v = w / norm (w, 1);
+ 
+ endfunction
+ 
+ %!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)
+ 
+ ## Yes, these test bounds are really loose.  There's
+ ## enough randomization to trigger odd cases with hilb().
+ 
+ %!test
+ %!  N = 6;
+ %!  A = hilb (N);
+ %!  cA = condest (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 = condest (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 = condest (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] = condest (A);
+ %!  x = A*v;
+ %!  assert (norm(x, inf), 0, eps);
*** ./scripts/linear-algebra/Makefile.in.orig23 2007-11-22 22:25:22.479513872 
+0100
--- ./scripts/linear-algebra/Makefile.in        2007-11-22 22:47:29.592593151 
+0100
***************
*** 33,41 ****
  INSTALL_PROGRAM = @INSTALL_PROGRAM@
  INSTALL_DATA = @INSTALL_DATA@
  
! SOURCES = __norm__.m commutation_matrix.m cond.m cross.m dmult.m \
!   dot.m duplication_matrix.m housh.m krylov.m krylovb.m logm.m \
!   null.m orth.m qzhess.m rank.m rref.m trace.m vec.m vech.m
  
  DISTFILES = $(addprefix $(srcdir)/, Makefile.in $(SOURCES))
  
--- 33,41 ----
  INSTALL_PROGRAM = @INSTALL_PROGRAM@
  INSTALL_DATA = @INSTALL_DATA@
  
! SOURCES = __norm__.m commutation_matrix.m cond.m condest.m cross.m \
!   dmult.m dot.m duplication_matrix.m housh.m krylov.m krylovb.m logm.m \
!   null.m onenormest.m orth.m qzhess.m rank.m rref.m trace.m vec.m vech.m
  
  DISTFILES = $(addprefix $(srcdir)/, Makefile.in $(SOURCES))
  
*** ./scripts/linear-algebra/onenormest.m.orig23        2007-11-22 
22:25:10.177078611 +0100
--- ./scripts/linear-algebra/onenormest.m       2007-11-22 22:59:22.923847897 
+0100
***************
*** 0 ****
--- 1,265 ----
+ ## Copyright (C) 2007, Regents of the University of California
+ ##
+ ## This file is part of Octave.
+ ##
+ ## Octave is free software; you can redistribute it and/or modify it
+ ## under the terms of the GNU General Public License as published by
+ ## the Free Software Foundation; either version 3 of the License, or (at
+ ## your option) any later version.
+ ##
+ ## Octave is distributed in the hope that it will be useful, but
+ ## WITHOUT ANY WARRANTY; without even the implied warranty of
+ ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ ## General Public License for more details.
+ ##
+ ## You should have received a copy of the GNU General Public License
+ ## along with Octave; see the file COPYING.  If not, see
+ ## <http://www.gnu.org/licenses/>.
+ 
+ ## -*- texinfo -*-
+ ## @deftypefn {Function File} address@hidden, @var{v}, @var{w}, @var{iter}] 
=} onenormest (@var{A}, @var{t}) 
+ ## @deftypefnx {Function File} address@hidden, @var{v}, @var{w}, @var{iter}] 
=} onenormest (@var{apply}, @var{apply_t}, @var{n}, @var{t})
+ ##
+ ## Apply Higham and Tisseur's randomized block 1-norm estimator to
+ ## matrix @var{A} using @var{t} test vectors. If @var{t} exceeds 5, then
+ ## only 5 test vectors are used.
+ ##
+ ## If the matrix is not explicit, e.g. when estimating the norm of 
+ ## @code{inv (@var{A})} given an LU factorization, @code{onenormest} 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.
+ ## @url{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." @url{http://citeseer.ist.psu.edu/223007.html}
+ ## @end itemize
+ ##
+ ## @seealso{condest, norm, cond}
+ ## @end deftypefn
+ 
+ ## Code originally licensed under
+ ##
+ ##  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.
+ ##
+ ## Relicensed to GPL for inclusion in Octave.
+ 
+ ## Author: Jason Riedy <address@hidden>
+ ## Keywords: linear-algebra norm estimation
+ ## Version: 0.2
+ 
+ function [est, v, w, iter] = onenormest (varargin)
+ 
+   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 ("onenormest: 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
+ 
+ %!demo
+ %!  N = 100;
+ %!  A = randn(N) + eye(N);
+ %!  [L,U,P] = lu(A);
+ %!  nm1inv = onenormest(@(x) U\(L\(P*x)), @(x) P'*(L'\(U'\x)), N, 30)
+ %!  norm(inv(A), 1)
+ 
+ %!test
+ %!  N = 10;
+ %!  A = ones (N);
+ %!  [nm1, v1, w1] = onenormest (A);
+ %!  [nminf, vinf, winf] = onenormest (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] = onenormest (@(x) A*x, @(x) A'*x, N, 3);
+ %!  [nminf, vinf, winf] = onenormest (@(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] = onenormest (A);
+ %!  [nminf, vinf, winf] = onenormest (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] = onenormest (A);
+ %!  [nminf, vinf, winf] = onenormest (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)
*** ./doc/interpreter/sparse.txi.orig23 2007-11-22 22:27:13.781404593 +0100
--- ./doc/interpreter/sparse.txi        2007-11-22 22:28:13.678655027 +0100
***************
*** 490,498 ****
  @item Linear algebra:
    @dfn{matrix_type}, @dfn{spchol}, @dfn{cpcholinv}, 
    @dfn{spchol2inv}, @dfn{spdet}, @dfn{spinv}, @dfn{spkron},
!   @dfn{splchol}, @dfn{splu}, @dfn{spqr}, @dfn{normest}, 
    @dfn{sprank}
! @c @dfn{condest}, @dfn{spaugment}
  @c @dfn{eigs}, @dfn{svds} but these are in octave-forge for now
  
  @item Iterative techniques:
--- 490,498 ----
  @item Linear algebra:
    @dfn{matrix_type}, @dfn{spchol}, @dfn{cpcholinv}, 
    @dfn{spchol2inv}, @dfn{spdet}, @dfn{spinv}, @dfn{spkron},
!   @dfn{splchol}, @dfn{splu}, @dfn{spqr}, @dfn{normest}, @dfn{condest},
    @dfn{sprank}
! @c @dfn{spaugment}
  @c @dfn{eigs}, @dfn{svds} but these are in octave-forge for now
  
  @item Iterative techniques:
***************
*** 828,833 ****
--- 828,835 ----
  
  @DOCSTRING(normest)
  
+ @DOCSTRING(condest)
+ 
  @DOCSTRING(spchol)
  
  @DOCSTRING(spcholinv)