Efficiently computing a few localized eigenvectors

Question

Let $H = \triangle + V(x) : \mathbb{R}^2 \rightarrow \mathbb{R}^2$. I am interested in domain decomposition for an eigenproblem involving $H$.

The lowest 1000 eigenfunctions of $H$, $ \psi_i $, can be partitioned using a region, $\Omega \subset \mathbb{R}^2$, such that each $\psi_i$ localizes either inside of $\Omega$ or outside of $\Omega$. $\Omega$ is not a subspace of $\mathbb{R}^2$ as it may be an oddly shaped region.

Label the inner eigenfunctions $\psi_i^{in}$ and the outer ones $\psi_i^{out}$. There's only about 10 $\psi_i^{in}$s. Given $\Omega$, my goal is to efficiently compute the $\psi_i^{in}$.

One way to find the $\psi_i^{in}$ would be to discretize, compute all 1000 $\psi_i$s, and then partition. This is what I do now (5-point stencil for $\triangle$ on a $10^3 \times 10^3$ grid). The problem is that this requires diagonalizing over a 1000 dimensional space in order to get 10 eigenvectors. It seems like there would be a cheaper way to compute the $\psi_i^{in}$.

Does anyone know an existing method that solves this problem? Also, this is a repost from https://mathoverflow.net/questions/88171/efficiently-computing-a-few-localized-eigenvectors

Edit I think I can solve this if I can at least figure a way to solve \begin{equation} \max \psi^T H \psi \text{ subject to } P\psi = \psi \text{ and } \psi^T \psi = 1 \end{equation} where $P$ is projection onto the space of functions localized over $\Omega$. If this is doable then something like inverse iteration should be doable which will give me what I want.

Answer 1

Discretize by choosing a basis of the subspace of localized functions where you expect the eigenvectors of interest to be. These define the columns of $P$. Choose another basis of a slightly bigger subspace also containing slightly less localized functions and slightly more wiggly functions. These define the columns of $Q$.

Now calculate the matrices $H':=Q^*HP$ and $K':=Q^*P$, and solve the overdetermined eigenvalue problem $\|H'\phi-EK'\phi\|=\min$, using OEIG from http://www.mat.univie.ac.at/~neum/software/oeig/ If the residiuals are not small enough, you can always increase the basis sizes.

Answer 2

As per your edit you should be able to do the following: $$ G=H-\|H\|I $$ Now solve $$PGP\psi=\lambda \psi$$ for the first $n$ eigenpairs, where $\psi\in L^2(\Omega)$.

The shift ensures that $G$ only has non-positive Eigenvalues. The projector then maps all exterior functions to 0. Hence the lowest eigenvectors are your requested eigenvectors with respective eigenvalues $\lambda+\|H\|$.

EDIT:

Rewrite the problem: $$\max_{\lambda,\psi} \|P\psi\|^2 s.t. \|\psi\|^2=1, \|H\psi-\lambda\psi\|^2=0$$ Now you have a constrained non-linear problem in which all your functions are convex (in each variable).

Now you can use a non-linear optimizer (e.g., IPOPT). If you're using a mesh on $\mathbb{R}^2$, all involved operators should be sparse.

Answer 3

EDIT: I misunderstood the question. This answers a different question of spectral locality, not locality in the spatial domain.

I don't know a way to guarantee that you have found all eigenpairs in a particular region, but you can certainly ask for those with certain properties. The table below is from the SLEPc manual. Eigenvalues in the interior are more difficult to compute, so you will normally need to use methods like harmonic extension (available with some eigensolver methods, does not change the way the subspace is built) or shift-and-invert (more powerful, uses a transformed problem). Robust, scalable implementations of these methods are available within SLEPc. Note that MATLAB eigs() can be used for smaller problems. For finding a large number of interior eigenpairs, you might consider a method like Zhang, Smith, Sternberg, and Zapol (2007).