# Efficiently computing a few localized eigenvectors

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 http://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 $$\max \psi^T H \psi \text{ subject to } P\psi = \psi \text{ and } \psi^T \psi = 1$$ 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.

-
Is $\Omega$ a subspace? –  Geoff Oxberry Feb 23 '12 at 5:43
I'm thinking of this as domain decomposition for $H$ where $\Omega$ is one of your domains. So yes. If you discretize a region where $H$ acts using a $10^3 \times 10^3$ grid and then label points in $\Omega$ you'll get a subspace of $\mathbb{R}^{10^6}$. –  dranxo Feb 23 '12 at 5:54
I guess the other way to think of the subspace question is to note that the sum/scale of two functions localized over $\Omega$ is again localized over $\Omega$. –  dranxo Feb 23 '12 at 6:04
I'm just wondering if $\Omega$ is sufficiently low-dimensional that you could take advantage of linear transformations to reduce the dimensionality of your problem even further. –  Geoff Oxberry Feb 23 '12 at 6:15
I think I see how my wording was confusing. I edited the problem statement for clarity. –  dranxo Feb 23 '12 at 6:28

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.

-
Any problem that needs to be done 'efficiently' should not be solved with Matlab. –  Mr. F Mar 26 '12 at 5:14
@EMS: It is not difficult to port the Matlab code to another language. But the expensive part is in the formation of $H'$ and $K'$, not in the part done with Matlab. –  Arnold Neumaier Mar 26 '12 at 18:05
As a practitioner, I disagree very much that it is easy to translate Matlab to, say, Python. Sure, there are tools that work for very basic Matlab code, such as OMPC, but this is not the typical case. I work in academia where many peers write lazy, bad code in Matlab because it provides more instant gratification. They use a ton of built-in functions that abstract away tons of critical implementation choices, and often their research conclusions are affected by this in ways not reported. –  Mr. F Mar 26 '12 at 18:38
One can translate matlab automatically to C. –  Arnold Neumaier Mar 26 '12 at 18:41
And it is completely inscrutable C code. Almost completely useless for any reasonably-sized program. –  Mr. F Mar 26 '12 at 18:42

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.

-
Does this behave well when there are degenerate eigenvalues? –  qubyte Feb 23 '12 at 18:17
The shift doesn't affect degeneracy and $\|H\|$ is larger than any eigenvalue. The projection introduces degeneracies at 0. The shift ensures that the sought eigenvalues are relatively far from 0. If your eigensolver has no issues with degeneracies this shouldn't be a problem. –  Deathbreath Feb 23 '12 at 18:30
Cool, always nice to have some extra tricks up my sleeve! –  qubyte Feb 23 '12 at 18:33
I can see that eigenvectors of $PGP$ give me eigenvectors of $PG$. I don't see how they give me eigenvectors of $G$. –  dranxo Feb 23 '12 at 22:15
@rcompton They don't unless $H$ and $PGP$ commute. Your remark implies to me that you really want to solve $$\max \frac{\|P\psi\|}{\|\psi\|} s.t. H\psi=\lambda\psi$$ –  Deathbreath Feb 24 '12 at 1:58

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

-
For eigenpairs in a given spectral region we are using sciencedirect.com/science/article/pii/S0021999110003451 . This problem is different as we would like to reduce computation further by localizing in spatial as well as spectral regions. –  dranxo Feb 23 '12 at 21:45
I thought that region refers to part of the domain of the differential equation, not to a region in the spectral plane. –  Arnold Neumaier Mar 26 '12 at 18:06