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.



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

| cite | improve this answer | |
  • $\begingroup$ Any idea how to generalize this to constraints with non-zero right-hand sides? $Bx = b$? $\endgroup$ – Alec Jacobson Oct 16 '14 at 18:49
  • 2
    $\begingroup$ Adding the constraint $Bx = b$ no longer makes the solution an eigenvalue problem. However, you can still tackle it using optimization methods for quadratic programs with equality constraints. For example sequential quadratic programming methods would be one option. If you have trouble finding references, repost your question. $\endgroup$ – user2457602 Oct 22 '14 at 14:09
  • $\begingroup$ For low frequency modes, I need the action of multiplying by $(Z^TAZ)^{-1}$. But $Z$ is rectangular. Is there a way to approach this without wrecking sparsity? $\endgroup$ – Alec Jacobson May 20 at 13:38

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$.

| cite | improve this answer | |

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.