2.1 Eigenvalues and eigenvectors

To develop our intuition about the eigenvalue structure of NEPs, it will be instructive to first think of a scalar function $F:\unicode[STIX]{x1D6FA}\rightarrow \mathbb{C}$ . The associated NEP $F(\unicode[STIX]{x1D706})v=0$ , $v\neq 0$ , is just a scalar root-finding problem. It becomes immediately apparent that an NEP on $\unicode[STIX]{x1D6FA}=\mathbb{C}$ may have

• ∙ no solutions at all, e.g. if $F(z)=\exp (z)$ ,

• ∙ finitely many solutions, e.g. if $F(z)=z^{3}-1$ ,

• ∙ infinitely many solutions, e.g. if $F(z)=\cos (z)$ .

The fact that even an NEP of size $1\times 1$ can have more than one eigenvalue tells us that the eigenvectors belonging to distinct eigenvalues need not be linearly independent: indeed the ‘vector’ $v=[\,1\,]$ is an eigenvector for all eigenvalues of an $1\times 1$ NEP. An example of size $2\times 2$ is

which is a singular matrix exactly at the points $z\in \mathbb{C}$ where $\text{e}^{\text{i}z^{2}}=1$ . Hence the eigenvalues of $F$ are $\unicode[STIX]{x1D706}_{k}=\pm \sqrt{2\unicode[STIX]{x1D70B}k}$ ( $k=0,\pm 1,\pm 2,\ldots$ ) and $[{1\atop -1}]$ is a right and left eigenvector for all of them. It will be useful to keep the function $F$ defined by (2.1) in mind as we will often refer back to it for illustration.

In the following discussions we will assume that $F:\unicode[STIX]{x1D6FA}\rightarrow \mathbb{C}^{n\times n}$ is holomorphic in a connected open set $\unicode[STIX]{x1D6FA}\subseteq \mathbb{C}$ (the domain), denoted by $F\in H(\unicode[STIX]{x1D6FA},\mathbb{C}^{n\times n})$ . This assumption is satisfied for most NEPs encountered in practice, where the elements of $F(z)$ are usually polynomial, rational, algebraic or exponential functions of $z\in \unicode[STIX]{x1D6FA}$ , or a combination thereof. As the definition of a matrix determinant involves only finite sums and products of the matrix entries, $F\in H(\unicode[STIX]{x1D6FA},\mathbb{C}^{n\times n})$ also implies that $\det F(z)\in H(\unicode[STIX]{x1D6FA},\mathbb{C})$ . Hence the eigenvalues of $F$ are the roots of a scalar holomorphic function in the domain $\unicode[STIX]{x1D6FA}$ . The following discreteness result is a consequence of that and can be found, for example, in Mennicken and Möller (Reference Mennicken and Möller2003, Theorem 1.3.1).

The condition that the resolvent set be non-empty is equivalent to saying that $\det F(z)$ does not vanish identically on $\unicode[STIX]{x1D6FA}$ . In this case we say that $F$ is regular. While a regular $F$ can have an infinite number of eigenvalues in $\unicode[STIX]{x1D6FA}$ , the set of eigenvalues $\unicode[STIX]{x1D6EC}(F)$ does not have accumulation points in $\unicode[STIX]{x1D6FA}$ . An illustrating example is the scalar function $F(z)=\sin (1/z)$ , which is holomorphic in $\unicode[STIX]{x1D6FA}=\mathbb{C}\setminus \{0\}$ with an infinite number of roots accumulating at $0\not \in \unicode[STIX]{x1D6FA}$ .

The algebraic multiplicity of an eigenvalue $\unicode[STIX]{x1D706}$ is defined as the multiplicity of the root of $\det F(z)$ at $z=\unicode[STIX]{x1D706}$ , that is, the smallest integer $j\geq 1$ such that

The algebraic multiplicity of an isolated eigenvalue must be finite, but in contrast to linear eigenvalue problems, in the nonlinear case the algebraic multiplicity of an eigenvalue is not necessarily bounded by the problem size  $n$ . This is illustrated by the scalar example $F(z)=z^{n+1}$ , which has the eigenvalue $\unicode[STIX]{x1D706}=0$ of algebraic multiplicity $n+1$ .

By Definition 1.1, the right eigenvectors associated with an eigenvalue  $\unicode[STIX]{x1D706}$ are the nonzero vectors in the null space of $F(\unicode[STIX]{x1D706})$ , denoted by $\operatorname{null}F(\unicode[STIX]{x1D706})$ . The dimension of this null space, $\dim (\operatorname{null}F(\unicode[STIX]{x1D706}))$ , is called the geometric multiplicity of $\unicode[STIX]{x1D706}$ . An eigenvalue $\unicode[STIX]{x1D706}$ is called semisimple if its geometric multiplicity is equal to its algebraic multiplicity, and it is called simple if its algebraic multiplicity equals one.

2.2 Root functions and generalized eigenvectors

An elegant definition of generalized eigenvectors is based on root functions, a concept introduced by Trofimov (Reference Trofimov1968). To motivate the idea of root functions, it is helpful to first consider a linear eigenvalue problem $Av=\unicode[STIX]{x1D706}v$ for $A\in \mathbb{C}^{n\times n}$ and $0\neq v\in \mathbb{C}^{n}$ , and to recall from matrix analysis that a sequence $(v_{0},v_{1},\ldots ,v_{m-1})$ of nonzero vectors is called a Jordan chain for the eigenvalue $\unicode[STIX]{x1D706}$ if

see, for example, the monograph by Horn and Johnson (Reference Horn and Johnson1985, Chapter 3). The vector $v_{0}$ is a right eigenvector of $A$ , or equivalently, a nonzero right null vector of  $A-\unicode[STIX]{x1D706}I$ , and the $v_{0},v_{1},\ldots ,v_{m-1}$ are generalized eigenvectors of  $A$ . Generalized eigenvectors of a matrix can be easily seen to be linearly independent. The maximal length $m$ of a Jordan chain starting with $v_{0}$ is called the rank of  $v_{0}$ .

Interestingly, by defining the functions $v(z)=\sum _{j=0}^{m-1}(z-\unicode[STIX]{x1D706})^{j}v_{j}$ and $F(z)=A-zI$ , the $m$ equations in (2.2) can be rewritten as

where $F^{(k)}(z)={\textstyle \frac{\text{d}^{k}}{\text{d}z^{k}}}F(z)$ . In other words, $F(z)v(z)$ has a root of multiplicity at least $m$ at $z=\unicode[STIX]{x1D706}$ . While (2.2) is specific to linear eigenvalue problems, (2.3) extends straightforwardly to the nonlinear case. The main difference for NEPs is that generalized eigenvectors need not necessarily be linearly independent, and even the zero vector is admitted as a generalized eigenvector. We make this precise in the following definition, which has been adopted from Mennicken and Möller (Reference Mennicken and Möller2003, Definition 1.6.1).

It can be shown (Mennicken and Möller Reference Mennicken and Möller2003, Proposition 1.6.4) that a complete system of Jordan chains always exists, which essentially requires the verification that for a regular NEP the rank $r(v_{0})$ of any eigenvector  $v_{0}$ is finite. The numbers $m_{j}$ satisfy $m_{1}\geq m_{2}\geq \cdots \geq m_{d}$ by definition and are called the partial multiplicities of $\unicode[STIX]{x1D706}$ . It can be shown that $\sum _{j=1}^{d}m_{j}$ corresponds to the algebraic multiplicity of $\unicode[STIX]{x1D706}$ (Mennicken and Möller Reference Mennicken and Möller2003, Proposition 1.8.5). The number $m_{1}$ is called the index of $\unicode[STIX]{x1D706}$ . If $m_{1}=\cdots =m_{d}=1$ , then $\unicode[STIX]{x1D706}$ is a semisimple eigenvalue; if in addition $d=1$ , then $\unicode[STIX]{x1D706}$ is a simple eigenvalue.

So far we have only considered right eigenvectors of $F$ , but clearly left eigenvectors of $F$ can be obtained as the right eigenvectors of $F^{\ast }$ . It turns out that there is a natural way to complement a complete system of Jordan chains for $F$ at $\unicode[STIX]{x1D706}$ with a complete system of Jordan chains of $F^{\ast }$ at $\unicode[STIX]{x1D706}$ , and the latter is uniquely determined if certain normalization conditions are enforced. The following result makes this precise (Mennicken and Möller Reference Mennicken and Möller2003, Theorem 1.6.5).

While the normalization conditions in Theorem 2.5 look rather messy in general, they simplify if the eigenvalue $\unicode[STIX]{x1D706}$ is semisimple. In this case we have $m_{1}=\cdots =m_{d}=1$ so that (2.4) becomes

that is, the left and right eigenvectors can be chosen $F^{\prime }(\unicode[STIX]{x1D706})$ -biorthonormal.

2.3 Factorizations: Smith form

It is often useful to have a factorization of an NEP about one or multiple eigenvalues, giving rise to local or global factorizations, respectively. A very useful factorization is the so-called Smith form, which we state in its global version; see Leiterer (Reference Leiterer1978), Gohberg and Rodman (Reference Gohberg and Rodman1981) and Kozlov and Maz’ja (Reference Kozlov and Maz’ja1999, Chapter A.6). It uses the notion of a unimodular matrix-valued function $P\in H(\unicode[STIX]{x1D6FA},\mathbb{C}^{n\times n})$ , which means that $\det P(z)$ is a nonzero constant over $\unicode[STIX]{x1D6FA}$ .

From the Smith form it is easy to identify eigenvectors of $F$ for an eigenvalue $\unicode[STIX]{x1D706}_{i}$ . Let $Q$ be partitioned columnwise as $Q(z)=[q_{1}(z),q_{2}(z),\ldots ,q_{n}(z)].$ By the diagonal structure of $D$ , the canonical unit vectors $e_{1},e_{2},\ldots ,e_{d_{i}}$ form a basis of $\operatorname{null}D(\unicode[STIX]{x1D706}_{i})$ , and therefore also a basis of $\operatorname{null}(P(\unicode[STIX]{x1D706}_{i})F(\unicode[STIX]{x1D706}_{i})Q(\unicode[STIX]{x1D706}_{i}))$ . Hence we have found a basis of $d_{i}$ right eigenvectors of $F$ for $\unicode[STIX]{x1D706}$ , namely,

By the same argument one can identify the first $d_{i}$ rows of

with left eigenvectors of $F$ for $z=\unicode[STIX]{x1D706}_{i}$ .

Another application of the Smith form is for the characterization of the singularities of the resolvent given by

Here the diagonal entries $\unicode[STIX]{x1D6FF}_{j}(z)^{-1}$ of $D(z)^{-1}$ are scalar meromorphic functions with poles of multiplicities $m_{i,j}$ at the eigenvalues $\unicode[STIX]{x1D706}_{i}$ ,

2.4 Expansions: Keldysh’s theorem

Let us consider an eigenvalue $z=\unicode[STIX]{x1D706}_{i}$ of $F$ . Then each term $\unicode[STIX]{x1D6FF}_{j}(z)^{-1}q_{j}(z)p_{j}(z)^{\ast }$ in (2.6) can be expanded into a matrix-valued Laurent series (Ahlfors Reference Ahlfors1953, Chapter 5.1.3) that is valid in some neighbourhood ${\mathcal{U}}\subseteq \unicode[STIX]{x1D6FA}$ about  $\unicode[STIX]{x1D706}_{i}$ ,

with coefficient matrices $S_{i,j,k}\in \mathbb{C}^{n\times n}$ , $S_{i,j,m_{i,j}}\neq O$ , and the remainder function $R_{j}\in H({\mathcal{U}},\mathbb{C}^{n\times n})$ . Combining these series for every term in (2.6) into a single expansion at $z=\unicode[STIX]{x1D706}_{i}$ , we arrive at

with a remainder function $R\in H({\mathcal{U}},\mathbb{C}^{n\times n})$ . This shows that the resolvent $F(z)^{-1}$ is meromorphic in ${\mathcal{U}}$ with a pole of multiplicity $\max _{j=1,\ldots ,d_{i}}\{m_{i,j}\}=m_{i,1}$ at $\unicode[STIX]{x1D706}_{i}$ .

Let us now assume that $F$ has only finitely many eigenvalues $\unicode[STIX]{x1D706}_{1},\unicode[STIX]{x1D706}_{2},\ldots ,\unicode[STIX]{x1D706}_{s}$ in the domain $\unicode[STIX]{x1D6FA}$ . Then we can apply the expansion (2.7) recursively at all the eigenvalues as follows: starting with (2.7) about $z=\unicode[STIX]{x1D706}_{1}$ , we can expand the remainder $R$ into another Laurent series about $z=\unicode[STIX]{x1D706}_{2}$ , and so forth, giving rise to the global expansion valid on the whole of $\unicode[STIX]{x1D6FA}$ ,

with $\widetilde{R}\in H(\unicode[STIX]{x1D6FA},\mathbb{C}^{n\times n})$ .

It is possible to characterize the matrices $S_{i,j,k}$ in (2.8) in terms of generalized eigenvectors of $F$ , a result that is known as Keldysh’s theorem. More precisely, for each distinct eigenvalue $\unicode[STIX]{x1D706}_{i}\in \unicode[STIX]{x1D6FA}$ ( $i=1,2,\ldots ,s$ ), let

be pairs of complete systems of Jordan chains for $F$ and $F^{\ast }$ at $\unicode[STIX]{x1D706}_{i}$ as defined by Theorem 2.5. Then the matrices $S_{i,j,k}$ in (2.8) are given as

Keldysh’s theorem can be found in this form, for example, in Mennicken and Möller (Reference Mennicken and Möller2003, Theorem 1.6.5) and Kozlov and Maz’ja (Reference Kozlov and Maz’ja1999, Theorem A.10.2). A more compact representation can be obtained by defining Jordan blocks

and arranging the generalized eigenvectors in columns of matrices

(Note the descending order of the left generalized eigenvectors in $W_{ij}$ .) By the Jordan-form definition of a matrix function (Higham Reference Higham2008, Definition 1.2), we have

and hence using (2.9),

We can now rewrite (2.8) and state Keldysh’s theorem in the following form.

Let us briefly consider two special cases of Theorem 2.8. Assume that all eigenvalues $\unicode[STIX]{x1D706}_{i}$ of $F$ are semisimple, that is, $m_{i,1}=\cdots =m_{i,d_{i}}=1$ . Then the matrix $J$ will be diagonal with each eigenvalue $\unicode[STIX]{x1D706}_{i}$ appearing exactly $d_{i}$ times. Further, if all $\unicode[STIX]{x1D706}_{i}$ of $F$ are simple (i.e.  $d_{1}=\cdots =d_{s}=1$ ), then $J$ is an $s\times s$ diagonal matrix with distinct eigenvalues, and we have

where $v_{0}^{i}$ and $w_{0}^{i}$ are right and left eigenvectors, respectively, satisfying $w_{0}^{i\ast }F^{\prime }(\unicode[STIX]{x1D706}_{i})v_{0}^{i}=1$ .

Theorem 2.8 shows that, up to a holomorphic additive term, the behaviour of the resolvent $F^{-1}(z)$ at the eigenvalues $\unicode[STIX]{x1D706}_{i}$ is captured by the inverse of a shifted Jordan matrix. In the linear case $F(z)=zI-A$ , we can set $\widetilde{R}\equiv 0$ in (2.12) and choose $V$ as a basis of $n=\overline{m}$ generalized eigenvectors and $W^{\ast }=V^{-1}$ , resulting in the standard Jordan decomposition $A=VJV^{-1}$ . In the nonlinear case, however, $V$ and $W$ can have an arbitrary (but identical) number $\overline{m}$ of columns which need not be linearly independent.

2.5 Invariant pairs

Invariant pairs play an important role in the analysis of NEPs and the derivation of algorithms. Assume again that $F\in H(\unicode[STIX]{x1D6FA},\mathbb{C}^{n\times n})$ on a non-empty domain $\unicode[STIX]{x1D6FA}\subseteq \mathbb{C}$ . Then by the Cauchy integral formula (Ahlfors Reference Ahlfors1953, Chapter 4.2) we have the representation

where $\unicode[STIX]{x1D6E4}\subset \unicode[STIX]{x1D6FA}$ is a contour containing the point $\unicode[STIX]{x1D707}$ in its interior. This formula is the basis for the concept of invariant pairs as defined by Beyn, Effenberger and Kressner (Reference Beyn, Effenberger and Kressner2011) and further advocated in Effenberger (Reference Effenberger2013a ).

Note that $F\in H(\unicode[STIX]{x1D6FA},\mathbb{C}^{n\times n})$ can always be written in the ‘split form’

with at most $\ell =n^{2}$ coefficient matrices $C_{j}$ and scalar functions $f_{j}\in H(\unicode[STIX]{x1D6FA},\mathbb{C})$ (one for each matrix entry). In many cases, $F$ is naturally given in this form and $\ell \ll n^{2}$ . Clearly, an invariant pair $(V,M)$ for $F$ satisfies

where each

is a matrix function. Hence invariant pairs can be seen as a generalization to NEPs of invariant subspaces for a single matrix, and standard pairs for matrix polynomials (Gohberg, Lancaster and Rodman Reference Gohberg, Lancaster and Rodman2009). The next result shows that invariant pairs can easily be constructed from root functions; see for example Gohberg, Kaashoek and van Schagen (Reference Gohberg, Kaashoek and van Schagen1993).

By the block structure of the matrices $V$ and $J$ defined in Theorem 2.8, we immediately find that $(V,J)$ is an invariant pair. An invariant pair $(V,M)$ as defined in Definition 2.10 may contain redundant information as, for example, $([V,V],\operatorname{diag}(M,M))$ is an invariant pair of at least twice the necessary size. Also, any $(V,M)$ with $V=O$ is a trivial invariant pair. The following concept of minimal pairs prevents such anomalies.

For the matrix-valued function $F\in H(\unicode[STIX]{x1D6FA},\mathbb{C}^{2\times 2})$ in (2.1), the pair $(V,J)$ with $J,V$ as in (2.13)–(2.14) is a minimal invariant pair with minimality index  $6$ . Note that the minimality index of an invariant pair $(V,M)$ is generically equal to one when $m\leq n$ . The following result gives a simple lower and upper bound on the minimality index. The lower bound appears in Beyn (Reference Beyn2012, Lemma 5.1) and the upper bound in Kressner (Reference Kressner2009).

Let $(V,M)$ be a minimal invariant pair for $F$ with minimality index $p$ , and write $F$ in the form (2.15). If $M=UTU^{\ast }$ is the Schur decomposition of  $M$ , then on using $f_{j}(M)=Uf_{j}(T)U^{\ast }$ we have that (2.16) is equivalent to

that is, $(VU,T)$ is an invariant pair for $F$ and it is also minimal since ${\mathcal{V}}_{p}(VU,T)={\mathcal{V}}_{p}(V,M)U$ . If we let $v=VUe_{1}$ , then $v\neq 0$ because

since ${\mathcal{V}}_{p}(VU,T)$ is of full column rank. Now, $t_{11}=\unicode[STIX]{x1D706}$ is an eigenvalue of $M$ and $VUf_{j}(\unicode[STIX]{x1D706})e_{1}=f_{j}(\unicode[STIX]{x1D706})VUe_{1}=f_{j}(\unicode[STIX]{x1D706})v$ so that

showing that $\unicode[STIX]{x1D706}$ is also an eigenvalue of $F$ .

The following theorem summarizes the correspondence between (minimal) invariant pairs and Jordan chains for $F$ at $\unicode[STIX]{x1D706}$ . Results of this type have been derived by Gohberg et al. (Reference Gohberg, Kaashoek and van Schagen1993, Lemma 2.1), Beyn et al. (Reference Beyn, Effenberger and Kressner2011, Proposition 2.4) and Effenberger (Reference Effenberger2013a , Theorem 3.1.13).

Note that we have not made any explicit assumptions on the integers $\widehat{d}_{i}$ and $\widehat{m}_{ij}$ in Theorem 2.15. If we assume that $\widehat{d}_{i}$ equals the geometric multiplicity of each $\unicode[STIX]{x1D706}_{i}$ , then if $(V,J)$ is a minimal pair the columns of each  $V_{i}$ in Theorem 2.15 form a complete system of Jordan chains. In this case the algebraic multiplicities of the eigenvalues of $J$ coincide with the algebraic multiplicities of the corresponding eigenvalues of $F$ . This gives rise to the following definition.

Complete invariant pairs have been introduced by Kressner (Reference Kressner2009) under the name simple invariant pairs. We prefer to use the term complete as it better aligns with the terminology of a complete system of Jordan chains in Definition 2.3.

2.6 Nonlinear Rayleigh functionals

Nonlinear Rayleigh functionals are the generalization of Rayleigh quotients for matrices. Duffin (Reference Duffin1955) introduces them for symmetric overdamped quadratic eigenvalue problems. They are generalized to larger classes of Hermitian NEPs in Rogers (Reference Rogers1964), Hadeler (Reference Hadeler1967) and Voss (Reference Voss2009), and to general NEPs in Schreiber (Reference Schreiber2008) and Schwetlick and Schreiber (Reference Schwetlick and Schreiber2012).

Let $\unicode[STIX]{x1D706}$ be an eigenvalue of the matrix-valued function $F\in H(\unicode[STIX]{x1D6FA},\mathbb{C}^{n\times n})$ with corresponding right and left eigenvectors $v$ and $w$ . For some $\unicode[STIX]{x1D70C}>0$ and $\unicode[STIX]{x1D700}<\unicode[STIX]{x1D70B}/2$ , let

be neighbourhoods of $\unicode[STIX]{x1D706}$ and $v$ , respectively. Here $\angle (x,v)$ denotes the angle between the nonzero vectors $x$ and $v$ , that is,

A functional

is a two-sided nonlinear Rayleigh functional if the following conditions hold:

The condition (2.20) restricts the functional to vectors close to eigenvectors corresponding to a simple eigenvalue.

The next result, which appears in Schwetlick and Schreiber (Reference Schwetlick and Schreiber2012, Theorem 5), concerns the local uniqueness of $p(x,y)$ and bounds the distance of $p(x,y)$ to the exact eigenvalue in terms of the angles between eigenvectors and approximations to eigenvectors.

It then follows from Theorem 2.17 that $p(v,w)=\unicode[STIX]{x1D706}$ , that is, $p(v,w)$ is an eigenvalue of $F$ .

Schreiber (Reference Schreiber2008) shows that the nonlinear Rayleigh functional $p$ is stationary at the eigenvectors $(v,w)$ , that is,

where $\unicode[STIX]{x1D706}=p(v,w)$ . In other words, the first-order terms in a perturbation expansion of $p$ at the eigenvectors $(v,w)$ vanish identically.

2.7 Some perturbation results

It is often of interest, for example in sensitivity analysis and also in numerical computing, to give simple regions in which the eigenvalues of $F$ are located. Probably the most famous result for linear eigenvalue problems is Gershgorin’s theorem; see Horn and Johnson (Reference Horn and Johnson1991, Chapter 6). Bindel and Hood (Reference Bindel and Hood2013, Theorem 3.1) provide a generalization for NEPs.

Also of interest is the sensitivity of an eigenvalue $\unicode[STIX]{x1D706}$ of $F\in H(\unicode[STIX]{x1D6FA},\mathbb{C}^{n\times n})$ under small perturbations $\unicode[STIX]{x1D6E5}F\in H(\unicode[STIX]{x1D6FA},\mathbb{C}^{n\times n})$ . To address this point, we assume that $\unicode[STIX]{x1D706}$ is a nonzero simple eigenvalue with right eigenvector $v$ and left eigenvector $w$ . We consider $F$ expressed as in (2.15) with the perturbation

To measure the sensitivity of $\unicode[STIX]{x1D706}$ , we can use the normwise condition number

where the $\unicode[STIX]{x1D6FC}_{j}$ are nonnegative parameters that allow freedom in how perturbations are measured – for example, in an absolute sense ( $\unicode[STIX]{x1D6FC}_{j}=1$ ) or in a relative sense ( $\unicode[STIX]{x1D6FC}_{j}=\Vert C_{j}\Vert _{2}$ ). By setting $\unicode[STIX]{x1D6FC}_{j}=0$ we can force $\unicode[STIX]{x1D6E5}C_{j}=0$ and thus keep $C_{j}$ unperturbed.

Pseudospectra are another established tool for gaining insight into the sensitivity of the eigenvalues of a matrix to perturbations; see Trefethen and Embree (Reference Trefethen and Embree2005) and the references therein. The $\unicode[STIX]{x1D700}$ -pseudospectrum of $F\in H(\unicode[STIX]{x1D6FA},\mathbb{C}^{n\times n})$ is the set

where, for example,

If $F$ is expressed as in (2.15), then we can consider perturbations $\unicode[STIX]{x1D6E5}F$ of the form (2.21) and measure them using

where the $\unicode[STIX]{x1D6FC}_{j}$ are nonnegative parameters defined as for (2.22). Other measures of the perturbations can be used too; see for example Tisseur and Higham (Reference Tisseur and Higham2001), Higham and Tisseur (Reference Higham and Tisseur2002) and Michiels, Green, Wagenknecht and Niculescu (Reference Michiels, Green, Wagenknecht and Niculescu2006).

The following characterizations of $\unicode[STIX]{x1D6EC}_{\unicode[STIX]{x1D700}}(F)$ (see Bindel and Hood Reference Bindel and Hood2013, Proposition 4.1, and Michiels et al. Reference Michiels, Green, Wagenknecht and Niculescu2006, Theorem 1) provide an easy way to check whether a point $z\in \unicode[STIX]{x1D6FA}$ is in the $\unicode[STIX]{x1D700}$ -pseudospectrum of $F$ or not.

The $\unicode[STIX]{x1D700}$ -pseudospectrum is connected to the backward error of an approximate eigenvalue $\widetilde{\unicode[STIX]{x1D706}}$ of $F$ defined by

where

with perturbations measured using, for example, (2.24) or (2.25). Then, by comparing the definitions (2.23) and (2.26), we have that

For an approximate eigenpair $(\widetilde{\unicode[STIX]{x1D706}},\widetilde{v})$ of $F$ and for perturbations measured using (2.25), a straightforward generalization of Tisseur (Reference Tisseur2000, Theorem 1 and Lemma 3) leads to explicit expressions for $\unicode[STIX]{x1D702}_{F}(\widetilde{\unicode[STIX]{x1D706}},\widetilde{v})$ and $\unicode[STIX]{x1D702}_{F}(\widetilde{\unicode[STIX]{x1D706}})$ ,

where $f(\widetilde{\unicode[STIX]{x1D706}})=\sum _{j=1}^{\ell }\unicode[STIX]{x1D6FC}_{j}|f_{j}(\widetilde{\unicode[STIX]{x1D706}})|$ .

Figure 2.22.

Spectral portrait of the matrix-valued function $F$ defined in (2.1). The bar shows the mapping of colours to $\log _{10}$ of the resolvent norm, that is, the mapping $z\mapsto \log _{10}\Vert F(z)^{-1}\Vert _{2}$ .

The largest possible real part of points in the $\unicode[STIX]{x1D700}$ -pseudospectrum of $F$ is called the $\unicode[STIX]{x1D700}$ -pseudospectral abscissa of $F$ , that is,

This quantity is of particular interest for NEPs associated with systems of differential equations such as the DDE (1.1), since perturbations of size $\unicode[STIX]{x1D700}$ may cause the system to become unstable whenever $\unicode[STIX]{x1D6FC}_{\unicode[STIX]{x1D700}}(F)\geq 0$ . The distance to instability of a stable system associated with $F$ can then be defined as

Algorithms to compute $\unicode[STIX]{x1D6FC}_{\unicode[STIX]{x1D700}}(F)$ can be found in Michiels and Guglielmi (Reference Michiels and Guglielmi2012) and Verhees et al. (Reference Verhees, Van Beeumen, Meerbergen, Guglielmi and Michiels2014).

4.1 Newton’s method for scalar functions

For a given initial guess $\unicode[STIX]{x1D706}^{(0)}$ , Newton’s method for finding the roots of a scalar equation $f(z)=0$ is given by

This iteration has local quadratic convergence to simple roots (Wilkinson Reference Wilkinson1965, Section 7.25). In order to apply it to the solution of an NEP $F(\unicode[STIX]{x1D706})v=0$ , we only need a scalar function  $f$ whose roots are the eigenvalues of $F$ , and its derivative $f^{\prime }$ . Different methods result from different choices of $f$ .

Probably the most obvious choice is

The trace theorem (e.g. Lancaster Reference Lancaster1966, Theorem 5.1), states that if the entries of $F(z)$ are differentiable functions of $z$ , then for any $z\in \mathbb{C}$ for which $f(z)=\det F(z)\neq 0$ we have

and so the Newton iteration (4.1) can be rewritten as

We refer to (4.3) as the Newton-trace iteration (Lancaster Reference Lancaster1966, Section 5.5); see also Khazanov and Kublanovskaya (Reference Khazanov and Kublanovskaya1988). Note that we only need the diagonal entries of $F^{-1}(\unicode[STIX]{x1D706}^{(k)})F^{\prime }(\unicode[STIX]{x1D706}^{(k)})$ . Nevertheless, with a straightforward implementation each iteration requires the factorization of an $n\times n$ matrix $F(\unicode[STIX]{x1D706}^{(k)})$ and $2n$ triangular solves, making it a rather expensive method. A basic MATLAB implementation of the Newton-trace iteration (4.2) is provided in Figure 4.1, where we use a simple stopping criterion based on the relative residual $|f(\unicode[STIX]{x1D706}^{(k)})|/\Vert F(\unicode[STIX]{x1D706}^{(k)})\Vert _{F}$ .

Figure 4.1.

Basic MATLAB implementation of the Newton-trace iteration (4.3) for problem (2.1). The NEP parameters $F$ and $F^{\prime }$ are specified in lines 2–3 and the method’s parameters in line 4. Upon convergence, lambda is the eigenvalue and nbr_iter the number of iterations.

Kublanovskaya (Reference Kublanovskaya1970) proposes applying Newton’s method to

where $r_{nn}(z)$ is the bottom-right entry of the $R$ factor in a rank-revealing QR decomposition of $F(z)$ ,

with $\unicode[STIX]{x1D6F1}(z)$ being a permutation matrix which ensures that $r_{nn}(z)$ becomes zero before any other diagonal entry does; see also Kublanovskaja (Reference Kublanovskaja1969). Since $r_{nn}(z)=0$ is equivalent to $\det F(z)=\det R(z)=0$ , we know that the roots of $f(z)$ are exactly the roots of $\det F(z)=0$ .

The permutation $\unicode[STIX]{x1D6F1}(z)$ in (4.4) is not a continuous function of $z$ , but in a small neighbourhood of an approximate eigenvalue $\unicode[STIX]{x1D706}^{(k)}$ of $F$ the permutation can be kept constant. So if we let

then in a small neighbourhood of $\unicode[STIX]{x1D706}_{k}^{(k)}$ we can use

with $Q(\unicode[STIX]{x1D706}^{(k)})=Q_{k}$ and $R(\unicode[STIX]{x1D706}^{(k)})=R_{k}$ . Garrett et al. (Reference Garrett, Bai and Li2016) show that

where $e_{n}=[0,\ldots ,0,1]^{T}$ and $p$ solves $R_{11}^{(k)}p=r_{12}^{(k)}$ , and therefore

This leads to the Newton-QR iteration for a root of $r_{nn}(z)$ ,

At convergence, we can take

as approximations for the right and left eigenvectors of the converged approximate eigenvalue, respectively.

Figure 4.2.

Basic MATLAB implementation of the Newton-QR iteration for problem (2.1). The NEP parameters $F$ and $F^{\prime }$ are specified in lines 2–3 and the method’s parameters in line 4. Upon convergence, lambda is the eigenvalue and nbr_iter the number of iterations.

A basic MATLAB implementation of the Newton-QR iteration (4.6) is given in Figure 4.2. It requires the repeated computation of a rank-revealing QR factorization of $F(\unicode[STIX]{x1D706}^{(k)})$ and hence is only feasible for problems of moderate size (unless $F$ has some structure that can be exploited; see Section 4.5). However, the iteration (4.6) can be useful in the context of iterative refinement. Yang (Reference Yang1983) and Wobst (Reference Wobst1987) consider an approach similar to that of Kublanovskaya by using an LU factorization with column pivoting in place of a rank-revealing QR factorization.

To avoid working with nearly singular matrices $F(\unicode[STIX]{x1D706}^{(k)})$ , Andrew, Chu and Lancaster (Reference Andrew, Chu and Lancaster1995) propose using in place of $F(z)$ the bordered matrix

where $b,c\in \mathbb{C}^{n}$ are such that the matrix $G(\unicode[STIX]{x1D706})$ is nonsingular and well-conditioned at a simple eigenvalue $\unicode[STIX]{x1D706}$ . Now if $c^{T}v\neq 0$ and $w^{T}b\neq 0$ , where $v$ and $w$ are the right and left eigenvectors of $F$ associated with the eigenvalue  $\unicode[STIX]{x1D706}$ , then $G(\unicode[STIX]{x1D706})$ is nonsingular (Andrew, Chu and Lancaster Reference Andrew, Chu and Lancaster1993, Theorem 10.1). For a vector $x$ such that $c^{T}x=1$ , the authors introduce the linear system

For $z$ near $\unicode[STIX]{x1D706}$ , the matrix $G(z)$ is nonsingular and hence $x=x(z)$ and $f=f(z)$ are smooth functions of $z$ . By Cramer’s rule we have

and therefore $f(z)=0$ if and only if $\det F(z)=0$ . Differentiating (4.8) with respect to $z$ leads to

The BDS method of Andrew et al. (Reference Andrew, Chu and Lancaster1995) and the implicit determinant method of Spence and Poulton (Reference Spence and Poulton2005) consist of applying a Newton iteration to  $f$ in (4.9). The latter method is detailed in Algorithm 4.3, with a basic MATLAB implementation given in Figure 4.4.

At convergence, $x^{(k)}$ provides an approximation to the right eigenvector $v$ . Note that the factorization used to solve (4.8) can be re-used to solve (4.10). Andrew et al. (Reference Andrew, Chu and Lancaster1995) propose different choices for the vectors $b$ and $c$ leading to the row/column deletion method ( $b=e_{i}$ , $c=e_{j}$ , where $e_{k}$ denotes the $k$ th column of the identity matrix), the column substitution method ( $b\approx v$ , $c=e_{j}$ ), giving the name BDS (bordered, deletion, substitution) for this class of methods.

Figure 4.4.

Basic MATLAB implementation of the implicit determinant method for problem (2.1). The NEP parameters $F$ and $F^{\prime }$ are specified in lines 2–3 and the method’s parameters in line 4. Upon convergence, lambda is the eigenvalue and nbr_iter the number of iterations.

Deflation may be necessary when computing several nearby eigenvalues in order to avoid the Newton iteration converging to an already computed eigenvalue. Suppose we have already computed $m$ roots $\unicode[STIX]{x1D706}_{\ell }$ ( $\ell =1,\ldots ,m$ ) of $f$ near $\unicode[STIX]{x1D706}_{0}$ , each of multiplicity $m_{\ell }$ . As suggested by Wilkinson (Reference Wilkinson1965, Section 7.48), we can apply Newton’s method to

which leads to the iteration

This strategy has been shown to work well in practice (Garrett et al. Reference Garrett, Bai and Li2016).

When the derivative of $F$ is not available, it can be replaced by a finite difference approximation, leading to a quasi-Newton method. In particular, the secant method is obtained by using the approximation

Figure 4.6.

Basins of convergence for various Newton-type methods discussed in Example 4.5 for problem (2.1). Each initial point in the square region  $\unicode[STIX]{x1D6FA}$ is coloured according to the eigenvalue a method is converging to, with the five exact eigenvalues of the NEP in $\unicode[STIX]{x1D6FA}$ indicated by white $\star$ symbols. White areas indicate that no convergence is achieved within $20$ iterations or the method converged to an eigenvalue outside  $\unicode[STIX]{x1D6FA}$ .

4.2 Newton’s method for the vector equation

Newton’s method can be applied directly to the NEP $F(\unicode[STIX]{x1D706})v=0$ rather than to a scalar equation $f(z)=0$ whose roots are eigenvalues of $F$ . For this, a normalization condition on the eigenvector of the form $u^{\ast }v=1$ for some nonzero vector $u$ is added to $F(\unicode[STIX]{x1D706})v=0$ so as to have $n+1$ equations for the $n+1$ unknowns in $(\unicode[STIX]{x1D706},v)$ . Note that $u^{\ast }v=1$ ensures that $v$ is not the zero vector. This approach is briefly mentioned in Unger (Reference Unger1950) and discussed further in Ruhe (Reference Ruhe1973). In order to apply Newton’s method to ${\mathcal{N}}[{v\atop \unicode[STIX]{x1D706}}]=0$ , where

we need the Jacobian matrix of the $(n+1)$ -dimensional operator ${\mathcal{N}}$ , which is readily calculated as

Newton’s iteration is now given as

If we assume that the approximate eigenvector $v^{(k)}$ is normalized such that $u^{\ast }v_{}^{(k)}=1$ for some nonzero vector $u$ , then (4.14) can be rewritten as

The Newton iteration (4.15)–(4.16) is equivalent to the nonlinear inverse iteration (Unger Reference Unger1950), which is given in Algorithm 4.7 with a basic MATLAB implementation given in Figure 4.8.

The normalization of $v^{(k+1)}$ is necessary to avoid numerical overflow or underflow and any vector norm can be used. When $\unicode[STIX]{x1D706}^{(k)}$ is close to an eigenvalue, the matrix $F(\unicode[STIX]{x1D706}^{(k)})$ of the linear system to be solved is nearly singular. However, standard theory of inverse iteration (see e.g. Ipsen Reference Ipsen1997, Section 6.3, or Peters and Wilkinson Reference Peters and Wilkinson1979, Section 2) shows that the error in the computed vector $\widetilde{v}^{(k+1)}$ will be almost parallel to the exact solution, that is, the inaccuracy is concentrated in the length of $\widetilde{v}^{(k+1)}$ and not its direction. As it is derived from a Newton iteration, the nonlinear inverse iteration converges locally and quadratically for any simple eigenpair. The vector $u$ in (4.16) and in Algorithm 4.7 can be chosen in a number of ways as discussed below.

• ∙ A simple choice is to take $u$ to be the $i$ th unit vector $e_{i}$ , which corresponds to keeping the $i$ th entry of the vectors $v^{(k)}$ constant.

• ∙ To prevent the iteration from converging to previously computed eigenpairs, $u$ can be chosen orthogonal to already computed eigenvectors (Anselone and Rall Reference Anselone and Rall1968).

• ∙ Choosing $u=F(\unicode[STIX]{x1D706}^{(k)})^{\ast }w^{(k)}$ , where $w^{(k)}$ is an approximate left eigenvector, allows us to rewrite the eigenvalue update as

  $$\begin{eqnarray}\unicode[STIX]{x1D706}^{(k+1)}=\unicode[STIX]{x1D706}^{(k)}-\frac{w^{(k)\ast }F(\unicode[STIX]{x1D706}^{(k)})v^{(k)}}{w^{(k)\ast }F^{\prime }(\unicode[STIX]{x1D706}^{(k)})v^{(k)}},\end{eqnarray}$$
  which is the generalized Rayleigh quotient iteration that Lancaster (Reference Lancaster1966) derived for the polynomial eigenvalue problem.

Figure 4.8.

Basic MATLAB implementation of the nonlinear inverse iteration for problem (2.1). The NEP parameters $F$ and $F^{\prime }$ are specified in lines 2–3 and the method’s parameters in lines 4–5. Upon convergence, lambda is the eigenvalue and nbr_iter the number of iterations.

As mentioned in the third point above, the nonlinear inverse iteration can be combined with the computation of a left eigenvector, and a two-sided Rayleigh functional as defined in Section 2.6, to yield the two-sided Rayleigh functional iteration (Schreiber Reference Schreiber2008), given in Algorithm 4.9.

Schreiber (Reference Schreiber2008) shows that Algorithm 4.9 achieves local cubic convergence for a simple eigentriple $(\unicode[STIX]{x1D706},v,w)$ . Note that at each iteration, an LU factorization of $F(\unicode[STIX]{x1D706}^{(k)})$ used to solve the first linear system can be re-used to solve the second linear system involving $F^{\ast }(\unicode[STIX]{x1D706}^{(k)})$ .

Each step of the nonlinear inverse iteration, Algorithm 4.7, requires the factorization of $F(\unicode[STIX]{x1D706}^{(k)})$ to solve the linear system, which cannot be re-used at the next iteration since $\unicode[STIX]{x1D706}^{(k)}$ varies. We could replace $\unicode[STIX]{x1D706}^{(k)}$ by a fixed value  $\unicode[STIX]{x1D70E}$ , say $\unicode[STIX]{x1D70E}=\unicode[STIX]{x1D706}^{(0)}$ , but then the resulting iteration would converge to an eigenpair $(\unicode[STIX]{x1D707},v)$ of the linear problem $Av=\unicode[STIX]{x1D707}Bv$ , where $A=F(\unicode[STIX]{x1D70E})$ and $B=F^{\prime }(\unicode[STIX]{x1D706}_{\ast })$ for some $\unicode[STIX]{x1D706}_{\ast }\in \mathbb{C}$ . If $F(z)$ is linear in $z$ so that $F^{\prime }(z)$ is a constant matrix, we could easily recover an eigenpair $(\unicode[STIX]{x1D706},v)$ of $F(\unicode[STIX]{x1D706})$ from $(\unicode[STIX]{x1D707},v)$ , since $F(\unicode[STIX]{x1D706})v=0$ and $Av=\unicode[STIX]{x1D707}Bv$ are directly related eigenvalue problems. This is not the case if $F(z)$ is truly nonlinear in $z$ .

Neumaier (Reference Neumaier1985) shows that this difficulty can be avoided by considering a variant of the nonlinear inverse iteration based on the use of the residual. If $F$ is twice continuously differentiable, iteration (4.15) can be rewritten as

Ignoring the second-order term leads to an iteration

where the new approximant $\unicode[STIX]{x1D706}^{(k+1)}$ for the eigenvalue $\unicode[STIX]{x1D706}$ needs to be determined beforehand. Neumaier (Reference Neumaier1985) shows that replacing $\unicode[STIX]{x1D706}^{(k)}$ by a fixed shift $\unicode[STIX]{x1D70E}$ in (4.17) does not destroy the convergence of the iteration to the wanted eigenpair. This leads to the residual inverse iteration, the pseudocode of which is given in Algorithm 4.10.

The initial vector $v^{(0)}$ can be computed by solving

for some nonzero vector $b$ . As suggested by Wilkinson (Reference Wilkinson1965) for the standard inverse iteration, we can choose $u=e$ , where $e$ is the vector of all ones, and $b=Le$ , where $F(\unicode[STIX]{x1D70E})=LU$ is an LU factorization. Then $\widetilde{v}^{(0)}=U^{-1}e$ . Note that compared with the nonlinear inverse iteration, Algorithm 4.7, in the residual inverse iteration the eigenvector is updated by solving a linear system whose right-hand side is the NEP residual, giving its name to the iteration. Due to the fixed shift $\unicode[STIX]{x1D70E}$ , the residual inverse iteration converges only linearly, and in practice, it is advisable to update the shift every now and then, in particular when the convergence is slow.

4.3 Deflation of computed eigenvalues and eigenpairs

The methods we have described so far are directed towards computing one eigenvalue or one eigenpair only. In principle, several runs of the same iteration with different initial guesses could return several eigenpairs, but special care needs to be taken to prevent the iteration from converging to already computed eigenpairs. This is what deflation aims to achieve.

A standard deflation technique consists of applying a non-equivalence transformation to the matrix-valued function $F$ that maps the already computed eigenvalues to infinity. This can be done as follows. Suppose we have computed $\ell$ simple eigenvalues of $F$ , $\unicode[STIX]{x1D706}_{1},\ldots ,\unicode[STIX]{x1D706}_{\ell }$ , and let $x_{i},y_{i}\in \mathbb{C}^{n}$ be such that $y_{i}^{\ast }x_{i}^{}=1$ , $i=1,\ldots ,\ell$ . Consider

Then it is not difficult to see that

Indeed, since $y_{i}^{\ast }x_{i}^{}=1$ ,

so that $\unicode[STIX]{x1D706}_{i}$ is not a root of $\det \widetilde{F}(z)$ since it is a simple root of $\det F(z)$ . With the exception of the eigenvalues $\unicode[STIX]{x1D706}_{i}$ , $i=1,\ldots ,\ell$ , $\widetilde{F}$  and  $F$ have the same finite eigenvalues. Now suppose that $F$ does not have an eigenvalue at infinity, that is, $G(z)=F(1/z)$ does not have an eigenvalue at zero (see Remark 2.2), and let $\widetilde{G}(z)=\widetilde{F}(1/z)$ . Then from

we see that $\widetilde{G}$ has $\ell$ eigenvalues at zero, and therefore $\widetilde{F}$ has $\ell$ eigenvalues at infinity. It follows from (4.19) that if $\widetilde{v}$ is an eigenvector of $\widetilde{F}$ with eigenvalue  $\unicode[STIX]{x1D706}$ , then

is an eigenvector of $F$ associated with the eigenvalue $\unicode[STIX]{x1D706}$ . To compute the next eigenpair of $F$ , we apply any Newton-based method to $\widetilde{F}(z)$ and recover the eigenvector using (4.20).

The choice of the vectors $x_{i}$ and $y_{i}$ affects the basins of convergence of Newton’s method. A choice that seems to work well in practice is to set $x_{i}$ and $y_{i}$ to be (approximate) right and left eigenvectors corresponding to $\unicode[STIX]{x1D706}_{i}$ . Note that similar deflation strategies have been used in Ferng, Lin, Pierce and Wang (Reference Ferng, Lin, Pierce and Wang2001) and Huang, Lin and Mehrmann (Reference Huang, Lin and Mehrmann2016).

Another approach proposed by Effenberger (Reference Effenberger2013b ) makes use of minimal invariant pairs as defined in Section 2.5. Effenberger’s deflation strategy works as follows. Let $F\in H(\unicode[STIX]{x1D6FA},\mathbb{C}^{n\times n})$ and suppose that we have already computed a minimal invariant pair $(V,M)\in \mathbb{C}^{n\times m}\times \mathbb{C}^{m\times m}$ for $F$ with minimality index  $p$ . We want to extend $(V,M)$ into another minimal invariant pair

of one size larger. By Definition 2.10, the pair $(\widehat{V},\widehat{M})$ is invariant if and only if ${\mathcal{F}}(\widehat{V},\widehat{M})=O$ . Since

we have

with

and $\unicode[STIX]{x1D6E4}$ a contour enclosing the eigenvalues of $M$ and $\unicode[STIX]{x1D706}$ . Since the pair $(V,M)$ is invariant, ${\mathcal{F}}(V,M)=O$ , and so $(\widehat{V},\widehat{M})$ is an invariant pair if and only if

The condition that $(\widehat{V},\widehat{M})$ be minimal is more involved. Effenberger (Reference Effenberger2013b , Lemma 3.4) shows that $(\widehat{V},\widehat{M})$ is minimal with minimality index not exceeding $p+1$ if

where the $m\times n$ matrix-valued function $A(\unicode[STIX]{x1D706})$ and the $m\times m$ matrix-valued function $B(\unicode[STIX]{x1D706})$ are given by

with

It then follows from (4.23) and (4.24) that $(\unicode[STIX]{x1D706},[{x\atop b}])$ is an eigenpair of the $(n+m)\times (n+m)$ NEP

where

The matrix-valued function $\widetilde{F}(\unicode[STIX]{x1D706})$ is holomorphic since $U(\unicode[STIX]{x1D706})$ is holomorphic (Effenberger Reference Effenberger2013b , Lemma 3.2), and $A(\unicode[STIX]{x1D706})$ and $B(\unicode[STIX]{x1D706})$ are matrix polynomials. If $F(\unicode[STIX]{x1D706})$ is regular, then so is $\widetilde{F}(\unicode[STIX]{x1D706})$ . Moreover, if $(\unicode[STIX]{x1D706},[{x\atop b}])$ is an eigenpair of $\widetilde{F}(\unicode[STIX]{x1D706})$ , then $([Vx],[{M\atop 0}{b\atop \unicode[STIX]{x1D706}}])$ is a minimal invariant pair for $F(\unicode[STIX]{x1D706})$ (Effenberger Reference Effenberger2013b , Theorem 3.6). This shows that the pair $(V,M)$ is deflated from the computation when solving the NEP (4.25) in place of $F(\unicode[STIX]{x1D706})v=0$ . An eigenpair for the latter can be computed using any of the methods described in Sections 4.1–4.2. But for these, we need to evaluate $U(\unicode[STIX]{x1D706})$ in (4.22) and its derivative at a scalar $\unicode[STIX]{x1D706}$ . This can be done via numerical integration techniques similar to those described in Section 5.2. Note also that, by using $\unicode[STIX]{x1D706}I-M=(zI-M)-(z-\unicode[STIX]{x1D706})I$ , (4.22), and the definition of invariant pairs, we have

see Effenberger (Reference Effenberger2013b , Lemma 4.2). Hence, if $\unicode[STIX]{x1D706}$ is not an eigenvalue of $M$ , then $U(\unicode[STIX]{x1D706})=F(\unicode[STIX]{x1D706})V(\unicode[STIX]{x1D706}I-M)^{-1}$ . For large NEPs, Effenberger (Reference Effenberger2013b ) combines this deflation strategy with a Jacobi–Davidson method (see Section 4.5) and uses a contour integration method as described in Section 5 to solve the projected NEP.

4.4 Block Newton method

Kressner (Reference Kressner2009) presents a block analogue of the Newton iteration on vectors which computes a complete invariant pair $(V,M)\in \mathbb{C}^{n\times m}\times \mathbb{C}^{m\times m}$ for $F\in H(\unicode[STIX]{x1D6FA},\mathbb{C}^{n\times n})$ with minimality index $p$ . Recall from Definition 2.16 that an invariant pair $(V,M)$ is complete if it is minimal and the algebraic multiplicities of the eigenvalues of $M$ are the same as the algebraic multiplicities of the eigenvalues of $F$ .

Let $(V,M)$ be an invariant pair for $F$ so that

where ${\mathcal{F}}$ is as in Definition 2.10, and assume that $(V,M)$ satisfies the normalization condition

where ${\mathcal{N}}(V,M)=U^{\ast }{\mathcal{V}}_{p}(V,M)-I_{m}$ with $U\in \mathbb{C}^{pn\times m}$ a fixed matrix of full column rank, and ${\mathcal{V}}_{p}(V,M)$ is as in (2.17). To apply Newton’s method to (4.26)–(4.27), we need the Fréchet derivatives of ${\mathcal{F}}$ and ${\mathcal{N}}$ at $(V,M)$ , which are given by

where $U_{}^{\ast }=[U_{0}^{\ast },\ldots ,U_{p-1}^{\ast }]$ with $U_{j}\in \mathbb{C}^{n\times m}$ . Kressner (Reference Kressner2009, Theorem 10) shows that the invariant pair $(V,M)$ is complete if and only if the linear matrix operator

corresponding to the Jacobian of (4.26)–(4.27) at $(V,M)$ is invertible. Then, by the implicit function theorem for holomorphic functions (Krantz Reference Krantz1982), we have that the entries of a complete invariant pair $(V,M)$ vary analytically under analytic changes of $F$ . Hence, complete invariant pairs are well-posed and Newton’s method for solving (4.26)–(4.27) converges locally quadratically to a complete invariant pair. Now assuming that $U^{\ast }{\mathcal{V}}_{p}(V,M)=I_{m}$ , the Newton correction $(\unicode[STIX]{x1D6E5}V,\unicode[STIX]{x1D6E5}M)$ satisfies

The block Newton method for computing an invariant pair $(V,M)\in \mathbb{C}^{n\times m}\times \mathbb{C}^{m\times m}$ of $F$ is given by the pseudocode in Algorithm 4.13.

The initial pair $(\widetilde{V}^{(0)},\widetilde{M}^{(0)})$ can be constructed with a block variant of inverse iteration (Kressner Reference Kressner2009, Algorithm 2). The solution to the linear matrix equation (4.28) is the most expensive part of the iteration. Kressner (Reference Kressner2009) shows how to do that efficiently using ideas from Beyn and Thümmler (Reference Beyn and Thümmler2009), and provides a MATLAB function nlevp_newtonstep that implements these ideas and returns the pair $(\unicode[STIX]{x1D6E5}V,\unicode[STIX]{x1D6E5}M)$ .

4.5 Large sparse NEPs

Most of the methods we have described so far require the solution of linear systems. For large sparse matrices, there are efficient direct methods implementing Gaussian elimination with some form of pivoting that make clever use of the sparsity structure to avoid fill-in and save storage. These include HSL (2016), MUMPS (2016), PARDISO (Schenk and Gärtner Reference Schenk and Gärtner2004) and UMFPACK (Davis Reference Davis2004). In place of direct methods to solve the sparse linear systems in Algorithms 4.7–4.10, we can also use iterative methods, for example Krylov subspace methods such as the generalized minimal residual method (GMRES) of Saad and Schultz (Reference Saad and Schultz1986), the biconjugate gradient method in its stabilized form (BICGSTAB) by van der Vorst (Reference van der Vorst1992), or the quasi-minimal residual method (QMR) of Freund and Nachtigal (Reference Freund and Nachtigal1996). Solving the linear systems iteratively leads to inexact versions of nonlinear inverse iteration, residual inverse iteration, and Rayleigh functional iteration. Let us consider in particular the inexact nonlinear inverse iteration, Algorithm 4.7, with the exact solve replaced by an inexact solve. Let $\unicode[STIX]{x1D70F}^{(k)}$ be a tolerance such that the approximate solution $\widetilde{v}^{(k+1)}$ returned by the iterative solver at step $k$ satisfies

and let

denote the error at iteration $k$ . Szyld and Xue (Reference Szyld and Xue2013a , Theorem 6) show that if, at each iteration,

is ensured with a constant $c$ independent of $k$ , then the inexact nonlinear inverse iteration converges at least quadratically. In a similar way, they show that for inexact versions of the residual inverse iteration and the two-sided Rayleigh functional iteration, the same order of convergence as for the exact iterations can be achieved if an appropriate sequence of tolerances is used for the inner solves.

Garrett et al. (Reference Garrett, Bai and Li2016) propose a variation of Kublanovskaya’s Newton-QR method for either large banded NEPs with narrow bandwidths relative to the size $n$ of the NEP, or large sparse NEPs that can be reduced to such banded NEPs. They make use of a special data structure for storing the matrices to keep memory and computational costs low. They replace the rank-revealing QR factorization of $F(\unicode[STIX]{x1D706}^{(k)})$ in (4.5) by a banded QR factorization followed by a rank-revealing QR factorization of the $R$ factor. MATLAB and C++ implementations of this approach including the deflation strategy (4.11) are publicly available (Garrett and Li Reference Garrett and Li2013).

For large sparse NEPs, the iterations described in Sections 4.1–4.2 can be used as part of an iterative projection method in the following way. Suppose we are given a matrix $U\in \mathbb{C}^{n\times k}$ with $k\ll n$ having orthonormal columns that span a subspace (the search space) containing an eigenvector of interest. Then instead of solving $F(\unicode[STIX]{x1D706})v=0$ we can solve the $k\times k$ projected NEP

where $Q\in \mathbb{C}^{n\times k}$ is some matrix with orthonormal columns spanning the test space. Now let $(\unicode[STIX]{x1D717},x)$ be an eigenpair of the projected problem (4.30), selected so that $\unicode[STIX]{x1D717}$ is as close as possible to the target eigenvalue. The pair $(\unicode[STIX]{x1D717},v)$ with $v=Ux$ is the corresponding Ritz eigenpair for $F$ , and if the residual $\Vert F(\unicode[STIX]{x1D717})v\Vert$ is small enough, we can accept $(\unicode[STIX]{x1D717},v)$ as an approximate eigenpair for $F$ . If the residual $\Vert F(\unicode[STIX]{x1D717})v\Vert$ is too large, then the search space $\operatorname{span}\{U\}$ can be extended by one step of the Newton iteration with initial guess $(\unicode[STIX]{x1D717},v)$ and $v$ normalized to have unit 2-norm. The Newton step in (4.14) is rewritten as

where $\unicode[STIX]{x1D6E5}\unicode[STIX]{x1D717}$ and $\unicode[STIX]{x1D6E5}v$ define the Newton correction (we omit the index $k$ ). Since we only need $\unicode[STIX]{x1D6E5}v$ to extend the search space, we premultiply the first equation in (4.31) by the oblique projector $I_{n}-F^{\prime }(\unicode[STIX]{x1D717})vq^{\ast }$ , where $q\in \mathbb{C}^{n}$ is such that $q^{\ast }F(\unicode[STIX]{x1D717})v=0$ and normalized such that $q^{\ast }F^{\prime }(\unicode[STIX]{x1D717})v=1$ . This yields

Because of the orthogonality condition $v^{\ast }\unicode[STIX]{x1D6E5}v=0$ , we can rewrite the first equation in (4.32) as

which has the form of a Jacobi–Davidson correction equation. We can use $\operatorname{span}\{U,\unicode[STIX]{x1D6E5}v\}$ as the new search space and $\operatorname{span}\{Q,F(\unicode[STIX]{x1D717})v\}$ as the new test space. This process is repeated until we obtain a Ritz eigenpair with a small residual. We can expect quadratic convergence of this process if the correction equation (4.33) is solved exactly. As for linear eigenproblems (Sleijpen and van der Vorst Reference Sleijpen and van der Vorst1996), the correction equation does not need to be solved accurately to preserve convergence.

Pseudocode is presented in Algorithm 4.19. Usually only a few steps of a preconditioned Krylov subspace method such as GMRES (Saad and Schultz Reference Saad and Schultz1986) are necessary to solve the correction equation. If $P$ is a preconditioner for $F(\unicode[STIX]{x1D717})$ so that $Py=r$ is easy to solve, the preconditioner $P$ should be modified to

for the correction equation (4.33). Davidson (Reference Davidson1975) proposes using $P=\operatorname{diag}F(\unicode[STIX]{x1D717})$ for linear problems.

Variations of the Jacobi–Davidson approach are proposed by Sleijpen et al. (Reference Sleijpen, Booten, Fokkema and van der Vorst1996) for polynomial eigenvalue problems, Betcke and Voss (Reference Betcke and Voss2004) for symmetric NEPs, and Voss (Reference Voss2007) and Effenberger (Reference Effenberger2013b ) for nonsymmetric NEPs. Variants of a two-sided Jacobi–Davidson method are discussed in Hochstenbach and Sleijpen (Reference Hochstenbach and Sleijpen2003). The nonlinear Arnoldi method of Voss (Reference Voss2004a ) corresponds to using one step of residual inverse iteration in place of the Jacobi–Davidson correction equation, that is, (4.33) is replaced by $F(\unicode[STIX]{x1D70E})\unicode[STIX]{x1D6E5}v=F(\unicode[STIX]{x1D717})v$ . The linear system is then solved approximately, $\unicode[STIX]{x1D6E5}v\approx MF(\unicode[STIX]{x1D717})v$ , using an approximation $M$ to $F(\unicode[STIX]{x1D70E})^{-1}$ .

4.6 Hermitian NEPs

Assume that the matrix-valued function $F(z)$ is Hermitian, that is, $F(\bar{z})=F(z)^{\ast }$ for all $z\in \mathbb{C}$ . We can take advantage of this property in the algorithms described in Sections 4.1–4.2 when factorizing the Hermitian matrix $F(\unicode[STIX]{x1D706}^{(k)})$ using a block $\text{LDL}^{\ast }$ factorization. Also, when $\unicode[STIX]{x1D706}$ is real, the solution to the nonlinear scalar equation in (4.18) in the first step of the residual inverse iteration can be replaced by finding $z$ (i.e. the Rayleigh functional) such that $v^{(k)\ast }F(z)v^{(k)}=0$ . The resulting residual inverse iteration converges quadratically (Neumaier Reference Neumaier1985).

Special algorithms can be designed when the Hermitian matrix-valued function $F$ satisfies properties (A1)–(A3) in Section 3 on some real interval  $\mathbb{I}$ . It follows from Theorem 3.1 that the eigenvalues of $F$ can be characterized as minmax and maxmin values of the Rayleigh functional $p(x)$ , which is the only root of $x^{\ast }F(p(x))x=0$ in $\mathbb{I}$ . Indeed, Theorem 3.1 asserts that if $\unicode[STIX]{x1D706}_{k}$ is a $k$ th eigenvalue of $F(z)$ (that is, $\unicode[STIX]{x1D707}=0$ is the $k$ th largest eigenvalue of the Hermitian matrix $F(\unicode[STIX]{x1D706}_{k})$ ), then

where $\mathbb{S}_{k}$ denotes the set of all subspaces of $\mathbb{C}^{n}$ of dimension $k$ , and $\mathbb{K}_{p}$ is a subspace of $\mathbb{C}^{n}$ onto which the Rayleigh functional $p$ is defined. The minimum in (4.34) is attained for an invariant subspace of the Hermitian matrix $F(\unicode[STIX]{x1D706}_{k})$ corresponding to its $k$ largest eigenvalues and the maximum is attained for some $x\in \operatorname{null}F(\unicode[STIX]{x1D706}_{k})$ . This suggests the method in Algorithm 4.20 for computing the $j$ th eigenvalue of $F$ , called safeguarded iteration (Werner Reference Werner1970).

Niendorf and Voss (Reference Niendorf and Voss2010) show that if $\unicode[STIX]{x1D706}_{j}$ is a simple eigenvalue of $F$ , then the safeguarded iteration converges locally and quadratically to $\unicode[STIX]{x1D706}_{j}$ . For large sparse problems, the safeguarded iteration can be embedded into an iterative projection method such as the Jacobi–Davidson scheme or the nonlinear Arnoldi scheme described in Section 4.5; see Voss (Reference Voss and Neittaanmäki2004b ) and Niendorf and Voss (Reference Niendorf and Voss2010).

For large sparse NEPs satisfying properties (A1)–(A3) defined in Section 3, Szyld and Xue (Reference Szyld and Xue2015) describe several variants of a preconditioned conjugate gradient method, which make use of the nonlinear variational principle (Theorem 3.1) and the nonlinear Cauchy interlacing theorem (Theorem 3.2).

Figure 4.22.

Basic MATLAB implementation of the safeguarded iteration to compute the five smallest eigenvalues of the loaded_string problem (1.3). Lines 2–5 define the NEP $F$ and line 6 specifies the parameters for the iteration.