<- previous index next ->
Eigenvalue and Eigenvector computation may be the most prolific for
special case numerical computation. Considering the size and speed
of modern computers, I use a numerical solution for a general
complex matrix. Thus, I do not have to worry about the constraints
on the matrix (is it numerically positive definite?)
This type of numerical algorithm, you do not want to develop yourself.
The technique is to find a working numerical code, test it yourself,
then adapt the code to your needs, possibly converting the code to
a different language.
Eigenvalues have application in the solution of some physics problems.
They are also used in some solutions of differential equations.
My interest comes from the vast types of testing that can, and should,
be performed on any code that claims to compute eigenvalues and
eigenvectors. Thorough testing of any numeric code you are going
to use is essential.
Information about eigenvalues, e (no lambda in plain ASCII)
and eigenvectors, v, for arbitrary n by n complex matrix A.
There are exactly n eigenvalues (some may have multiplicity greater than 1)
For every eigenvalue there is a corresponding eigenvector.
For eigenvalues with multiplicity greater than 1,
each has a unique eigenvector.
The set of n eigenvectors form a basis, they are all mutually orthogonal.
(The dot product of any pair of eigenvectors is zero.)
View this page in a fixed width font, else the matrices are shambles.
| 1 0 0 |
det|A-eI| = 0 defines e, where I is an n by n identity matrix | 0 1 0 |
zero except 1 on the diagonal | 0 0 1 |
For a 3 by 3 matrix A:
| a11-e a12 a13 |
det| a21 a22-e a23 | = (a11-e)*(a22-e)*(a33-e)+a12*a23*a31+a13*a32*a21-
| a31 a32 a33-e | a31*(a22-e)*a13-a21*a12*(a33-e)-a32*a23*(a11-e)
Writing out the above determinant gives the "Characteristic Equation"
of the matrix A.
Combining terms gives c_n * e^n + ... + c_2 * e^2 + c_1 * e + c_0 = 0
Dividing through by c_n gives exactly n coefficients.
There are exactly n roots for an nth order polynomial and the
n roots of the characteristic equation are the n eigenvalues.
The relation between each eigenvalue and its corresponding eigenvector is
Av = ev where |v| is non zero.
Given a matrix A and a non singular matrix P and P inverse p^-1
B = P A P^-1 the matrix B has the same eigenvalues as matrix A.
B is a similarity transform of A.
The diagonal elements of a diagonal matrix are the eigenvalues of the matrix.
| a11 0 0 0 |
| 0 a22 0 0 |
| 0 0 a33 0 | has eigenvalues a11, a22, a33, a44 and
| 0 0 0 a44 |
corresponding eigenvectors
v1=|1 0 0 0| v2=|0 1 0 0| v3=|0 0 1 0| v4=|0 0 0 1|
Notice that the eigenvectors are not necessarily unique and may be
scaled by an arbitrary, non zero, constant. Normalizing the length
of each eigenvector to 1.0 is common.
The eigenvalues of a 2 by 2 matrix are easily computed as the roots
of a second order equation.
det| a11-e a12 |=0 or (a11-e)*(a22-e) - a12*a21 = 0 or
| a21 a22-e |
e^2 - (a11+a22) e + (a11*a22-a12*a21) = 0
Let a=1, b=-(a11+a22), c= (a11*a22-a12*a21)
then e = (-b +/- sqrt(b^2-4*a*c)/2*a computes the two eigenvalues.
Note that a matrix with all real coefficients may have complex eigenvalues.
Computing the characteristic equation is usually not a good way
to compute eigenvalues for n greater than 4 or 5.
It becomes difficult to compute the coefficients of the characteristic
equation accurately and it is also difficult to compute the roots accurately.
Note that given a high order polynomial, a matrix can be set up from
the coefficients such that the eigenvalues of the matrix are the roots
of the polynomial. x^4 + c_3 x^3 + c_2 x^2 + c_1 x + c_0 = 0
| -c_3 -c_2 -c_1 -c_0 | where c_4 of the polynomial is 1.0
| 1 0 0 0 |
| 0 1 0 0 |
| 0 0 1 0 |
The maximum value of the sum of the absolute values of each row and column
is a upper bound on the absolute value of the largest eigenvalue.
This maximum value is typically called the L1 norm of the matrix.
Scaling a matrix by multiplying every element by a constant causes
every eigenvalue to be multiplied by that constant. The constant may
be less than one, thus dividing works the same.
The eigenvalues of the inverse of a matrix are the reciprocals
of the eigenvalues of the original matrix.
A matrix with all elements equal has all zero eigenvalues.
The eigenvectors are typically chosen as the unit vectors.
| 0 0 0 |
A = | 0 0 0 | has three eigenvalues, all equal to zero
| 0 0 0 |
| 1 1 1 | | 2 2 2 2 |
| 1 1 1 | | 2 2 2 2 | 1 has eigen values 0, 0, 3
| 1 1 1 | | 2 2 2 2 | 2 has eigen values 0, 0, 0, 8
| 2 2 2 2 | in general n-1 zeros and a row sum
Each row that increases the singularity of a matrix, increases
the multiplicity of some eigenvalue.
The trace of a matrix is the sum of the diagonal elements.
The trace of a matrix is equal to the sum of the eigenvalues.
In order to keep the same eigenvalues, interchanging two rows
of a matrix, then requires interchanging the corresponding two columns.
The eigenvectors will probably be different.
Testing a program that claims to compute eigenvalues and eigenvectors
is interesting because there are many possible tests. All should be
used.
Given A is an n by n complex matrix (that may have all real elements),
using IEEE 64-bit floating point and good algorithms:
1) Evaluate the determinant det|A-eI| for each eigenvalue, e.
The result should be near zero. 10^-9 or smaller can be expected
when A is small and eigenvalues are about the same magnitude.
2) Evaluate each eigenvalue with its corresponding eigenvector.
Av-ev should be a vector with all elements near zero.
Typically check the magnitude of the largest element.
10^-9 or smaller can be expected when A is small.
3) Invert the matrix of eigenvectors with code that computes
the degree of singularity. An n by n matrix should have
n-1 singularities, only one independent row.
Alternately, compute the dot product of every pair of eigenvectors
and check for near zero.
4) Compute the trace of A and subtract the sum of the eigenvalues.
The result should be near zero.
5) Compute the maximum of the sum of the absolute values of each
row and column of A. Check that the absolute value of every
eigenvalue is less that or equal this maximum.
Generating test matrices to be used for testing.
1) Try matrices with n=1,2,3 first.
All zero matrix, all eigenvalues zero and eigenvectors should
be the unit basis vectors. If the length of the eigenvectors
is not 1.0, then you have to normalize them.
2) Try diagonal matrices with n=1,2,3,4
Typically put 1, 2, 3, 4 on the diagonal to make it easy to
check the values of the computed eigenvalues.
3) Generate a random n by n matrix, P, with real and imaginary values.
Compute P inverse, P^-1.
Compute matrix B = P A P^-1 for the A matrices in 2)
The eigenvalues should be the same yet the eigenvectors
may be different.
4) Randomly interchange some rows of B in 3) and randomly interchange
some columns of the same B. The eigenvalues should be the same
yet the eigenvectors may be different.
5) Choose a set of values, typically complex values e1, e2, ..., en.
Compute the polynomial that has those roots
(x-e1)*(x-e2)*...*(x-en) and convert to the form
x^n + c_n-1 x^n-1 + ... c_2 x^2 + c_1 x + c_0
Create the matrix n by n with the first row being negative c's
and the subdiagonal being 1's.
| -c_n-1 ... -c_2 -c_1 -c_0 |
| 1 0 0 0 |
...
| 0 0 0 0 |
| 0 1 0 0 |
| 0 0 1 0 |
The eigenvalues should be e1, e2, ..., en
Then do a similarity transform with a random matrix
and check that the same eigenvalues are computed.
Yes, this is a lot of testing, yet once coded, it can be used
for that "newer and better" version the used software salesman
is trying to sell you.
Now, look at some code that computes eigenvalues and eigenvectors
and the associated test code.
First, MatLab
% eigen.m demonstrate eigenvalues in MatLab
format compact
e1=1; % desired eigenvalues
e2=2;
e3=3;
e4=4;
P=poly([e1 e2 e3 e4]) % make polynomial from roots
r=roots(P) % check that roots come back
A=zeros(4);
A(1,1)=-P(2); % build matrix A
A(1,2)=-P(3);
A(1,3)=-P(4);
A(1,4)=-P(5);
A(2,1)=1;
A(3,2)=1;
A(4,3)=1
[v e]=eig(A) % compute eigenvectors and eigenvalues
I=eye(4); % identity matrix
for i=1:4
z1=det(A-I.*e(i,i)) % should be near zero
end
for i=1:4
z2=A*v(:,i)-e(i,i)*v(:,i) % note columns, near zero
end
z3=trace(A)-trace(e) % should be near zero
% annotated output
%
%P =
% 1 -10 35 -50 24
%r =
% 4.0000 these should be eigenvalues
% 3.0000
% 2.0000
% 1.0000
%A =
% 10 -35 50 -24 polynomial coefficients
% 1 0 0 0
% 0 1 0 0
% 0 0 1 0
%v =
% 0.9683 0.9429 0.8677 0.5000 eigenvectors
% 0.2421 0.3143 0.4339 0.5000 are
% 0.0605 0.1048 0.2169 0.5000 columns
% 0.0151 0.0349 0.1085 0.5000
%e =
% 4.0000 0 0 0 eigenvalues
% 0 3.0000 0 0 as diagonal
% 0 0 2.0000 0 of the matrix
% 0 0 0 1.0000
%z1 =
% 1.1280e-13 i=1 first eigenvalue
%z1 =
% 5.0626e-14
%z1 =
% 1.3101e-14
%z1 =
% 2.6645e-15
%z2 =
% 1.0e-14 * note multiplier i=1
% -0.4441
% -0.1110
% -0.0333
% -0.0035
%z2 =
% 1.0e-14 *
% -0.4441
% -0.1554
% -0.0444
% -0.0042
%z2 =
% 1.0e-14 *
% -0.7994
% -0.2776
% -0.0833
% -0.0028
%z2 =
% 1.0e-13 *
% -0.3103
% -0.0600
% -0.0122
% -0.0033
%z3 =
% -7.1054e-15 trace check
Now very similar code in MatLab using complex eigenvalues.
A similarity transform is applied and scaling is applied.
One eigenvalue check is now accurate to about 10^-7.
eigen2.m
eigen2_m.out
Using complex similarity transform:
eigen3.m
eigen3_m.out
Note that the MatLab help on "eig" says they use the
LAPACK routines. The next lecture covers some LAPACK.
A Fortran program to compute eigenvalues and eigenvectors from
TOMS, ACM Transactions on Mathematical Software, algorithm 535:
535.for
535.dat
535.out
535_roots.out
535_d.out
A Java program to compute eigenvalues and eigenvectors is:
Eigen2.java
TestEigen2.java
TestEigen2.out
An Ada program to compute eigenvalues and eigenvectors is:
generic_complex_eigenvalues.ada
test_generic_complex_eigenvalues.ada
generic_complex_eigen_check.ada
A rather limited use eigenvalue computation method is the
power method. It works sometimes and may require hundreds
of iterations. The following code shows that is can work
for finding the largest eigenvalue of a small real matrix.
eigen_power.c
eigen_power_c.out
<- previous index next ->