# A FEAST Algorithm for the Linear Response Eigenvalue Problem
## 3 The FEAST Algorithm for LREP
### The Main Algorithm
,\quad\text{where }G=R^{-\rm{T}}U^{\rm{T}}KUR^{-1}.\]  
Let the eigenvalues of \(G\) be \(\rho_{j}^{2}\) in ascending order and the associated eigenvectors be \(q_{j}\), i.e., \(Gq_{j}=\rho_{j}^{2}q_{j}\) for \(1\leq j\leq s\). Then,  
\[H_{\rm{sa}}\hat{z}_{j}=\rho_{j}\hat{z}_{j},\quad\hat{z}_{j}= \begin{bmatrix} \rho_{j}q_{j}\\ q_{j}\end{bmatrix} .\]  
Approximating eigenpairs of \(H\) is then taken to be \((\rho_{j},\hat{z}_{j})\) where:  
\[\hat{z}_{j}= \begin{bmatrix}\rho_{j}\tilde{\psi}_{j}\\ \tilde{\psi}_{j}\end{bmatrix} ,\quad\tilde{\psi}_{j}=UR^{-1}q_{j}\quad\text{ and}\quad\tilde{\psi}_{j}=VR^{-1}q_{j}. \tag{9}\]  
We summarize the pseudo-code of the FEAST algorithm for LREP in Algorithm 1. A few remarks regarding Algorithm 1 are in order:  
**Remark 1**. :  
_(a) In practice, we do not need to save_ \(V_{\omega}\)_,_ \(U_{\omega}\) _and_ \(G_{\omega}\) _for every_ \(i\) _in Algorithm_ 1_. In fact, the subscript_ \"_\"_\" is just convenience for the convergence analysis in the next subsection._  
_(b) Since the trapezoidal rule yields much faster decay outside_ \(\mathcal{C}\) _than the Gauss-Legendre quadrature_ _[_23_]__, our numerical examples use the trapezoidal rule to evaluate_ \(f(KM)Y\)_, i.e., in_ (6)_,_  
\[\mu_{i} =c+re^{i\pi t_{i}},\quad t_{i}=\frac{i-1}{q-1}, \text{for }1\leq i\leq q,\] \[\omega_{1} =\frac{\pi}{2(q-1)},\quad\omega_{q}=\frac{\pi}{2(q-1)},\quad \omega_{i}=\frac{\pi}{q-1}, \text{for }2\leq i\leq q-1.\]  
_(c) In Algorithm 1, it usually is required to know the eigenvalue counts_ \(s\) _inside the contour_ \(\mathcal{C}\) _in advance. In practice,_ \(s\) _is unknown a priori, and an estimator has to be used instead. Some available methods have been proposed in_ _[_27; 28_]_ _to estimate eigenvalue counts, based on stochastic evaluations of:_  
\[\frac{1}{2\pi t}\oint_{\mathcal{C}}\operatorname{trace}\left((\mu I-KM)^{-1} \right)d\mu.\]  
_(d) In our numerical implementation of Algorithm 1, we monitor the convergence of a computed eigenpair_ \((\rho_{j},\hat{z}_{j})\) _by its normalized relative residual norm:_  
\[r(\rho_{j})=\frac{\|H\hat{z}_{j}-\rho_{j}\hat{z}_{j}\|_{1}}{(\|H\|_{1}+|\rho_{ j}|)\|\hat{z}_{j}\|_{1}}. \tag{10}\]  
_The approximate eigenpair_ \((\rho_{j},\hat{z}_{j})\) _is considered as converged when:_\[\rho_{j}\in(\underline{\omega},\,\overline{\omega})\quad\text{and}\quad r(\rho_{j})<\varepsilon,\]  
_where \(\underline{\omega}=c-r,\,\overline{\omega}=c+r,\) and \(\varepsilon\) is a preset tolerance._  
```
0: Given an initial block \(Y_{{}_{(0)}}\in\mathbb{R}^{N\times s}\). 0: Converged approximated eigenpairs \((\rho_{j},\,z_{j})\). 1:for\(i=1,2,\ldots,\) until convergence do
2: Compute \(V_{{}_{(0)}}=\tilde{f}(KM)Y_{{}_{(0)}}\) by (6), and \(U_{{}_{(0)}}=MV_{{}_{(0)}}\). 3: Compute \(V_{{}_{(0)}}^{\mathrm{T}}U_{0}=R_{{}_{(0)}}^{\mathrm{T}}R_{{}_{(0)}}\), \(\tilde{U}_{{}_{(0)}}=U_{{}_{(0)}}R_{{}_{(0)}}^{-1}\), \(\tilde{V}_{{}_{(0)}}=V_{{}_{(0)}}R_{{}_{(0)}}^{-1}\), and \(G_{{}_{(0)}}=\tilde{U}_{{}_{(0)}}^{\mathrm{T}}K\tilde{U}_{{}_{(0)}}\). 4: Compute the spectral decomposition \(G_{{}_{(0)}}=Q_{{}_{(0)}}\Omega_{{}_{(0)}}Q_{{}_{(0)}}^{ \mathrm{T}}\) and approximate eigenpairs \((\rho_{j_{{}_{(0)}}},\tilde{z}_{j_{{}_{(0)}}})\) where \(\tilde{z}_{j_{{}_{(0)}}}=\left[\rho_{j_{{}_{(0)}}}\tilde{V}_{{}_{(0)}}^{ \mathrm{T}}\tilde{V}_{{}_{(0)}}^{\mathrm{T}}\right]^{\mathrm{T}}\) by (9) for \(j=1,\ldots,s\). 5: If convergence is not reached then go to Step 2, with \(Y_{{}_{(0)}}=\tilde{V}_{{}_{(0)}}Q_{{}_{(0)}}\). 6:endfor
```  
**Algorithm 1** The FEAST algorithm for LREP.