Generalized eigenvalue problem for symmetric, low rank matrix

I'd like to solve a generalized eigenvalue problem of the form:

$$\mathrm{A}x = \lambda \mathrm{B}x$$ $$s.t. x_i^T\mathrm{B}x_i=1.$$

Where $\mathrm{A}$ and $\mathrm{B}$ are symmetric but low-rank matrices and $x_i$ is the eigenvector associated with the largest eigenvalue $\lambda_i$.

I've been playing around with the geigen package in R, which provides a nice wrapper to LAPACK routines, and it seems that the constraint is satisfied only when either $\mathrm{A}$ or $\mathrm{B}$ is positive definite (I'm assuming it is implementing a Cholesky decomposition).

In other cases, geigen default's to LAPACK's dggev routine which results in eigenvectors that are orthogonal to $\mathrm{B}$. That is, $x_i\mathrm{B}x_i=0$, and so scaling/normalizing the eigenvectors have no effect.

I don't know enough to determine if this is expected behavior or if there another implementation to this generalized eigenvalue problem which satisfies this constraint?

thanks.

Answering my own question here thanks to this paper and some careful Googling, I get the desired result thanks to simultaneous diagonalization, with two applications of SVD. I've adapted the notation to be consistent with my original question.

$$\mathrm{X}^T\mathrm{A}\mathrm{X}=\Lambda$$ $$\mathrm{X}^T\mathrm{B}\mathrm{X}=\mathrm{I}$$

The goal is to find the transformation matrix $\mathrm{X}$ and diagonal matrix $\Lambda$, so that the columns contain the eigenvectors and diagonal entries contain eigenvalues of the generalized eigenvalue problem, respectively.

First perform a singular value decomposition of $\mathrm{B}$:

$$\mathrm{B}=\mathrm{V}_r\Sigma_r\mathrm{V}_r^T\equiv \mathrm{X}'(\mathrm{X}')^T$$

where $\mathrm{X}' \equiv \mathrm{V}_r\Sigma_r^{-\frac{1}{2}}$ and $r=\mathrm{rank}(\mathrm{B})$, so we're taking the first $r$ singular values and vectors to define matrix $\mathrm{X}'$.

Next, define a matrix $\mathrm{C}=(\mathrm{X}')^T\mathrm{A}\mathrm{X}'$ and perform an svd to obtain:

$$\mathrm{C}=(\mathrm{X}'')^T\Sigma' \mathrm{X}''$$

then $\mathrm{X} = \mathrm{X}'\mathrm{X}''$ and the above equations are satisfied if $\mathrm{I}$ and $\Lambda$ are of dimension $r\times r$.

• not sure if this implementation is the most efficient though, since it involves several matrix multiplications and 2 rounds of SVD. – zzk Jan 16 '15 at 16:54