Using IDL: Mathematics

Eigenvalues and Eigenvectors

Consider a system of equations that satisfies the array-vector relationship Ax = lx, where A is an n-by-n array, x is an n-element vector, and l is a scalar. A scalar l and nonzero vector x that simultaneously satisfy this relationship are referred to as an eigenvalue and an eigenvector of the array A, respectively. The set of all eigenvectors of the array A is then referred to as the eigenspace of A. Ideally, the eigenspace will consist of n linearly-independent eigenvectors, although this is not always the case.

IDL computes the eigenvalues and eigenvectors of a real symmetric n-by-n array using Householder transformations and the QL algorithm with implicit shifts. The eigenvalues of a real, n-by-n nonsymmetric array are computed from the upper Hessenberg form of the array using the QR algorithm. Eigenvectors are computed using inverse subspace iteration.

Although it is not practical for numerical computation, the problem of computing eigenvalues and eigenvectors can also be defined in terms of the determinant function. The eigenvalues of an n-by-n array A are the roots of the polynomial defined by det(A - lI), where I is the identity matrix (an array with 1s on the main diagonal and 0s elsewhere) with the same dimensions as A. By expressing eigenvalues as the roots of a polynomial, we see that they can be either real or complex. If an eigenvalue is complex, its corresponding eigenvectors are also complex.

The following examples demonstrate how to use IDL to compute the eigenvalues and eigenvectors of real, symmetric and nonsymmetric n-by-n arrays. Note that it is possible to check the accuracy of the computed eigenvalues and eigenvectors by algebraically manipulating the definition given above to read Ax - lx = 0; in this case 0 denotes an n-element vector, all elements of which are zero.

Symmetric Array with n Distinct Real Eigenvalues

To compute eigenvalues and eigenvectors of a real, symmetric, n-by-n array, begin with a symmetric array A.

A = [[ 3.0,  1.0, -4.0], $  
     [ 1.0,  3.0, -4.0], $  
     [-4.0, -4.0,  8.0]]  
; Compute the tridiagonal form of A:  
TRIRED, A, D, E  
; Compute the eigenvalues (returned in vector D) and  
; the eigenvectors (returned in the rows of the array A):  
TRIQL, D, E, A  
; Print eigenvalues:  
PRINT, D  

IDL prints:

  2.00000  4.76837e-07  12.0000  

The exact values are: [2.0, 0.0, 12.0].

;Print the eigenvectors, which are returned as row vectors in A:  
PRINT, A  

IDL prints:

  0.707107  -0.707107   0.00000  
 -0.577350  -0.577350  -0.577350  
 -0.408248  -0.408248   0.816497  

The exact eigenvectors are:

Nonsymmetric Array with n Distinct Real and Complex Eigenvalues

To compute the eigenvalues and eigenvectors of a real, nonsymmetric n-by-n array, begin with an array A. In this example, there are n distinct eigenvalues and n linearly-independent eigenvectors.

A = [[ 1.0, 0.0,  2.0], $  
     [ 0.0, 1.0, -1.0], $  
     [-1.0, 1.0,  1.0]]  
; Reduce to upper Hessenberg format:  
hes = ELMHES(A)  
; Compute the eigenvalues:  
evals = HQR(hes)  
; Print the eigenvalues:  
PRINT, evals	  

IDL prints:

(  1.00000, -1.73205)(  1.00000,  1.73205)  
(  1.00000,  0.00000)  

; Initialize a variable to contain the residual:  
residual = 1  
; Compute the eigenvectors and the residual for each  
; eigenvalue/eigenvector pair, using double-precision arithmetic:  
evecs = EIGENVEC(A, evals, /DOUBLE, RESIDUAL=residual)  
; Print the eigenvectors, which are returned as  
; row vectors in evecs:  
PRINT, evecs[*,0]  

IDL prints:

(  0.68168704,  0.18789033)( -0.34084352, -0.093945164)  
(  0.16271780, -0.59035830)  
PRINT, evecs[*,1]   

IDL prints:

(  0.18789033,  0.68168704)( -0.093945164, -0.34084352)  
( -0.59035830,  0.16271780)  
PRINT, evecs[*,2]   

IDL prints:

(  0.70710678,  0.0000000)(  0.70710678,  0.0000000)  
( -2.3570226e-21, 0.0000000)  

We can check the accuracy of these results using the relation Ax - lx = 0. The array contained in the variable specified by the RESIDUAL keyword contains the result of this computation.

PRINT, residual  

IDL prints:

( -1.2021898e-07,  1.1893681e-07)(  6.0109490e-08, -5.9468404e-08)  
(  1.0300230e-07,  1.0411269e-07)  
(  1.1893681e-07, -1.2021898e-07)( -5.9468404e-08,  6.0109490e-08)  
(  1.0411269e-07,  1.0300230e-07)  
(      0.0000000,      0.0000000)(      0.0000000,      0.0000000)  

The results are all zero to within machine precision.

Repeated Eigenvalues

To compute the eigenvalues and eigenvectors of a real, nonsymmetric n-by-n array, begin with an array A. In this example, there are fewer than n distinct eigenvalues, but n independent eigenvectors are available.

A = [[8.0, 0.0, 3.0], $  
     [2.0, 2.0, 1.0], $  
     [2.0, 0.0, 3.0]]   
; Reduce A to upper Hessenberg form and compute the eigenvalues.  
; Note that both operations can be combined into a single command.  
evals = HQR(ELMHES(A))  
; Print the eigenvalues:  
PRINT, evals  

IDL prints:

(  9.00000,  0.00000) (  2.00000,  0.00000)   
(  2.00000,  0.00000)   

; Initialize a variable to contain the residual:  
residual = 1  
; Compute the eigenvectors and residual, using  
; double-precision arithmetic:  
evecs = EIGENVEC(A, evals, /DOUBLE, RESIDUAL=residual)  
; Print the eigenvectors:  
PRINT, evecs[*,0]   

IDL prints:

(  0.90453403,  0.0000000)(  0.30151134,  0.0000000)  
(  0.30151134,  0.0000000)  
PRINT, evecs[*,1]   

IDL prints:

( -0.27907279,  0.0000000)( -0.78140380,  0.0000000)  
(  0.55814557,  0.0000000)  
PRINT, evecs[*,2]   

IDL prints:

( -0.27907279,  0.0000000)( -0.78140380,  0.0000000)  
(  0.55814557,  0.0000000)  

We can compute an independent eigenvector for the repeated eigenvalue (2.0) by perturbing it slightly, allowing the algorithm EIGENVEC to recognize the eigenvalue as distinct and to compute a linearly-independent eigenvector.

newresidual = 1  
evecs[*,2] = EIGENVEC(A, evals[2]+1.0e-6, /DOUBLE, $  
   RESIDUAL = newresidual)   
PRINT, evecs[*,2]  

IDL prints:

( -0.33333333,  0.0000000)(  0.66666667,  0.0000000)  
(  0.66666667,  0.0000000)  

Once again, we can check the accuracy of these results by checking that each element in the residuals -for both the original eigenvectors and the perturbed eigenvector- is zero to within machine precision.

The So-called Defective Case

In the so-called defective case, there are fewer than n distinct eigenvalues and fewer than n linearly-independent eigenvectors. Begin with an array A:

A = [[2.0, -1.0], $  
     [1.0,  0.0]]   
; Reduce A to upper Hessenberg form and compute the eigenvalues.  
; Note that both operations can be combined into a single command.  
evals = HQR(ELMHES(A))  
; Print the eigenvalues:  
PRINT, evals  

IDL prints:

(  1.00000,  0.00000)(  1.00000,  0.00000)  

; Compute the eigenvectors, using double-precision arithmetic:  
evecs = EIGENVEC(A, evals, /DOUBLE)  
; Print the eigenvectors:  
PRINT, evecs[*,0]   

IDL prints:

(  0.70710678,  0.0000000)(  0.70710678,  0.0000000)  
PRINT, evecs[*,1]   

IDL prints:

(  0.70710678,  0.0000000)(  0.70710678,  0.0000000)  

We attempt to compute an independent eigenvector using the method described in the previous example:

evecs[*,1] = EIGENVEC(A, evals[1]+1.0e-6, /DOUBLE)   
PRINT, evecs[1,*]  

IDL prints:

(  0.70710678,  0.0000000)(  0.70710678,  0.0000000)  

In this example, n independent eigenvectors do not exist. This situation is termed the defective case and cannot be resolved analytically or numerically.

Routines for Computing Eigenvalues and Eigenvectors

See "Eigenvalues and Eigenvectors" (in the functional category Mathematics) for a brief description of IDL routines for computing eigenvalues and eigenvectors.Detailed information is available in the IDL Reference Guide.

  IDL Online Help (June 16, 2005)