octave_packages/linear-algebra-2.2.0/condeig.m

   1 ## Copyright (C) 2006, 2007 Arno Onken <asnelt@asnelt.org>
   2 ##
   3 ## This program is free software; you can redistribute it and/or modify it under
   4 ## the terms of the GNU General Public License as published by the Free Software
   5 ## Foundation; either version 3 of the License, or (at your option) any later
   6 ## version.
   7 ##
   8 ## This program is distributed in the hope that it will be useful, but WITHOUT
   9 ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  10 ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
  11 ## details.
  12 ##
  13 ## You should have received a copy of the GNU General Public License along with
  14 ## this program; if not, see <http://www.gnu.org/licenses/>.
  15
  16 ## -*- texinfo -*-
  17 ## @deftypefn {Function File} {@var{c} =} condeig (@var{a})
  18 ## @deftypefnx {Function File} {[@var{v}, @var{lambda}, @var{c}] =} condeig (@var{a})
  19 ## Compute condition numbers of the eigenvalues of a matrix. The
  20 ## condition numbers are the reciprocals of the cosines of the angles
  21 ## between the left and right eigenvectors.
  22 ##
  23 ## @subheading Arguments
  24 ##
  25 ## @itemize @bullet
  26 ## @item
  27 ## @var{a} must be a square numeric matrix.
  28 ## @end itemize
  29 ##
  30 ## @subheading Return values
  31 ##
  32 ## @itemize @bullet
  33 ## @item
  34 ## @var{c} is a vector of condition numbers of the eigenvalue of
  35 ## @var{a}.
  36 ##
  37 ## @item
  38 ## @var{v} is the matrix of right eigenvectors of @var{a}. The result is
  39 ## the same as for @code{[v, lambda] = eig (a)}.
  40 ##
  41 ## @item
  42 ## @var{lambda} is the diagonal matrix of eigenvalues of @var{a}. The
  43 ## result is the same as for @code{[v, lambda] = eig (a)}.
  44 ## @end itemize
  45 ##
  46 ## @subheading Example
  47 ##
  48 ## @example
  49 ## @group
  50 ## a = [1, 2; 3, 4];
  51 ## c = condeig (a)
  52 ## @result{} [1.0150; 1.0150]
  53 ## @end group
  54 ## @end example
  55 ## @end deftypefn
  56
  57 function [v, lambda, c] = condeig (a)
  58
  59   # Check arguments
  60   if (nargin != 1 || nargout > 3)
  61     print_usage ();
  62   endif
  63
  64   if (! isempty (a) && ! ismatrix (a))
  65     error ("condeig: a must be a numeric matrix");
  66   endif
  67
  68   if (columns (a) != rows (a))
  69     error ("condeig: a must be a square matrix");
  70   endif
  71
  72   if (issparse (a) && (nargout == 0 || nargout == 1) && exist ("svds", "file"))
  73     ## Try to use svds to calculate the condition as it will typically be much
  74     ## faster than calling eig as only the smallest and largest eigenvalue are
  75     ## calculated.
  76     try
  77       s0 = svds (a, 1, 0);
  78       v = svds (a, 1) / s0;
  79     catch
  80       ## Caught an error as there is a singular value exactly at Zero!!
  81       v = Inf;
  82     end_try_catch
  83     return;
  84   endif
  85
  86   # Right eigenvectors
  87   [v, lambda] = eig (a);
  88
  89   if (isempty (a))
  90     c = lambda;
  91   else
  92     # Corresponding left eigenvectors
  93     vl = inv (v)';
  94     # Normalize vectors
  95     vl = vl ./ repmat (sqrt (sum (abs (vl .^ 2))), rows (vl), 1);
  96
  97     # Condition numbers
  98     # cos (angle) = (norm (v1) * norm (v2)) / dot (v1, v2)
  99     # Norm of the eigenvectors is 1 => norm (v1) * norm (v2) = 1
 100     c = abs (1 ./ dot (vl, v)');
 101   endif
 102
 103   if (nargout == 0 || nargout == 1)
 104     v = c;
 105   endif
 106
 107 endfunction
 108
 109 %!test
 110 %! a = [1, 2; 3, 4];
 111 %! c = condeig (a);
 112 %! expected_c = [1.0150; 1.0150];
 113 %! assert (c, expected_c, 0.001);
 114
 115 %!test
 116 %! a = [1, 3; 5, 8];
 117 %! [v, lambda, c] = condeig (a);
 118 %! [expected_v, expected_lambda] = eig (a);
 119 %! expected_c = [1.0182; 1.0182];
 120 %! assert (v, expected_v, 0.001);
 121 %! assert (lambda, expected_lambda, 0.001);
 122 %! assert (c, expected_c, 0.001);