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

Edit: I reposted to http://scicomp.stackexchange.com/questions/1396/efficiently-computing-a-few-localized-eigenvectors#comment2200_1396 and hopefully clarified the problem statement.

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$. My guess is that this will end up looking like power iterations with a projection step built in between matrix applies. If this is doable then something like inverse iteration should be doable which will give me what I want.

-
From what I know, if you are willing to compute 10% of the spectrum, you might as well go all the way. Besides, 1000 x 1000 is quite small by today's standards---or am I missing something obvious? –  Suvrit Feb 11 '12 at 4:48
If no help here may be ask at scicomp.stackexchange.com –  Alexander Chervov Feb 11 '12 at 12:34
Suvrit, I've never heard of that 10% rule before. A full decomposition of a 1000x1000 matrix is probably possible on cellphones, but I need to compute these eigenvectors thousands of times as part of a larger computation so any speedup is useful. –  dranxo Feb 13 '12 at 17:42
Also, I was using "1000","100", and "10" simply because I thought it made the exposition clearer than "n","m","k" etc. The actual dimensions change depending on the simulation. –  dranxo Feb 13 '12 at 17:49
@rcompton: that 10% is a ballpark number---dense matrix computations build on BLAS3 kernels, which can give them an advantage. I'd say try out both ARPACK and LAPACK for the kinds of matrices you have (including embedding the eigendecomp into an inner loop), and then decide. –  Suvrit Feb 13 '12 at 19:41

To build on Federico's answer, why not run a restarted Lanczos iteration but compute harmonic Ritz vectors to get approximations to the interior eigenpairs? For something so small, though, why not just diagonalize, as has already be stated.

-
Yes Ritz vectors should improve over the Lanczos method. I suppose if the Lanczos idea can work then this will work better. –  dranxo Feb 13 '12 at 18:05
@rcompton, just to be clear, I am not talking about regular Ritz vectors. I am talking about harmonic Ritz vectors which yield approximations to eigenvectors associated to eigenvalues near the origin (the so-called interior eigenvalues). –  Kirk S. Feb 16 '12 at 6:48
@Kirk I've never worked with those before. Is there a standard reference? –  dranxo Feb 16 '12 at 23:36
@rcompton Absolutely! Here is a reference by Ronald Morgan that may be a good starting point: Computing interior eigenvalues of large matrices. Linear Algebra Appl. 154/156 (1991), 289–309. –  Kirk S. Feb 17 '12 at 3:08

Just a random idea:

The standard method for getting a small part of the spectrum in large and sparse symmetric problems is the restarted Lanczos method. Essentially, you run some iterations of the Lanczos method, then you check the eigenpairs that have been computed, keep some of them and throw away the rest. Typically the pairs to drop and keep are selected based on the eigenvalues (you want the ones with smallest or largest modulus, for instance), but in this case you could try to modify the method and keep the ones that are "localized" in your region of interest.

Problems:

1. I cannot tell you for sure that this would work --- as far as I know there is no easy way to tell to which eigenvalues Lanczos will converge, and it is known that it has a tendency to converge to those at the border of the spectrum, so your efforts to steer it away from selected eigenpairs could be useless.
2. As far as I know, hooks for the selection strategy are not present in the usual Lanczos implementations, so you may have to code it yourself.

In any case I agree with Suvrit's comment --- your dimensions are kind of borderline for Lanczos to be more effective than full diagonalization.

-
Actually, a friend of mine who does eigenvectors for a living insisted, that if we are willing to settle for 10% of the spectrum, we might as well compute the whole spectrum for $n$ even as large as 50,000. So I guess, unless absolutely short on storage, for a tiny $1000 \times 1000$ matrix, I would go ahead and compute the entire spectrum. It takes only 3 seconds on my 5 year old laptop to compute the full spectrum of a 1000 x 1000 matrix! –  Suvrit Feb 11 '12 at 19:36
As for the size, I'd like to work with larger matrices in the future. The Lanczos idea is interesting. I don't have a good intuition for how the eigenvectors converge so I'm hesitant to throw out candidates that are not localizing early on. –  dranxo Feb 13 '12 at 18:02