# A FEAST Algorithm for the Linear Response Eigenvalue Problem
## 3 The FEAST Algorithm for LREP  
### The Main Algorithm  
For LREP, suppose \(\lambda_{\ell_{i}}\) for \(i=1,\ldots,s\) are all the wanted eigenvalues whose square are inside the circle \(\mathcal{C}\) with center at \(c\) and radius \(r\) on the complex plane, as illustrated in the following Figure 1.  
Given a random matrix \(Y\in\mathbb{R}^{N\times s}\), for any filter function \(f(\cdot)\) that serves the purpose of filtering out the unwanted eigen-directions, we have, by Lemma 2,  
\[f(KM)Y=\Psi f(\Lambda^{2})\Phi^{\mathrm{T}}Y.\]  
As we known in [10], the ideal filter for computing \(\lambda_{\ell_{i}}\) for \(i=1,\ldots,s\) should map all \(s\) wanted eigenvalues to one and all unwanted ones to zero. For such a purpose, we also consider the filter function applied in [17], i.e.,  
Figure 1: Schematic representation of the wanted eigenvalues. \[f(\alpha)=\frac{1}{2\pi t}\oint_{C}(\mu-\alpha)^{-1}d\mu.\]  
By [PERSON]'s residue theorem [26] in complex analysis, we have:  
\[f(KM)Y =\Psi f(\Lambda^{2})\Phi^{\mathrm{T}}Y\] \[=\Psi\left(\frac{1}{2\pi t}\oint_{C}(\mu I-\Lambda^{2})^{-1}d\mu \right)\Phi^{\mathrm{T}}Y\] \[=\Psi_{\langle t_{1}t_{2}\rangle}\Phi_{\langle t_{2}t_{2}\rangle }^{\mathrm{T}}Y. \tag{5}\]  
As stated in [24] [Lemma 1], the matrix \(\Phi_{\langle t_{1}t_{2}\rangle}^{\mathrm{T}}Y\), with probability one, is nonsingular if the entries of \(Y\) are random numbers from a continuous distribution and they are independent and identically distributed. Therefore, naturally, the columns of \(f(KM)Y\) span the invariant subspace of \(KM\) corresponding to \(\lambda_{\ell_{i}}^{2}\) for \(i=1,\ldots,s\). However, it needs to calculate the contour integral. In practice, this is done by using numerical quadratures to compute this filter \(f(KM)\) approximately. It is noted that \(\lambda_{\ell_{i}}^{2}\) for \(i=1,\ldots,s\) are all real. Let \(\mu=c+re^{it}\) for \(0\leq t\leq 2\pi\). Then, we have:  
\[\frac{1}{2\pi t}\oint_{C}(\mu I-\Lambda^{2})^{-1}d\mu =\frac{r}{2\pi}\int_{0}^{2\pi}e^{it}\left((c+re^{it})I-\Lambda^{2 }\right)^{-1}dt\] \[=\frac{r}{2\pi}\left(\int_{0}^{\pi}e^{it}(\mu I-\Lambda^{2})^{-1 }dt+\int_{\pi}^{2\pi}e^{it}(\mu I-\Lambda^{2})^{-1}dt\right)\] \[=\frac{r}{2\pi}\int_{0}^{\pi}\left(e^{it}(\mu I-\Lambda^{2})^{-1 }+\overline{e^{it}}(\overline{\mu}I-\Lambda^{2})^{-1}\right)dt\] \[f(KM)Y=\Psi f(\Lambda^{2})\Phi^{\mathrm{T}}Y.\]0  
where \(\mathrm{Re}()\) stands for the real part. A \(q\)-point interpolatory quadrature rule, such as the trapezoidal rule or Gauss-Legendre quadrature, leads to:  
\[f(KM)Y=\Psi f(\Lambda^{2})\Phi^{\mathrm{T}}Y.\]1  
where \(t_{i}\) and \(\omega_{i}\) are the quadrature nodes and weights, respectively, and \(\mu_{i}=c+re^{it_{i}}\) for \(1\leq i\leq q\). Let \(V=\tilde{f}(KM)Y=\frac{r}{\pi}\sum_{i=1}^{q}\omega_{i}\,\mathrm{Re}\left(e^{t _{i}}(\mu_{i}I-KM)^{-1}Y\right)\). Then, by (6), \(V\) is computed by solving shifted linear systems of the form:  
\[f(KM)Y=\Psi f(\Lambda^{2})\Phi^{\mathrm{T}}Y.\]2  
Furthermore, for any interpolatory quadrature rule, we have by [18] that:  
\[f(KM)Y=\Psi f(\Lambda^{2})\Phi^{\mathrm{T}}Y.\]3  
It follows by Lemma 3 that if the columns of \(V\) span the \(y\)-component of an approximate invariant subspace of \(H\), then the columns of \(U=MV\) should span the \(x\)-component of the same approximate invariant subspace. Therefore, \(\{\mathcal{R}(U)\), \(\mathcal{R}(V)\}\) is a pair of approximated deflating subspaces of \(H\). By [3], the best approximations of \(\lambda_{\ell_{i}}\) for \(i=1,\ldots,s\) with the pair of approximate deflating subspaces \(\{\mathcal{R}(U)\), \(\mathcal{R}(V)\}\) are eigenvalues of:\[H_{\rm{sa}}= \begin{bmatrix}0&W_{1}^{-T}U^{\rm{T}}KUW_{1}^{-1}\\ W_{2}^{-T}V^{\rm{T}}MVW_{2}^{-1}&0\end{bmatrix} ,\]  
where \(W_{1}\) and \(W_{2}\) are two nonsingular factors of \(V^{\rm{T}}U=W_{1}^{\rm{T}}W_{2}\). In particular, we choose \(W_{1}=W_{2}=R\) where \(R^{\rm{T}}R\) is \(V^{\rm{T}}U^{\prime}\)s Cholesky decomposition. This leads to:  
\[H_{\rm{sa}}= \begin{bmatrix}0&G\\ I&0\end{bmatrix}