# A FEAST Algorithm for the Linear Response Eigenvalue Problem
## 1 Introduction  
In computational quantum chemistry and physics, the random phase approximation (RPA) or the Bethe-Salpeter (BS) equation describe the excitation states and absorption spectra for molecules or the surfaces of solids [1; 2]. One important question in the RPA or BS equation is how to compute a few eigenpairs associated with several of the smallest positive eigenvalues of the following eigenvalue problem:  
\[\mathcal{H}w= \begin{bmatrix}A&B\\ -B&-A\end{bmatrix}\begin{bmatrix}u\\ v\end{bmatrix} =\lambda \begin{bmatrix}u\\ v\end{bmatrix} =\lambda w, \tag{1}\]  
where \(A,B\in\mathbb{R}^{N\times N}\) are both symmetric matrices and \( \begin{bmatrix}A&B\\ B&A\end{bmatrix} \) is positive definite [3; 4; 5]. Through a similarity transformation, the eigenvalue problem (1) can be equivalently transformed into:  
\[Hz= \begin{bmatrix}0&K\\ M&0\end{bmatrix}\begin{bmatrix}y\\ x\end{bmatrix} =\lambda \begin{bmatrix}y\\ x\end{bmatrix} =\lambda z, \tag{2}\]  
where \(K=A-B\) and \(M=A+B\). The eigenvalue problem (2) was still referred to as the linear response eigenvalue problem (LREP) [3; 6] and will be so in this paper, as well. The condition imposed upon \(A\) and \(B\) in (1) implies that both \(K\) and \(M\) are \(N\times N\) real symmetric positive definite matrices. However, there are cases where one of them may be indefinite [7]. Therefore, to be consistence with [8; 9; 10], throughout the rest of the paper, we relax the condition on \(K\) and \(M\) to that they are symmetric and one of them is positive definite, while the other may be indefinite. There has been immense recent interest in developing new theories, efficient numerical algorithms of LREP, and the associated excitation response calculations of molecules for materials' design in energy science [11; 12; 13; 14; 15]. From (2), we have \(Kx=\lambda y\) and \(My=\lambda x\), and they together lead to:  
\[KMy =\lambda^{2}y,\] \[MKx =\lambda^{2}x.\]  
Since \(KM=KM^{1/2}M^{1/2}\) has the same eigenvalues as \(M^{1/2}KM^{1/2}\), which is symmetric, all eigenvalues of \(KM\) are real. Denote these eigenvalues by \(\omega_{i}\) (\(1\leq i\leq N\)) in ascending order, i.e.,  
\[\omega_{1}\leq\omega_{2}\leq\cdots\leq\omega_{N}. \tag{3}\]  
The eigenvalues of \(MK\) are \(\omega_{i}\) (\(1\leq i\leq N\)), as well. Let \(\iota=\sqrt{-1}\), the imaginary unit, and:  
\[\lambda_{i}= \begin{cases}\sqrt{\omega_{i}},&\text{if }\omega_{i}\geq 0,\\ \iota\sqrt{-\omega_{i}},&\text{if }\omega_{i}<0.\end{cases}  \tag{4}\]  
The eigenvalues of \(H\) are:  
\[\pm\lambda_{i}\quad\text{for }i=1,2,\ldots,N.\]  
This practice of enumerating the eigenvalues of \(H\) will be used later for the much smaller projection of \(H\), as well. Since the dimension \(N\) is usually very large for systems of practical interest, solving LREP (2) is a very challenging problem. The iterative methods are recent efforts for computing the partial spectrum of (2), such as the locally-optimal block preconditioned 4D conjugate gradient method (LOBP4 DCG) [6] and its space-extended variation [16], the block Chebyshev-Davidson method [10], as well as the generalized Lanczos-type methods [8; 9]. These algorithms are all based on the so-called pair of deflating subspaces, which is a generalization of the concept of the invariant subspace in the standard eigenvalue problems. Thus, these methods can be regarded as extensions of the associated classical methods for the standard eigenvalue problem. Recently, in [17], a density-matrix-based algorithm, called FEAST, was proposed by [PERSON] for the calculation of a segment of eigenvalues and their associated eigenvectors of a real symmetric or Hermitian matrix. The algorithm promises high potential for parallelism by dividing the spectrum of a matrix into an arbitrary number of nonintersecting intervals. The eigenvalues in each interval and the associated eigenvectors can be solved independently of those in the other intervals.