changeset 11649:196a759fc7e8 release-3-0-x

new norm arg for cond
author John W. Eaton <jwe@octave.org>
date Fri, 15 Feb 2008 16:37:32 -0500
parents 7f1f0b6c2b30
children a8e0f7184a59
files scripts/ChangeLog scripts/linear-algebra/cond.m
diffstat 2 files changed, 54 insertions(+), 23 deletions(-) [+]
line wrap: on
line diff
--- a/scripts/ChangeLog
+++ b/scripts/ChangeLog
@@ -1,3 +1,8 @@
+2008-02-15  Rolf Fabian  <r.fabian@jacobs-university.de>
+
+	* linear-algebra/cond.m: New optional second argument to
+	specify 1-norm, inf-norm, or frobenius-norm. 
+
 2008-02-08  Kostas Poulios  <poulios.konstantinos@googlemail.com>
 
 	* plot/__quiver__.m: make arrow head be in z-plane of the arrow
--- a/scripts/linear-algebra/cond.m
+++ b/scripts/linear-algebra/cond.m
@@ -1,5 +1,5 @@
 ## Copyright (C) 1993, 1994, 1995, 1996, 1997, 1999, 2000, 2003, 2004,
-##               2005, 2006, 2007 John W. Eaton
+##               2005, 2006, 2007, 2008 John W. Eaton
 ##
 ## This file is part of Octave.
 ##
@@ -18,40 +18,66 @@
 ## <http://www.gnu.org/licenses/>.
 
 ## -*- texinfo -*-
-## @deftypefn {Function File} {} cond (@var{a})
-## Compute the (two-norm) condition number of a matrix. @code{cond (a)} is
-## defined as @code{norm (a) * norm (inv (a))}, and is computed via a
-## singular value decomposition.
-## @seealso{norm, svd, rank}
+## @deftypefn {Function File} {} cond (@var{a},@var{p})
+## Compute the @var{p}-norm condition number of a matrix. @code{cond (@var{a})} is
+## defined as @code{norm (@var{a}, @var{p}) * norm (inv (@var{a}), @var{p})}.
+## By default @code{@var{p}=2} is used which implies a (relatively slow)
+## singular value decomposition. Other possible selections are 
+## @code{@var{p}= 1, Inf, inf, 'Inf', 'fro'} which are generally faster.
+## @seealso{norm, inv, det, svd, rank}
 ## @end deftypefn
 
 ## Author: jwe
 
-function retval = cond (a)
+function retval = cond (a, p)
 
-  if (nargin == 1)
+  if (nargin && nargin < 3)
     if (ndims (a) > 2)
-      error ("cond: Only valid on 2-D objects")
+      error ("cond: only valid on 2-D objects")
+    endif
+
+    if (nargin <2)
+      p = 2;
     endif
 
-    [nr, nc] = size (a);
-    if (nr == 0 || nc == 0)
-      retval = 0.0;
-    endif
-    if (any (any (isinf (a) | isnan (a))))
-      error ("cond: argument must not contain Inf or NaN values");
+    if (! isstr (p) && p == 2)
+      [nr, nc] = size (a);
+      if (nr == 0 || nc == 0)
+        retval = 0.0;
+      elseif (any (any (isinf (a) | isnan (a))))
+        error ("cond: argument must not contain Inf or NaN values");
+      else
+        sigma   = svd (a);
+        sigma_1 = sigma(1);
+        sigma_n = sigma(end);
+        if (sigma_1 == 0 || sigma_n == 0)
+          retval = Inf;
+        else
+          retval = sigma_1 / sigma_n;
+        endif
+      endif
     else
-      sigma = svd (a);
-      sigma_1 = sigma(1);
-      sigma_n = sigma(length (sigma));
-      if (sigma_1 == 0 || sigma_n == 0)
-        retval = Inf;
-      else
-        retval = sigma_1 / sigma_n;
-      endif
+      retval = norm (a, p) * norm (inv (a), p);  
     endif
   else
     print_usage ();
   endif
 
 endfunction
+
+%!test
+%! y= [7, 2, 3; 1, 3, 4; 6, 4, 5];
+%! tol = 1e-6;
+%! type = {1, 2, 'fro', 'inf', inf};
+%! for n = 1:numel(type)
+%!   rcondition(n) = 1 / cond (y, type{n});
+%! endfor
+%! assert (rcondition, [0.017460, 0.019597, 0.018714, 0.012022, 0.012022], tol);
+
+%!assert (abs (cond ([1, 2; 2, 1]) - 3) < sqrt (eps));
+
+%!assert (cond ([1, 2, 3; 4, 5, 6; 7, 8, 9]) > 1.0e+16);
+
+%!error cond ();
+
+%!error cond (1, 2, 3);