As computational power moves towards the exascale era, the complexity of scientific simulations continues to increase. A key challenge facing many scientific applications is the efficient and accurate solution of large-scale eigenvalue problems. One such effort to address this challenge is our ISOLV-BSE innovation study. Led by José E. Román from the Universitat Politècnica de València, ISOLV-BSE targets an important class of problems involving structured pseudo-Hermitian matrices in the context of the Bethe-Salpeter Equation (BSE) for material science simulations.

In an informative interview within the Inno4scale project, Román provided insights into the motivations, methods, challenges, and future potential of ISOLV-BSE. ISOLV-BSE aims to combine the material science code Yambo with the SLEPc library to develop a highly efficient solver for structured pseudo-Hermitian matrices, which can be deployed in exascale computing environments.

**A Collaboration of Expertise: Yambo and SLEPc**

The ISOLV-BSE project is the result of a collaboration between two institutions: the Universitat Politècnica de València and the Italian National Research Council. While the latter provides expertise in material science and serves as the application expert, Universitat Politècnica de València contributes its deep knowledge in high-performance computing (HPC) and the development of numerical solvers, particularly the SLEPc library.

At the heart of the project is the integration of Yambo, a well-established first-principles material science code, and SLEPc, a parallel library for solving large-scale eigenvalue problems. Yambo specializes in the study of optical properties of materials and quasiparticle energies using advanced techniques such as Many-Body Perturbation Theory (MBPT), where the Bethe-Salpeter Equation (BSE) is used to model electron-hole interactions. As Román explains, “Yambo is different from other density functional theory (DFT) codes because it goes beyond DFT, which typically focuses on ground-state properties. Instead, Yambo is used to calculate the optical properties and excited states of materials, which requires solving a different, and much more computationally demanding, set of equations.”

The computational complexity in Yambo arises from the Bethe-Salpeter equation (BSE), which generates large matrices that are both non-Hermitian and dense. These matrices must be solved to compute the eigenvalues and eigenvectors that describe the excitonic properties of materials. In contrast to the Hermitian matrices that arise in traditional DFT approaches, the BSE matrix are not hermitian and require much larger memory usage. Accordingly specialized numerical methods are needed, as standard Hermitian solvers cannot efficiently handle the problem, which are able, at the same time to handle large memory consumption.

SLEPc, developed by Román’s team over the past 25 years, offers a robust framework for addressing such large-scale eigenvalue problems. It supports various problem types, including standard, generalized, and nonlinear eigenvalue problems, and is designed to work efficiently on parallel systems, including GPUs and MPI clusters. However, in the context of the BSE, additional innovation was required to handle the dense, non-Hermitian nature of the matrix.

### Pseudo-Hermitian Matrices and Their Challenges

A central focus of ISOLV-BSE is the development of solvers for structured pseudo-Hermitian matrices, which exhibit unique algebraic properties that can be exploited for computational efficiency. Structured pseudo-Hermitian matrices arise in the Bethe-Salpeter framework as a result of the mathematical structure of the excitonic Hamiltonian, which describes the interaction between electron-hole pairs. The matrix typically takes the form:

Where R is a complex Hermitian matrix, and C is complex symmetric. This matrix is pseudo-Hermitian because it is self-adjoint with respect to a signature matrix, which introduces certain symmetries in the eigenvalues that can be exploited. In particular, the eigenvalues of a structured pseudo-Hermitian matrix appear in pairs, such that for every eigenvalue λ, there exists a corresponding −λ. This symmetry reduces the number of eigenvalues that need to be computed, thus significantly lowering the computational cost.

Despite these useful properties, solving structured pseudo-Hermitian matrices is far from trivial. Unlike Hermitian matrices, where the eigenvalues are guaranteed to be real and the eigenvectors are orthogonal, non-Hermitian matrices can have complex eigenvalues, and their eigenvectors are not orthogonal. This poses additional challenges in terms of computation and accuracy. Fortunately, the Bethe-Salpeter matrix has real eigenvalues and their eigenvectors satisfy a certain special orthogonality condition.

Moreover, the Bethe-Salpeter matrix is dense, meaning it requires a substantial amount of memory and computational resources to handle, especially when dealing with large-scale simulations on exascale supercomputers.

### The Need for Specialized Solvers

The integration of SLEPc into Yambo had been previously attempted but faced several limitations due to the dense, non-Hermitian nature of the BSE matrix. “Yambo needed a more specialized solution because it requires not only right eigenvectors, but also left eigenvectors,” Román explains. This led to the development of the ISOLV-BSE project, which aims to create a specific solver for structured pseudo-Hermitian matrices , reducing the computational load and improving the efficiency of BSE simulations.

Traditional non-Hermitian solvers treat the problem as general, without exploiting the specific structure of the matrix. In contrast, the algorithm being developed by ISOLV-BSE takes advantage of the block structure of the Bethe-Salpeter matrix to minimize the amount of computation required. “We designed our solver to only operate on the upper two blocks of the matrix, using the symmetry of the eigenvalues to avoid unnecessary calculations,” Román elaborates. This approach drastically reduces the memory requirements as well as the computational cost.

Another critical aspect of the new algorithm is its ability to avoid the computation of left eigenvectors directly. In structured pseudo-Hermitian matrices, the left eigenvectors can be derived from the right eigenvectors due to the symmetry of the matrix, a feature that Román’s team has capitalized on to save computational resources. This results in an order-of-magnitude improvement in speed compared to standard non-Hermitian solvers, making it possible to solve problems that were previously computationally intractable.

**Algorithmic Techniques**

The team considered several algorithmic strategies to improve convergence. One of these is the development of a Lanczos-type iterative solver for structured pseudo-Hermitian matrices. The Lanczos method is particularly well suited for solving large, sparse Hermitian problems, but it required significant modifications to handle the block structure and non-Hermitian nature of the BSE matrix. This method, based on the Krylov subspace technique, constructs a sequence of vectors that approximates the eigenvectors of the matrix.

The original plan was to use shift-and-invert spectral transform, a technique that accelerates the convergence of eigenvalues near a specified point. However, during the project, they discovered that the new algorithm naturally avoids the negative eigenvalues, making the shift-and-invert technique unnecessary and very costly for a dense matrix because of the need to inverse the matrix which needs a lot of computation

Instead, the team used other approaches such as polynomial filtering and preconditioned eigensolvers, both of which aim to improve convergence without the need for a matrix factorization, which is less costly than the shift-and-invert approach is. “We are still testing these approaches. They are ways of doing the same thing, or the similar thing as the shift-and-invert technique, but cheaper,” adds Román.

### Efficient Simulation of Optical Properties in CrI₃ Monolayers

CrI₃ monolayers were utilized as a benchmark to validate the performance of the newly developed solver for structured pseudo-Hermitian matrices CrI₃, a two-dimensional magnetic material, exhibits unique optical properties due to its layered structure of chromium and iodine atoms. These properties were modeled using the Bethe-Salpeter equation, generating dense, non-Hermitian matrices.

By integrating the specialized structured pseudo-Hermitian matrix solver into the SLEPc library, the ISOLV-BSE project achieved significant reductions in computational cost while maintaining high accuracy. This approach allowed for more efficient calculation of optical absorption spectra and excitonic interactions in CrI₃, proving the solver’s potential for scaling to exascale computing environments.

**Addressing the Challenges of Exascale Computing**

Scaling the algorithm to run efficiently on exascale systems was a primary challenge for ISOLV-BSE due to the goals set by the Inno4scale project. Dense matrices, like those produced by BSE, require a substantial amount of memory and processing power, making them difficult to scale efficiently across many nodes. However, Román’s team has made significant progress in optimizing the memory usage of the algorithm. “By operating only on the upper two blocks of the matrix and avoiding the storage of left eigenvectors, we’ve managed to reduce the memory requirements and ensure that the algorithm scales well as the number of processors increases,” Román explains.

The project’s focus on parallel scalability is critical, especially as simulations become larger and more complex. The integration of GPU support in SLEPc is another key feature that allows the solver to take advantage of modern HPC architectures. “We’re now able to use GPUs in addition to MPI, which opens up new possibilities for even larger-scale simulations,” Román mentions.

**Long-Term Impact on Science and Computing**

The impact of ISOLV-BSE extends beyond its immediate applications in material science. The innovations developed in the project could influence a wide range of fields that rely on solving large-scale eigenvalue problems. Moreover, the project’s contribution to numerical linear algebra is significant, as it introduces a new class of solvers that can be adapted to other types of structured eigenproblems.

In terms of HPC, the project underscores the importance of investing in algorithmic improvements. “A better algorithm can provide much larger performance gains than simply investing in faster hardware,” Román emphasized. “With ISOLV-BSE, we’ve shown that by optimizing the algorithm, we can achieve a tenfold improvement in performance without needing to upgrade the hardware.”

The SLEPc library, which already has a global user base, will include the new solver developed through ISOLV-BSE, making it available to a wide range of scientific communities. “This is one of the most important results of the project. By releasing the solver in SLEPc, we’re enabling researchers from many fields to benefit from our work,” concludes Román.