Eigenvalue algorithm

Eigenvalue algorithm

In linear algebra, one of the most important problems is designing efficient and stable algorithms for finding the eigenvalues of a matrix. These eigenvalue algorithms may also find eigenvectors.


Characteristic polynomial

Given a square matrix A, an eigenvalue λ and its associated eigenvector v are, by definition, a pair obeying the relation

A{\bold v} = \lambda{\bold v} \,\!,

where v is nonzero. Equivalently, (A−λI)v  =  0 (where I is the identity matrix), implying det(A−λI)  =  0. This determinant is a polynomial in λ, known as the characteristic polynomial of A. One common method for determining the eigenvalues of a small matrix is by finding roots of the characteristic polynomial.

Unfortunately, this method has some limitations. A general polynomial of order n > 4 cannot be solved by a finite sequence of arithmetic operations and radicals (see Abel–Ruffini theorem). There do exist efficient root-finding algorithms for higher order polynomials. However, finding the roots of the characteristic polynomial may be an ill-conditioned problem even when the underlying eigenvalue problem is well-conditioned. For this reason, this method is rarely used.

The above discussion implies a restriction on all eigenvalue algorithms. It can be shown that for any polynomial, there exists a matrix (see companion matrix) having that polynomial as its characteristic polynomial (actually, there are infinitely many). If there did exist a finite sequence of arithmetic operations for exactly finding the eigenvalues of a general matrix, this would provide a corresponding finite sequence for general polynomials, in contradiction of the Abel–Ruffini theorem. Therefore, general eigenvalue algorithms are expected to be iterative.

Power iteration

Let the eigenvalues of A be \lambda_1,\lambda_2, \dots, \lambda_k. Assume that λ1 has absolute value strictly larger than that of \lambda_2, \dots, \lambda_k. This is an essential restriction: for a matrix with real coefficients, if the eigenvalue with highest absolute value is not real, its complex conjugate is always an eigenvalue with the same absolute value.

With the preceding assumptions in mind, the idea of the method is to choose an (arbitrary) unit length vector b and then repeatedly multiply it by the matrix A and re-scale. One carries out the computation b_1:=\frac{A b}{||Ab||}, b_2:=\frac{A b_1}{||Ab_1||}, b_3:=\frac{A b_2}{||Ab_2||}, .... Let the vector b have a generalized eigenspace decomposition b=v_1+v_2+\dots+v_k, where v1 belongs to the generalized eigenspace corresponding to eigenvalue λ1, v2 belongs to the generalized eigenspace corresponding to eigenvalue λ2, etc. Each iteration of the algorithm will decrease the "contribution" of the components in the eigenspaces of eigenvalues \lambda_2,\dots, \lambda_k relative to the contribution of λ1 and therefore the vector bi will converge to a unit eigenvector of eigenvalue λ1.

The power iteration algorithm for finding the (largest) eigenvalue is simple to implement but is otherwise not very useful in practice. Its convergence is slow except for special cases of matrices. Without modification, it can only find the dominant eigenvalue (and the corresponding eigenvector), provided they exist.

A few of the more advanced eigenvalue algorithms are variations of power iteration. In addition, some of the better algorithms for the generalized eigenvalue problem are based on power iteration.

Matrix eigenvalues

In mathematics, and in particular in linear algebra, an important tool for describing eigenvalues of square matrices is the characteristic polynomial: saying that λ is an eigenvalue of A is equivalent to stating that the system of linear equations (A - λI) v = 0 (where I is the identity matrix) has a non-zero solution v (namely an eigenvector), and so it is equivalent to the determinant det (A - λI) being zero. The function p(λ) = det (A - λI) is a polynomial in λ since determinants are defined as sums of products. This is the characteristic polynomial of A: the eigenvalues of a matrix are the zeros of its characteristic polynomial.

It follows that we can compute all the eigenvalues of a matrix A by solving the equation pA(λ) = 0. If A is an n-by-n matrix, then pA has degree n and A can therefore have at most n eigenvalues. Conversely, the fundamental theorem of algebra says that this equation has exactly n roots (zeroes), counted with multiplicity. All real polynomials of odd degree have a real number as a root, so for odd n, every real matrix has at least one real eigenvalue. In the case of a real matrix, for even and odd n, the non-real eigenvalues come in conjugate pairs.

An example of a matrix with no real eigenvalues is the 90-degree rotation

\begin{bmatrix}0 & -1\\ 1 & 0\end{bmatrix}

whose characteristic polynomial is x2 + 1 and so its eigenvalues are the pair of complex conjugates i, -i.

The Cayley–Hamilton theorem states that every square matrix satisfies its own characteristic polynomial, that is, pA(A) = 0.


Eigenvalues of 2×2 matrices

An analytic solution for the eigenvalues of 2×2 matrices can be obtained directly from the quadratic formula: if

A = \begin{bmatrix} a  & b \\ c & d \end{bmatrix}

then the characteristic polynomial is

\rm det \begin{bmatrix} a-\lambda & b \\ c & d-\lambda \end{bmatrix}=(a-\lambda)(d-\lambda)-bc=\lambda^2-(a+d)\lambda+(ad-bc)

so the solutions are

 \lambda = \frac{a + d}{2}  \pm \sqrt{\frac{(a + d)^2}{4} + bc - ad} = \frac{a + d}{2}  \pm \frac{\sqrt{4bc + (a - d)^2  }}{2}

Notice that the characteristic polynomial of a 2×2 matrix can be written in terms of the trace tr(A) = a + d and determinant det(A) = adbc as

{\rm det} \begin{bmatrix} a-\lambda & b \\ c & d-\lambda \end{bmatrix}
  = {\rm det} \left[ A - \lambda I_{2}\right] 
  = \lambda^2- \lambda {\rm tr}(A)+ {\rm det}(A)

where I2 is the 2×2 identity matrix. The solutions for the eigenvalues of a 2×2 matrix can thus be written as

\lambda = \frac{1}{2} \left({\rm tr}(A) \pm \sqrt{{\rm tr}^2 (A) - 4 {\rm det}(A)} \right)

Thus, for the very special case where the 2×2 matrix has zero determinant, but non-zero trace, the eigenvalues are zero and the trace (corresponding to the negative and positive roots, respectively). For example, the eigenvalues of the following matrix are 0 and (a2 + b2):

\begin{bmatrix} a^2  & a b \\ a b & b^2 \end{bmatrix}

This formula holds for only a 2×2 matrix.

Eigenvalues of 3×3 matrices


A = \begin{bmatrix} a & b & c \\ d & e & f \\ g & h & i \end{bmatrix}

then the characteristic polynomial of A is

\det \begin{bmatrix} a-\lambda & b & c \\ d & e-\lambda & f \\ g & h & i-\lambda \end{bmatrix}= -\lambda^3 + \lambda^2 ( a + e + i ) + \lambda ( db + gc + fh - ae - ai - ei ) + ( aei - afh - dbi + dch + gbf - gce ).

Alternatively the characteristic polynomial of a 3×3 matrix can be written in terms of the trace tr(A) and determinant det(A) as

{\rm det} \left[ A - \lambda I_{3}\right]
= -\lambda^3 
  + \lambda^2 {\rm tr}(A) 
  + \lambda \frac{1}{2}\left[ {\rm tr}(A^2) - {\rm tr}^2(A) \right]
  + {\rm det}(A)

where I3 is the 3×3 identity matrix.

The eigenvalues of the matrix are the roots of this polynomial, which can be found using the method for solving cubic equations.

A formula for the eigenvalues of a 4×4 matrix could be derived in an analogous way, using the formulae for the solutions of the quartic equation.

A programmatical approach to finding eigenvalues

In case you don’t have access to scientific software that has an eigenvalue function, below are a set of steps you could use as pseudo code to help you create a program to solve problems.

For a 2×2:

//we know there are two eigenvalues, so the equations from above can be hard-coded in with minimal effort
Eig11 = (a+d)/2 + sqrt(((a+d)*(a+d))/4 + bc – ad)
Eig12 = (a+d)/2 - sqrt(((a+d)*(a+d))/4 + bc – ad)
//eig11 is the value found using + and eig12 is the value using -
//or we get the same two eigenvalues with the modified equation
Eig21 = (a+d)/2 + sqrt(4bc + (a-d)(a-d))/2
Eig22 = (a+d)/2 - sqrt(4bc + (a-d)(a-d))/2

For a 3×3:

//a 3×3 is more complicated and requires several helper equations to accomplish due to the λ³ term.
//These steps should help you calculate eigen values for a matrix that has ***3 REAL eigen values***.
Define x,y,z
//Use the equation from above to get your cubic equation and combine all constant terms possible to
//give you a reduced equation we will use a, b, c and d to denote the coefficients of this equation.
//Eqn = aλ³ + bλ² + cλ + d = 0
x = ((3c/a)(/))/3
y = ((2b³/)(9bc/) + (27d/a))/27
z =/4 +/27
Define I, j, k, m, n, p (so equations are not so cluttered)
i = sqrt(/4 - z)
j = -i^(1/3)
k = arccos(-(y/2i))
m = cos(k/3)
n = sqrt(3)*sin(k/3)
p = b/3a
Define Eig1, Eig2, Eig3
Eig1 = 2j*m + p
Eig2 = -j *(m + n) + p
Eig3 = -j*(m - n) + p

Eigenvalues of a Symmetric 3x3 Matrix

(Reference: Oliver K. Smith: Eigenvalues of a symmetric 3 × 3 matrix. Commun. ACM 4(4): 168 (1961) ) Note: this method does not work for singular matrices (matrices with one or more zero eigenvalues). This is a matlab version:

% Given symmetric 3x3 matrix M, compute the eigenvalues
m = trace(M)/3;
K = M-m*eye(3);
q = det(K)/2;
p = 0
for i=1:3
    for j=1:3
        p = p + K(i,j)^2;
p = p/6;
phi = 1/3*acos(q/p^(3/2));
% NOTE: the follow formula assume accurate computation and therefor q/p^(3/2) should be in range of [1,-1], 
% but in real code, because of numerical errors, it must be checked. Thus, in case abs(q) >= abs(p^(3/2)), set phi = 0;
if(abs(q) >= abs(p^(3/2)))
    phi = 0;
eig1 = m + 2*sqrt(p)*cos(phi)
eig2 = m - sqrt(p)*(cos(phi) + sqrt(3)*sin(phi))
eig3 = m - sqrt(p)*(cos(phi) - sqrt(3)*sin(phi))

Or here is a version in python that uses numpy:

def smith_eigenvals(mat):
   """Returns the eigenvalues of a 3x3 symmetric matrix.                       
    Or can use numpy.linalg.eigvals function.
    import numpy
    m = numpy.trace(mat) / 3.
    K = mat - m*numpy.eye(3)
    q = numpy.linalg.det(K) / 2.
    p = numpy.sum(K*K) / 6.
    # NB in Smith's paper he uses phi = (1./3.)*arctan(sqrt(p**3 - q**2)/q), which is equivalent to below:
    phi = (1./3.)*numpy.arccos(q/(p**(3./2.)))
    print "m", m
    print "p", p
    print "q", q
    print "phi", phi
    if abs(q) >= abs(p**(3./2.)):
        phi = 0.
    if phi < 0.:
        phi = phi + numpy.pi/3.
    eig1 = m + 2.*numpy.sqrt(p)*numpy.cos(phi)
    eig2 = m - numpy.sqrt(p)*(numpy.cos(phi) + numpy.sqrt(3.)*numpy.sin(phi))
    eig3 = m - numpy.sqrt(p)*(numpy.cos(phi) - numpy.sqrt(3.)*numpy.sin(phi))
    return eig1, eig2, eig3

Eigenvalues and eigenvectors of special matrices

For matrices satisfying A2 = α one can write explicit formulas for the possible eigenvalues and the projectors on the corresponding eigenspaces.



AP_+=\sqrt{\alpha}P_+ \quad AP_-=-\sqrt{\alpha}P_-


P_+P_+=P_+ \quad P_-P_-=P_- \quad P_+P_-=P_-P_+=0

This provides the following resolution of identity

I=P_+ + P_-= \frac{1}{2}\left(I+\frac{A}{\sqrt{\alpha}}\right) 
+ \frac{1}{2}\left(I-\frac{A}{\sqrt{\alpha}}\right)

The multiplicity of the possible eigenvalues is given by the rank of the projectors.

Example computation

The computation of eigenvalue/eigenvector can be realized with the following algorithm.

Consider an n-square matrix A

1. Find the roots of the characteristic polynomial of A. These are the eigenvalues.
  • If n different roots are found, then the matrix can be diagonalized.
2. Find a basis for the kernel of the matrix given by A − λnI. For each of the eigenvalues. These are the eigenvectors
  • The eigenvectors given from different eigenvalues are linearly independent.
  • The eigenvectors given from a root-multiplicity are also linearly independent.

Let us determine the eigenvalues of the matrix

A = \begin{bmatrix}
     0 & 1 & -1 \\
     1 & 1 &  0 \\
    -1 & 0 &  1 

which represents a linear operator R³ → R³.

Identifying eigenvalues

We first compute the characteristic polynomial of A:

p(\lambda) = \det( A - \lambda I) =
\det \begin{bmatrix}
    -\lambda &    1      &   -1 \\
        1    & 1-\lambda &    0     \\
       -1    &    0      & 1-\lambda
= -\lambda^3 + 2\lambda^2 +\lambda - 2.

This polynomial factors to p(λ) = − (λ − 2)(λ − 1)(λ + 1). Therefore, the eigenvalues of A are 2, 1 and −1.

Identifying eigenvectors

With the eigenvalues in hand, we can solve sets of simultaneous linear equations to determine the corresponding eigenvectors. Since we are solving for the system (A - \lambda I)v = 0\,, if λ = 2 then,

       -2    & 1   &   -1 \\
        1    & -1  &    0 \\
       -1    & 0   &   -1
       v_1 \\
       v_2 \\
\end{bmatrix} = 0.

Now, reducing (A - 2I)\, to row echelon form:

       -2    & 1   &   -1 \\
        1    & -1  &    0 \\
       -1    & 0   &   -1
\end{bmatrix} \to
        1 & 0 & 1 \\
        0 & 1 & 1 \\
        0 & 0 & 0

allows us to solve easily for the eigenspace E2:

        1 & 0 & 1 \\
        0 & 1 & 1 \\
        0 & 0 & 0
       v_1 \\
       v_2 \\
\end{bmatrix} = 0 \to
    v_1 + v_3 = 0 \\
    v_2 + v_3 = 0
\to E_2 = {\rm span}\begin{bmatrix}1 \\ 1 \\ -1\end{bmatrix}.

We can confirm that a simple example vector chosen from eigenspace E2 is a valid eigenvector with eigenvalue λ = 2:

A \begin{bmatrix} \; 1 \\ \; 1 \\ -1 \end{bmatrix} 
= \begin{bmatrix} \; 2 \\ \; 2 \\ -2 \end{bmatrix} 
= 2 \begin{bmatrix} \; 1 \\ \; 1 \\ -1 \end{bmatrix}.

Note that we can determine the degrees of freedom of the solution by the number of pivots.

If A is a real matrix, the characteristic polynomial will have real coefficients, but its roots will not necessarily all be real. The complex eigenvalues come in pairs which are conjugates. For a real matrix, the eigenvectors of a non-real eigenvalue z , which are the solutions of (AzI)v = 0, cannot be real.

If v1, ..., vm are eigenvectors with different eigenvalues λ1, ..., λm, then the vectors v1, ..., vm are necessarily linearly independent.

The spectral theorem for symmetric matrices states that if A is a real symmetric n-by-n matrix, then all its eigenvalues are real, and there exist n linearly independent eigenvectors for A which are mutually orthogonal. Symmetric matrices are commonly encountered in engineering.

Our example matrix from above is symmetric, and three mutually orthogonal eigenvectors of A are

v_1 = \begin{bmatrix}\; 1  \\ \;1 \\   -1 \end{bmatrix},\quad v_2 = \begin{bmatrix}\; 0\;\\   1 \\    1 \end{bmatrix},\quad v_3 = \begin{bmatrix}\; 2  \\  -1 \\ \; 1 \end{bmatrix}.

These three vectors form a basis of R³. With respect to this basis, the linear map represented by A takes a particularly simple form: every vector x in R³ can be written uniquely as

x = x1v1 + x2v2 + x3v3

and then we have

Ax = 2x1v1 + x2v2x3v3.

Advanced methods

A popular method for finding eigenvalues is the QR algorithm, which is based on the QR decomposition. Besides, combining Householder transformation with LU decomposition can get better convergence than QR algorithm. [1] Other advanced methods include:

Most eigenvalue algorithms rely on first reducing the matrix A to Hessenberg or tridiagonal form. This is usually accomplished via reflections.

See also


Wikimedia Foundation. 2010.

Игры ⚽ Нужно сделать НИР?

Look at other dictionaries:

  • Jacobi eigenvalue algorithm — The Jacobi eigenvalue algorithm is a numerical procedure for the calculation of all eigenvalues and eigenvectors of a real symmetric matrix. Description Let varphi in mathbb{R}, , 1 le k < l le n and let J(varphi, k, l) denote the n imes n matrix …   Wikipedia

  • Divide-and-conquer eigenvalue algorithm — Divide and conquer eigenvalue algorithms are a class of eigenvalue algorithms for Hermitian or real symmetric matrices that have recently (circa 1990s) become competitive in terms of stability and efficiency with more traditional algorithms such… …   Wikipedia

  • Eigenvalue, eigenvector and eigenspace — In mathematics, given a linear transformation, an Audio|De eigenvector.ogg|eigenvector of that linear transformation is a nonzero vector which, when that transformation is applied to it, changes in length, but not direction. For each eigenvector… …   Wikipedia

  • QR algorithm — In numerical linear algebra, the QR algorithm is an eigenvalue algorithm; that is, a procedure to calculate the eigenvalues and eigenvectors of a matrix. The QR transformation was developed in 1961 by John G.F. Francis (England) and by Vera N.… …   Wikipedia

  • Root-finding algorithm — A root finding algorithm is a numerical method, or algorithm, for finding a value x such that f(x) = 0, for a given function f. Such an x is called a root of the function f. This article is concerned with finding scalar, real or complex roots,… …   Wikipedia

  • Lanczos algorithm — The Lanczos algorithm is an iterative algorithm invented by Cornelius Lanczos that is an adaptation of power methods to find eigenvalues and eigenvectors of a square matrix or the singular value decomposition of a rectangular matrix. It is… …   Wikipedia

  • Grover's algorithm — is a quantum algorithm for searching an unsorted database with N entries in O(N1/2) time and using O( log N) storage space (see big O notation). It was invented by Lov Grover in 1996.Classically, searching an unsorted database requires a linear… …   Wikipedia

  • Block Lanczos algorithm for nullspace of a matrix over a finite field — The Block Lanczos algorithm for nullspace of a matrix over a finite field is a procedure for finding the nullspace of a matrix using only multiplication of the matrix by long, thin matrices. These long, thin matrices are considered as vectors of… …   Wikipedia

  • Eigenvalues and eigenvectors — For more specific information regarding the eigenvalues and eigenvectors of matrices, see Eigendecomposition of a matrix. In this shear mapping the red arrow changes direction but the blue arrow does not. Therefore the blue arrow is an… …   Wikipedia

  • Singular value decomposition — Visualization of the SVD of a 2 dimensional, real shearing matrix M. First, we see the unit disc in blue together with the two canonical unit vectors. We then see the action of M, which distorts the disk to an ellipse. The SVD decomposes M into… …   Wikipedia

Share the article and excerpts

Direct link
Do a right-click on the link above
and select “Copy Link”