Chebyshev Accelerated Subspace Eigensolver for Pseudo-hermitian Hamiltonians

This paper extends the Chebyshev Accelerated Subspace iteration Eigensolver (ChASE) to efficiently compute thousands of the smallest positive eigenpairs of pseudo-hermitian Hamiltonians for excitonic materials by introducing an oblique Rayleigh-Ritz projection with quadratic convergence and a parallel implementation of the Chebyshev filter optimized for exascale systems.

The Gist

Imagine you are trying to find the most important "notes" in a massive, chaotic symphony played by a material (like silicon or molybdenum disulfide). These notes determine how the material interacts with light—whether it glows, conducts electricity, or absorbs energy.

In the world of physics, this symphony is represented by a giant mathematical object called a Hamiltonian. For many materials, this object has a special, tricky property: it's Pseudo-Hermitian.

Think of a standard "Hermitian" matrix like a perfectly balanced seesaw. If you push down on one side, the other side goes up in a predictable, symmetrical way. But a Pseudo-Hermitian matrix is like a seesaw where one side is painted black and the other white, and the rules of balance are slightly different. It has a "positive" side and a "negative" side that are mirror images of each other.

The problem? Scientists need to find the smallest positive notes (the lowest energy levels) in this symphony. But because of the tricky "black and white" nature of the matrix, standard tools get confused. They might accidentally listen to the "negative" notes or get stuck in the middle of the spectrum, unable to hear the quiet, important notes at the bottom.

Here is how the authors of this paper fixed the problem, using a clever new version of a tool called ChASE.

1. The Problem: The "Middle of the Spectrum" Trap

Imagine the notes of the symphony are arranged on a long line from very low (negative) to very high (positive).

• Standard tools are great at finding the very loudest notes at the far ends.
• The Challenge: The scientists need the quietest positive notes, which are hiding right in the middle of the line, surrounded by a sea of negative notes.
• The Trap: If you try to use a standard filter to pick out the positive notes, the "negative" mirror images get in the way, making the search slow and messy.

2. The Solution: The "Squaring Trick"

The authors realized that if you square the music (mathematically speaking), something magical happens.

• If you take a negative number (like -5) and square it, it becomes positive (25).
• If you take a positive number (like 5) and square it, it stays positive (25).

By applying this "squaring" trick to their math, they turned the entire symphony into a purely positive one. Now, the "negative" notes and the "positive" notes both look like positive numbers. This allows them to use a powerful filter (called a Chebyshev Filter) to easily isolate the smallest numbers they care about, without getting confused by the negative ones.

3. The Shortcut: The "Magic Mirror"

Even with the squaring trick, the math is still twice as big as it needs to be because of that "black and white" mirror structure. Doing all the calculations for the whole thing would be like trying to clean a house by cleaning every room twice.

The authors discovered a Magic Mirror shortcut.

• Because the negative notes are just a mirror reflection of the positive ones, they don't need to calculate the negative side from scratch.
• Once they calculate the "positive" side, they can just look in the mirror (using a simple sign-flip operation) to instantly know what the "negative" side looks like.
• Result: They do half the work, but get the full picture. This saves massive amounts of time and computer power.

4. The "Oblique" Lens: Seeing Clearly

In the final step, they need to extract the exact notes from the filtered sound. Standard tools use a "straight-on" lens (Orthogonal Projection). But because of the weird "black and white" rules of this material, a straight-on lens blurs the image.

The authors invented a new "Oblique" lens.

• Imagine looking at a reflection in a funhouse mirror. A straight-on view makes it look distorted. But if you tilt your head at a specific angle (an oblique angle), the reflection suddenly snaps into perfect focus.
• This new mathematical lens allows them to separate the positive notes from the negative ones perfectly, ensuring the computer converges (finds the answer) incredibly fast.

5. The Result: Super-Fast on Super-Computers

They tested this new method on the JUPITER supercomputer, which uses thousands of powerful graphics cards (GPUs) working together.

• Speed: They could find thousands of the smallest energy notes in just a few seconds.
• Scale: They solved problems with matrices so large they would crash standard computers.
• Efficiency: The method scales beautifully. Whether you use 4 computers or 256, it keeps working efficiently, just like a well-organized assembly line.

The Big Picture

Before this paper, simulating how new materials interact with light was slow and often required making rough approximations that missed important details.

This new method is like giving scientists a high-speed, noise-canceling headset that can instantly tune into the specific, quiet frequencies of a material's electronic structure. It allows them to design better solar cells, faster computer chips, and more efficient batteries by simulating them with extreme accuracy and speed.

In short: They took a messy, confusing math problem, used a "squaring" trick to simplify it, a "mirror" trick to cut the work in half, and a "tilted lens" to see the answer clearly, resulting in a tool that is fast, scalable, and ready for the next generation of materials science.

Technical Summary

1. Problem Statement

The paper addresses the computational challenge of solving the Bethe-Salpeter Equation (BSE) for excitonic materials. Numerically, the BSE is formulated as an eigenvalue problem involving a pseudo-Hermitian Hamiltonian ($H$) of size $n \times n$ (where $n=2m$).

• Structure: The Hamiltonian has a specific block structure:
  $$H = \begin{bmatrix} A & B \\ -\bar{B} & -\bar{A} \end{bmatrix}$$
  where $A$ is resonant, $B$ is the coupling term, and the matrix satisfies $SH = H^*S$ with $S = \text{diag}(I, -I)$.
• The Challenge:
  • Spectral Location: Unlike standard Hermitian problems where target eigenvalues are at the spectrum's edge, the smallest positive eigenvalues of $H$ lie in the middle of the spectrum (due to positive-negative symmetry).
  • Scale: Materials science applications require computing thousands of the smallest positive eigenpairs.
  • Limitations of Existing Methods:
    • Direct Solvers (e.g., ELPA): Too expensive ($O(n^3)$) for large systems when only a fraction of the spectrum is needed.
    • Standard Iterative Solvers (e.g., Lanczos, SLEPc): Struggle with the middle-of-spectrum location and the non-Hermitian nature of $H$. They often fail to scale efficiently to thousands of eigenpairs due to reorthogonalization costs.
    • Tamm-Dancoff Approximation (TDA): Neglects the coupling term $B$ to reduce $H$ to a Hermitian block $A$. While faster, this approximation often yields inaccurate optical properties.

2. Methodology: Extending ChASE

The authors extend the Chebyshev Accelerated Subspace iteration Eigensolver (ChASE), originally designed for Hermitian matrices, to handle pseudo-Hermitian Hamiltonians. The extension involves five key algorithmic modifications:

A. Spectral Folding via $H^2$

Since the target eigenvalues are near zero (in the middle of the spectrum), standard Chebyshev filtering is inefficient.

• Strategy: The algorithm implicitly filters using $H^2$.
• Effect: This folds the spectrum at zero, mapping both positive and negative eigenvalues to the same positive range. The smallest positive eigenvalues of $H$ become the smallest eigenvalues of $H^2$, allowing standard Chebyshev filtering to target them effectively.

B. Symmetry Exploitation and Search Space Reduction

The pseudo-Hermitian structure implies that if $v$ is a right eigenvector for $\lambda$, then $K\bar{v}$ is an eigenvector for $-\lambda$ (where $K$ is a specific permutation matrix).

• Optimization: Instead of filtering the full search space $W = [W_+, W_-]$, the algorithm filters only the positive half ($W_+$).
• Recovery: The negative half ($W_-$) is recovered analytically via the relation $W_- = K W_+$. This reduces the computational cost of the filtering step by roughly half compared to a naive approach.
• Communication Avoidance: The matrix-vector product $H^2 W_+$ is implemented using the relation $H = SH^*S$. This allows processes to operate on local data with minimal global communication (only sign-flipping operations are needed), leveraging the massive parallelism of GPUs.

C. Oblique Rayleigh-Ritz Projection

Standard orthogonal Rayleigh-Ritz fails for non-Hermitian matrices because the right eigenvectors are not orthogonal.

• Dual Basis Construction: The authors propose an oblique variant using a dual basis $Q_L = S Q (Q^* S Q)^{-1}$.
• Hermitian Equivalence: This construction transforms the reduced problem into a Hermitian eigenvalue problem ($G Z = Z \tilde{\Lambda}$) where $G = (Q^* S Q)^{-1} Q^* S H Q$.
• Benefit: This allows the use of highly optimized Hermitian solvers (like HEEVD) for the reduced problem and preserves the positive-negative spectral symmetry.

D. Quadratic Convergence Proof

A critical theoretical contribution is the proof that this oblique variant achieves quadratic convergence for the Ritz values (similar to the Hermitian case).

• The proof relies on showing that the chosen dual basis $Q_L$ ensures the projection errors for left and right eigenvectors are balanced ($O(\sigma_u) \approx O(\sigma_v)$) and that the overlap factor $|\tilde{\delta}|^{-1}$ remains bounded.

E. Subspace Initialization

To ensure numerical stability, the initial random subspace is constrained to the S-positive manifold (where $v^* S v > 0$). This is achieved by bounding the ratio of the norms of the upper and lower blocks of the vectors, preventing the Rayleigh-Ritz procedure from encountering singularities.

3. Key Contributions

1. Algorithmic Extension: Successfully generalized the ChASE solver from Hermitian to pseudo-Hermitian matrices, enabling the calculation of thousands of eigenpairs.
2. Efficiency via Symmetry: Developed a strategy to filter only half the search space by exploiting the $K$-conjugate symmetry, significantly reducing computational load and memory bandwidth requirements.
3. Communication Avoidance: Adapted the parallel matrix-matrix multiplication for the pseudo-Hermitian structure to minimize global communication on distributed GPU clusters.
4. Theoretical Guarantee: Proved that the proposed oblique Rayleigh-Ritz scheme preserves quadratic convergence of Ritz values, a property usually lost in non-Hermitian settings.
5. Robust Implementation: Created a production-ready implementation on distributed GPU systems (JUPITER supercomputer) handling matrices up to $n \approx 10^5$.

4. Results

The authors tested the solver on six large-scale pseudo-Hermitian Hamiltonians (Silicon and Molybdenum Disulfide systems) ranging from $n \approx 3,000$ to $n \approx 105,000$.

• Convergence: The solver converged in fewer than 25 iterations for all test cases (often <10), even when targeting 3% of the spectrum.
• Performance (Strong Scaling):
  • Tested on the JUPITER supercomputer (NVIDIA Grace-Hopper GH200 nodes).
  • Throughput: Achieved up to 4.6 PFLOP/s on 256 GPUs.
  • Time-to-Solution: Computed 3,144 eigenpairs for a matrix of size $n=104,832$ in 37.3 seconds.
  • Comparison: Significantly outperformed SLEPc (Lanczos-based) and ELPA (direct solver) for this specific regime (thousands of eigenpairs, partial spectrum).
• Scalability: Demonstrated strong scaling efficiency up to 256 GPUs, maintaining ~30% parallel efficiency even at large scales, with effective FLOP/s exceeding 60% of the theoretical peak on single nodes.

5. Significance

This work bridges a critical gap in computational materials science. By enabling the efficient, scalable solution of the full pseudo-Hermitian Bethe-Salpeter equation without the Tamm-Dancoff approximation, it allows researchers to simulate optoelectronic properties of complex materials with higher accuracy. The ability to compute thousands of eigenpairs on modern exascale GPU architectures makes previously intractable large-scale simulations feasible, offering a robust alternative to both expensive direct diagonalization and unstable iterative methods.