{'fraction_non_alphanumeric': 0.031203376102988192, 'fraction_numerical': 0.010060104863804936, 'mean_word_length': 4.769798327594688, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Structure Preserving Parallel Algorithms for Solving the Bethe-Salpeter Eigenvalue Problem

The Bethe-Salpeter eigenvalue problem is a dense structured eigenvalue problem arising from discretized Bethe-Salpeter equation in the context of computing exciton energies and states. A computational challenge is that at least half of the eigenvalues and the associated eigenvectors are desired in practice. We establish the equivalence between Bethe-Salpeter eigenvalue problems and real Hamiltonian eigenvalue problems. Based on theoretical analysis, structure preserving algorithms for a class of Bethe-Salpeter eigenvalue problems are proposed. We also show that for this class of problems all eigenvalues obtained from the Tamm-Dancoff approximation are overestimated. In order to solve large scale problems of practical interest, we discuss parallel implementations of our algorithms targeting distributed memory systems. Several numerical examples are presented to demonstrate the efficiency and accuracy of our algorithms.


Introduction
The absorption of a photon by a molecular system or solid can promote an electron in an occupied single-particle state (or orbital) to an unoccupied state while keeping the charge neutrality. In the physics community, this process is often described as the simultaneous creation of a negatively charged quasielectron (or simply electron) and a positively charged quasihole (or hole) in the material that was originally in the lowest energy electronic configuration (the ground state). Upon absorbing a photon, the entire molecular or extended system is in an excited state that contains a correlated electron-hole pair, which is referred to as an exciton. The amount of energy required to trigger this excitation gives an important characterization of the material. In many-body physics, a two-particle collective excitation is often described by a two-particle Green's function, with the excitation energy level being a pole of this function. It has been shown that the two-particle Green's function satisfies an equation often known as the Bethe-Salpeter equation (BSE) [22].
The poles of the two-particle Green's function can be obtained by computing the eigenvalues of a Hamiltonian operator H associated with this Green's function. It can be shown that, with an appropriate discretization scheme, the finite dimensional representation of the Bethe-Salpeter Hamiltonian has the following block structure where A, B ∈ C n×n , with Here we use A H to denote the conjugate transpose of A and B T to denote the transpose of B. We will refer to an eigenvalue problem of the form (1) with the additional symmetry given by Equation (2) as a Bethe-Salpeter eigenvalue problem.
In principle, we are interested in all possible excitation energies, although some excitations are more likely to occur than others. Such likelihood can often be measured in term of what is known as the spectral density or density of states of H, which is defined to be the number of eigenvalues per unit energy interval. To obtain a highly accurate representation of the spectral density, the computation of all eigenvalues is required.
In general, both A and B are dense. The dimension of these matrices is proportional to the product of the number of occupied and unoccupied states, both of which are proportional to the number of electrons n e in the system. Hence it can become quite large for large systems that contain many atoms (and electrons).
We are interested in efficient and reliable parallel algorithms for computing all eigenvalues and the corresponding eigenvectors of (1). Because H is a non-Hermitian matrix, we need to compute both the left and the right eigenvectors.
Although it is possible to treat H as a general non-Hermitian matrix and use existing parallel algorithms [9,12] implemented in ScaLAPACK to solve such a problem, this approach does not take advantage of the special structure of the Bethe-Salpeter Hamiltonian and is thus not efficient. Nor does this approach preserve some desirable properties of the eigenvalues and eigenvectors. Moreover, the current release of ScaLAPACK only has a subroutine that performs a Schur decomposition of a complex non-Hermitian matrix [9].
In the following section, we show that H belongs to a class of matrices known as complex Hamiltonian matrices whose eigenvalues satisfy a special symmetry property. Although several algorithms have been developed for computing the eigenvalues and eigenvectors of this class of matrices, efficient parallel implementations of these algorithms are not available, and are not easy to develop.
In this paper, we develop a special algorithm that can leverage existing parallel computational kernels available in the ScaLAPACK software package. Our algorithm is based on the observation that computing the eigenpairs of (1) is equivalent to computing the eigenpairs of a real Hamiltonian matrix of the form where Re(A) and Im(A) denote the real and imaginary parts of A, respectively. Furthermore, when A and B satisfy the property which often holds in practice [28], it can be shown that all eigenvalues of iH r and thus of H are real, and they come in positive and negative pairs [13,14]. We present an efficient parallel algorithm for computing the positive eigenvalues and the corresponding eigenvectors. The algorithm makes use of existing kernel in ScaLAPACK as well as a new kernel we developed for computing eigenpairs of a skew symmetric tridiagonal matrix. A simple transformation can be used to obtain the eigenvectors associated with the negative eigenvalues.
When H is real, which is the case for systems with real-space inversion symmetry, the Bethe-Salpeter eigenvalue problem can be transformed into a product eigenvalue problem. We propose an efficient and accurate parallel algorithm for solving the product eigenvalue problem.
When facing the challenge of computing the eigenpairs of the non-Hermitian matrix (1), many researchers in the physics community chose to drop the off-diagonal blocks, B and −B, and compute eigenpairs of the Hermitian matrix A only. This approach is often known as the Tamm-Dancoff approximation (TDA) [7,21,23]. We show that when the condition (4) holds, each eigenvalue of A is an upper bound of the corresponding positive eigenvalue of H when all eigenvalues of A and H are sorted. Our numerical experiment shows that TDA can introduce a non-negligible shift of the spectral density of H, which is consistent with what has been reported in the physics literature.
The rest of the paper is organized as follows. In Section 2, we discuss some basic properties of the Bethe-Salpeter eigenvalue problem. Then in Section 3, we develop structure preserving algorithms built on these properties as well as the additional assumption (4). Finally, we demonstrate the efficiency and accuracy of our proposed algorithms in Section 4 by several examples from physical models.
with X 12 and X 21 complex symmetric. It can also be verified that this Hamiltonian structure is preserved under what is called a symplectic similarity transformation where S ∈ C 2n×2n is a complex symplectic matrix that satisfies S T J n S = J n [17]. These properties are key to the development of several numerical algorithms for computing the eigenvalues of dense complex Hamiltonian matrices. Examples of these algorithms include the Hamiltonian-Jacobi algorithms [6,18] and SR-like algorithms [10].
The Bethe-Salpeter Hamiltonian matrix defined in (1) is clearly a complex Hamiltonian matrix because Although the algorithms for complex Hamitonian matrices can be used to solve the Bethe-Salpeter eigenvalue problems, they do not take advantage of the additional symmetry relationship between the (1, 2) and (2, 1) blocks in H. This symmetry leads to the symmetry of Λ(H) as stated in the following theorem. Unfortunately, complex symplectic transformations in general do not preserve the structure of Bethe-Salpeter Hamiltonian matrices. To seek a class of fully structure preserving similarity transformations, we establish the equivalence between Bethe-Salpeter eigenvalue problems and real Hamiltonian eigenvalue problems in the following theorem. This is the main theoretical result of this paper.
Theorem 2. A Bethe-Salpeter eigenvalue problem can be recuded to a real Hamiltonian eigenvalue problem, and vice versa. Then is a real symmetric matrix. Therefore, any Bethe-Salpeter eigenvalue problem can be converted to a real Hamiltonian eigenvalue problem. On the other hand, let H r be a real 2n × 2n Hamiltonian matrix of the form It can be verified that A H = A, B T = B, and Therefore, we convert a real Hamiltonian eigenvalue problem to a Bethe-Salpeter eigenvalue problem. This completes the proof.
Theorem 2 fully characterizes the spectral properties of general Bethe-Salpeter eigenvalue problems. Theorem 1 becomes a direct consequence of Theorem 2. As a result, several existing algorithms (see [3]) for solving Hamiltonian eigenvalue problems can be applied to the matrix H r defined in (3).
A primary interest in this paper is drawn to the case when the property (4) holds. In this case, the Bethe-Salpeter eigenvalue problem is essentially a symmetric eigenvalue problem because is the product of two Hermitian matrices, and in addition, one of them is positive definite [13]. The property of Λ(H) is summarized in the following theorem.
Theorem 3. Let H be of the form (1) satisfying (2) and (4). Then H is diagonalizable. In addition, the eigenvalues of H are real and appear in pairs {λ, −λ}.
The proof of the above theorem is relatively straightforward. We now consider a relatively simple case in which H is a real matrix. In this case (1) simplifies to where both A and B are n × n real symmetric matrices. By definition H is a real Hamiltonian matrix. By performing a symplectic orthogonal similarity transformation, we obtain a block cyclic form 1 √ 2 This suggests that the eigenvalues of H are the square roots of the eigenvalues of (A + B)(A − B).
Since both A and B are real, (4) simplifies to that or equivalently, both A + B and A − B are positive definite. Under this condition, the real Bethe-Salpeter eigenvalue problem is also known as a linear response eigenvalue problem which recently has attracted a lot of attention (see, for example, [1,2]). In contrast to many recent developments [2,20] in linear response eigenvalue problems that focus on large sparse eigensolvers, we develop dense eigensolvers for this eigenvalue problem in the next section.

Tamm-Dancoff approximation
When the off-diagonal blocks of the Bethe-Salpeter Hamiltonian are set to zero, a practice known in the physics community as the Tamm-Dancoff approximation [7,21,23], the Bethe-Salpeter eigenvalue problem reduces to a Hermitian eigenvalue problem. One can use efficient algorithms available in ScaLAPACK to compute eigenpairs of A. In many cases, the results are found to be sufficiently close to the eigenvalues of the full Bethe-Salpeter Hamitonian and explain experiment. However, in general, this simplification can lead to noticable difference in the computed spectrum.
In this subsection, we show that Tamm-Dancoff approximation consistently overestimate the positive eigenvalues when the property (4) holds. More precisely, for j = 1, 2, . . ., n, where λ j (·) denote the jth largest eigenvalue. That is, every positive eigenvalue obtained by TDA is greater than or equal to the corresponding one of H. This theoretical result is consistent with several computational experiments made in [13,19].
Our proof makes use of the following lemma.
With the help of this arithmetic-geometric inequality on eigenvalues, we prove the claim we have made in the following theorem.
indicating that (4) is equivalent to A + B 0. Notice that Therefore, A − B is also positive definite. In the proof of Theorem 2, we have shown that H is unitarily similar to −iJ n ( A + B). Since the eigenvalues of H are real and appear in pairs Applying Lemma 1 yields , that is, the eigenvalues of A are the same as those of A with doubled multiplicity, we arrive at The theorem is thus proved.

Algorithms and implementations
In this section, we present structure preserving algorithms for solving the Bethe-Salpeter eigenvalue problem. As we have shown in Theorem 2, Bethe-Salpeter eigenvalue problems are equivalent to real Hamiltonian eigenvalue problems. Thus any Hamiltonian eigensolver (see, for example, [3]) can always be used to solve this type of eigenvalue problem. However, when (4) holds, a more efficient algorithm can be used to solve the Bethe-Salpeter eigenvalue problem. Throughout this section the condition (4) is assumed to be satisfied. In most cases, only the positive eigenvalues and the corresponding eigenvectors are required for studying the properties of materials. We will demonstrate how to compute them efficiently. We will also show that the negative eigenvalues and the associated eigenvectors, if needed, can be constructed by symmetry with negligible computational cost.

Complex Bethe-Salpeter eigenvalue problems
From (7), we can write A straightforward approach is to feed this problem to a generalized Hermitian-definite eigensolver (for example, ZHEGV in LAPACK). However, there are several drawbacks to this simple approach. For example, the computed spectrum is not guaranteed to be symmetric. Recovering the left eigenvectors also requires extra manipulations after calling the eigensolver. It is possible to first compute the Cholesky factorization and then apply a standard Hermitian eigensolver to the transformed problem However, the structure of H is still destroyed in the sense that the backward error in the computed Cholesky factorization cannot be interpreted as a structured backward error of H.
To develop a structure preserving approach, we make use of the observation we made in Section 2. We have observed that Q H HQ = −iJ n M , where Notice that both J n and M are real matrices. Thus by working with these matrices, we can avoid the use of complex arithmetic. The basic steps of our algorithm are outlined in Algorithm 1, in which we only compute the positive eigenvalues and the corresponding eigenvectors. After the matrix M is constructed, we perform a Cholesky factorization, M = LL T to transform the non-Hermitian matrix, −iJ n LL T , into the Hermitian matrix, −iL T J n L. This step can be accomplished using real arithmetic only. Notice that the backward error in the computed Cholesky factorization of M can be converted to a structured backward error of H because all entries in the lower (upper) triangular part of M are independent. Then we observe that L T J n L is a real skew-symmetric matrix whose spectrum is symmetric with respect to the origin. Applying a skew-symmetric eigensolver (for example, the one in [24,25]) preserves the structure of Λ(H) and avoids the explicit use of complex arithmetic. We will discuss in detail about the skew-symmetric eigensolver in Section 3.3. An alternative approach is to compute an SVD-like decomposition of L and avoids the explicit formulation of L T J n L. We refer to [26,27] for such an approach. Once the eigenvalues and eigenvectors of L T J n L have been computed, the eigenvectors of H can be recovered by simply accumulating all similarity transformations. Complex arithmetic can also be avoided in this step by carefully manipulating the real and imaginary parts of the eigenvectors. The left and right eigenvectors can be easily converted to each other through These eigenvectors are in general normalized to satisfy that X H 1 X 1 − X H 2 X 2 = I n . Finally, we remark that the eigenvectors correspond to the negative eigenvalues of H, if needed, can also be constructed with negligible computational cost from the outputs of Algorithm 1. This is based on the fact that

Real Bethe-Salpeter eigenvalue problems
As we have seen in Section 2, the real Bethe-Salpeter eigenvalue problem can be reduced to a product eigenvalue problem. Directly feeding A + B and A − B to a symmetric-definite eigensolver squares the eigenvalues and can potentially spoil the accuracy of the eigenvalues and eigenvectors. When the accuracy requirement is high, an alternative approach is needed.
Notice that the spectral decomposition of (L T 2 L 1 )(L T 2 L 1 ) T can be obtained from the singular value decomposition of L T 2 L 1 . Theorem 5 summarizes this observation. We omit the proof since it merely consists of simple algebraic manipulations.
Then the spectral decomposition of H is given by Based on this theorem, we propose Algorithm 2 as a real Bethe-Salpeter eigensolver. We remark that such an eigensolver can also be used as a dense kernel within structure preserving projection methods for linear response eigenvalue problems [2,15]. In this case negative eigenvalues and the corresponding eigenvectors may be needed sometimes. Just like the complex case, these eigenvectors can be easily constructed from the outputs of Algorithm 2 according to Theorem 5.

Parallel implementations
To solve large scale Bethe-Salpeter eigenvalue problems arising from quantum physics, parallelization of Algorithms 1 and 2 on distributed memory systems is required in practical computations. All O(n 3 ) operations in Algorithm 2 consist of basic linear algebra operations, and can be accomplished by calling linear algebra libraries such as BLAS/LAPACK. There also exist ScaLAPACK subroutines for these operations, which allow us to parallelize the algorithm directly. So we will mainly discuss implementation issues for Algorithm 1.
The main obstacle to efficiently implementing Algorithm 1 is the lack of a skew-symmetric eigensolver in LAPACK/ScaLAPACK. The algorithm described in [24,25] is based on level 1 BLAS Algorithm 2 Algorithm for the real Bethe-Salpeter eigenvalue problem Input: A = A T , B = B T ∈ R n×n such that A + B 0 and A − B 0. Output: X + , Y + ∈ R 2n×n and Λ + = diag {λ 1 , . . . , λ n } satisfying and λ i > 0 for i = 1, . . ., n.

1: Compute the Cholesky factorizations
operations, hence is not efficient on modern architectures with memory hierarchy. Therefore, we have to implement step 4 in Algorithm 1 by ourselves. To make use of ScaLAPACK as much as possible, we propose the following strategy shown in Algorithm 3. Several remarks on the implementation are in order: 1. Tridiagonalizing a skew-symmetric matrix can be accomplished by applying a sequence of Householder reflections. This is in fact slightly simpler compared to the procedure of tridiagonalizing a symmetric matrix. We implemented a modified version of ScaLAPACK's PDSYTRD to achieve this goal. Several BLAS/PBLAS-like subroutines for skew-symmetric matrix operations are also implemented.

Suppose that
is a real symmetric tridiagonal matrix, whose spectral decomposition can be easily computed by calling ScaLAPACK. This technique is essentially the same as the one in [25]. We then use the bisection method (PDSTEBZ/PDSTEIN) to compute the positive eigenvalues and the corresponding eigenvectors. The eigenvectors are reorthogonalized when the accuracy requirement is high.

3.
Step 4 in Algorithm 3 is not performed by PBLAS. Because Q and D are both known, the application of these unitary transformations can be accomplished with O(n 2 ) operations.
4. In the last step of Algorithm 3, we separate the computation for real and imaginary parts: Re(X + ) = Re(X)V + , Im(X + ) = Im(X)V + . This is based on the fact that V + is real, and one PZGEMM call is about twice as expensive as two PDGEMM calls.
Finally, we remark that our parallel algorithms/implementations are just proof-of-concept ones. There is certainly room of improvements (for example, two-stage tridiagonalization techniques [5], alternative tridiagonal eigensolvers [8], fine tuning in BLAS/PBLAS-like subroutines, etc.). But investigation in this direction exceeds the scope of this paper.

Numerical experiments
In this section, we present numerical examples for three matrices obtained from discretized Bethe-Salpeter equations. The numerical experiments are performed on the Cray XE6 machine, Hopper, at National Energy Research Scientific Computing Center (NERSC). 3 Each Hopper node consists of two twelve-core AMD 'MagnyCours' 2.1-GHz processors, and has 32 GB DDR3 1333-MHz memory. Each core has its own L1 and L2 caches, with 64 KB and 512 KB, respectively. Hopper's compute nodes are connected via a 3D torus Cray Gemini network with a maximum bandwidth of 8.3 gigabytes per second. The internode latency ranges from 1.27 microseconds to 1.38 microseconds. The latency between cores is less than 1 microsecond.
The examples we use here correspond to discretized Bethe-Salpeter Hamiltonians associated with naphthalene, gallium arsenide (GaAs), and boron nitride (BN), respectively. The dimensions of these matrices are 64×64, 256×256, and 4608×4608, respectively. We implemented Algorithm 1 in Fortran 90, using the message passing interface (MPI) for parallelization. We used ScaLAPACK to perform some basic parallel matrix computations. No multithreading features such as OpenMP or threaded linear algebra libraries are used. To make fair comparisons with other algorithms, all eigenvalues and eigenvectors are calculated in our implementation.
We first compare our implementation of Algorithm 1 with LAPACK's non-Hermitian eigensolver ZGEEV. Test runs are performed using a single core so that both solvers are sequential. From Table 1 we see that both solvers produce solutions with small residuals, and Algorithm 1 is in general more accurate. We also compute the eigenvalues of A using LAPACK's Hermitian eigensolver ZHEEV. In Figure 1, we plot the spectral density of the computed eigenvalues with and without TDA. The figure illustrates the difference between Λ(A) and Λ(H) for the naphthalene example. We observe up to 12% relative differences in eigenvalues in this example when TDA is used. We also see from the figure that Λ(A) is always to the right of Λ(H). For the other two examples, the error introduced by TDA are up to 1.2% and 0.93%, respectively. Computational results confirm that all eigenvalues obtained by TDA are larger than the true eigenvalues as predicted by Theorem 4. We remark that ZGEEV produces quite accurate eigenvalues in these examples, despite the fact that all computed eigenvalues are nonreal. Table 2 contains the execution time of different approaches. Algorithm 1 is about five times faster than the non-Hermitian solver ZGEEV. An interesting observation is that TDA is not much    faster than our full-BSE solver, especially when n gets large. This is not a surprising result. In fact, the major cost of Algorithm 1 is diagonalizing a real 2n × 2n skew-symmetric matrix, which is comparable to the expense of diagonalizing a complex n × n Hermitian matrix. Finally, we perform a simple scalability test on our parallel implementation of Algorithm 1. We use the BN example since it is of moderate size. Test runs are performed on 1 × 1, 2 × 2, . . ., 9 × 9, 10 × 10 processor grids, with block factor n b = 64 for the 2D block cyclic data layout. In Figure 2, we illustrate the overall execution time, as well as the performance profile for each component of the algorithm. We can see from the figure that all components scale similarly in this example. This indicates that Algorithm 1 scales reasonably well as the overall parallel scalability is close to that of the PDGEMM component.