Eigenvalue algorithm

The adjoint M^* of a complex matrix M is the transpose of the conjugate of M: M ^* = M ^T. A square matrix A is called normal if it commutes with its adjoint: A^*A = AA^*. It is called Hermitian if it is equal to its adjoint: A^* = A. All Hermitian matrices are normal. If A has only real elements, then the adjoint is just the transpose, and A is Hermitian if and only if it is symmetric. When applied to column vectors, the adjoint can be used to define the canonical inner product on Cⁿ: 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 matrix A, Cⁿ has an orthonormal basis consisting of eigenvectors of A. The corresponding matrix of eigenvectors is unitary.
• The eigenvalues of a Hermitian matrix are real, since (λ − λ)v = (A^* − A)v = (A − A)v = 0 for a non-zero eigenvector v.
• If A is real, there is an orthonormal basis for Rⁿ consisting of eigenvectors of A if and only if A is symmetric.

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

Since the determinant of a triangular matrix is the product of its diagonal entries, if T is triangular, then ${\textstyle \det(\lambda I-T)=\prod _{i}(\lambda -T_{ii})}$. Thus the eigenvalues of T are its diagonal entries.

If p is any polynomial and p(A) = 0, then the eigenvalues of A also satisfy the same equation. If p happens to have a known factorization, then the eigenvalues of A lie among its roots.

For example, a projection is a square matrix P satisfying P² = 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 the nullity of P, while the multiplicity of 1 is the rank of P.

Another example is a matrix A that satisfies A² = α²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 AP_{+}=\alpha P_{+}\quad AP_{-}=-\alpha P_{-}}$

and

${\displaystyle P_{+}P_{+}=P_{+}\quad P_{-}P_{-}=P_{-}\quad P_{+}P_{-}=P_{-}P_{+}=0.}$

The column spaces of P₊ and P₋ are the eigenspaces of A 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 the root formulas makes this approach less attractive.

For the 2×2 matrix

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

the characteristic polynomial is

${\displaystyle \det {\begin{bmatrix}\lambda -a&-b\\-c&\lambda -d\end{bmatrix}}=\lambda ^{2}\,-\,\left(a+d\right)\lambda \,+\,\left(ad-bc\right)=\lambda ^{2}\,-\,\lambda \,{\rm {tr}}(A)\,+\,\det(A).}$

Thus the eigenvalues can be found by using the quadratic formula:

${\displaystyle \lambda ={\frac {{\rm {tr}}(A)\pm {\sqrt {{\rm {tr}}^{2}(A)-4\det(A)}}}{2}}.}$

Defining ${\textstyle {\rm {gap}}\left(A\right)={\sqrt {{\rm {tr}}^{2}(A)-4\det(A)}}}$ to be the distance between the two eigenvalues, it is straightforward to calculate

${\displaystyle {\frac {\partial \lambda }{\partial a}}={\frac {1}{2}}\left(1\pm {\frac {a-d}{{\rm {gap}}(A)}}\right),\qquad {\frac {\partial \lambda }{\partial b}}={\frac {\pm c}{{\rm {gap}}(A)}}}$

with similar formulas for c and d. From this it follows that the calculation is well-conditioned if the eigenvalues are isolated.

Eigenvectors can be found by exploiting the Cayley–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, then A is a multiple of the identity and any non-zero vector is an eigenvector.)

For example, suppose

${\displaystyle A={\begin{bmatrix}4&3\\-2&-3\end{bmatrix}},}$

then tr(A) = 4 − 3 = 1 and det(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,

${\displaystyle A-3I={\begin{bmatrix}1&3\\-2&-6\end{bmatrix}},\qquad A+2I={\begin{bmatrix}6&3\\-2&-1\end{bmatrix}}.}$

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 by A.

The characteristic equation of a symmetric 3×3 matrix A is:

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

This equation may be solved using the methods of Cardano or Lagrange, but an affine change to A will simplify the expression considerably, and lead directly to a trigonometric solution. If A = pB + qI, then A and B have the same eigenvectors, and β is an eigenvalue of B if and only if α = pβ + q is an eigenvalue of A. Letting ${\textstyle q={\rm {tr}}(A)/3}$ and ${\textstyle p=\left({\rm {tr}}\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 identity cos 3θ = 4cos³ θ − 3cos θ reduces the equation to cos 3θ = det(B) / 2. Thus

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

If det(B) is complex or is greater than 2 in absolute value, the arccosine should be taken along the same branch for all three values of k. This issue doesn't arise when A 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 radians

p1 = A(1,2)^2 + A(1,3)^2 + A(2,3)^2
if (p1 == 0) 
   % A is diagonal.
   eig1 = A(1,1)
   eig2 = A(2,2)
   eig3 = A(3,3)
else
   q = trace(A)/3               % trace(A) is the sum of all diagonal values
   p2 = (A(1,1) - q)^2 + (A(2,2) - q)^2 + (A(3,3) - q)^2 + 2 * p1
   p = sqrt(p2 / 6)
   B = (1 / p) * (A - q * I)    % I is the identity matrix
   r = 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 / 3
   elseif (r >= 1)
      phi = 0
   else
      phi = acos(r) / 3
   end

   % the eigenvalues satisfy eig3 <= eig2 <= eig1
   eig1 = q + 2 * p * cos(phi)
   eig3 = q + 2 * p * cos(phi + (2*pi/3))
   eig2 = 3 * q - eig1 - eig3     % since trace(A) = eig1 + eig2 + eig3
end

Once again, the eigenvectors of A can be obtained by recourse to the Cayley–Hamilton theorem. If α₁, α₂, α₃ are distinct eigenvalues of A, 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 the generalized eigenspace of α₁ is spanned by the columns of A − α₂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

${\displaystyle A={\begin{bmatrix}3&2&6\\2&2&5\\-2&-1&-4\end{bmatrix}}.}$

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,

${\displaystyle A-I={\begin{bmatrix}2&2&6\\2&1&5\\-2&-1&-5\end{bmatrix}},\qquad A+I={\begin{bmatrix}4&2&6\\2&3&5\\-2&-1&-3\end{bmatrix}}}$

and

${\displaystyle (A-I)^{2}={\begin{bmatrix}-4&0&-8\\-4&0&-8\\4&0&8\end{bmatrix}},\qquad (A-I)(A+I)={\begin{bmatrix}0&4&4\\0&2&2\\0&-2&-2\end{bmatrix}}}$

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 of A. 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. The cross 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 not 0, 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.