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,$
where
$\sigma^z_1 = kron(\sigma^z, I, I),$
$\sigma^x_2 = kron(I, \sigma^x, I),$
etc.
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.
