6
$\begingroup$

Assume that we have generalized eigenvalue problem:

$B^HB\textbf{x} = \lambda A\textbf{x}$

where $A$ is an nxn Hermitian sparse matrix (n is very large, so we do not have $A^{-1}$ but can solve using iterative methods) and full-rank, and $B$ is a 2xn matrix such that $B^HB$ is also nxn but only rank 2. Thus, we know that this problem can only have 2 non-zero eigenvalues. Is there any simple way for finding the two eigenpairs corresponding to nonzero eigenvalues by taking advantage of the very low rank of $B^HB$? Assume that we have the two eigenvectors of $B$.

If I am only interested in the eigenvector corresponding to the largest eigenvalue, is there a faster way of finding it than using simple power iteration on the transformed standard eigenvalue problem: $A^{-1}B^HB\textbf{x} = \lambda\textbf{x}$?

Thanks!

$\endgroup$
4
$\begingroup$

This answer is essentially a fix of the approach suggested by @WolfgangBangerth, as there is not enough space in the comments.

Starting from $$ B^H B x = \lambda A x, $$ if we are interested in eigenpairs corresponding to nonzero eigenvalues, then we must have that $B^H B x$ lies in the range of $A$, and $Ax$ lies in the range of $B^H B$, which is to say that, since $A$ is invertible, $$ B^H B x \in \mathrm{Range}(A) = \mathbf{C}^n, $$ and $$ Ax \in \mathrm{Range}(B^H B) = \mathrm{span}(B^H). $$ Now, the first constraint is trivially satisfied, but we must ensure that $Ax \in \mathrm{span}(B^H)$, which is equivalent to the constraint $$ x \in \mathrm{span}(A^{-1} B^H). $$ Then if the columns of a unitary matrix $Q$ span the columns of $A^{-1}B^{H}$, we have that $$ x = Q Q^H x $$ for any eigenvector corresponding to a nonzero eigenvalue.

We are now ready to use the mechanism from Wolfgang's approach:

  1. Compute $W := A^{-1} B^H$ through two (preconditioned) Krylov solves
  2. Compute $[Q,R]=\mathrm{qr}(W)$
  3. Form $K := (B Q)^H (B Q)$ and $M := Q^H (A Q)$
  4. Solve the $2 \times 2$ eigenvalue problem $K U = M U \Lambda$
  5. Form the interesting global eigenvectors, $Z := Q U$.
$\endgroup$
  • $\begingroup$ Thanks, Jack! Are you sure you need to do QR? I don't think the two vectors comprising the W matrix need to necessarily be orthonormal? (A quick test shows that it works even if you solve $W^HB^HBW\textbf{y} = \lambda W^HAW\textbf{y}$ ) $\endgroup$ – Costis May 22 '12 at 2:08
  • $\begingroup$ The QR decomposition for an $m \times n$ matrix, $m \ge n$, is $O(mn^2)$. In this case, $n=2$, so the cost is linear and should be dominated by the Krylov solves. I am skeptical of how this would work without a QR decomposition. $\endgroup$ – Jack Poulson May 22 '12 at 2:26
  • $\begingroup$ $W=QR$, so substituting: $R^HQ^HB^BQR\textbf{y}=\lambda R^HQ^HAQR\textbf{y}$. Multiply both sides by $R^{-H}$ and substitute $\textbf{x}=R\textbf{y}$. You end up with $Q^HB^BBQ\textbf{x}=\lambda Q^HAQ\textbf{x}$ which has the same eigenvalues as if you just used W instead of Q. $\endgroup$ – Costis May 22 '12 at 2:34
  • $\begingroup$ I can get the eigenvector by just doing: $\textbf{w}=W\textbf{y}$ which implicitly multiplies by R since $W=QR$. Just tried a quick test case and it seems to work, although as you said I think it would be trivial do QR as compared to doing the Krylov solves. $\endgroup$ – Costis May 22 '12 at 2:45
  • $\begingroup$ Ah, good point! I would still rather work with $Q$ though, as the cost of computing it is insignificant, and it will be more numerically stable. $\endgroup$ – Jack Poulson May 22 '12 at 2:48
3
$\begingroup$

If $B$ is $2\times n$, then the only two non-trivial eigenvectors (i.e. the eigenvectors corresponding to the two non-zero eigenvalues) can be written as linear combinations of the vectors that form the two rows of $B$. Let's call these two vectors $b_1, b_2$ so that $B=\left[\begin{matrix}b_1^T\\b_2^T\end{matrix}\right]$.

Now, let $P \in {\mathbb R}^{2\times n}$ be the projector from ${\mathbb R}^n$ onto the two-dimensional space spanned by $b_1,b_2$. Since we are only interested in vectors in this space, we know that the two non-trivial eigenvectors must satisfy $x = P^TPx$. The eigenvalue problem can then be written as $$ B^H B P^T P x = \lambda A P^T P x. $$ Even though this linear system has $n$ rows, it is really only a two-dimensional problem since we can only determine only two components of $x$. The remainder of the linear system is over-determined, but we can select the two independent equations by projecting onto the non-trivial subspace: $$ P B^H B P^T P x = \lambda P A P^T P x. $$

In other words, you only have to solve the $2 \times 2$ eigenvalue problem $$ P B^H B P^T y = \lambda (P A P^T) y. $$ This is easy to solve since the matrices involved are only $2\times 2$ and the matrix on the right can easily be computed using just two matrix-vector and two vector-vector products.

$\endgroup$
  • $\begingroup$ I don't think your assumption that $x=P^T Px$ is valid for non-trivial $A$. $\endgroup$ – Jack Poulson May 21 '12 at 20:14
  • $\begingroup$ I think that $P$ will need to be modified to span a space including the columns of $B^H$ and $A^{-1} B^H$, which is at most rank 4, and only requires two solves with $A$ to set up. $\endgroup$ – Jack Poulson May 21 '12 at 21:37
  • $\begingroup$ Hmmm.. I think you only need the columns of $A^{-1}B^H$ actually, so if you take $P=A^{-1}B^H$ Wolfgang's approach will work. $\endgroup$ – Costis May 22 '12 at 1:14
  • $\begingroup$ @Costis: I think you are right. I have a nice explanation of why which I will post as an answer, as there is not enough space here. $\endgroup$ – Jack Poulson May 22 '12 at 1:31
  • $\begingroup$ I feel like the nomenclature is slightly unclear. $P$ cannot be a projector, because $P$ is neither square, nor idempotent. For the sake of clarity, the relevant (orthogonal) projector appears to be $P^{T}P$. Your explanation also seems to implicitly rely on knowing that the eigenvectors of $B^{H}B$ form an orthonormal basis, which is why an orthogonal projector is appropriate. $\endgroup$ – Geoff Oxberry May 22 '12 at 2:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.