I'm interested in numerically diagonalizing a class of structured, symmetric eigenvalue problems with potentially extremely small eigenvalue gaps. The question I have is how to design a numerically stable way to get the eigenvalues and eigenvectors.

While usually small gaps destroy precision in numerical eigenproblems, for reasons I will explain below the structure of these problems allow for accurate solutions in principal, i.e. the eigenvalues and eigenvectors don't change dramatically as a function of perturbations of the non-zero matrix elements. An explicit example will follow, but the reason for the stability is the following: the structure of the small parameters is such that they can always be brought to the block off-diagonal form

$H = \left( \begin{array}{c|c} D_0 & \epsilon \\ \epsilon^T & D_1 \end{array} \right),$

where $D_0$, $D_1$ are diagonal and $\epsilon$ is small compared to the level spacing in $D_0$ and $D_1$. In this case, perturbation theory says that the eigenvalues $\lambda_i$ are close to those of $D_0$ or $D_1$, i.e.

$|\lambda_i - d_i | < \frac{||\epsilon||^2}{d_i - d_j} $.

The eigenproblems of interest come from Hamiltonians of very short quantum spin chains, for example

$H = -h_1 \sigma^z_1 - h_2\sigma^z_2 - h_3\sigma^z_3 - J_1 \sigma^x_1 \sigma^x_2 - J_2 \sigma^x_2 \sigma^x_3,$


$\sigma^z_1 = kron(\sigma^z, I, I),$

$\sigma^x_2 = kron(I, \sigma^x, I),$


Thus the full matrix takes a form like

$H = \left( \begin{array}{cccccccc} -d_0 & 0 & 0 & -J_1 & 0 & 0 & -J_2 & 0 \\ 0 & -d_1 & -J_1 & 0 & 0 & 0 & 0 & -J_2 \\ 0 & -J_1 & -d_2 & 0 & -J_2 & 0 & 0 & 0 \\ -J_1 & 0 & 0 & -d_3 & 0 & -J_2 & 0 & 0 \\ 0 & 0 & -J_2 & 0 & d_3 & 0 & 0 & -J_1 \\ 0 & 0 & 0 & -J_2 & 0 & d_2 & -J_1 & 0 \\ -J_2 & 0 & 0 & 0 & 0 & -J_1 & d_1 & 0 \\ 0 & -J_2 & 0 & 0 & -J_1 & 0 & 0 & d_0 \\ \end{array} \right),$

where $d_0 = h_1 + h_2 + h_3$, $d_1 = -h_1 + h_2 + h_3$, $d_2 = h_1 - h_2 + h_3$, $d_3 = h_1 + h_2 - h_3$.

The ratios of the parameters $h_1, h_2, h_3, J_1, J_2$ are potentially much smaller than machine precision. This leads to very small gaps in the eigenvalues and trouble for the typical numerical eigensolvers working at finite precision. However, I do not think it should be impossible to determine the eigenvectors, because when the ratios are large we have precisely the situation where perturbation theory can be used to approximately determine the eigenvalues and eigenvectors in the following manner:

For example, if $h_2$ is the largest scale, then perturbation theory in the small parameters $J_i/h_2$ shows that the eigenvalues are very close to the eigenvalues of these smaller matrices:

$ \tilde{H}_{\pm} = \pm h_2 I - h_1 \sigma^z_1 - h_3 \sigma^z_2 \pm \frac{J_1 J_2}{h_2} \sigma^x_1 \sigma^x_2 + \mathcal{O}(\frac{J_i^4}{h_2^3}).$

This argument can be repeated again if one of the remaining coefficients $h_1, h_3, \frac{J_1 J_2}{h_2}$ is large compared to the others. For example, if $\frac{J_1 J_2}{h_2}$ is the largest remaining scale, the eigenvalues will be very close to

$ \lambda_{\pm, \pm, \pm} \approx \pm h_2 \pm \frac{J_1 J_2}{h_2} \pm \frac{h_1 h_2 h_3}{J_1 J_2}.$ The eigenvectors are also accurately determined by perturbation theory.

My goal is to accurately determine a single eigenvector of these matrices. When using numerical eigensolvers, the issue is as follows: when the small gaps (in the example $2\frac{h_1 h_2 h_3}{J_1 J_2}$) are less than machine precision, a numerical eigensolver will mix the eigenstates. I could just use perturbation theory in this case. In the situations where the ratios of the parameters are not so small, perturbation theory fails to be accurate but numerical eigensolvers have no problem. Is it possible to make an eigensolver that captures both cases accurately?

P.S. For the specific test problem chosen here, a mapping to free fermions shows that the exact eigenvalues are given by the formula $ \lambda_{\pm, \pm, \pm} = \pm \sigma_1 \pm \sigma_2 \pm \sigma_3$ where $\sigma_i$ are the singular values of the matrix

$M = \left( \begin{array}{ccc} h_1 & J_1 & 0 \\ 0 & h_2 & J_2 \\ 0 & 0 & h_3 \\ \end{array} \right).$

This isn't true for the general case I am interested in but can serve as a useful check that the eigensolver is actually accurate.

  • 4
    $\begingroup$ What about using high (arbitrary precision) arithmetic? Do you have computation time or memory restrictions which would preclude that? For example (not an endorsement) advanpix.com/2011/10/12/… . $\endgroup$ Commented Jun 5, 2018 at 18:41

1 Answer 1


Is it possible to make an eigensolver that captures both cases accurately?

In general no, for the following reason. When a typical eigensolver produces approximate eigenvalues of a matrix $A$, the answer is usually backward stable: the output eigenvalues are the exact eigenvalues of a nearby matrix $A+\delta A$ where $\|\delta A\|\leq \epsilon_{\mathrm{mach}}\|A\|$. This is also what causes them to produce inaccurate almost-double eigenvalues.

In your case, you have a special matrix $H$ that can only be perturbed in specific ways: as far as an eigensolver is concerned, the nearby matrix can contain 64 independent perturbation, one per matrix element, but in your matrix only the 5 parameters can be perturbed and the space of potential perturbations $\delta H$ is much smaller than what is allowed by $\|\delta H\| \leq \epsilon_{\mathrm{mach}}\|H\|$.

This means that if you did have an eigensolver that could handle this, it would have to be designed in advance with knowledge of the special structure of your matrix. This is in fact what you've done already: you've described an eigensolver that dispatches either to a standard solver or to perturbation theory. I think this is probably the best you can do anyway, short of having some closed form solution or something.


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.