Dual application of Chebyshev polynomial for efficiently
computing thousands of central eigenvalues in many-spin systems

Energy level statistics provides an essential characterization of quantum chaos Haake (2010); Brody et al. (1981). Integrable systems often imply level clustering and a Poisson distribution of energy level spacings Brody et al. (1981), while chaotic systems exhibit strong level repulsion and a Wigner-Dyson distribution Bohigas et al. (1984). Other useful statistical tools like the $\delta_{3}$ statistic Mehta (1991) and the power spectrum of the $\delta_{n}$ statistic also depend on the level distribution Relaño et al. (2002). Numerical observation of eigenvalues is important to exactly characterize these level statistics, because direct theoretical derivation is not available. In addition, the need to access individual eigenstates in the middle of the spectrum, which correspond to the “infinite temperature” limit, is emphasized in studying the phenomenon of many-body localization (MBL) Basko et al. (2006); Nandkishore and Huse (2015); Pietracaprina et al. (2018). Numerical simulations remain instrumental in understanding quantitatively many aspects of the MBL problem. As the many-body problems of interest often involve a huge Hilbert space dimension growing exponentially in terms of system size, it is rather challenging to conduct a full diagonalization of the Hamiltonian or to solve the time-independent Schrödinger equation. Seeking to resolve the eigenvalue problem in a small part of the spectrum is thus an unavoidable and desirable substitution.

To obtain the interior eigenstates, a hybrid strategy of matrix spectroscopy is often invoked Ericsson and Ruhe (1980); Wyatt (1995); Minehardt et al. (1997). It aims at computing eigenstates in selected regions of the spectrum and couples the ground state solvers with a spectral filter, where the filter is designed to cast the selected interior regions to the edges of the spectrum. Among these filters, the Green function $(\mathcal{H}-\lambda I)^{-1}$ is an excellent one. After applying the Green function filter, the cluster of eigenvalues near the test energy $\lambda$ is mapped to very large positive and negative values. The level spacings near $\lambda$ are greatly amplified, which thus improves the convergence. The Lanczos method Lanczos (1950) has been coupled with such a filter Ericsson and Ruhe (1980); Wyatt (1995), and the Chebyshev polynomial expansion of the Green function is implemented Mandelshtam and Taylor (1995). In particular, the shift-invert method Bai et al. (2000) essentially utilizes this spectral transformation and has been widely used in quantum many-spin systems Pietracaprina et al. (2018); Luitz et al. (2015); Sierant and Zakrzewski (2020); Hopjan and Heidrich-Meisner (2020); Šuntajs et al. (2020). Moreover, it was even proved to be the most efficient one for the MBL problem Pietracaprina et al. (2018). However, for large systems this method suffers rapid increases in both computation time and memory consumption, due to the factorization of $(\mathcal{H}-\lambda I)^{-1}$. Other spectral filters have been also proposed, including the Dirac delta function filter $\delta(\mathcal{H}-\lambda I)$ Li et al. (2016); Fang and Saad (2012); Sierant et al. (2020); Guan and Zhang (2020) and the rectangular window function filter Pieper et al. (2016), both of which are expanded by the Chebyshev polynomial. Note that all the methods mentioned above are iterative ones, where the computation time is usually proportional to the required number of eigenvalues, i.e., more eigenvalues means more filtrations and reorthogonalizations Pietracaprina et al. (2018); Fang and Saad (2012); Sierant et al. (2020).

In this paper, we propose an exact numerical method, DACP, to calculate thousands of eigenvalues in the middle of the energy band, which is already helpful to reveal the level statistics Shklovskii et al. (1993); Georgeot and Shepelyansky (1998). For spin systems, such a region usually locates at the peak of the density of states where the energy levels are so close that they are nearly degenerate. It is extremely challenging to distinguish these near-degenerate eigenvalues without factorizing $\mathcal{H}^{-1}$. In the DACP method, we construct an exponential of semicircle (exp-semicircle) filter to quickly damp the unwanted part of the spectrum by employing the Chebyshev polynomial for the first time. The second application of the Chebyshev polynomial is to fast search a set of states to span the specific subspace consisting of all the eigenstates corresponding to the desired eigenvalues. The DACP method transforms the original high dimensional eigenvalue problem to a low dimensional one. For practical problems in many-spin systems, the DACP method is very efficient, due to its full exploration of several excellent properties of the Chebyshev polynomial while other methods usually utilize only one of them. For a large class of many-spin systems, the DACP method exhibits a significant increase in the computation speed, up to a factor of $30$, in comparison with the shift-invert method. The memory saving is more drastic, up to a factor of $100$. Moreover, the DACP method distinguishes itself from those iterative filtering ones, as its convergence time varies only slightly when the required number of eigenvalues is vastly changed.

The paper is organized as follows. The detailed formalism of the DACP method, including the exp-semicircle filtration, Chebyshev evolution, and subspace diagonalization, are described in Sec. II. Each of the processes relies on a property of the Chebyshev polynomial. In Sec. III we discuss details of two many-spin systems, the Ising model and the spin glass shards, and present the computational results. We conclude in Sec. IV.

To access central eigenvalues of large spin systems, we are restricted by the matrix-free mode, i.e., the matrix of the Hamiltonian must not be explicitly expressed/stored. Instead, we treat the Hamiltonian $\mathcal{H}$ as a function whose input and output are both states/vectors. Therefore, we shall only operate with the quantum states $\left|\psi\right>,\mathcal{H}\left|\psi\right>,\cdots,\mathcal{H}^{k}\left|\psi\right>$, where $k$ is a positive integer (not too large), in the original Hilbert space of dimension up to $5\times 10^{5}$. In this paper, we set our goal to find $5,000$ eigenvalues in the middle of the spectrum for many-spin systems as an illustration.

The DACP idea is fairly straightforward. We first transform a randomly initialized state into a wave packet in the subspace spanned essentially by at least $5,000$ central eigenstates. This transformation is realized by an exp-semicircle filter implemented by the Chebyshev polynomial. The Chebyshev polynomial is used to perform an exponential decay. With this particular state in hand, we then generate a large amount of linearly independent states, as large as possible, to approximately span the subspace consisting of the required eigenstates. The Chebyshev polynomial is used again to oscillate (complex exponential) the states, producing approximately linearly independent states. Once the generating set for the desired eigenstates is known, one may explicitly calculate $H$, the reduced representation of $\mathcal{H}$ in this subspace. The remaining operations are restricted in the subspace of dimension around $10^{4}$. Finally, direct diagonalization of $H$, which is of size $10^{4}\times 10^{4}$, gives the required eigenvalues. The detailed procedures and discussions are given in the subsections below.

We apply the DACP method to the quantum spin-$1/2$ systems with two-body interactions. Such systems are good models for investigating a large class of important problems in quantum computing, solid state theory, and quantum statistics Nielsen and Chuang (2011); Bogolubov and N. N. Bogolubov (2009); Sachdev (2011); Dutta et al. (2015). A large amount of exact eigenvalues helps us to obtain the statistical properties, to distinguish quantum chaos from integrability, and serves as a benchmark to evaluate other approximate methods as well.

Generally speaking, the DACP method can deal with the spin systems consisting of couplings between any pair of the $N$ spins. Each of the Pauli matrix $\sigma^{\alpha}$ or the two coupling Pauli matrices $\sigma^{\alpha}\otimes\sigma^{\beta}$, where $\alpha,\beta=x,y,z$, is properly represented by a specific function.

We specify the spin model for two physical systems. One is the disordered one-dimensional transverse field Ising model Fisher (1995), where the Hamiltonian is

with $\sigma_{i}$ the Pauli matrices for the spin $i$. This system is exactly solvable by Jordan-Wigner transformation Sachdev (2011), making it an ideal correctness checker for the DACP method. The nearest neighbor exchange interaction constants $J_{i,i+1}$ are random numbers that uniformly distributed in $[-J/\sqrt{N},J/\sqrt{N}]$ with $J=10$. The local random magnetic fields are represented by $\Gamma^{z}_{i}$, which are random numbers that uniformly distributed in the interval $[0,\Gamma]$ with $\Gamma=1$.

Another system is the spin glass shards Georgeot and Shepelyansky (1998), which represents a class of global-range interacting systems that require relatively large bond dimensions to be tackled by the DMRG methods Schollwöck (2005). The Hamiltonian describing the system is

All symbols and parameters are the same as that of the above Ising model, except that the first summation runs over all possible spin pairs. This system is interesting because it presents two crossovers from integrability to quantum chaos and back to integrability again. In the limit $J/\Gamma\rightarrow 0$, the ground state is paramagnetic with all spins in the local field direction and the system is integrable Georgeot and Shepelyansky (1998). In the opposite limit $J/\Gamma\rightarrow\infty$, the ground state is spin glass and the system is also integrable since there are $N$ operators ($\sigma_{i}^{x}$) commuting with the Hamiltonian. A quantum chaos region exists between these two limits. $J=10\,\Gamma$ is approximately the border from the quantum chaos to the integrable (the spin glass side) when $N=20$ Georgeot and Shepelyansky (1998).

By employing the upper-bound-estimator, which costs little extra computation and bounds up the largest absolute eigenvalue $E_{0}$, one may estimate $E_{\rm max}=E_{0}$ and $E_{\rm min}=-E_{0}$ Zhou and Li (2011). For this setting we have utilized the symmetry of the density of states (DOS), a bell-shape profile centered at zero, in the many-spin systems. Since we require $5,000$ central eigenvalues, we may set $n=4,000$, corresponding to a dimension $8,001$ and being adequate to span the whole subspace $\mathbb{L}$. The overlap matrix $S$ is generally singular. The approximate distribution of DOS $\rho(E)$ may be efficiently calculated through the Fourier transformation of a time evolved wave function or through a better estimation method given in Hams and De Raedt (2000). The parameter $a$ is appropriately chosen to ensure that the number of eigenstates contained in $[-a,a]$ is a little less than $8,000$ (as illustrated in Fig. 2, the precision of some converged eigenvalues may be lower than required).

In practice, sometimes there may exist highly near-degenerate eigenvalues, with level spacings as small as $10^{-7}\Gamma$ while the average spacing is $10^{-5}\Gamma$. It is still hard (two magnitudes longer time) for the Chebyshev evolution to discriminate such close pair of eigenvalues. To circumvent this challenge, we employ the block filter technique Zhou (2010), which means a block of states is filtered or evolved “simultaneously”, in programming of the DACP method. The idea is that two or several random states in the degenerate subspace are usually linearly independent. For each numerical test, a block of $5$ initial trial states is randomly generated and employed with the parameter $n$ being adjusted to $n=800$ accordingly.

By these settings, we perform numerical tests on the above two systems to show the high exactness and efficiency of DACP method. For this work, we consider only the eigenvalues computations. All the timing information reported in this manuscript is obtained from calculations on the Intel(R) Xeon(R) CPU E7-8880 v4, using sequential mode.

In Fig. 2 we present the relative error $\eta$ in logarithmic scale versus the system energy $E_{\rm exact}$, for the Ising model with $N=19$ and the spin glass shards with $N=17$. We have defined the relative error $\eta$ of the computed central eigenvalues $E$ as

Exact eigenvalues of both systems have been obtained by other reliable methods. For the Ising model, we make use of the famous Jordan-Wigner transformation to reduce the original $2^{N}\times 2^{N}$ matrix to a $2N\times 2N$ one, and restore the full spectrum of the original Hamiltonian Sachdev (2011). For the spin glass shards we simply utilize the function $eigs$ of MATLAB, to find $5,000$ eigenvalues closest to $E=0$. As for our numerical tests, the parameter $a$ in Fig. 2 is $0.036\Gamma$ and $0.16\Gamma$ for (a) and (b), respectively. In computing the Ising model, the number of eigenvalues satisfying $\eta<10^{-6}$ is not enough (less than $5,000$) by the settings mentioned above. By adjusting the block size to $10$ and the parameter $n$ to $500$, we then collect enough eigenvalues. The number of converged eigenvalues, i.e., computed eigenvalues satisfying the condition $\eta<10^{-6}$ is $5,385$ for (a) and $5,000$ (all the exact eigenvalues we have) for (b), while the total number of computed eigenvalues is $6,232$ and $5,910$, respectively.

The spike around $E=0$ for both figures is due to the smallness of the denominator $E_{\rm exact}$. The smallest absolute eigenvalue is about $4.4\times 10^{-6}\Gamma$ for (a) and $2.9\times 10^{-5}\Gamma$ for (b). Besides, there is a flat plateau in the middle of the figures, indicating that for those eigenvalues around $E=0$ we encounter the numerical error, i.e., the absolute error reaches the limit of the double precision representation. In Fig. 2(a) one may observe other spikes, which does not appear in (b). Each of those spikes corresponds to a cluster of eigenvalues lying extremely close to each other. The integrability of the Ising model implies the level clustering and a Poissonian distribution of energy level spacings. We speculate that the Chebyshev evolution does not function well for closely lied eigenvalues, as they give quite similar phase contributions to the final states. To significantly amplify the phase difference $\exp({\rm i}\Delta Et)$, it requires $t\gtrsim 1/\Delta E$, where $\Delta E$ is the energy difference. An alternative solution is to increase the block size to match up with the dimension of near-degenerate subspace. Ignoring the central plateau and the spikes, the distribution of $\eta$ shows an inverted shape of the exp-semicircle filter.

In Fig. 3 we compare the computation time $T$ (CPU time in seconds) of the DACP method with that of the shift-invert approach. The latter is widely used in computing eigenpairs at the middle of the spectrum for quantum spin systems, to name a few, like in Pietracaprina et al. (2018); Luitz et al. (2015); Sierant and Zakrzewski (2020); Hopjan and Heidrich-Meisner (2020). In our tests, it is implemented by the $eigs$ function of Matlab R2019b, which employs the implicitly restarted Arnoldi method (ARPACK) Lehoucq et al. (1998). Each of the numerical tests finds around $5,000$ eigenvalues for either of the two spin systems. For the DACP method, we have discarded the time consumed during the subspace diagonalization, as it is a constant and does not affect the scaling behavior. Specifically, the subspace diagonalization spends roughly $700$ CPU seconds in each test. Such a constant is completely negligible for $N\geq 17$ systems.

As shown in Fig. 3, for spin systems with $13\leq N\leq 17$, the DACP method is about $30$ times faster than the shift-invert approach for the Ising model and $8$ times faster for the spin glass shards. Since the shift-invert method essentially finds the ground energy of $\mathcal{H}^{-1}$, which is usually not a sparse matrix, its computation time $T$ is not affected by the sparsity of the Hamiltonian and the two dashed lines nearly coincide. On the contrary, the DACP method employs the polynomial combination of $\mathcal{H}$ acting on the states, its computation time $T$ is heavily influenced by the number of Pauli operators of the specific Hamiltonian. For example, when $N=18$ the computation time $T$ for the spin glass shards is $5.5$ times that for the Ising model, while the number of Pauli operators is $6.2$ times. Considering this effect, the DACP method is still advantageous over the shift-invert approach in the worst case where the exchange interactions run over all possible spin pairs.

Furthermore, the scaling among these two methods is comparable. To compare quantitatively, we extract the scaling constants $\alpha$ of the four lines by fitting the numerical results, where $\alpha$ is defined by $T=T_{0}\exp(\alpha N)$. The values of $\alpha$ are shown in Table 1. Indeed, the four scaling constants are quite close. Suppose that $\alpha$ remains unchanged for a bigger $N$, then the DACP method’s advantage in time cost may keep on for rather large systems. In addition, as illustrated by the numerical tests of the shift-invert method in Pietracaprina et al. (2018), the factorization time for finding $\mathcal{H}^{-1}$ dominates other computation steps. Considering the factorization part only, the execution time in Pietracaprina et al. (2018) exhibits a scaling constant $\alpha\simeq 1.66$, indicating a worse scaling behavior for systems with large spins. The efficiency of the DACP method is thus confirmed.

In Fig. 4 we show the computation time versus the number of converged eigenvalues, for the two systems of size $N=15$. The horizontal axis represents the number of computed eigenvalues satisfying the condition $\eta(E)<10^{-6}$ in each test. Recall that the DACP method is divided into three parts: the exp-semicircle filtration, Chebyshev evolution, and subspace diagonalization. We denote their time consumption as $T_{F}$, $T_{E}$, and $T_{D}$, respectively, while solving time is the full computation time $T=T_{F}+T_{E}+T_{D}$. One immediately notice that for the two systems, computation time $T$ decreases on the left side, reaching a plateau at the middle region, while quickly increases on the right side. Yet the dashed lines $T_{E}$ present a quite small increasing rate for the whole region, as the horizontal axis is in logarithmic scale. In the middle region, where the DACP method performs the best, the solving time $T$ and the evolution time $T_{E}$ roughly coincide with each other.

These features are easily understood if one considers the following four facts. (i) Notice that the DOS $\rho(E)$ is usually a bell-shape profile peaked at zero in spin systems and that $a$ is typically a tiny parameter ($a/E_{\rm max}\sim 10^{-2}$ for $N=19$), thus one may safely take $\rho(E)$ as a constant $\bar{\rho}$ in $[-a,a]$. Consequently, the number of required eigenvalues is $R\simeq 2\bar{\rho}a$. (ii) Setting the action of the Hamiltonian on the state, $\mathcal{H}\left|\psi\right>$, as a basic step, and denoting the computation time of the basic step as $\tau$, we may count the filtration time as $\tau$ times twice the cut-off order $K=18E_{\rm max}/a$, i.e., $T_{F}\simeq 36E_{\rm max}\tau/a\propto R^{-1}$. (iii) Similarly, we find the Chebyshev evolution time is $T_{E}\simeq\lfloor R\pi E_{\rm max}/a\rfloor\tau\propto R^{0}$ since $R\propto a$. (iv) It is well known that the full diagonalization time is proportional to the cube of the matrix size, i.e., $T_{D}\propto R^{3}$. The combination of these computation times clearly explains the behavior shown in Fig. 4. For the left side, the exp-semicircle filtration dominates, as $T_{F}$ is inversely proportional to the parameter $R$. Whereas, when it comes to the intermediate region where the number of eigenvalues is several hundreds to thousands, the evolution time $T_{E}$ consists of the majority of the computation time while the other two are negligible. As shown in Fig. 4, the performance of DACP method is approximately the same in finding $100$ to $3,000$ eigenvalues. As $R$ keeps increasing, the subspace diagonalization time $T_{D}$ eventually consumes the most computation time.

We remark that the plateau in the middle region in Fig. 4 distinguishes the DACP method from those iterative filtering ones Wyatt (1995); Guan and Zhang (2020); Pieper et al. (2016); Zhou and Saad (2007); Zhou (2010); Fang and Saad (2012), although all of them take use of the Chebyshev polynomials. The iterative methods usually require a larger amount of filtration and reorthogonalization, thus a longer convergence time, in finding more eigenvalues Pietracaprina et al. (2018); Fang and Saad (2012); Sierant et al. (2020). In fact, as illustrated in Ref. Pietracaprina et al. (2018), one finds that the computation time is roughly proportional to the number of eigenvalues required for the shift-invert approach. The DACP method is thus highly desirable for large scale eigenvalues computations.

Finally, we note that the memory requirements by DACP method is rather small, as it works in a matrix-vector product mode which avoids an explicit matrix representation of the Hamiltonian. Moreover, as shown in Append. C, the whole set of states in Eq. (8) is not preserved in memory. The major memory consumption is the storage of elements $H_{ij}$ and $S_{ij}$, thus the memory requirements of the DACP method relates only to the dimension $d$ of subspace $\mathbb{L}$. By the settings of this work, it occupies around $5.6$ GB of memory for appropriately calculating $5,000$ eigenvalues in $N\leq 20$ systems. On the contrary, for the shift-invert method, the matrix size of $\mathcal{H}^{-1}$ grows rather fast as the system size $N$ increases (proportional to $4^{N}$), demanding a large amount of memory. For example, in our tests it consumes $573$ GB of memory for $N=17$ systems, which is already a factor of $100$ more than that of the DACP method.

We propose the DACP method to efficiently calculate a large scale (at least $5,000$) of central eigenvalues with high precision, by employing twice the Chebyshev polynomials. The proposed method is not a general eigensolver, as the generality gives way to the efficiency. It explores the excellent properties of Chebyshev polynomial to efficiently filter out those non-central parts of spectrum and to construct the appropriate basis states for the central spectrum. The numerical tests for the Ising model and the spin glass shards confirm the exactness and efficiency of the DACP method. Compared to the widely used shift-invert method, the DACP method gives a considerable increase in the speed of computations, for the Ising model up to a factor of $30$ while for the spin glass shards the increase in speed is less but still considerable (a factor of $8$). The memory requirement is drastically decreased, up to a factor of $100$ for $N=17$ spin systems and even more for larger systems. Moreover, as shown in Append. D, the DACP method is both more stable and more efficient than the Lanczos method applied with the exp-semicircle filter. As a powerful tool for central eigenvalues calculation, the DACP method may find potential applications in many physical problems, such as many-body localization in condensed matter physics Nandkishore and Huse (2015); Abanin et al. (2019); Luitz et al. (2015); Hopjan and Heidrich-Meisner (2020), and level statistics in quantum chaos Haake (2010); Georgeot and Shepelyansky (1998); Sierant and Zakrzewski (2020); Mondal et al. (2020).