6.2: Numerical Eigensolvers

As discussed above, Abel's impossibility theory tells us that there is no general algebraic formula for calculating the eigenvalues of an \(N\times N\) matrix, for \(N\ge 5\). In practice, however, there exist numerical methods, called eigensolvers, which can compute eigenvalues (and eigenvectors) even for very large matrices, with hundreds of rows/columns, or larger. How could this be?

The answer is that numerical eigensolvers are approximate, not exact. But even though their results are not exact, they are very precise—they can approach the exact eigenvalues to within the fundamental precision limits of floating-point arithmetic. Since we're limited by those floating-point precision limits anyway, that's good enough!

6.2.1 Sketch of the Eigensolver Method

We will not go into detail about how numerical eigensolvers work, as that would involve a rather long digression. For those interested, the following paper provides a good pedagogical presentation of the subject:

• Bruno Lang, Direct Solvers for Symmetric Eigenvalue Problems. Modern Methods and Algorithms of Quantum Chemistry, Proceedings, Second Ed. (2000). PDF download link

Here is a very brief sketch of the basic method. Similar to Gaussian elimination, the algorithm contains two phases, a relatively costly/slow initial phase and a relatively fast second phase. The first phase, which is called Householder reduction, applies a carefully-chosen set of similarity transformations to the input matrix \(\mathbf{A}_0\):

\[\mathbf{A}_0 \;\rightarrow \;\mathbf{A}_1 = \mathbf{X}_1^{-1} \mathbf{A}_0 \mathbf{X}_1 \;\rightarrow \; \mathbf{A}_2 = \mathbf{X}_2^{-1} \mathbf{A}_1 \mathbf{X}_2, \;\; \mathrm{etc.}\]

The end result is a matrix \(\mathbf{A}_k\) which is in Hessenberg form: the elements below the first subdiagonal are all zero (the elements immediately below the main diagonal, i.e. along the first subdiagonal, are allowed to be nonzero). The entire Householder reduction phase requires \(O(N^{3})\) arithmetic operations, where \(N\) is the size of the matrix.

The second phase of the algorithm is called QR iteration. Using a different type of similarity transformation, the elements along the subdiagonal of the Hessenberg matrix are reduced in magnitude. When these elements become negligible, the matrix becomes upper-triangular; in that case, the eigenvalues are simply the elements along the diagonal.

The QR process is iterative, in that it progressively reduces the magnitude of the matrix elements along the subdiagonal. Formally, an infinite number of iterations would be required to reduce these elements to zero—that's why Abel's impossibility theorem isn't violated! In practice, however, QR iteration converges extremely quickly, so this phase ends up taking only \(O(N^{2})\) time.

Hence, the overall runtime for finding the eigenvalues of a matrix scales as \(O(N^{3})\). The eigenvectors can also be computed as a side-effect of the algorithm, with no extra increase in the runtime scaling.

6.2.2 Python Implementation

There are four main numerical eigensolvers implemented in Scipy, which are all found in the scipy.linalg package:

• scipy.linalg.eig returns the eigenvalues and eigenvectors of a matrix.
• scipy.linalg.eigvals returns the eigenvalues (only) of a matrix.
• scipy.linalg.eigh returns the eigenvalues and eigenvectors of a Hermitian matrix.
• scipy.linalg.eigvalsh returns the eigenvalues (only) of a Hermitian matrix.

The reason for having four separate functions is efficiency. The runtimes of all four functions scale as \(O(N^{3})\), but for each \(N\) the actual runtimes of eigvals and eigvalsh will be shorter than eig and eigh, because the eigensolver is only asked to find the eigenvalues and need not construct the eigenvectors. Furthermore, eigvalsh is faster than eigvals, and eigh is faster than eig, because the eigensolver algorithm can make use of certain numerical shortcuts which are valid only for Hermitian matrices.

If you pass eigvalsh or eigh a matrix that is not actually Hermitian, the results are unpredictable; the function may return a wrong value without signaling any error. Therefore, you should only use these functions if you are sure that the input matrix is definitely Hermitian (which is usually because you constructed the matrix that way); if the mtrix is Hermitian, eigvalsh or eigh are certainly preferable to use, because they run faster than their non-Hermitian counterparts.

Here is a short program that uses eigvals to find the eigenvalues of a \(3\times 3\) matrix:

from scipy import *
import scipy.linalg as lin

A = array([[1,3,1],[1, 3, 4],[2, 4, 2]])
lambd = lin.eigvals(A)

print(lambd)

Running the program outputs:

[ 7.45031849+0.j         -0.72515925+0.52865751j -0.72515925-0.52865751j]

The return value of eigvals is a 1D array of complex numbers, storing the eigenvalues of the input. The eigvalsh function behaves similarly, except that a real array is returned (since Hermitian matrices have real eigenvalues). Note: we cannot use lambda as the name of a variable, because lambda is reserved as a special keyword in Python.

Here is an example of using eig:

>>> A = array([[1,3,1],[1, 3, 4],[2, 4, 2]])
>>> lambd, Q = lin.eig(A)
>>> lambd
array([ 7.45031849+0.j    , -0.72515925+0.52865751j,   -0.72515925-0.52865751j])
>>> Q
array([[ 0.40501343+0.j        ,  0.73795979+0.j        ,  0.73795979-0.j        ],
       [ 0.65985810+0.j        , -0.51208724+0.22130102j, -0.51208724-0.22130102j],
       [ 0.63289132+0.j        ,  0.26316357-0.27377508j,  0.26316357+0.27377508j]])

The eig function returns a pair of values; the first is a 1D array of eigenvalues (which we name lambd in the above example), and the second is a 2D array containing the corresponding eigenvectors in each column (which we name Q). For example, the first eigenvector can be accessed with Q[:,0]. We can verify that this is indeed an eigenvector:

>>> dot(A, Q[:,0])
array([ 3.01747903+0.j,  4.91615298+0.j,  4.71524187+0.j])
>>> lambd[0] * Q[:,0]
array([ 3.01747903+0.j,  4.91615298+0.j,  4.71524187+0.j])

The eigh function behaves similarly, except that the 1D array of eigenvalues is real.

6.2.3 Generalized Eigenvalue Problem

Sometimes, you might also come across generalized eigenvalue problems, which have the form

\[\mathbf{A}\, \vec{x} = \lambda\, \mathbf{B}\, \vec{x},\]

for known equal-sized square matrices \(\mathbf{A}\) and \(\mathbf{B}\). We call \(\lambda\) a "generalized eigenvalue", and \(\vec{x}\) a "generalized eigenvector", for the pair of matrices \((\mathbf{A},\mathbf{B})\). The generalized eigenvalue problem reduces to the ordinary eigenvalue problem when \(\mathbf{B}\) is the identity matrix.

The naive way to solve the generalized eigenvalue problem would be to compute the inverse of \(\mathbf{B}^{-1}\), and then solve the eigenvalue problem for \(\mathbf{B}^{-1}\mathbf{A}\). However, it turns out that the generalized eigenvalue problem can be solved directly with only slight modifications to the usual numerical eigensolver algorithm. In fact, the Scipy eigensolvers described in the previous section will solve the generalized eigenvalue problem if you pass a 2D array \(\mathbf{B}\) as the second input (if that second input is omitted, the eigensolvers solve the ordinary eigenvalue problem, as described above).

Here is an example program for solving a generalized eigenvalue problem:

from scipy import *
import scipy.linalg as lin

A = array([[1,3,1],[1, 3, 4],[2, 4, 2]])
B = array([[0,2,1], [0, 1, 1], [2, 0, 1]])

lambd, Q = lin.eig(A, B)

## Verify the solution for first generalized eigenvector:
lhs = dot(A,Q[:,0])             # A . x
rhs = lambd[0] * dot(B, Q[:,0]) # lambda B . x

print(lhs)
print(rhs)

Running the above program prints:

[-0.16078694+0.j -0.07726949+0.j  0.42268561+0.j]
[-0.16078694+0.j -0.07726949+0.j  0.42268561+0.j]

The Hermitian eigensolvers, eigh and eigvalsh, can be used to solve the generalized eigenvalue problem only if both the \(\mathbf{A}\) and \(\mathbf{B}\) matrices are Hermitian.