# HG changeset patch # User Rik # Date 1468851307 25200 # Node ID e859caa53399519cecb4f694faad3387313e17bc # Parent faae09d7ed5b50996be1b68765603e2bcb7430d8 condeig.m: New function imported from Octave-Forge linear-algebra package. * scripts/linear-algebra/condeig.m: New function. * scripts/linear-algebra/module.mk: Add to build system. * linalg.txi: Add to manual. * cond.m: Rephrase to state that cond calculates a condition number with respect to inversion. Mention condeig in the seealso block. diff -r faae09d7ed5b -r e859caa53399 doc/interpreter/linalg.txi --- a/doc/interpreter/linalg.txi Sun Jul 17 20:40:06 2016 -0400 +++ b/doc/interpreter/linalg.txi Mon Jul 18 07:15:07 2016 -0700 @@ -100,6 +100,8 @@ @DOCSTRING(cond) +@DOCSTRING(condeig) + @DOCSTRING(det) @DOCSTRING(eig) diff -r faae09d7ed5b -r e859caa53399 scripts/linear-algebra/cond.m --- a/scripts/linear-algebra/cond.m Sun Jul 17 20:40:06 2016 -0400 +++ b/scripts/linear-algebra/cond.m Mon Jul 18 07:15:07 2016 -0700 @@ -19,7 +19,8 @@ ## -*- texinfo -*- ## @deftypefn {} {} cond (@var{A}) ## @deftypefnx {} {} cond (@var{A}, @var{p}) -## Compute the @var{p}-norm condition number of a matrix. +## Compute the @var{p}-norm condition number of a matrix with respect to +## inversion. ## ## @code{cond (@var{A})} is defined as ## @tex @@ -40,7 +41,7 @@ ## indicates small changes (such as underflow or round-off error) will produce ## large changes in the resulting output. In such cases the solution results ## from numerical computing are not likely to be accurate. -## @seealso{condest, rcond, norm, svd} +## @seealso{condest, rcond, condeig, norm, svd} ## @end deftypefn ## Author: jwe diff -r faae09d7ed5b -r e859caa53399 scripts/linear-algebra/condeig.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/linear-algebra/condeig.m Mon Jul 18 07:15:07 2016 -0700 @@ -0,0 +1,135 @@ +## Copyright (C) 2006, 2007 Arno Onken +## +## 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 +## . + +## -*- texinfo -*- +## @deftypefn {} {@var{c} =} condeig (@var{a}) +## @deftypefnx {} {[@var{v}, @var{lambda}, @var{c}] =} condeig (@var{a}) +## Compute condition numbers of a matrix with respect to eigenvalues. +## +## The condition numbers are the reciprocals of the cosines of the angles +## between the left and right eigenvectors; Large values indicate that the +## matrix has multiple distinct eigenvalues. +## +## The input @var{a} must be a square numeric matrix. +## +## The outputs are: +## +## @itemize @bullet +## @item +## @var{c} is a vector of condition numbers for the eigenvalues of +## @var{a}. +## +## @item +## @var{v} is the matrix of right eigenvectors of @var{a}. The result is +## equivalent to calling @code{[@var{v}, @var{lambda}] = eig (@var{a)}}. +## +## @item +## @var{lambda} is the diagonal matrix of eigenvalues of @var{a}. The +## result is equivalent to calling +## @code{[@var{v}, @var{lambda}] = eig (@var{a})}. +## @end itemize +## +## Example +## +## @example +## @group +## a = [1, 2; 3, 4]; +## c = condeig (a) +## @result{} [1.0150; 1.0150] +## @end group +## @end example +## @seealso{eig, cond, balance} +## @end deftypefn + +function [v, lambda, c] = condeig (a) + + if (nargin != 1) + print_usage (); + endif + + if (! (isnumeric (a) && issquare (a))) + error ("condeig: A must be a square numeric matrix"); + endif + + if (issparse (a) && nargout <= 1) + ## Try to use svds to calculate the condition number as it will typically + ## be much faster than calling eig as only the smallest and largest + ## eigenvalue are calculated. + + ## FIXME: This calculates one condition number for the entire matrix. + ## In the full case, separate condition numbers are calculated for every + ## eigenvalue. + try + s0 = svds (a, 1, 0); # min eigenvalue + v = svds (a, 1) / s0; # max eigenvalue + catch + ## Caught an error as there is a singular value exactly at zero!! + v = Inf; + end_try_catch + return; + endif + + ## Right eigenvectors + [v, lambda] = eig (a); + + if (isempty (a)) + c = []; + else + ## Corresponding left eigenvectors + vl = inv (v)'; + ## Normalize vectors + vl ./= repmat (sqrt (sum (abs (vl .^ 2))), rows (vl), 1); + + ## Condition numbers + ## Definition: cos (angle) = (norm (v1) * norm (v2)) / dot (v1, v2) + ## Eigenvectors have been normalized so `norm (v1) * norm (v2)' = 1 + c = abs (1 ./ dot (vl, v)'); + endif + + if (nargout <= 1) + v = c; + endif + +endfunction + + +%!test +%! a = [1, 2; 3, 4]; +%! c = condeig (a); +%! expected_c = [1.0150; 1.0150]; +%! assert (c, expected_c, 0.001); + +%!test +%! a = [1, 3; 5, 8]; +%! [v, lambda, c] = condeig (a); +%! [expected_v, expected_lambda] = eig (a); +%! expected_c = [1.0182; 1.0182]; +%! assert (v, expected_v, 0.001); +%! assert (lambda, expected_lambda, 0.001); +%! assert (c, expected_c, 0.001); + +# Test empty input +%!assert (condeig ([]), []) + +## Test input validation +%!error condeig () +%!error condeig (1,2) +%!error condeig ({1}) +%!error condeig (ones (3,2)) +%!error condeig (ones (2,2,2)) + diff -r faae09d7ed5b -r e859caa53399 scripts/linear-algebra/module.mk --- a/scripts/linear-algebra/module.mk Sun Jul 17 20:40:06 2016 -0400 +++ b/scripts/linear-algebra/module.mk Mon Jul 18 07:15:07 2016 -0700 @@ -5,6 +5,7 @@ scripts/linear-algebra/commutation_matrix.m \ scripts/linear-algebra/cond.m \ scripts/linear-algebra/condest.m \ + scripts/linear-algebra/condeig.m \ scripts/linear-algebra/cross.m \ scripts/linear-algebra/duplication_matrix.m \ scripts/linear-algebra/expm.m \