1 ## Copyright (C) 2006, 2007 Arno Onken <asnelt@asnelt.org>
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
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
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/>.
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.
23 ## @subheading Arguments
27 ## @var{a} must be a square numeric matrix.
30 ## @subheading Return values
34 ## @var{c} is a vector of condition numbers of the eigenvalue of
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)}.
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)}.
46 ## @subheading Example
52 ## @result{} [1.0150; 1.0150]
57 function [v, lambda, c] = condeig (a)
60 if (nargin != 1 || nargout > 3)
64 if (! isempty (a) && ! ismatrix (a))
65 error ("condeig: a must be a numeric matrix");
68 if (columns (a) != rows (a))
69 error ("condeig: a must be a square matrix");
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
80 ## Caught an error as there is a singular value exactly at Zero!!
87 [v, lambda] = eig (a);
92 # Corresponding left eigenvectors
95 vl = vl ./ repmat (sqrt (sum (abs (vl .^ 2))), rows (vl), 1);
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)');
103 if (nargout == 0 || nargout == 1)
112 %! expected_c = [1.0150; 1.0150];
113 %! assert (c, expected_c, 0.001);
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);