Sparse smallest eigenvalue problem on a linear subspace?

Question

I am interested in methods for solving the optimization problem

$$ \begin{array}{rl} \arg\min_x & x^T A x \\ \mathrm{s.t.} & x^T x = 1 \\ & Bx = 0 \end{array} $$

where $A$ is symmetric and full rank with well-separated eigenvalues, and $B$ has a nontrivial null space. In other words, I want to restrict $A$ to the null space of $B$ and find the eigenvector corresponding to the smallest eigenvalue. In the absence of the linear constraint, this problem is a standard eigenvalue problem and can be solved using power iterations. Ideally I would like to find a similar (iterative) procedure, since in my case I can gain a lot of performance by first computing a sparse factorization of $A$ (and potentially $B$ as well).

I am not interested in methods for dense problems, nor answers that simply amount to "just use software package X"; I am interested in implementing the method myself. For this reason, simple methods are preferred, even over methods that are perhaps slightly better/faster/more robust.

I am aware of the paper "Some Modified Matrix Eigenvalue Problems" by Gene Golub, but do not see a way to efficiently incorporate the generalized inverse in the sparse case. I am also aware of this question, but do not want to rely on a sophisticated solver for generalized eigenvalue problems; I'm really looking for something "not so different" from the power method.

Thanks!

Answer 1

Your main concern is not destroying sparsity - a transformation does not necessarily destroy it. Consider the transformed eigenvalue problem, that eliminates the constraint $x^Tx = 1$

$$ \min_x \frac{x^TAx}{x^Tx} \qquad s.t.\quad Bx = 0$$

Consider a matrix $Z$ that is a basis for the nullspace of $B$, i.e. $BZ = 0$. Then, restricting $x$ to be in the nullspace of $B$ as $x = Zy$, the eigenvalue problem now becomes an unconstrained optimization problem

$$ \min_{y} \frac{y^TZ^TAZy}{y^TZ^TZy}$$

The minimizer of this problem is the eigenvector corresponding to the smallest eigenvalue of the generalized eigenvalue problem

$$ Z^TAZ y = \lambda Z^T Z y$$

To solve this eigenvalue problem one does not have to form $Z^TAZ$ and destroy sparsity (which is your concern). One only needs to form matrix-vector products with $Z^TAZ$, this can be done in turn. How to form $Zx$? One possibility is to consider a partitioning using a permutation matrix $P$ such that $B_p$ is well-conditioned of rank $m$, where $B \in \mathbb{R}^{m\times n}$

$$BP = [B_p, B_n] \qquad Z = P \begin{bmatrix} -B_p^{-1}B_n \\ I \end{bmatrix}$$ In practice one doesn't invert $B_p$ but computes a sparse LU factorization. Once you know how to form $Zx$, any Krylov based eigenvalue solver for generalized eigenvalues can be used (their convergence is often better than power method). Other possibilities are discussed in section 6 (Null space methods) in this paper http://mathcs.emory.edu/~benzi/Web_papers/bgl05.pdf

Here is a good source for templates for implementing eigensolvers - you can pick your favorite http://web.eecs.utk.edu/~dongarra/etemplates/node155.html

Answer 2

Nobody has answered this question, so let me take a stab at an attempt. I don't quite know how to solve the problem exactly, but here's a method that should work.

Let us assume that you have a basis of the null space and that it forms the columns of a matrix $Z$, i.e., $BZ=0$. Then consider the eigenvalue problem $$ \begin{array}{rl} \arg\min_{x_\varepsilon} & {x_\varepsilon}^T A {x_\varepsilon} + \frac{1}{\varepsilon} {x_\varepsilon}^T Z^T Z {x_\varepsilon} \\ \mathrm{s.t.} & {x_\varepsilon}^T {x_\varepsilon} = 1 \end{array} $$ If you consider the limit $\varepsilon\rightarrow 0$, then any vector not in the null space of $B$ produces an infinite objective function, and you will only get those eigenvalues you care about.

Beyond this, let me point out one other idea (though I do not know enough about the numerical realization of eigenspectrum solvers to really say much about this). Namely, when you compute the second, third, etc, $L$th eigenvalue of a matrix $A$, then this is often written as $$ \begin{array}{rl} \arg\min_{x} & {x}^T A {x_\varepsilon} \\ \mathrm{s.t.} & {x_\varepsilon}^T {x_\varepsilon} = 1 \\ &x^T v_1 = 0 \\ & \vdots \\ &x^T v_{L-1} = 0 \end{array} $$ In other words, you try to minimize in a subspace. This is no different than what you want to do, except that your subspace isn't given by the previous eigenvectors $v_1, \ldots, v_{L-1}$ but by the basis of your null space $Z$.