Advanced Linear Algebra: Foundations to Frontiers

Subsection 10.1.1 Subspace iteration with a Hermitian matrix

The idea behind subspace iteration is to perform the Power Method with more than one vector in order to converge to (a subspace spanned by) the eigenvectors associated with a set of eigenvalues.

We continue our discussion by restricting ourselves to the case where \(A \in \Cmxm\) is Hermitian. Why? Because the eigenvectors associated with distinct eigenvalues of a Hermitian matrix are mutually orthogonal (and can be chosen to be orthonormal), which will simplify our discussion. Here we repeat the Power Method:

In previous discussion, we used \(v^{(k)} \) for the current approximation to the eigenvector. We now add the subscript to it, \(v_0^{(k)} \text{,}\) because we will shortly start iterating with multiple vectors.

[ lambda_0, v0 ] = PowerMethod( A, x, maxiters, illustrate, delay )

Recall that when we analyzed the convergence of the Power Method, we commented on the fact that the method converges to an eigenvector associated with the largest eigenvalue (in magnitude) if two conditions are met:

• \(\vert \lambda_0 \vert \gt \vert \lambda_1 \vert \text{.}\)

• \(v_0^{(0)} \) has a component in the direction of the eigenvector, \(x_0 \text{,}\) associated with \(\lambda_0 \text{.}\)

A second initial vector, \(v_1^{(0)} \text{,}\) does not have a component in the direction of \(x_0 \) if it is orthogonal to \(x_0 \text{.}\) So, if we know \(x_0 \text{,}\) then we can pick a random vector, subtract out the component in the direction of \(x_0 \text{,}\) and make this our vector \(v_1^{(0)} \) with which we should be able to execute the Power Method to find an eigenvector, \(x_1 \text{,}\) associated with the eigenvalue that has the second largest magnitude, \(\lambda_1 \) . If we then start the Power Method with this new vector (and don't introduce roundoff error in a way that introduces a component in the direction of \(x_0 \)), then the iteration will home in on a vector associated with \(\lambda_1 \) (provided \(A \) is Hermitian, \(\vert \lambda_1 \vert \gt \vert \lambda_2 \vert \text{,}\) and \(v_1^{(0)} \) has a component in the direction of \(x_1 \text{.}\)) This iteration would look like

[ lambda_1, v1 ] = PowerMethodLambda1( A, x, x0, maxiters, illustrate, delay )

Because we should be concerned about the introduction of a component in the direction of \(x_0 \) due to roundoff error, we may want to reorthogonalize with respect to \(x_0 \) in each iteration:

[ lambda_1, v1 ] = PowerMethodLambda1Reorth( A, x, v0, maxiters, illustrate, delay  );

We now observe that the steps that normalize \(x_0 \) to have unit length and then subtract out the component of \(v_1 \) in the direction of \(x_0 \text{,}\) normalizing the result, are exactly those performed by the Gram-Schmidt process. More generally, it is is equivalent to computing the QR factorization of the matrix \(\left( \begin{array}{c|c} x_0 \amp v_1 \end{array} \right).\) This suggests the algorithm

Obviously, this redundantly normalizes \(x_0 \text{.}\) It puts us on the path of a practical algorithm for computing the eigenvectors associated with \(\lambda_0 \) and \(\lambda_1 \text{.}\)

The problem is that we typically don't know \(x_0 \) up front. Rather than first using the power method to compute it, we can instead iterate with two random vectors, where the first converges to a vector associated with \(\lambda_0 \) and the second to one associated with \(\lambda_1 \text{:}\)

We observe:

• If \(\vert \lambda_0 \vert \gt \vert \lambda_1 \vert \text{,}\) the vectors \(v_0^{(k)} \) will converge linearly to a vector in the direction of \(x_0 \) at a rate dictated by the ratio \(\vert \lambda_1 \vert / \vert \lambda_0 \vert \text{.}\)

• If \(\vert \lambda_0 \vert \gt \vert \lambda_1 \vert \gt \vert \lambda_2 \vert, \text{,}\) the vectors \(v_1^{(k)} \) will converge linearly to a vector in the direction of \(x_1 \) at a rate dictated by the ratio \(\vert \lambda_2 \vert / \vert \lambda_1 \vert \text{.}\)

• If \(\vert \lambda_0 \vert \geq \vert \lambda_1 \vert \gt \vert \lambda_2 \vert\) then \(\Span( \{ v_0^{(k)}, v_1^{(k)} \} ) \) will eventually start approximating the subspace \(\Span( \{ x_0, x_1 \} ) \text{.}\)

What we have described is a special case of subspace iteration. The associated eigenvalue can be approximated via the Rayleigh quotient:

Alternatively,

if \(A \) is Hermitian, \(\vert \lambda_1 \vert \gt \vert \lambda_2 \vert \text{,}\) and \(v^{(0)} \) and \(v^{(1)} \) have components in the directions of \(x_0 \) and \(x_1 \text{,}\) respectively.

The natural extention of these observations is to iterate with \(n \) vectors:

By extending the reasoning given so far in this unit, if

• \(A \) is Hermitian,

• \(\vert \lambda_0 \vert \gt \vert \lambda_1 \vert \gt \cdots \gt \vert \lambda_{n-1} \vert \gt \vert \lambda_{n} \vert \text{,}\) and

• each \(v_j \) has a component in the direction of \(x_j \text{,}\) an eigenvector associated with \(\lambda_j \text{,}\)

then each \(v_j^{(j)} \) will converge to a vector in the direction \(x_j \text{.}\) The rate with which the component in the direction of \(x_p \text{,}\) \(0 \leq p \lt n \text{,}\) is removed from \(v_j^{(k)} \text{,}\) \(n \leq j \lt m \text{,}\) is dictated by the ratio \(\vert \lambda_p \vert / \vert \lambda_j \vert \text{.}\)

If some of the eigenvalues have equal magnitude, then the corresponding columns of \(\widehat V^{(k)} \) will eventually form a basis for the subspace spanned by the eigenvectors associated with those eigenvalues.

[ Lambda, V ] = SubspaceIteration( A, V, maxiters, illustrate, delay  );