Innumerical analysis, one of the most important problems is designing efficient andstablealgorithms for finding theeigenvalues of amatrix. Theseeigenvalue algorithms may also find eigenvectors.

Given ann ×nsquare matrixA ofreal orcomplex numbers, aneigenvalueλ and its associatedgeneralized eigenvectorv are a pair obeying the relation^[1]

${\displaystyle {\left(A-\lambda I\right)}^k\mathbf{v}=0,}$

wherev is a nonzeron × 1 column vector,I is then ×nidentity matrix,k is a positive integer, and bothλ andv are allowed to be complex even whenA is real. Whenk = 1, the vector is called simply aneigenvector, and the pair is called aneigenpair. In this case,Av =λv. Any eigenvalueλ ofA has ordinary^{[note 1]} eigenvectors associated to it, for ifk is the smallest integer such that(A −λI)^kv = 0 for a generalized eigenvectorv, then(A −λI)^k−1v is an ordinary eigenvector. The valuek can always be taken as less than or equal ton. In particular,(A −λI)ⁿv = 0 for all generalized eigenvectorsv associated withλ.

For each eigenvalueλ ofA, thekernelker(A −λI) consists of all eigenvectors associated withλ (along with 0), called theeigenspace ofλ, while the vector spaceker((A −λI)ⁿ) consists of all generalized eigenvectors, and is called thegeneralized eigenspace. Thegeometric multiplicity ofλ is the dimension of its eigenspace. Thealgebraic multiplicity ofλ is the dimension of its generalized eigenspace. The latter terminology is justified by the equation

${\displaystyle p_A\left(z\right)=det\left(zI-A\right)=\prod _{i=1}^k(z-{\lambda }_i{)}^{{\alpha }_i},}$

wheredet is thedeterminant function, theλ_i are all the distinct eigenvalues ofA and theα_i are the corresponding algebraic multiplicities. The functionp_A(z) is thecharacteristic polynomial ofA. So the algebraic multiplicity is the multiplicity of the eigenvalue as azero of the characteristic polynomial. Since any eigenvector is also a generalized eigenvector, the geometric multiplicity is less than or equal to the algebraic multiplicity. The algebraic multiplicities sum up ton, the degree of the characteristic polynomial. The equationp_A(z) = 0 is called thecharacteristic equation, as its roots are exactly the eigenvalues ofA. By theCayley–Hamilton theorem,A itself obeys the same equation:p_A(A) = 0.^{[note 2]} As a consequence, the columns of the matrix$\prod _{i\ne j}(A-{\lambda }_{i}I{)}^{{\alpha }_i}$ must be either 0 or generalized eigenvectors of the eigenvalueλ_j, since they are annihilated by${\displaystyle (A-{\lambda }_{j}I{)}^{{\alpha }_j}}$. In fact, thecolumn space is the generalized eigenspace ofλ_j.

Any collection of generalized eigenvectors of distinct eigenvalues is linearly independent, so a basis for all ofCⁿ can be chosen consisting of generalized eigenvectors. More particularly, this basis{v_i}ⁿ
_i=1 can be chosen and organized so that

• ifv_i andv_j have the same eigenvalue, then so doesv_k for eachk betweeni andj, and
• ifv_i is not an ordinary eigenvector, and ifλ_i is its eigenvalue, then(A −λ_iI)v_i =v_i−1 (in particular,v₁ must be an ordinary eigenvector).

If these basis vectors are placed as the column vectors of a matrixV = [v₁v₂ ⋯v_n], thenV can be used to convertA to itsJordan normal form:

where theλ_i are the eigenvalues,β_i = 1 if(A −λ_i+1)v_i+1 =v_i andβ_i = 0 otherwise.

More generally, ifW is any invertible matrix, andλ is an eigenvalue ofA with generalized eigenvectorv, then(W⁻¹AW −λI)^kW^−kv = 0. Thusλ is an eigenvalue ofW⁻¹AW with generalized eigenvectorW^−kv. That is,similar matrices have the same eigenvalues.

TheadjointM^* of a complex matrixM is the transpose of the conjugate ofM:M^* =M^T. A square matrixA is callednormal if it commutes with its adjoint:A^*A =AA^*. It is calledHermitian if it is equal to its adjoint:A^* =A. All Hermitian matrices are normal. IfA has only real elements, then the adjoint is just the transpose, andA is Hermitian if and only if it issymmetric. When applied to column vectors, the adjoint can be used to define the canonical inner product onCⁿ:w ⋅v =w^*v.^{[note 3]} Normal, Hermitian, and real-symmetric matrices have several useful properties:

• Every generalized eigenvector of a normal matrix is an ordinary eigenvector.
• Any normal matrix is similar to a diagonal matrix, since its Jordan normal form is diagonal.
• Eigenvectors of distinct eigenvalues of a normal matrix are orthogonal.
• The null space and the image (or column space) of a normal matrix are orthogonal to each other.
• For any normal matrixA,Cⁿ has an orthonormal basis consisting of eigenvectors ofA. The corresponding matrix of eigenvectors isunitary.
• The eigenvalues of a Hermitian matrix are real, since(λ −λ)v = (A^* −A)v = (A −A)v = 0 for a non-zero eigenvectorv.
• IfA is real, there is an orthonormal basis forRⁿ consisting of eigenvectors ofA if and only ifA is symmetric.

It is possible for a real or complex matrix to have all real eigenvalues without being Hermitian. For example, a realtriangular matrix has its eigenvalues along its diagonal, but in general is not symmetric.

Any problem of numeric calculation can be viewed as the evaluation of some functionf for some inputx. Thecondition numberκ(f,x) of the problem is the ratio of the relative error in the function's output to the relative error in the input, and varies with both the function and the input. The condition number describes how error grows during the calculation. Its base-10 logarithm tells how many fewer digits of accuracy exist in the result than existed in the input. The condition number is a best-case scenario. It reflects the instability built into the problem, regardless of how it is solved. No algorithm can ever produce more accurate results than indicated by the condition number, except by chance. However, a poorly designed algorithm may produce significantly worse results. For example, as mentioned below, the problem of finding eigenvalues for normal matrices is always well-conditioned. However, the problem of finding the roots of a polynomial can bevery ill-conditioned. Thus eigenvalue algorithms that work by finding the roots of the characteristic polynomial can be ill-conditioned even when the problem is not.

For the problem of solving the linear equationAv =b whereA is invertible, thematrix condition numberκ(A⁻¹,b) is given by||A||_op||A⁻¹||_op, where|| ||_op is theoperator norm subordinate to the normalEuclidean norm onCⁿ. Since this number is independent ofb and is the same forA andA⁻¹, it is usually just called the condition numberκ(A) of the matrixA. This valueκ(A) is also the absolute value of the ratio of the largestsingular value ofA to its smallest. IfA isunitary, then||A||_op = ||A⁻¹||_op = 1, soκ(A) = 1. For general matrices, the operator norm is often difficult to calculate. For this reason, othermatrix norms are commonly used to estimate the condition number.

For the eigenvalue problem,Bauer and Fike proved that ifλ is an eigenvalue for adiagonalizablen ×n matrixA witheigenvector matrixV, then the absolute error in calculatingλ is bounded by the product ofκ(V) and the absolute error inA.^[2]As a result, the condition number for findingλ isκ(λ,A) =κ(V) = ||V ||_op ||V⁻¹||_op. IfA is normal, thenV is unitary, andκ(λ,A) = 1. Thus the eigenvalue problem for all normal matrices is well-conditioned.

The condition number for the problem of finding the eigenspace of a normal matrixA corresponding to an eigenvalueλ has been shown to be inversely proportional to the minimum distance betweenλ and the other distinct eigenvalues ofA.^[3] In particular, the eigenspace problem for normal matrices is well-conditioned for isolated eigenvalues. When eigenvalues are not isolated, the best that can be hoped for is to identify the span of all eigenvectors of nearby eigenvalues.

The most reliable and most widely used algorithm for computing eigenvalues isJohn G. F. Francis' andVera N. Kublanovskaya'sQR algorithm, considered one of the top ten algorithms of 20th century.^[4]

Any monic polynomial is the characteristic polynomial of itscompanion matrix. Therefore, a general algorithm for finding eigenvalues could also be used to find the roots of polynomials. TheAbel–Ruffini theorem shows that any such algorithm for dimensions greater than 4 must either be infinite, or involve functions of greater complexity than elementary arithmetic operations and fractional powers. For this reason algorithms that exactly calculate eigenvalues in a finite number of steps only exist for a few special classes of matrices. For general matrices, algorithms areiterative, producing better approximate solutions with each iteration.

Some algorithms produce every eigenvalue, others will produce a few, or only one. However, even the latter algorithms can be used to find all eigenvalues. Once an eigenvalueλ of a matrixA has been identified, it can be used to either direct the algorithm towards a different solution next time, or to reduce the problem to one that no longer hasλ as a solution.

Redirection is usually accomplished by shifting: replacingA withA −μI for some constantμ. The eigenvalue found forA −μI must haveμ added back in to get an eigenvalue forA. For example, forpower iteration,μ =λ. Power iteration finds the largest eigenvalue in absolute value, so even whenλ is only an approximate eigenvalue, power iteration is unlikely to find it a second time. Conversely,inverse iteration based methods find the lowest eigenvalue, soμ is chosen well away fromλ and hopefully closer to some other eigenvalue.

Reduction can be accomplished by restrictingA to the column space of the matrixA −λI, whichA carries to itself. SinceA -λI is singular, the column space is of lesser dimension. The eigenvalue algorithm can then be applied to the restricted matrix. This process can be repeated until all eigenvalues are found.

If an eigenvalue algorithm does not produce eigenvectors, a common practice is to use an inverse iteration based algorithm withμ set to a close approximation to the eigenvalue. This will quickly converge to the eigenvector of the closest eigenvalue toμ. For small matrices, an alternative is to look at the column space of the product ofA −λ'I for each of the other eigenvaluesλ'.

A formula for the norm of unit eigenvector components of normal matrices was discovered by Robert Thompson in 1966 and rediscovered independently by several others.^{[5][6][7][8][9]}IfA is an$n\times n$ normal matrix with eigenvaluesλ_i(A) and corresponding unit eigenvectorsv_i whose component entries arev_i,j, letA_j be the$n-1\times n-1$ matrix obtained by removing thei-th row and column fromA, and letλ_k(A_j) be itsk-th eigenvalue. Then$${\displaystyle |v_{i,j}{|}^2\prod _{k=1,k\ne i}^n({\lambda }_i(A)-{\lambda }_k(A))=\prod _{k=1}^{n-1}({\lambda }_i(A)-{\lambda }_k(A_j))}$$

If${\displaystyle p,p_j}$ are the characteristic polynomials of${\displaystyle A}$ and${\displaystyle A_j}$, the formula can be re-written as$${\displaystyle |v_{i,j}{|}^2=\frac{p_j({\lambda }_i(A))}{p^{\prime }({\lambda }_i(A))}}$$assuming the derivative${\displaystyle p^{\prime }}$ is not zero at${\displaystyle {\lambda }_i(A)}$.

Because the eigenvalues of a triangular matrix are its diagonal elements, for general matrices there is no finite method likegaussian elimination to convert a matrix to triangular form while preserving eigenvalues. But it is possible to reach something close to triangular. Anupper Hessenberg matrix is a square matrix for which all entries below thesubdiagonal are zero. A lower Hessenberg matrix is one for which all entries above thesuperdiagonal are zero. Matrices that are both upper and lower Hessenberg aretridiagonal. Hessenberg and tridiagonal matrices are the starting points for many eigenvalue algorithms because the zero entries reduce the complexity of the problem. Several methods are commonly used to convert a general matrix into a Hessenberg matrix with the same eigenvalues. If the original matrix was symmetric or Hermitian, then the resulting matrix will be tridiagonal.

When only eigenvalues are needed, there is no need to calculate the similarity matrix, as the transformed matrix has the same eigenvalues. If eigenvectors are needed as well, the similarity matrix may be needed to transform the eigenvectors of the Hessenberg matrix back into eigenvectors of the original matrix.

For symmetric tridiagonal eigenvalue problems all eigenvalues (without eigenvectors) can be computed numerically in time O(n log(n)), using bisection on the characteristic polynomial.^[11]

Iterative algorithms solve the eigenvalue problem by producing sequences that converge to the eigenvalues. Some algorithms also produce sequences of vectors that converge to the eigenvectors. Most commonly, the eigenvalue sequences are expressed as sequences of similar matrices which converge to a triangular or diagonal form, allowing the eigenvalues to be read easily. The eigenvector sequences are expressed as the corresponding similarity matrices.

While there is no simple algorithm to directly calculate eigenvalues for general matrices, there are numerous special classes of matrices where eigenvalues can be directly calculated. These include:

Since the determinant of atriangular matrix is the product of its diagonal entries, ifT is triangular, then$det(\lambda I-T)=\prod _i(\lambda -T_{ii})$. Thus the eigenvalues ofT are its diagonal entries.

Ifp is any polynomial andp(A) = 0, then the eigenvalues ofA also satisfy the same equation. Ifp happens to have a known factorization, then the eigenvalues ofA lie among its roots.

For example, aprojection is a square matrixP satisfyingP² =P. The roots of the corresponding scalar polynomial equation,λ² =λ, are 0 and 1. Thus any projection has 0 and 1 for its eigenvalues. The multiplicity of 0 as an eigenvalue is thenullity ofP, while the multiplicity of 1 is the rank ofP.

Another example is a matrixA that satisfiesA² =α²I for some scalarα. The eigenvalues must be±α. The projection operators

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

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

satisfy

${\displaystyle A{P}_{+}=\alpha P_{+}A{P}_{-}=-\alpha P_{-}}$

and

${\displaystyle P_{+}{P}_{+}=P_{+}{P}_{-}{P}_{-}=P_{-}{P}_{+}{P}_{-}=P_{-}{P}_{+}=0.}$

Thecolumn spaces ofP₊ andP₋ are the eigenspaces ofA corresponding to+α and−α, respectively.

For dimensions 2 through 4, formulas involving radicals exist that can be used to find the eigenvalues. While a common practice for 2×2 and 3×3 matrices, for 4×4 matrices the increasing complexity of theroot formulas makes this approach less attractive.

For the 2×2 matrix

the characteristic polynomial is

Thus the eigenvalues can be found by using thequadratic formula:

${\displaystyle \lambda =\frac{\mathrm{t}\mathrm{r}(A)\pm \sqrt{{\mathrm{t}\mathrm{r}}^2(A)-4det(A)}}{2}.}$

Defining$\mathrm{g}\mathrm{a}\mathrm{p}\left(A\right)=\sqrt{{\mathrm{t}\mathrm{r}}^2(A)-4det(A)}$ to be the distance between the two eigenvalues, it is straightforward to calculate

${\displaystyle \frac{\mathrm{\partial }\lambda }{\mathrm{\partial }a}=\frac{1}{2}\left(1\pm \frac{a-d}{\mathrm{g}\mathrm{a}\mathrm{p}(A)}\right),\frac{\mathrm{\partial }\lambda }{\mathrm{\partial }b}=\frac{\pm c}{\mathrm{g}\mathrm{a}\mathrm{p}(A)}}$

with similar formulas forc andd. From this it follows that the calculation is well-conditioned if the eigenvalues are isolated.

Eigenvectors can be found by exploiting theCayley–Hamilton theorem. Ifλ₁,λ₂ are the eigenvalues, then(A −λ₁I)(A −λ₂I) = (A −λ₂I)(A −λ₁I) = 0, so the columns of(A −λ₂I) are annihilated by(A −λ₁I) and vice versa. Assuming neither matrix is zero, the columns of each must include eigenvectors for the other eigenvalue. (If either matrix is zero, thenA is a multiple of the identity and any non-zero vector is an eigenvector.)

For example, suppose

thentr(A) = 4 − 3 = 1 anddet(A) = 4(−3) − 3(−2) = −6, so the characteristic equation is

${\displaystyle 0={\lambda }^2-\lambda -6=(\lambda -3)(\lambda +2),}$

and the eigenvalues are 3 and -2. Now,

In both matrices, the columns are multiples of each other, so either column can be used. Thus,(1, −2) can be taken as an eigenvector associated with the eigenvalue -2, and(3, −1) as an eigenvector associated with the eigenvalue 3, as can be verified by multiplying them byA.

The characteristic equation of a symmetric 3×3 matrixA is:

${\displaystyle det\left(\alpha I-A\right)={\alpha }^3-{\alpha }^2\mathrm{t}\mathrm{r}(A)-\alpha \frac{1}{2}\left(\mathrm{t}\mathrm{r}(A^2)-{\mathrm{t}\mathrm{r}}^2(A)\right)-det(A)=0.}$

This equation may be solved using the methods ofCardano orLagrange, but an affine change toA will simplify the expression considerably, and lead directly to atrigonometric solution. IfA =pB +qI, thenA andB have the same eigenvectors, andβ is an eigenvalue ofB if and only ifα =pβ +q is an eigenvalue ofA. Letting$q=\mathrm{t}\mathrm{r}(A)/3$ and$p={\left(\mathrm{t}\mathrm{r}\left((A-qI{)}^2\right)/6\right)}^{1/2}$, gives

${\displaystyle det\left(\beta I-B\right)={\beta }^3-3\beta -det(B)=0.}$

The substitutionβ = 2cosθ and some simplification using the identitycos 3θ = 4cos³θ − 3cosθ reduces the equation tocos 3θ = det(B) / 2. Thus

${\displaystyle \beta =2\cos \left(\frac{1}{3}\arccos \left(det(B)/2\right)+\frac{2k\pi }{3}\right),k=0,1,2.}$

Ifdet(B) is complex or is greater than 2 in absolute value, the arccosine should be taken along the same branch for all three values ofk. This issue doesn't arise whenA is real and symmetric, resulting in a simple algorithm:^[17]

% Given a real symmetric 3x3 matrix A, compute the eigenvalues% Note that acos and cos operate on angles in radiansp1=A(1,2)^2+A(1,3)^2+A(2,3)^2if(p1==0)% A is diagonal.eig1=A(1,1)eig2=A(2,2)eig3=A(3,3)elseq=trace(A)/3% trace(A) is the sum of all diagonal valuesp2=(A(1,1)-q)^2+(A(2,2)-q)^2+(A(3,3)-q)^2+2*p1p=sqrt(p2/6)B=(1/p)*(A-q*I)% I is the identity matrixr=det(B)/2% In exact arithmetic for a symmetric matrix  -1 <= r <= 1% but computation error can leave it slightly outside this range.if(r<=-1)phi=pi/3elseif(r>=1)phi=0elsephi=acos(r)/3end% the eigenvalues satisfy eig3 <= eig2 <= eig1eig1=q+2*p*cos(phi)eig3=q+2*p*cos(phi+(2*pi/3))eig2=3*q-eig1-eig3% since trace(A) = eig1 + eig2 + eig3end

Once again, the eigenvectors ofA can be obtained by recourse to theCayley–Hamilton theorem. Ifα₁,α₂,α₃ are distinct eigenvalues ofA, then(A −α₁I)(A −α₂I)(A −α₃I) = 0. Thus the columns of the product of any two of these matrices will contain an eigenvector for the third eigenvalue. However, ifα₃ =α₁, then(A −α₁I)²(A −α₂I) = 0 and(A −α₂I)(A −α₁I)² = 0. Thus thegeneralized eigenspace ofα₁ is spanned by the columns ofA −α₂I while the ordinary eigenspace is spanned by the columns of(A −α₁I)(A −α₂I). The ordinary eigenspace ofα₂ is spanned by the columns of(A −α₁I)².

For example, let

The characteristic equation is

${\displaystyle 0={\lambda }^3-{\lambda }^2-\lambda +1=(\lambda -1{)}^2(\lambda +1),}$

with eigenvalues 1 (of multiplicity 2) and -1. Calculating,

and

Thus(−4, −4, 4) is an eigenvector for −1, and(4, 2, −2) is an eigenvector for 1.(2, 3, −1) and(6, 5, −3) are both generalized eigenvectors associated with 1, either one of which could be combined with(−4, −4, 4) and(4, 2, −2) to form a basis of generalized eigenvectors ofA. Once found, the eigenvectors can be normalized if needed.

If a 3×3 matrix${\displaystyle A}$ is normal, then the cross-product can be used to find eigenvectors. If${\displaystyle \lambda }$ is an eigenvalue of${\displaystyle A}$, then the null space of${\displaystyle A-\lambda I}$ is perpendicular to its column space. Thecross product of two independent columns of${\displaystyle A-\lambda I}$ will be in the null space. That is, it will be an eigenvector associated with${\displaystyle \lambda }$. Since the column space is two dimensional in this case, the eigenspace must be one dimensional, so any other eigenvector will be parallel to it.

If${\displaystyle A-\lambda I}$ does not contain two independent columns but is not0, the cross-product can still be used. In this case${\displaystyle \lambda }$ is an eigenvalue of multiplicity 2, so any vector perpendicular to the column space will be an eigenvector. Suppose${\displaystyle \mathbf{v}}$ is a non-zero column of${\displaystyle A-\lambda I}$. Choose an arbitrary vector${\displaystyle \mathbf{u}}$ not parallel to${\displaystyle \mathbf{v}}$. Then${\displaystyle \mathbf{v}\times \mathbf{u}}$ and${\displaystyle (\mathbf{v}\times \mathbf{u})\times \mathbf{v}}$ will be perpendicular to${\displaystyle \mathbf{v}}$ and thus will be eigenvectors of${\displaystyle \lambda }$.

This does not work when${\displaystyle A}$ is not normal, as the null space and column space do not need to be perpendicular for such matrices.

• List of eigenvalue algorithms

• Bojanczyk, Adam W.; Adam Lutoborski (Jan 1991)."Computation of the Euler angles of a symmetric 3X3 matrix".SIAM Journal on Matrix Analysis and Applications.12 (1):41–48.doi:10.1137/0612005.

• Numerical linear algebra

• Articles with short description
• Short description is different from Wikidata