<- 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. Some statistical analysis uses eigenvalues. Sometimes only the largest or second largest eigenvalues are important. 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. Vertical bars, | |, in this lecture is not absolute value, it means vector or matrix | 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 | (n e's for n by n A) 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 e^n + c_n-1 * e^n-1 + ... + c_2 * e^2 + c_1 * e + c_0 = 0 (Divide through by c_n, to have exactly n unknown coefficients 0..n-1). (Special case because of " = 0 " ) 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. Typically, we require the length of v to be 1. 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. The eigenvectors may be different. 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. The roots and thus the eigenvalues may be complex. Note that a matrix with all real coefficients may have complex eigenvalues and or complex eigenvectors: A= | 1 -1 | Lambda= 1 + i vec x= | i -1 | | 1 1 | Lambda= 1 - i vec x= |-i 1 | 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 | Thus, if you have code that can find eigenvalues, you have code that can find roots of a polynomial. Not the most efficient method. 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 one eigenvalue equal to the sum of a row and all other eigenvalues equal to zero. 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 eigenvalues 0, 0, 3 | 1 1 1 | | 2 2 2 2 | 2 has eigenvalues 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 code that claims to compute eigenvalues
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) 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. The trace of a matrix is the sum of the diagonal elements of the matrix. 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 of B should be the same as the eigenvalues of A, yet the eigenvectors may be different. 4) Randomly interchange some rows and corresponding columns of 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.Sample code in several languages
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. (Matrix initialized with complex values 4+1i, ...) 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 535_double.for double precision 535.dat 535_double_f.out 535_2.dat 535_2.out 535b.dat 535b.out 535_double.f90 just ! for C 535.rand 535_rand_f90.out 535.ranA 535_ranA_f90.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 The general eigen problem is A*v = e*B*v Checks are also det|V|=0 det|(B^-1*A)-eI|=0 Each eigenvalue e has an eigenvector v as a column of V. A rather quick and dirty translation of 535_double.for to Ada: test_eigen.adb cxhes.adb Hessenberg reduction cxval.adb compute eigenvalues cxvec.adb compute eigenvectors cxcheck.adb compute residual eigdet.adb check eigenvalues evalvec.adb check eigenvectors cxinverse.adb invert matrix complex_arrays.ads complex_arrays.adb real_arrays.ads real_arrays.adb test_eigen_ada.out Built command gnatmake test_eigen.adb Run command test_eigen > test_eigen_ada.out rmatin.adb read 535_double.f90 data test_535_eigen.adb test_535_eigen_ada.out 535.dat test_535b_eigen_ada.out 535b.dat test_535_2_eigen_ada.out 535_2.dat test_535_x12_eigen_ada.out 535_x12_1.dat In aircraft design, there are stability questions that can be answered using eigenvalue computations. In Matlab, you can draw objects, such as this crude airplane, as well as doing numerical computation.drawn by plane.m wing_2d.m rotx_.m roty_.m rotz_.m
OK, long lecture, extra topic: Generating random numbers, valuable for some testing.
First we generate uniformly distributed random numbers in the range 0.0 to 1.0. A uniform distribution should have roughly the same number of values in each range, bin. Gaussian or Normal distributions should have roughly a bell shaped curve, based on a given mean and sigma. Sigma is the standard deviation. Sigma squared is the variance. udrnrt.h three generators udrnrt.c three generators code plot_udrnrt.c test_program plot_udrnrt1.out udrnrt output and plot2.0*udrnrt()-1.0 output, -1.0 to 1.0 plot_udrnrt2.out gauss output and plot
plot_udrnrt3.out gauss output and plot
plot_udrnrt4.out sine wave output and plot
test_rand.py uniform and gauss plot
1,000,000 uniformly distributed samples
1,000,000 Gaussian distributed samples, mean=0.0, sigma=3.0
<- previous index next ->