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.

  • $\begingroup$ Is $\Omega$ a subspace? $\endgroup$ Feb 23, 2012 at 5:43
  • $\begingroup$ 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}$. $\endgroup$
    – dranxo
    Feb 23, 2012 at 5:54
  • $\begingroup$ 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$. $\endgroup$
    – dranxo
    Feb 23, 2012 at 6:04
  • $\begingroup$ 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. $\endgroup$ Feb 23, 2012 at 6:15
  • $\begingroup$ I think I see how my wording was confusing. I edited the problem statement for clarity. $\endgroup$
    – dranxo
    Feb 23, 2012 at 6:28

3 Answers 3


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.

  • $\begingroup$ Any problem that needs to be done 'efficiently' should not be solved with Matlab. $\endgroup$
    – user563
    Mar 26, 2012 at 5:14
  • $\begingroup$ @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. $\endgroup$ Mar 26, 2012 at 18:05
  • $\begingroup$ 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. $\endgroup$
    – user563
    Mar 26, 2012 at 18:38
  • $\begingroup$ One can translate matlab automatically to C. $\endgroup$ Mar 26, 2012 at 18:41
  • $\begingroup$ And it is completely inscrutable C code. Almost completely useless for any reasonably-sized program. $\endgroup$
    – user563
    Mar 26, 2012 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\|$.


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.

  • $\begingroup$ Does this behave well when there are degenerate eigenvalues? $\endgroup$
    – qubyte
    Feb 23, 2012 at 18:17
  • $\begingroup$ 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. $\endgroup$ Feb 23, 2012 at 18:30
  • $\begingroup$ Cool, always nice to have some extra tricks up my sleeve! $\endgroup$
    – qubyte
    Feb 23, 2012 at 18:33
  • $\begingroup$ I can see that eigenvectors of $PGP$ give me eigenvectors of $PG$. I don't see how they give me eigenvectors of $G$. $\endgroup$
    – dranxo
    Feb 23, 2012 at 22:15
  • $\begingroup$ @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$$ $\endgroup$ Feb 24, 2012 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).

SLEPc eigenvalue selection

  • $\begingroup$ 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. $\endgroup$
    – dranxo
    Feb 23, 2012 at 21:45
  • 1
    $\begingroup$ I thought that region refers to part of the domain of the differential equation, not to a region in the spectral plane. $\endgroup$ Mar 26, 2012 at 18:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.