## 3 The Contour Integral Method for Solving NEP  
### A Contour Integral Method  
In this paper, a typical NEP under consideration is as follows. Given a (bounded) domain \(\Omega\subset\mathbb{C}\) and a matrix-valued function \({\bf T}:\,\Omega\to\mathbb{C}^{p\times p}\) analytic in \(\Omega\), one wants to compute the values \(\lambda\in\Omega\) (eigenvalues) and \({\bf v}\in\mathbb{C}^{p}\), \({\bf v}\
eq 0\) (eigenvectors) such that  
\[{\bf T}(\lambda){\bf v}=0 \tag{14}\]  
A class of numerical methods for NEPs has recently been developed based on contour integration [20; 21; 22; 23; 24]. These methods solve NEPs using the Smith form and [PERSON]'s theorem for the analytic matrix-valued function, respectively. One of the latter, [PERSON]'s algorithm [22], is used as a basis for introducing the method. For notational simplicity when explaining the algorithm, it is assumed that the eigenvalues in which one is interested are simple. As for the more general case, the interested reader is referred to [22]. Additionally, throughout this paper the superscripts H and \(-1\) for a matrix represent the conjugate transpose and the inverse operation, respectively. [PERSON]'s algorithm is based on [PERSON]'s theorem as follows:  
**Theorem 1**. : _([PERSON]'s theorem [22; 24]). Let \(\Phi\subset\Omega\) be a compact subset. Let \(\lambda_{k}\) for \(k=1,\,\ldots,\,m\) denote the eigenvalues of \({\bf T}(\lambda)\) in \(\Phi\). Let \({\mathbf{v}}_{k}\) and \({\mathbf{w}}_{k}\) for \(k=1,\,\ldots,\,m\) denote their left and right eigenvectors, such that_  
\[{\bf T}(\lambda_{k}){\bf v}_{k}=0,\,\,{\bf w}_{k}^{\rm H}{\bf T}(\lambda_{k})= 0,\,\,{\bf w}_{k}^{\rm H}{\bf T}^{\prime}(\lambda_{k}){\bf v}_{k}=1. \tag{15}\]  
Then, there is a neighbourhood \(\Pi\) of \(\Phi\) in \(\Omega\) and an analytic matrix-valued function \({\bf L}\): \(\Pi\to\mathbb{C}^{p\times p}\) such that the resolvent \({\bf T}^{-1}(z)\) can be written as  
\[{\bf T}^{-1}(z)=\sum_{k=1}^{m}{\bf v}_{k}{\bf w}_{k}^{\rm H}(z-\lambda_{k})^{ -1}+{\bf L}(z) \tag{16}\]  
for all \(z\in\Pi\backslash\{\lambda_{1},\,\ldots,\lambda_{m}\}\). The core idea of [PERSON]'s method is to use formula (16) for the resolvent, and derive the following expression for the corresponding contour integral by utilizing the classical residue theorem. **Lemma 1** [22]. : _Suppose that \({\bf T}(z)\) has no eigenvalues on a closed contour \(\partial C\subset\Omega\) and has finitely many eigenvalues \(\{\lambda_{1},\ldots,\lambda_{m}\}\) lying inside this contour, then, for any function \(\varphi(z):\,\Omega\to\mathbb{C}\) that is analytic in \(\Omega\),_  
\[(2\pi{\rm i})^{-1}\oint_{\partial C}\varphi(z){\bf T}^{-1}(z){\rm d}z=\sum_{k= 1}^{m}\varphi(z){\bf v}_{k}{\bf w}_{k}^{\rm H}={\bf V}\varphi(\Lambda){\bf W} ^{\rm H} \tag{17}\]  
_where the columns of \({\bf V}\), \({\bf W}\in\mathbb{C}^{p\times m}\) are the right and left eigenvectors, respectively, normalized as in (15), and \(\Lambda\in\mathbb{C}^{m\times m}\) is the diagonal matrix of the eigenvalues._  
For the problem discussed in this paper, the number \(m\) of eigenvalues \(\{\lambda_{1},\,\ldots,\lambda_{m}\}\) lying inside \(\partial C\) is usually larger than the size \(p\) of \({\bf T}(z)\). In this case, one needs to use higher-order moments in (17), that is, employ different choices of function \(\varphi(z)\), e.g., the monomials \(\varphi(z)=z^{0},z^{1},z^{2}\), etc. Now, one multiplies (17) from the right by a constant random probing matrix \(\mathbf{R}\in\mathbb{C}^{p\times p}\). For any \(k=0\), \(1\),..., \(q\) (\(q>0\) is a given as an integer), one defines  
\[\mathbf{A}_{k}=\left(2\pi\mathrm{i}\right)^{-1}\oint_{\partial\mathcal{C}}\ z^{k} \mathbf{T}^{-1}(z)\mathbf{R}\mathrm{d}z=\mathbf{V}\Lambda^{k}\mathbf{W}^{ \mathrm{H}}\mathbf{R} \tag{18}\]  
and arrange these matrices in \(qp\times qp\) block-Hankel matrices as follows:  
\[\mathbf{B}_{0}=\left[ \begin{array}{ccc}\mathbf{A}_{0}&\cdots&\mathbf{A}_{q- 1}\\ \vdots&\ddots&\vdots\\ \mathbf{A}_{q-1}&\cdots&\mathbf{A}_{2q-2}\end{array}