Eigenvalue Problem

Content

• Preliminary

  • Definition

  • Perturbation theory

• Direct Eigensolver

  • Power method

  • Inverse power method

  • Orthogonal iteration

  • QR iteration

• Direct Symmetric Eigensolver

  • Divide-and-Conquer

  • Bisection and inverse iteration

  • Jacobi’s method

• Iterative Symmetric Eigensolver

  • Lanczos algorithm

Preliminary

Definition

The polynomial $p(\lambda) = \det(A-\lambda I)$ is called the characteristic polynomial of $A$. The roots of $p(\lambda)=0$ are the eigenvalues of $A$.

A nonzero vector $x$ satisfying $Ax=\lambda x$ is a (right) eigenvector for the eigenvalue $\lambda$.

A nonzero vector $y$ such that $y^\ast A = \lambda y^\ast$ is a left eigenvector.

A Jordan block has one right eigenvector $e_1=[1, 0, \ldots, 0]^T$, and one left eigenvector $e_n=[0, \ldots, 0, 1]^T$. Therefore, a matrix is diagonalizable if all the Jordan blocks are scalar, otherwise it is called defective. An $n$-by-$n$ defective matrix does not have $n$ eigenvectors.

For Jordan block $J_m(\lambda)$, if $m=1$, then $\lambda$ is called a simple eigenvalue.

Perturbation theory

Eigenvalues pf $A$ are continuous functions of $A$, even if they are not differentiable.

Condition number

Let $\lambda$ be a simple eigenvalue of $A$ with normalized right eigenvector $x$ and normalized left eigenvector $y$. Let $\lambda + \delta \lambda$ be the corresponding eigenvalue of $A+\delta A$. Then

where $\Theta(y, x)$ is the acute angle between $y$ and $x$. And $\sec \Theta(y, x)=\frac{1}{\vert y^*x \vert}$ is the condition number of the eigenvalue $\lambda$.

NOTE: The condition number of Jordan block is $1 / \vert e_n^\ast e_1 \vert = \infty$.

Pertubation theory for symmetric matrix

Eigenvalue:

Weyl’s Theorem: Let $A$ and $E$ be $n$-by-$n$ symmetric matrices. Let $\alpha_1 \geq \cdots \geq \alpha_n$ be the eigenvalues of $A$ and $\hat{\alpha}_1 \geq \cdots \geq \hat{\alpha}_n$ be the eigenvalues of $A + E$. Then $\vert \alpha_i - \hat{\alpha}_i \vert \leq \vert\vert E \vert\vert_2$.

This theorem could be proved using Courant-Fischer minimax theorem.

Eigenvector:

Let $A=Q\Lambda Q^T = Q\text{diag}(\alpha_i)Q^T$ be an eigendecomposition of $A$. Let $A+E=\hat{A}=\hat{Q} \hat{\Lambda} \hat{Q}^T$ be the pertured eigendecomposition. Let $\theta$ denote the acute between $q_i$ and $\hat{q}_i$.

Define eigen-gap as $\text{gap}(i,A)=\min_{j\neq i} \vert \alpha_j - \alpha_i \vert$. Then

Therefore, we can not hope to get the individual eigenvectors accurately if we have a group of $k$ tightly clustered eigenvalues. However, it is possible to compute the $k$-dimensional invariant subspace spanned by these vectors quite accurately.

Direct Eigensolver

We call a method direct if experience shows that it (nearly) never fails to converge in a fixed number of iterations.

Power method

Note that multiplying a scalar will not change the direction of a vector.

Thus the space of $x_i$ is the same as $A^i x_0$, further, assume $\lambda_1$ is the largest eigenvalue with eigen-gap strictly greater than $0$,

Thus, the direction of normalized vector will converge to $S e_1$, i.e., eigenvector corresponding to the largest eigenvalue.

Inverse power method

With a shift $\sigma$, the method will converge to the eigenvalue closest to $\sigma$,

Orthogonal Iteration

Note that if $col(X)=col(Y)$, then $col(AX)=col(AY)$.

Thus, we have

Furthermore,

Assume $\Lambda_1^k$ and $C_{01}$ are invertible, then we have

Note that right-multiplying an invertible matrix will not change the column space.

Finally, Let $\mu_i = \lambda_{p+i}$ (the eigenvalues in $\Lambda_2$), and $\nu_j = \lambda_j$ (the eigenvalues in $\Lambda_1$), then

Thus, $E_k \rightarrow 0$ when $k \rightarrow \infty$. This is the invariant subspace corresponding to the $p$ largest eigenvalues.

Q: what is the target of the algorithm ?

A: eigen-subspace comes from $Z_i$, however, the eigenvalue should be computed from another quantity $Z_i ^T A Z_i$. Note that if $A$ is not symmetric, this quantity will only converge to the Schur form of $A$ under some conditions.

Q: how to we obtain eigenvectors ?

A: The algorithm converges to the Schur form $Q^\ast A Q = T$. Then if $Tx=\lambda x$, we have $AQ=QTx=\lambda Qx$, $Qx$ is an eigenvector of $A$. Thus we only need to find eigenvectors of a block upper triangle matrix, i.e., solve some triangluar linear systems.

QR Iteration

To incorporate shifting and inverting, and simplify the computation of eigenvalues, reorganize the orthogonal iteration as

A simplified expression of iteration is

One could prove the under some condition $A_i$ is identical to the matrix $Z_i^T A Z_i$ implicitly computed by orthogonal iteration and thus converges to Schur form.

To make QR iteration practical, we need to solve the following issues:

• QR decomposition is too expensive. Can we seek some cheaper structure throughout iteration?

• How do we recognize convergence?

• How shall we choose shift $\sigma_i$ to accelerate convergence?

Hessenberg Reduction

Hessenber form: matrix $A$ is zero below the first subdiagonal. An upper Hessenberg matrix is unreduced if all subdiagonals are nonzero.

NOTE: Hessenberg form is the maximal zero structure preserved by QR iteration.

Fact 1: QR decomposition of a Hessenberg matrix is cheap ($O(n^2)$ instead of $O(n^3)$)

Fact 2: If you start with Hessenberg, QR iteration keeps it Hessenberg! This means the expensive step becomes cheap forever.

Fact 3 (Implicit Q Theorem): The orthogonal matrix $Q$ such that $Q^T A Q = H$ is unreduced upper Hessenberg is uniquely determined by its first column (up to signs).

If we initially reduce the matrix to uppper Hessenberg form, the QR iteration preserved the Hessenberg form.

Furthermore, since our target is only $A_{i+1}= Q_i^T A_i Q_i$, we don’t need to explicitly implement the QR decomposition and matrix multiplication, instead, the uniqueness of Hessenberg transformation tells us that we only need to design the first column and then find Givens rotations to make the resulted matrix upper Hessenberg.

The first column of $Q_i$ is proportional to be the first column of $A_i-\sigma_i I$ since this is the first column of Q in the QR decomposition of $A_i -\sigma_i I$.

Convergence condition

Since $A_i$ is related to $A$ by a similarity transformation by an orthogonal matrix, we expect $A_i$ to have roundoff errors of size $O(\varepsilon \vert\vert A \vert\vert)$.

Thus, nay subdiagonal entry of $A_i$ smaller than $O(\varepsilon \vert\vert A \vert\vert)$ in magnitude could be set to zero.

Implicit shift

AS QR iteration converges, $A_i$ becomes upper triangular, and diagonal entries approach eigenvalues. The last row will be zero except for the eigenvalue in $(n,n)$ position.

Types of shift

• bottom-right entry shift: $\sigma_i = A_i(n, n)$

• Wilkinson shift: choose the eigenvalue of the bottom-right 2-by-2 block that is closer to $A_i(n, n)$

• double shift: choose $\sigma_i$ as one eigenvalue of the bottom-right 2-by-2 block, perform QR step with $\sigma_i$ and its conjugate $\bar{\sigma}_i$, which could be done in an implicit manner

Deflation

If $\sigma_i$ is an exact eigenvalue of $A_i$, then the QR iteration converges in one step and reduce to smaller eigenproblem.

Once an eigenvalue converges at the bottom, we remove it and work on a smaller matrix.

NOTE: although the shift accelarate the convergence, there are tiny sets of matrices where the algorithm failed.

Direct Symmetric Eigensolver

For symmetric matrix, there are always some simplification for general eigensolver and some new algorithms utilizing this property to give better efficiency and accuracy.

Divide-and-Conquer

Given a tridiagonal symmetric matrix $T$, we try to reduce the problem into matrices in half of the size with a recursive manner.

Suppose we have obtained eigen-decomposition of $T_1$ and $T_2$ as $T_i = Q_i \Lambda_i Q_i^T$ for $i=1,2$. Then we have

where

Hence, we only need to consider eigendecomposition of $D+\rho uu^T$.

Eigenvalues

For eigenvalues, using $\det(I+xy^T)=1+y^Tx$, we have the secular equation

Here we assume that $D-\lambda I$ is nonsingular. This is because for $(D + \rho uu^T)q=\lambda q$ with $Dq = \lambda q$,

• if $u^T q=0$, then $u_j=0$,
• if $u^T q \neq 0$, then we have $\lambda=d_j$ and further $u_j=0$,

both cases gives deflation for the secular equation, i.e., the number of roots need to solved become less. If $u_i=0$, then the corresponding eigenvector is $e_i$.

For solving the secular equation, we use the Newton’s method. Since $f(\lambda)$ is not well approximated by a straight line $l(x)$, we approximate it by another simple function $h(x)$.

Given the approximated zero $\lambda_j$ (note that this is not the notation for exact eigenvalue), we construct

Let

and

then we ask

Eigenvectors

If $\alpha$ is an eigenvalue of $D+\rho uu^T$, then $(D - \alpha I)^{-1} u$ is its eigenvectors.

However, this simple formula is not numerically stable. If $\alpha_i$ and $\alpha_{i+1}$ are close, then $(D - \alpha_i I)^{-1} u$ and $(D - \alpha_{i+1} I)^{-1} u$ will lose orthogonality.

To maintain the orthogonality between vectors is a difficult issue (recall the Gram-Schmidt process), some alternative methods have been developed for stable computation of eigenvectors.

Bisection and Inverse Iteration

The Sylvester’s inertia theorem asserts that for nonsingular $X$,

By Gaussian elimination,

Thus, we can count the number of eigenvalues in different signs from the sign of diagonal elements in $D$.

Define

which is simple and turns out to be stale to compute for tridiagonal $A$.

Suppose $z_1 < z_2$ , then the number of eigenvalues in the interval $[z_1, z_2)$ could be computed by

Combined with bisection algorithm, we can find all eigenvalues of $A$ insides $[a,b)$ to a given tolerance.

Once we have computed the eigenvalues, we use inverse iteration with shifts to compute eigenvectors. For orthogonality, one may perform a QR decomposition to reorthogonalize the computed eigenvectors.

Jacobi’s method

Jabobi’s method works on the original dense matrix and is usually much slower than the previous method. But it can sometimes compute tiny eigenvalues and their eigenvectors with much higher accuracy than the previous methods and can be easily parallelized.

We define the Jacobi rotation $J_i=R(j,k,\theta)$ as Givens rotation embedding into identity matrix aligning with the dimension of the matrix $A$.

We will make $J^T A J$ nearly diagonal by iteratively choosing $J_i$ to make one pair of offdiagonal entries of $A_{i+1}=J_i^T A_i J_i$ zero at a time, i.e.,

For each rotation step fixing index pair $(j,k)$, we have

Note that we have zero out $a_{jk}$, thus

By the design of $J_i=R(j,k,\theta)$, we know that

Since orthogonal rotation will not change the Frobenius norm and the sum-of-square of the four cross-term are the same, we have

The Jacobi method will eventually converge to the diagonal form and gives the desired eigenpairs.

Iterative Symmetric Eigensolver

Lanczos Algorithm

Suppose $A=A^T$ is a symmetric matrix.

Starting with a given $x_1$, the $k-1$ iteration of either power method or inverse iteration produce a sequence of vectors $x_1, \ldots, x_k$. These vectors span a Krylov subspace, defined as

Rather than taking $x_k$ as our approximate eigenvector, it is natural to ask for the “best” approximate eigenvector in $\mathcal{K}_k$. These best approximations are then called Ritz values and Ritz vectors.

Rayleigh-Ritz method

Given a partial unitary matrix $Q_k \in \mathbb{R}^{n\times k}$ such that $Q_k^T Q_k = I_k$, we approximate the eigenvalues of matrix $A$ by the eigenvalues of

The approximations are called Ritz values.

Let $T_k = V \Lambda V^T$ be the eigendecomposition of $T_k$, the eigenvector approximations are given by $Q_k V$ and are called Ritz vectors.

The Ritz values and Ritz vectors are “optimal” approximation in the sense that it minimizes the residual

or equivalent, let $P_k$ ranges all partial orthogonal matrix such that $\text{span}(P_k) = \text{span}(Q_k)$, it minimizes

Lanczos Eigensolver

When applying the Rayleigh-Ritz method with $Q_k$ computed by the Lanczos algorithm, the resulted

is a tridiagonal matrix and hence easy to compute the approximated Ritz values, Ritz vectors from $T_k = Q_k^T A Q_k$ and residuals.

One of the possible issues is that the lose of orthogonality between eigenvectors. A careful implementation with more cost is called Lanczos with full reorthogonalization.