Text provided under a Creative Commons Attribution license, CC-BY. All code is made available under the FSF-approved MIT license. (c) Kyle T. Mandli

Note: This material largely follows the text “Numerical Linear Algebra” by Trefethen and Bau (SIAM, 1997) and is meant as a guide and supplement to the material presented there.

Eigenproblems¶

Overview¶

We will now consider eigenproblems of the form

where $A \in \mathbb C^{m \times m}$, $\mathbf{x} \in \mathbb C^m$ and $\lambda \in \mathbb C$. The vector $\mathbf{x}$ is known as the eigenvector and $\lambda$ the eigenvalue. The set of all eigenvalues is called the spectrum of $A$.

The basics¶

The eigenproblem

can be rewritten as

which implies that the eigenvectors are in the Null space of $A-\lambda I$.

However for this matrix to have a non-trivial Null space, requires that $A-\lambda I$ is singular.

Characteristic Polynomial¶

If $A-\lambda I$ is singular, it follows that

where ${\cal P}_A(\lambda)$ can be shown to be a $m$th order polynomial in $\lambda$ known as the characteristic polynomial of a matrix $A$

We can then state the following theorem regarding the zeros of $\mathcal{P}_A$ and the eigenvalues of $A$:

Theorem: $\lambda$ is an eigenvalue of $A$ if and only if $\mathcal{P}_A(\lambda) = 0$.

i.e. the eigenvalues are the roots of ${\cal P}_A(\lambda)$, and therefore there are exactly $m$ eigenvalues.

Proof:

Note that this theorem implies that even though $A \in \mathbb R^{m \times m}$ that $\lambda \in \mathbb C$.

Computing Eigenvalues¶

In basic linear algebra classes we usually find the eigenvalues by directly calculating the roots of ${\cal P}_A(\lambda)$ which can work for low-degree polynomials. Unfortunately the following theorem (due to Galois) suggests this is not a good way to compute eigenvalues:

Theorem: For an $m \geq 5$ there is a polynomial $\mathcal{P}(z)$ of degree $m$ with rational coefficients that has a real root $\mathcal{P}(z_0) = 0$ with the property that $z_0$ cannot be written using any expression involving rational numbers, addition, subtraction, multiplication, division, and $k$th roots.

I.e., there is no way to find the roots of a polynomial of degree $>4$ in a deterministic, fixed number of steps.

Not all is lost however!¶

We just must use an iterative approach where we construct a sequence that converges to the eigenvalues.

Some Questions

• How does this relate to how we found roots previously?

• Why will it still be difficult to use our rootfinding routines to find Eigenvalues?

We will return to how we actually find Eigenvalues (and roots of polynomials) after a bit more review

Eigenvalue Factorization and Diagonalization¶

Given that there are exactly $m$ (possibly repeated) eigenvalues for a $m\times m$ system, the eigenproblem is really more correctly written as

Or in Matrix form as

where $X$ is the matrix formed by the eigenvectors $x$ as its columns and $\Lambda$ is a diagonal matrix with the eigenvalues along its diagonal.

Expanded, $A X = X \Lambda$ looks like:

Here we note that the eigenpair $(\mathbf{x}_j, \lambda_j)$ are matched as the $j$th column of $X$ and the $j$th element of $\Lambda$ on the diagonal.

Diagonalization¶

Eigenproblems can always be written as

However, if there are a linearly independent set of eigenvectors then the matrix $X$ is invertible (why?). Under this condition we can transform $A$ into the diagonal matrix $\Lambda$ by

Factorization¶

or more usefully, we can rewrite $A$ as a product of more useful matrices

And use it in similar ways to how we use other factorization such as $A=QR$ and $A=LU$

The Rules for diagonalizability¶

Unfortunately, not all matrices can be diagonalized (i.e. don’t have a full linearly independent set of eigenvectors that span $\mathbb{R}^m$ or $\mathbb{C}^m$). Here are the rules.

A matrix can be factored as $A = X\Lambda X^{-1}$ if

1. all eigenvalues are distinct (never repeat). These are known as simple eigenvalues

2. Eigenvalues repeat, but for every repeated eigenvalue there can be found a linearly independent eigenvector

3. The matrix is Hermitian ($A^\ast = A$, or if real, then symmetric $A^T = A$). In this case the eigenvalues are always real and the eigenvectors can always be chosen orthonormal. In this special case $X=Q$ and

The issue of repeated eigenvalues can be made a bit more precise by defining

• Algebraic multiplicity: the number of times an eigenvalue is repeated

• Geometric multiplicity: the number of linearly independent eigenvectors corresponding to each eigenvalue.

If the algebraic multiplicity is equal to the geometric multiplicity for all $\lambda$ then we can say that there is a full eigenspace and the matrix is diagonalizable.

Example: Computing Multiplicities¶

Compute the geometric and algebraic multiplicities for the following matrices. What is the relationship between the algebraic and geometric multiplicities?

1. The characteristic polynomial of $A$ is

so the eigenvalues are all $\lambda = 2$ so we know the algebraic multiplicity is 3 of this eigenvalue. The geometric multiplicity is determined by the number of linearly independent eigenvectors. For this matrix we have three eigenvectors that are all linearly independent which happen to be the unit vectors in each direction (check!). This means that the geometric multiplicity is also 3.

1. The characteristic polynomial of $B$ is the same as $A$ so again we know $\lambda = 2$ but now we need to be a bit careful about the eigenvectors. In this case the only eigenvector is a scalar multiple of $e_1$ so the geometric multiplicity is 1.

Interpretations of the Eigenspace¶

One way to interpret the eigenproblem is that of one that tries to find the subspaces of $\mathbb C^m$ which act like scalar multiplication by $\lambda$. The eigenvectors associated with one eigenvalue then form a subspace of $S \subseteq \mathbb C^m$.

When an eigenvalue has algebraic multiplicity that equals its geometric then it is called non-defective and otherwise defective. This property is also inherited to the matrix so in the above example $A$ and $B$ are non-defective and defective matrices respectively.

Determinant and Trace¶

Two important properties of matrices have important relationships with their eigenvalues, namely the determinant and trace. The determinant we have seen, the trace is defined as the sum of the elements on the diagonal of a matrix, in other words

The relationship between the determinant and the eigenvalues is not difficult to guess due to the nature of the characteristic polynomial. The trace of a diagonal matrix is clear and provides another suggestion to the relationship.

Theorem: The determinant $\det(A)$ and trace $\text{trace}(A)$ are equal to the product and sum of the eigenvalues of $A$ respectively counting algebraic multiplicity.

Similarity Transformations¶

Generally, we say any two matrices $A$ and $B$ are similar if they can be related through an invertible matrix $M$ as

Example

a diagonalizable matrix $A$ is similar to the diagonal matrix $\Lambda$ through the invertible matrix $X$

Theorem: If $A$ and $B$ are similar matrices, they have the same eigenvalues and their eigenvectors are related through an invertible matrix $M$

Proof: Let

or

if $A\mathbf{x} = \lambda\mathbf{x}$ then

or

which shows that $\lambda$ is also an eigenvalue of $B$ with corresponding eigenvector $\mathbf{y} = M\mathbf{x}$

Schur Factorization¶

A Schur factorization of a matrix $A$ is defined as

where $Q$ is unitary and $T$ is upper-triangular. Because $Q^\ast=Q^{-1}$ (for square unitary matrices). It follows directly that $A$ and $T$ are similar.

• Good News! $T$ is upper triangular so its eigenvalues can just be read of the diagonal

• Bad News! There is no deterministic way to calculate $T$ as that would violate Galois theory of polynomials

Theorem: Every matrix $A \in \mathbb C^{m \times m}$ has a Schur factorization.

Note that the above results imply the following

• An eigen-decomposition $A = X \Lambda X^{-1}$ exists if and only if $A$ is non-defective (it has a complete set of eigenvectors)

• A unitary transformation $A = Q \Lambda Q^\ast$ exists if and only if $A$ is normal ($A^\ast A = A A^\ast$)

• A Schur factorization always exists

Note that each of these lead to a means for isolating the eigenvalues of a matrix and will be useful when considering algorithms for finding them.

Condition Number of a Simple Eigenvalue¶

Before we discuss a number of approaches to computing eigenvalues it good to consider what the condition number of a given eigenproblem is.

Let

define the eigenvalue problem in question. Here we will introduce a related problem

where $y$ is the left eigenvector and from before $x$ is the right eigenvector. These vectors also can be shown to have the relationship $y^\ast x \neq 0$ for a simple eigenvalue.

Now consider the perturbed problem

Expanding this and throwing out quadratic terms and removing the eigenproblem we have

Multiple both sides of the above by the left eigenvector and use $y^\ast x \neq 0$ to find

where we again use the slightly different definition of the eigenproblem. We can then solve for $\delta \lambda$ to find

meaning that the ratio between the dot-product of the left and right eigenvectors and the conjugate dot-product of the matrix $\delta A$ then form a form of bound on the expected error in the simple eigenvalue.

Computing Eigenvalues¶

Almost all useful approaches to computing eigenvalues do so through the computation of the Schur factorization. The Schur factorization, as we have seen, will preserve the eigenvalues. The steps to compute the Schur factorization are usually broken down into two steps

1. Directly transform $A$ into a Hessenberg matrix, a matrix that contains zeros below its first sub-diagonal, directly using Householder reflections. This is as close to triangular that you can get by direct similarity transformations.

2. Use an iterative method to change the sub-diagonal into all zeros

Hessenberg and Tridiagonal form¶

What we want to do is construct a sequence of similarity transformations matrices that turns $A$ into a Hessenberg matrix with the same eigenvalues as $A$. We use Householder reflections to do this with the important distinction that we only want to remove zeros below the first sub-diagonal.

so we have the sequence $H = Q^\ast A Q$ which has the same eigenvalues as the original matrix $A$.

Question? Why can’t we just use Householder to take $A\rightarrow T$ like we did for the $QR$?

One important special case of this sequence of transformations is that if the matrix $A$ is hermitian (the matrix is its own conjugate transpose, $A = A^\ast$, or symmetric in the real case) then the Hessenberg matrix is tridiagonal.

We now will focus on how to formulate the iteration step of the eigenproblem. We will also restrict our attention to symmetric, real matrices. This implies that all eigenvalues will be real and have a complete set of orthogonal eigenvectors. Generalizations can be made of many of the following algorithms but is beyond the scope of this class.

Rayleigh Quotient and Inverse Iteration¶

There are a number of classical approaches to computing the iterative step above which we will review here. Inverse power iteration in particular is today still the dominant means of finding the eigenvectors once the eigenvalues are known.

Rayleigh Quotient¶

The Rayleigh quotient of a vector $x \in \mathbb R^m$ is the scalar

The importance of the Rayleigh quotient is made clear when we evaluate $r(x)$ at an eigenvector. When this is the case the quotient evaluates to the corresponding eigenvalue.

The Rayleigh quotient can be motivated by asking the question, given an eigenvector $x$, what value $\alpha$ acts most like an eigenvalue in an $\ell_2$ sense:

This can be reformulated as a least-squares problem noting that $x$ is the “matrix”, $\alpha$ is the unknown vector (scalar) and $Ax$ is the right-hand side so we have

which can be solved so that

Power Iteration¶

Power iteration is a straightforward approach to finding the eigenvector of the largest eigenvalue of $A$. The basic idea is that the sequence

will converge (although very slowly) to the desired eigenvector.

We implement this method by initializing the algorithm with some vector $v$ with $||v|| = 1$. We then apply the sequence of multiplications.

The reason why this works can be seen by considering the initial vector $v$ as a linear combination of the orthonormal eigenvectors (which we have assumed exist) such that

Multiplying $v^{(0)}$ by $A$ then leads to

here $c_1$ is some constant due to the fact the eigenvectors are not uniquely specified.

wRepeating this $k$ times we have

Since $\lambda_1 > \lambda_i$ for all $i \neq 1$ then in the limit the terms $\lambda_2^k / \lambda_1^k$ will approach zero and on normalization $v^{(k)}/||v^{(k)}||\rightarrow \mathbf{q}_1$.

Inverse Iteration with shifts¶

Inverse iteration with shifts uses a similar approach with the difference being that we can use it to find any of the eigenvectors for the matrix $A$.

Some Preliminaries: inverse and shift rules of Eigenvalues

Show that if $\mathbf{x}$ is an eigenvector of $A$ with eigenvalue $\lambda$, then

• $\mathbf{x}$ is an eigenvector of $A^{-1}$ with eigenvalue $1/\lambda$

• $\mathbf{x}$ is an eigenvector of $A -\sigma I$ with eigenvalue $\lambda - \sigma$

So...

If we want to find the smallest eigenvalue we can consider the power method on $A^{-1}$,

But we really don’t want to find $A^{-1}$ which is expensive, instead we can do the equivalent iteration

x[0] = x0
for i in range(MAX_ITS):
    solve A w[i] = x[i]
        x[i+1] = w[i]/norm(w[i]) 

as $\mathbf{w}_i = A^{-1}\mathbf{x}_i$

and if we want to find the eigenvalue closest to some number $\mu$ we can apply the power method to

the eigenvectors of this matrix are

where $\lambda_j$ are the eigenvalues of $A$.

If $\mu$ is close to a particular $\lambda_j$, say $\lambda_J$, then

will be larger than any of the other $(\lambda_j - \mu)^{-1}$. In this way we effectively have picked out the eigenvalue we want to consider in the power iteration!

Rayleigh Quotient Iteration¶

By themselves the above approaches are not particularly useful but combining them we can iterate back and forth to find the eigenvalue, eigenvector pair:

1. Compute the Rayleigh quotient and find an estimate for $\lambda_j$

2. Compute one step of inverse iteration to approximate $x_j$

3. Repeat...

QR/RQ Algorithm¶

All of the above methods pick out at most a few eigenvalues at a time. However it turns out we can use the $QR$ algorithm, to iterate towards the Schur factorization and find all the eigenvalues simultaneously.

The simplest algorithm just iterates

    while not converged:
        Q, R = numpy.linalg.qr(A)
        A = R.dot(Q)        

calculating the $QR$ factorization of $A$, then forming a new $A=RQ$, This sequence will eventually converge to the Schur decomposition of the matrix $A$.

Code this up and see what happens.

So why does this work? The first step is to find the $QR$ factorization of $A^{(k-1)}$ or

which is equivalent to finding

and multiplying on the right leads to

In this way we can see that this is a similarity transformation of the matrix $A^{(k-1)}$ since the $Q^{(k)}$ is an orthogonal matrix ($Q^{-1} = Q^T$). This of course is not a great idea to do directly but works great in this case as we iterate to find the upper triangular matrix $R^{(k)}$ which is exactly where the eigenvalues appear.

In practice this basic algorithm is modified to include a few additions:

1. Before starting the iteration $A$ is reduced to Hessenberg form.

2. Motivated by the inverse power iteration we observed we instead consider a shifted matrix $A^{(k)} - \mu^{(k)} I$ for factoring. The $\mu$ picked is related to the estimate given by the Rayleigh quotient. Here we have

   $$\mu^{(k)} = \frac{(q_m^{(k)})^T A q_m^{(k)}}{(q_m^{(k)})^T q_m^{(k)}} = (q_m^{(k)})^T A q_m^{(k)}.$$

   (45)

3. Deflation is used to reduce the matrix $A^{(k)}$ into smaller matrices once (or when we are close to) finding an eigenvalue to simplify the problem.

This has been the standard approach until recently for finding eigenvalues of a matrix.

Application: Finding the roots of a polynomial¶

Numpy has a nice function called roots which returns the $n$ roots of a $n$th degree polynomial

described by a $n+1$ vector of coefficients $\mathbf{c}$

This routine, does not try and actually find the roots of a high-order polynomial, instead it actually calculates the eigenvalues of a companion matrix $C$ whose characteristic polynomial $P_C(\lambda)$ is the monic polynomial $p_m(x) = p(x)/c_n$.

It can be shown that this matrix can be constructed as (see e.g.)

Alternatives¶

Jacobi¶

Jacobi iteration employs the idea that we know the eigenvalues of a matrix of size equal to or less than 4 (we know the roots of the characteristic polynomial directly). Jacobi iteration therefore attempts to break the matrix down into at most 4 by 4 matrices along the diagonal via a series of similarity transformations based on only diagonalizing sub-matrices 4 by 4 or smaller.

Bisection¶

It turns out if you do not want all of the eigenvalues of a matrix that using a bisection method to find some subset of the eigenvalues is often the most efficient way to get these. This avoids the pitfall of trying to find the eigenvalues via other root-finding approaches by only needing evaluations of the function and if a suitable initial guess is provided can find the eigenvalue quickly that is closest to the initial bracket provided.

Divide-and-conquer¶

This algorithm is actually the one used most often used if both eigenvalues and eigenvectors are needed and performs up to twice as fast as the $QR$ approach. The basic idea is to split the matrix into two pieces at every iteration by introducing zeros on the appropriate off-diagonals which neatly divides the problem into two pieces.

Arnoldi and Lanczos Iteration¶

Krylov subspace methods (which we will unfortunately not cover) are another approach to finding eigenvalues of a matrix. These methods generally use some piece of the $QR$ approach outlined above and are extremely effective at finding the “extreme” eigenvalues of the matrix.