Disguised and new Quasi-Newton methods for nonlinear eigenvalue problems

TL;DR

This paper develops new quasi-Newton algorithms for solving nonlinear eigenvalue problems, providing convergence analysis and unifying existing methods under a common framework.

Contribution

It introduces novel quasi-Newton methods for NEPs, analyzes their convergence, and interprets existing algorithms as special cases of these methods.

Findings

New quasi-Newton algorithms for NEPs are proposed.

Convergence properties are established using Keldysh's theorem.

Existing methods like residual inverse iteration are shown to be quasi-Newton methods.

Abstract

In this paper we take a quasi-Newton approach to nonlinear eigenvalue problems (NEPs) of the type $M(\lambda)v=0$, where $M:\mathbb{C}\rightarrow\mathbb{C}^{n\times n}$ is a holomorphic function. We investigate which types of approximations of the Jacobian matrix lead to competitive algorithms, and provide convergence theory. The convergence analysis is based on theory for quasi-Newton methods and Keldysh's theorem for NEPs. We derive new algorithms and also show that several well-established methods for NEPs can be interpreted as quasi-Newton methods, and thereby provide insight to their convergence behavior. In particular, we establish quasi-Newton interpretations of Neumaier's residual inverse iteration and Ruhe's method of successive linear problems.

Equations

$$M(\lambda )v=0,v\ne 0$$

$$u^{H}M(\lambda )=0,u\ne 0.$$

$$F\left(\left[\begin{array}{c}v\\ \lambda \end{array}\right]\right):=\left[\begin{array}{c}M(\lambda )v\\ c^{H}v-1\end{array}\right]=0$$

$${\tilde{J}}_k\left[\begin{array}{c}{x}_{k+1}-x_k\\ {\mu }_{k+1}-{\mu }_k\end{array}\right]=-\left[\begin{array}{c}M({\mu }_k)x_k\\ c^H{x}_k-1\end{array}\right]=-F_k$$

$${\tilde{J}}_k\approx J_k=J\left(\left[\begin{array}{c}{x}_k\\ {\mu }_k\end{array}\right]\right):=\left[\begin{array}{cc}M({\mu }_k)& M^{\prime }({\mu }_k)x_k\\ c^H& 0\end{array}\right]$$

$${\tilde{J}}_{1,k}:=\left[\begin{array}{cc}M(\sigma )& M^{\prime }(\sigma )x_0\\ c^H& 0\end{array}\right],$$

$${\tilde{J}}_{2,k}:=\left[\begin{array}{cc}M(\sigma )& M^{\prime }({\mu }_k)x_k\\ c^H& 0\end{array}\right].$$

$${\tilde{J}}_{3,k}:=\left[\begin{array}{cc}M(\sigma )& M[{\mu }_{k+1},{\mu }_k]x_k\\ c^H& 0\end{array}\right],$$

$$M[\lambda ,\mu ]=\{\begin{array}{ll}\frac{M(\lambda )-M(\mu )}{\lambda -\mu }& \lambda \ne \mu \\ M^{\prime }(\lambda )& \lambda =\mu .\end{array}$$

$${\tilde{J}}_{4,k}:=\left[\begin{array}{cc}M({\mu }_k)& M^{\prime }({\mu }_k)x_{k+1}\\ c^H& 0\end{array}\right].$$

$$c^H{x}_1=c^H{x}_2=\cdots =1,$$

$$\Delta x_k+\Delta {\mu }_{k}M(\sigma {)}^{-1}{M}^{\prime }(\sigma )x_0=-M(\sigma {)}^{-1}M({\mu }_k)x_k=-y_k$$

$$y_k:=M(\sigma {)}^{-1}M({\mu }_k)x_k.$$

$${\mu }_{k+1}-{\mu }_k=\Delta {\mu }_k=-{\alpha }_0{c}^H{y}_k$$

$$q_0$$

$${\alpha }_0$$

$$x_{k+1}=x_k-y_k-\Delta {\mu }_k{q}_0.$$

$$w_{\sigma }^H:=c^{H}M(\sigma {)}^{-1}.$$

$$\Delta {\mu }_k$$

$$z_k$$

$${\mu }_{k+1}$$

$$x_{k+1}$$

$$c^H\Delta x_k+w_{\sigma }^{H}M[{\mu }_{k+1},{\mu }_k]x_k\Delta {\mu }_k$$

$$w_{\sigma }^H\frac{M({\mu }_{k+1})-M({\mu }_k)}{{\mu }_{k+1}-{\mu }_k}{x}_k\Delta {\mu }_k$$

$$w_{\sigma }^{H}M({\mu }_{k+1})x_k$$

$$\Delta x_k+M(\sigma {)}^{-1}(M({\mu }_{k+1})-M({\mu }_k))x_k$$

$$x_{k+1}$$

$${\mu }_{k+1}$$

$$v_{k+1}$$

$$M({\mu }_k)(x_{k+1}-x_k)+\Delta {\mu }_k{M}^{\prime }({\mu }_k)x_{k+1}$$

$$M({\mu }_k)x_{k+1}+\Delta {\mu }_k{M}^{\prime }({\mu }_k)x_{k+1}$$

$$\left[\begin{array}{c}{x}_{k+1}-v\\ {\mu }_{k+1}-\lambda \end{array}\right]=A\left[\begin{array}{c}{x}_k-v\\ {\mu }_k-\lambda \end{array}\right]+\mathcal{O}\left({\Vert \left[\begin{array}{c}{x}_k-v\\ {\mu }_k-\lambda \end{array}\right]\Vert }^2\right).$$

$$A_1=\left[\begin{array}{cc}(I-{\alpha }_0{q}_0{c}^H)M(\sigma {)}^{-1}\left(M(\sigma )-M(\lambda )\right)& (I-{\alpha }_0{q}_0{c}^H)M(\sigma {)}^{-1}{M}^{\prime }(\lambda )v\\ {\alpha }_0{c}^{H}M(\sigma {)}^{-1}\left(M(\sigma )-M(\lambda )\right)& {\alpha }_0{c}^{H}M(\sigma {)}^{-1}{M}^{\prime }(\lambda )v\end{array}\right].$$

$$\left[\begin{array}{c}{v}_{k+1}\\ {\lambda }_{k+1}\end{array}\right]=\phi \left(\left[\begin{array}{c}{v}_k\\ {\lambda }_k\end{array}\right]\right)=\left[\begin{array}{c}{v}_k\\ {\lambda }_k\end{array}\right]-{\tilde{J}}_1{\left(\left[\begin{array}{c}{v}_k\\ {\lambda }_k\end{array}\right]\right)}^{-1}F\left(\left[\begin{array}{c}{v}_k\\ {\lambda }_k\end{array}\right]\right).$$

$$J_{\ast }:=J\left(\left[\begin{array}{c}v\\ \lambda \end{array}\right]\right)=\left[\begin{array}{cc}M(\lambda )& M^{\prime }(\lambda )v\\ c^H& 0\end{array}\right].$$

Full text

1 Introduction

One of the most common techniques to improve the convergence or efficiency of Newton’s method for nonlinear systems of equations is to replace the Jacobian matrix with a different matrix. Among these quasi-Newton method constructions, sometimes called inexact Newton methods, the most common variation is to keep the Jacobian matrix constant. The factorization of this matrix can be precomputed before carrying out the iterations. This is beneficial, e.g., in situations where the problem stems from a discretization of a PDE, as the resulting system is often large and the Jacobian matrix is sparse with a structure allowing a sparse LU-factorization to be pre-computed.

In this paper we consider nonlinear eigenvalue problems (NEPs) of the type

[TABLE]

where $M:\Omega\rightarrow\mathbb{C}^{n\times n}$. There are various flavors of Newton’s method available in the literature (further discussed below) for this class of NEPs. Some of these methods do have the property that the matrix in the linear system to be solved in every iteration remains constant. However, these methods are in general not seen as Jacobian matrix modifications of Newton’s method, but are often derived from quite different viewpoints. In this paper we investigate methods resulting from modifying the Jacobian matrix in various ways, and illustrate differences, similarities and efficiency of the resulting methods. It turns out that several well-established approaches for NEPs can be viewed as quasi-Newton methods.

In the NEP-class that we consider in this paper $\Omega\subset\mathbb{C}$ is a closed set, $M$ is analytic in $\Omega$ and we suppose $\lambda\in\Omega$. We call the vector $v$ a (right) eigenvector if it satisfies (1) and $u$ the left eigenvector if it satisfies

[TABLE]

We call $(\lambda,v,u)$ an eigentriplet of (1). Without loss of generality we phrase the NEP as a system of equations

[TABLE]

which is equivalent to (1) if $c\in\mathbb{C}^{n}$ is not orthogonal to the eigenvector $v$. This condition is not problematic in practice since $c$ can be chosen freely and thus will generically not be orthogonal to any eigenvector. The quasi-Newton approach to (3) consists of generating sequences of approximations $(\mu_{1},x_{1}),(\mu_{2},x_{2}),\ldots$ from the relation

[TABLE]

where $\tilde{J}_{k}$ is an approximation of the Jacobian matrix

[TABLE]

The eigenvector and eigenvalue updates will be denoted $\Delta x_{k}=x_{k+1}-x_{k}$ and $\Delta\mu_{k}=\mu_{k+1}-\mu_{k}$.

We consider four specific modifications of the Jacobian matrix, briefly justified as follows. From a quasi-Newton perspective, the most common approach consists of keeping the Jacobian matrix constant, i.e., setting

[TABLE]

where $\sigma=\lambda_{0}$ is the starting value for the eigenvalue and $x_{0}$ the starting vector for the eigenvector approximation. In this manuscript we refer to this as Quasi-Newton 1 (QN1). In Section 2.1 we show how (4) with Jacobian matrix approximation (6) can be reformulated such that in every iteration we need to solve one linear system associated with the matrix $M(\sigma)$.

We shall later show (in Section 2.1) that a more accurate approximation of the Jacobian matrix leads to an algorithm which has the same computational cost per iteration as QN1, i.e., it involves the solution of one linear system with the matrix $M(\sigma)$ per iteration. More precisely, we keep only the (1,1)-block constant by setting

[TABLE]

The quasi-Newton method (4) with Jacobian matrix approximation (7) will be refered to as Quasi-Newton 2 (QN2).

We also investigate a method (which we call Quasi-Newton 3) corresponding to keeping the $(1,1)$-block constant as in QN2, but also replace the derivative in the (1,2)-block. We replace the derivative with a finite difference involving the future eigenvalue approximation $\lambda_{k+1}$, i.e., we set

[TABLE]

where we use the standard notation for divided differences

[TABLE]

Note that the modification of the Jacobian matrix in (8) makes the iteration (4) implicit in the sense that the formula (4) for the new approximation $(\lambda_{k+1},v_{k+1})$ involves $\lambda_{k+1}$ in a nonlinear way. Many implicit variations of Newton’s method have been considered in the literature, see e.g., [12, 18] and references therein. It turns out that certain implicit variations of Newton’s method improve the convergence and sometimes even increase the convergence order. In contrast to many other implicit Newton methods, our choice of the Jacobian matrix is done with the goal of having a method whose iterates can be computed without solving a (computationally demanding) nonlinear system of equations. This is possible for the specific choice (8), as we shall illustrate in Section 2.

We consider one more modified Jacobian matrix which also leads to an implicit method, but now implicit in the eigenvector. The vector $x_{k}$ in the (1,2) block is replaced by the future vector $x_{k+1}$, such that we obtain

[TABLE]

Our study has the following conclusions and contributions:

• •

QN1 and QN2 can be phrased as algorithms only involving one linear solve with $M(\sigma)$ per iteration (Section 2.1-2.2)

• •

QN3 is equivalent to Neumaier’s residual inverse iteration [25] (Section 2.3)

• •

QN4 is equivalent to Ruhe’s method of successive linear problems [29] (Section 2.4)

• •

We provide exact characterizations of the convergence of factor for QN1 and QN2, and establish that the convergence factor of QN2 and QN3 are identical (Section 3.1-3.2)

• •

We show how to adapt fundamental theory for inexact Newton method [5] to study QN4 (Section 3.3)

• •

We provide generalizations of convergence rate dependence on eigenvalue clustering analogous to methods for linear eigenvalue methods (Section 4)

The studied properties of the methods are illustrated in numerical simulations in Section 5.

Newton’s method for linear and nonlinear eigenvalue problems has been studied for decades, and the field is still under active development, as can be observed in the summaries in [22, 40]. The technique of deriving methods using an augmented system as in (3) was investigated already in 1950’s by Unger [38], and was key to characterizing the relationship of inverse iteration (for linear eigenvalue problems) as described by Peters and Wilkinson in [27]. Newton’s method based on solving the nonlinear equation $\beta(\lambda)=0$ where $M(\lambda)v=\beta(\lambda)e_{p}$ was presented in [26] and independently by an essentially equivalent procedure by Lancaster [20]. In 1970’s Ruhe [29] also pointed out the relevance of Newton-type methods and how they relate to the inverse iteration. He used the augmented systems of equations to derive variations of the inverse iteration for NEPs as well as the method of successive linear problems.

Newton-type methods have also more recently received considerable attention, e.g., in the PhD thesis of Schreiber [30], where two-sided generalizations of inverse iteration methods as well as Jacobi-Davidson methods are developed. The recent results by Effenberger and Kressner [19, 7] contain a generalization of Newton-type methods that allows the computation of several eigenvalues simultaneously. This block Newton approach has been successful in the setting of continuation methods [4]. Variants of Newton methods where the linear system associated with $M(\sigma)$ is only solved to some accuracy have been studied in [34]. A recent variant of the rational Krylov method can also be interpreted in a Newton-setting [1]. There are several convergence results for the residual inverse iteration and other Newton type methods [34, 35, 17, 37] which are mostly presented in a separated fashion without using quasi-Newton interpretations and results for quasi-Newton methods.

2 Explicit reformulations of the quasi-Newton methods

The formulations of the QN-iterations above are not practical in general. In fact, it is even questionable to call the formulations of QN3 and QN4 iterative, since the Jacobian matrix depends on quantities in an implicit way. Nevertheless, it turns out that certain reformulations of QN3 and QN4 allow us to explicitly compute sequences of approximations, which satisfy (4). Also QN1 and QN2 have to be reformulated in order to become practical. We show how to carry out this reformulation to obtain algorithms which do not require solving many linear systems for different matrices, but only for $M(\sigma)$.

We first make an observation in common for all the considered methods. The last row of the correction equation (4) is the same for all methods and for all the choices of $\tilde{J}_{k}$ as in (6)-(10) we see that $c^{H}(x_{k+1}-x_{k})=-c^{H}x_{k}+1$. By induction, this implies that

[TABLE]

i.e., all iterates (except possibly the first iterate) are normalized.

2.1 The $n$-dimensional form of QN1

The iteration (4) with the Jacobian matrix approximation (6) can be reformulated as follows. By multiplying the first block equation in (4) with $M(\sigma)^{-1}$ from the left yields

[TABLE]

where we define

[TABLE]

Moreover, by multiplication of (12) from the left with $c^{H}$ and using the fact that $c^{H}\Delta x_{k}=0$ due to (11) we have

[TABLE]

with

[TABLE]

The above equations can be combined into an algorithm. As a precomputation we compute $\alpha_{0}$ in (15b) and $q_{0}$ in (15a), and in the iteration we compute $y_{k}$ from (13), $\Delta\mu_{k}$ from (14) and subsequently update

[TABLE]

The algorithm is summarized in Algorithm 1.

Remark 1** (Properties of Algorithm 1).**

*The advantage of Algorithm 1 over the original formulation (4) consists in the fact that the matrix in the linear system to be solved in every step is $M(\sigma)$. Therefore, a pre-factorization can be carried out in the original problem size. This is of advantage when the problem stems from a partial differential equation and the augmented matrix $\tilde{J}_{1,k}$ may be more difficult to factorize. Although the Jacobian matrix approximation $\tilde{J}_{1,k}\approx J_{k}$ is the most common in the context of quasi-Newton approximations, it does not appear very competitive in this setting. The other Jacobian matrix approximations, which are more accurate, lead to better methods in terms of convergence and do not in general require more computation. *

2.2 The $n$-dimensional formulation of QN2

A reformulation of (4) with the Jacobian matrix approximation (7) follows similar steps as the derivation in the previous section. To simplify the notation we set

[TABLE]

The derivation leads up to the formulas

[TABLE]

More precisely, formula (18a) stems from left-multiplying the first block row in (4) with $w_{\sigma}^{H}=c^{H}M(\sigma)^{-1}$ and (18b) and (18d) stems from left-multiplying the first block row in (4) with $M(\sigma)^{-1}$. We summarize the resulting method in Algorithm 2.

Remark 2** (Properties of Algorithm 2).**

Note that similar to Algorithm 1, Algorithm 2 requires only one linear solve with the matrix $M(\sigma)$ per iteration. Since Algorithm 2 corresponds to a more accurate approximation of the Jacobian matrix, it is expected to converge faster than Algorithm 1. This difference characterized theoretically and computationally in Section 3.1 and Section 5.

Algorithm 2 involves the vector $w_{\sigma}$, which can be computed as in (17), i.e., it would require one additional linear solve with $M(\sigma)^{H}$. This extra linear solve can however be avoided by treating $w_{\sigma}$ as a fixed vector (chosen by the user) and then using that $c$ is arbitrary such that we can chose it as $c^{H}=w_{\sigma}^{T}M(\sigma)$. This works rather well in practice, but fixing $w_{\sigma}$ instead of $c$ may make the convergence factor larger if $\sigma$ is close to the eigenvalue, as we shall further illustrate in Section 4.

2.3 The explicit formulation of QN3 is residual inverse iteration

Although the Jacobi approximation $J_{k}\approx\tilde{J}_{3,k}$ in (8) involves eigenvalue approximations not yet computed, we can proceed with an elimination procedure similar to Section 2.1 and Section 2.2.

We multiply the first block in equation (4) from the left with $w_{\sigma}^{H}=c^{H}M(\sigma)^{-1}$. This gives an equation which we can simplify as follows:

[TABLE]

Note that (19c) is a scalar-valued equation, in one unknown variable $\mu_{k+1}$, since $x_{k}$ can be viewed as a known vector. In fact, if we treat $\lambda$ as a function of $x$, this is the inverse function of $w^{H}M(\lambda)x=0$ is, which is commonly known as the Rayleigh functional or generalized Rayleigh quotient [39, 42, 41]. This function generally exists, at least in a neighborhood of a simple eigenvalue [17, Proposition 2.1], and it is a computable quantity for many problems.

By multiplying the first block row in (4) from the left by $M(\sigma)^{-1}$, we obtain

[TABLE]

Under the assumption that the Rayleigh functional in (19c) is computable, the relations (19c) and (20b) form an explicit algorithm. In fact, this algorithm is already extensively used in current research, where it is commonly known as residual inverse iteration and it was first introduced by Neumaier in [25]. Residual inverse iteration also forms the basis of some recent state-of-the-art algorithms for NEPs, most importantly the nonlinear Arnoldi method [39].

Theorem 3**.**

The Quasi-Newton method (4) with the modified Jacobian matrix (8) is equivalent to residual inverse iteration as described in [25].

Remark 4** (Relation between quasi Newton variant 2 and residual inverse iteration).**

*Due to the analyticity of $M(\lambda)$, the residual inverse iteration (19c) in (20b) can also be expressed as *

[TABLE]

From these formulas we see directly that for linear eigenvalue problems where $M(\lambda)=A-\lambda I$, QN2 and residual inverse iteration are equivalent. Hence, they are both generalizations of the standard inverse iteration method. This is consistent with the convergence analysis in Section 3 which shows that QN2 and QN3 have the same convergence factor.

2.4 The explicit formulation of QN4 is the method of successive linear problems

In the previous subsection we saw that iterates satisfying (4) with the modified Jacobian matrix (8) can be computed in practice and the resulting algorithm is in fact equivalent to a well-known method. The Jacobi approximation $J_{k}\approx\tilde{J}_{4,k}$ in (10) also involves a quantity which we do not have access to at iteration $k$, the vector $x_{k+1}$. We now show that similar to QN3, we can carry out an elimination such that the update can be computed in an explicit way. This algorithm also turns out to be equivalent to a well-established method.

The first block row in (4) with approximation (10) simplifies to

[TABLE]

Since we know that the iterates $x_{1},x_{2},\ldots$ are normalized, we directly identify (21b) as a (linear) generalized eigenvalue problem where $\Delta\mu_{k}$ is the eigenvalue. Hence, we can construct an iteration satisfying (4) with Jacobi approximation (10) by repeatedly solving the generalized eigenvalue problem (21b) and updating the eigenvalue $\mu_{k+1}=\mu_{k}+\Delta\mu_{k}$. This method is known as the method of successive linear problems and was studied and used by Ruhe in [29], where it was justified directly from a Taylor expansion of $M(\lambda)$.

Theorem 5**.**

The Quasi-Newton method (4) with the modified Jacobian matrix (10) is equivalent to the method of successive linear problems [29].

3 Local convergence analysis

3.1 Convergence factor analysis of QN1 and QN2

In order to characterize the convergence of QN1 and QN2 we will derive a first-order result. More precisely, we will describe the local behavior, if $(x_{k},\mu_{k})$ is close to $(\lambda,v)$, we describe the behavior with a matrix $A\in\mathbb{C}^{(n+1)\times(n+1)}$ such that

[TABLE]

In general, we have linear convergence, with a local convergence factor given by the spectral radius of $A$. The explicit form of $A$ for our first two quasi-Newton methods are given in following theorems.

Theorem 6** (Local convergence Algorithm 1).**

Suppose the sequence $(\mu_{1},x_{1})$, $(\mu_{2},x_{2}),\ldots$ is generated by Algorithm 1 started with $(\mu_{0},x_{0})$ and suppose the sequence converges to the eigenpair $(\lambda,v)$. Then, the sequence satisfies (22) with $A=A_{1}$ where

[TABLE]

Proof.

For notational convenience let $\tilde{J}_{1}:\mathbb{C}^{n+1}\rightarrow\mathbb{C}^{n+1}$ denote the function corresponding to $\tilde{J}_{1,k}$ in (6). In this fixed-point setting, our quasi-Newton method can be expressed as

[TABLE]

The $A$-matrix in (22), corresponding to a fixed point map, is given by the Jacobian of $\varphi$. In our case this can be explicitly expressed,

by using the structure of the Jacobian matrix evaluated in the eigenpair

[TABLE]

More precisely, we have

[TABLE]

∎

Theorem 7** (Local convergence Algorithm 2).**

Suppose the sequence $(\mu_{1},x_{1})$, $(\mu_{2},x_{2}),\ldots$ is generated by Algorithm 1 started with $(\mu_{0},x_{0})$ and suppose the sequence converges to the eigenpair $(\lambda,v)$. Then, the sequence satisfies (22) with $A=A_{2}$ where

[TABLE]

where $\alpha:=(w_{\sigma}^{H}M^{\prime}(\lambda)v)^{-1}$ and $q:=M(\sigma)^{-1}M^{\prime}(\lambda)v$.

Proof.

The proof follows the same reasoning as in the proof of Theorem 6, except that the Jacobian matrix of the fixed point map $\varphi$ in (26a) in this case becomes

[TABLE]

The application of the Schur complement formula directly leads to (27). ∎

3.2 Local convergence of QN3

Note that QN3 is not a fixed point iteration in the formulation (8). However, since residual inverse iteration and QN3 are equivalent in the sense of Theorem 3, we already have a convergence factor available in [17]. More surprisingly, the convergence factor for residual inverse iteration (given in [17]) is identical to the convergence factor of QN2.

Corollary 8** (Convergence factor equivalence QN2 and QN3).**

The non-zero eigenvalues of $A_{2}$ given in (27) are the same as the non-zero eigenvalues of

[TABLE]

and the convergence factors of QN2 and QN3 are the same.

Proof.

Since QN3 is equivalent to residual inverse iteration according to Theorem 3 we can directly use the convergence characterization in [17]. More precisely, [17, Theorem 3.1] states that the convergence factor of residual inverse iteration is the largest eigenvalue of the matrix $B$ in (28).

It remains to show that the non-zero eigenvalues of $B$ are the same as the non-zero eigenvalues of $A_{2}$ given in (27). Clearly, the non-zero eigenvalues of $A_{2}$ are the non-zero eigenvalues of the (1,1)-block of $A_{2}$. The equivalence is based on the general property that if a matrix $C$ and vectors $c$ and $v$ satisfy, $c^{H}C=0$ and $c^{H}v=1$, then $C$ and $C(I-vc^{H})$ have same non-zero eigenvalues. This can be seen from the fact that if $\gamma_{i}z_{i}=Cz_{i}$, then $\gamma_{i}c^{H}z_{i}=c^{H}Cz_{i}=0$ and $c^{H}z_{i}=0$ if $\gamma_{i}\neq 0$, i.e., $C(I-vc^{H})z_{i}=Cz_{i}-vc^{H}z_{i}=\gamma_{i}z_{i}$. By using this general property where $C$ is the (1,1)-block of (27), we obtain $B=C(I-vc^{H})$ where $B$ is given in (28). ∎

3.3 Local convergence of QN4

The quasi-Newton method corresponding to (10) is different in character in comparison to the other quasi-Newton methods we have considered. This methos does not involve the computation of a linear system for a constant matrix. Although some convergence results for the method of successive linear problems (and therefore also QN4) are available in the literature [15], it is natural in our quasi-Newton approach to characterize the convergence using general results for quasi-Newton methods. It turns out that (10) is a very accurate approximation of the Jacobian matrix and we can apply results for quasi-Newton methods given in [5].

The characterization in [5] is mainly based on a quantity which describes the inexactness the quasi-Newton method by comparing it with a step involving the exact Jacobian matrix. More precisely, we consider the vector $r_{k}$ in [5], which in our context becomes

[TABLE]

Various results are given in [5] in terms of the norm of $r_{k}$. In particular, the result [5, Theorem 3.3] demonstrates that we have a local quadratic convergence if $\|r_{k}\|\leq\mathcal{O}(\|F_{k}\|^{2})$. Due to the fact that all iterates are normalized as in (11), this condition can be simplified in our case and stated as

[TABLE]

In order to use this result, we now suppose that the Jacobian matrix in the eigenpair denoted $J_{*}$ and defined by (25) is invertible. The invertability condition is satisfied if the eigenpair $(\lambda,v)$ is such that $c^{H}v\neq 0$ and the eigenvalue $\lambda$ is simple.

This implies that the function $F:\mathbb{C}^{n+1}\rightarrow\mathbb{C}^{n+1}$ is invertible in a neighborhood of $(\lambda,v)$. If $(\mu_{k},x_{k})$ is in this neighborhood, we have from the implicit function theorem that

[TABLE]

Theorem 9** (Convergence QN4).**

Suppose $(\mu_{1},x_{1}),(\mu_{2},x_{2}),\ldots$ are iterates satisfyng (10) and suppose they converge to the eigenpair $(\lambda,v)$. If $(\lambda,v)$ is a simple or semi-simple eigenpair, then $(\mu_{k},x_{k})$ converges at least quadratically.

Proof.

From (29) and properties of the two-norm we have

[TABLE]

In the last step we used that, due to the assumption that $(\mu_{1},x_{1}),(\mu_{2},x_{2}),\ldots$ converges to an eigenpair, we have

[TABLE]

for sufficiently large $k$. By using the expansion (30) we find directly that $\|r_{k}\|\leq\mathcal{O}(\|M(\mu_{k})x_{k}\|^{2})$, which by [5, Theorem 3.3] implies quadratic convergence. ∎

4 Interpretation of convergence factors

In order to provide further insight to the linearly convergent quasi-Newton methods, we now make characterizations of the matrix $A$. To this end, we use results associated with (what is commonly referred to as) Keldysh’s theorem; see the general formulation in [23], and the more recent descriptions in the context of NEPs, e.g. [3, 34]. For simple eigenvalues, Keldysh’s theorem implies that there exists a function $R_{1}(\sigma)$, analytic in a neighborhood of the eigenvalue, such that

[TABLE]

where $(\lambda_{1},v_{1},u_{1})$ is the eigentriplet of a simple eigenvalue.

We first observe that the convergence factor of QN2 and QN3 can be directly analyzed with Keldysh’s theorem, since

[TABLE]

In the last inequality we used that

[TABLE]

The relationship (32) indicates that the convergence factor depends linearly on the shift eigenvalue distance, and linearly in $M^{\prime}(\lambda_{1})$.

Remark 10** (Double non-semisimple eigenvalues).**

In case the eigenvalue $\lambda_{1}$ is a double non-semisimple eigenvalue, there exist so called generalized eigenvectors $\widetilde{v}_{1}$ and $\widetilde{u}_{1}$ such that

[TABLE]

where $u$ and $v$ are the left and right eigenvectors corresponding to $\lambda_{1}$. According to [3, Theorem 2.6] with $L=1$ and $m_{1}=2$, there is a neighborhood $U$ of $\lambda_{1}$ where we have the expansion

[TABLE]

where $R(\sigma)$ is analytic in $U$. Then, instead of (27) the iteration matrix $A_{2}$ is of the form

[TABLE]

where $R(\sigma)$ contains the contribution from all the eigenvalues other than $\lambda_{1}$. Thus the iteration matrix $A-vc^{T}$ contains the factor $\frac{1}{\sigma-\lambda_{1}}$ and unlike in the case of a simple eigenvalue, the convergence factor is not asymptotically proportional to $|\lambda_{1}-\sigma|$.

4.1 Eigenvalue clustering and condition number

The convergence factor bound (32) provides insight on how the convergence depends on the shift if the shift-eigenvalue distance is small (consistent with what was pointed out in [17]). We now show that different insight can be provided by using a more general form of Keldysh’s theorem. This applies to the situation when the shift-eigenvalue distance is not necessarily small.

Inverse iteration for linear eigenvalue problems (with normalization $c^{H}v=1$) has the following property for diagonalizable matrices. The convergence factor for the eigentriplet $(\lambda_{1},v_{1},u_{1})$ can be bounded in terms of reciprocal eigenvalue-shift distances weighted with the condition number

[TABLE]

where $\kappa_{i}=\|u_{i}\|\|v_{i}\|/|u_{i}^{H}v_{i}|$ is the eigenvalue condition number (following the standard definition [21]) and $P_{1}$ is the projector $P_{1}=I-v_{1}c^{H}$.

In order to generalize this property, we use a more general form of Keldysh’s theorem. We let $\Gamma\subset\Omega$ be a simple, closed, piecewise-smooth curve and denote the eigenvalues in its interior by $\lambda_{1},\cdots,\lambda_{k}\in\operatorname{int}(\Gamma)$. Then, Keldysh’s theorem states that

[TABLE]

where $R_{\Gamma}$ is analytic in $\operatorname{int}(\Gamma)$. The following result provides an analogue of the eigenvalue clustering property (34), under the assumption that $R_{\Gamma}(\sigma)$ is small.

Corollary 11** (Eigenvalue clustering).**

*Suppose that $\Omega$ is a closed simply connected domain with boundary $\Gamma$. Suppose that $M$ is analytic in this domain and that all the eigenvalues are simple. Denote the corresponding eigentriplets by $(\lambda_{1},v_{1},u_{1})$,$\dots$,

$(\lambda_{k},v_{k},u_{k})$, with normalization $c^{H}v_{1}=\cdots=c^{H}v_{k}=1$. Then, the convergence factor for QN2 and QN3 are bounded by

[TABLE]

where $\kappa_{i}$ is the eigenvalue condition number for NEPs,

[TABLE]

and $R_{\Gamma}(\sigma)$ is the remainder-term in Keldysh’s theorem in (35).

Proof.

The result follows the steps in (32) but instead using the form (35) and matrix norm inequalities. ∎

4.2 Characterization of Keldysh’s remainder term $R_{\Gamma}(\sigma)$

The analysis in the previous subsection indicates that a dependence on the eigenvalue clustering similar to the linear case can be expected under the condition that $R_{\Gamma}$ is small. Keldysh’s theorem is, in a certain sense, a matrix version of the partial fraction expansion of an analytic function (also known as Mittag-Leffler’s theorem) and has been characterized in [10]. See [45] for a partial fraction expansion characterization for delay eigenvalue problems. A precise characterization of $R_{\Gamma}$ and of its norm for the general case using [10] is somewhat involved. We take a less ambitious approach and point out certain situations where it is small or vanishes. Although $R_{\Gamma}$ does not always vanish in the limit when $\Gamma$ encloses $\mathbb{C}$, it does vanish under certain assumptions. We characterize several of such situations next, and note that these results are general for nonlinear eigenvalue problems. Therefore they may be of interest also outside the scope of quasi-Newton methods.

First we need an explicit formulation of the remainder term $R_{\Gamma}$.

Lemma 12**.**

Let $\Gamma\subset\Omega$ be a simple, closed, piecewise-smooth curve and denote the eigenvalues of (1) in the interior of $\Gamma$ by $\lambda_{1},\cdots,\lambda_{k}\in\operatorname{int}(\Gamma)$ and suppose $\lambda_{1},\ldots,\lambda_{k}$ are simple. Then, the analytic function $R_{\Gamma}(z)$, i.e., the analytic part of $M(z)^{-1}$ in $\mathrm{int}(\Gamma)/\{\lambda_{1},\ldots,\lambda_{k}\}$, given by the Keldysh theorem, has the integral representation

[TABLE]

Proof.

We can apply the Cauchy integral formula, since $R_{\Gamma}$ is analytic, and use equation (35)

[TABLE]

The Cauchy residue theorem implies that $\int_{\Gamma}\frac{1}{(\lambda-\lambda_{i})(\lambda-z)}\,\mathrm{d}\lambda=0$, which proves (36). ∎

This leads directly to a sufficient condition for vanishing $R_{\Gamma}$, involving $M(\lambda)^{-1}$ in the limit $\lambda\rightarrow\infty$.

Lemma 13**.**

Suppose $M(\lambda)$ is analytic in $\mathbb{C}$ and suppose all eigenvalues are simple. Moreover, suppose

[TABLE]

Then, the set of eigenvalues is finite and

[TABLE]

for all $\lambda\in\mathbb{C}\setminus\{\lambda_{1},\lambda_{2},\ldots,\lambda_{k}\}$, where $\lambda_{i}$, $1\leq i\leq k$, denote the eigenvalues of $M$.

Proof.

Since (37) implies that for every $\varepsilon>0$ there exists an $R$ such that $\sup_{|z|>R}\|M(z)^{-1}\|<\varepsilon$, we have in particular that there exists $r$ such that $M(z)^{-1}$ has no poles outside a disk of radius $r$. The eigenvalues of (1) are roots of the analytic function $\det(M(\lambda))$. An analytic function only has a finite number of roots in a compact subset of the complex plane, and we therefore only have a finite number of eigenvalues in the disk of radius $r$. Take $\Gamma$ to be a circle of radius $r$. Then, using the representation given by Lemma 12, we get the bound

[TABLE]

Letting $r\rightarrow\infty$, the claim follows. ∎

Remark 14** (Generalizations to higher order multiplicities).**

Lemma 12 and Lemma 13 have the assumption that the eigenvalues are simple. A generalization to higher algebraic and/or geometric multiplicities seems feasible but more involved. Then, the expression (36) follows from [3, Corollary 2.8], and the fact that when the contour $\Gamma$ encircles $\lambda_{i}$ and $z$, and $\ell\geq 1$, it holds that

[TABLE]

by the residue theorem. Then the condition $||M(\lambda)^{-1}||\rightarrow 0$ as $|\lambda|\rightarrow\infty$, and the inequality (39) imply that $R_{\Gamma}(z)=0$.

From the above lemmas we conclude the following result which states that the $R_{\Gamma}$ vanishes if the NEP is the sum of a polynomial with leading non-singular coefficient and a term which decays sufficiently fast.

Theorem 15**.**

Suppose $M$ is analytic in $\mathbb{C}$ and suppose all eigenvalues are simple. Moreover, suppose $M(\lambda)$ is of the form

[TABLE]

where

$P(\lambda)=\sum_{i=0}^{N}A_{i}\lambda^{i}$* for some matrices $A_{0},\ldots,A_{N}\in\mathbb{C}^{n\times n}$ such that that $A_{N}$ is non-singular, and* 2. 2.

$\lambda^{-N}F(\lambda)\rightarrow 0$* as $\left|\lambda\right|\rightarrow\infty$.*

Then, $R_{\Gamma}(\lambda)=0$ and the representation (38) holds for $M(\lambda)^{-1}$.

Proof.

We see that

[TABLE]

where B(\lambda)=\big{(}\sum\limits_{\ell=0}^{N-1}\lambda^{\ell-N}A_{\ell}+\lambda^{-N}F(\lambda)\big{)}A_{N}^{-1}. By the assumption 2 above, it holds that

[TABLE]

Using the bound $||(I-A)^{-1}||\leq(1-||A||)^{-1}$ for $||A||<1$ (see [13, pp. 351]), and choosing $\lambda$ such that $||B(\lambda)||<1$, we see from (40) that

[TABLE]

and therefore $||M(\lambda)^{-1}||\rightarrow 0,\quad\textrm{as}\quad\left|\lambda\right|\rightarrow\infty$. ∎

In addition to the polynomial eigenvalue problem with invertible leading coefficient matrix (eventually perturbed in the sense of Theorem 15) the conditions of Lemma 13 (and the representation (38), subsequently) hold for several problems in the literature. For instance the following class of rational eigenvalue problems are often encountered in practice [2]. Let $A$, $B$ and $C_{i}$, $1\leq i\leq k$, be square matrices such that $B$ is invertible, and let

[TABLE]

where $\sigma_{i}$ are given poles.

These results hold also in the context of certain modifications of the symmetric eigenvalue problem [14] and  [33]. Let $A,B\in\mathbb{C}^{n\times n}$ be such that $B$ is invertible,

and let $s:\mathbb{C}\rightarrow\mathbb{C}$ be a function such that

$s(\lambda)\rightarrow C$, $C$ constant, as $\left|\lambda\right|\rightarrow\infty$. Let

[TABLE]

The condition on the limit behaviour of $s(\lambda)$ holds in particular for [14, Example 2].

Remark 16** (A counterexample).**

Although the above theory shows that $R_{\Gamma}$ does vanish in many situations similar to the linear case, it is not always zero, as can be seen from the example

[TABLE]

where $f$ is a scalar analytic function. Then,

[TABLE]

and when we select $\Gamma$ such that encloses the two eigenvalues $\lambda_{1}=1$ and $\lambda_{2}=2$, we have

[TABLE]

5 Numerical simulations

5.1 Rational eigenvalue problem

We illustrate several properties of the results with numerical simulations111The simulations are publicly available online: http://www.math.kth.se/~eliasj/src/qnewton/. We consider the problem ’loaded string’ from the NLEVP collection [2]. The problem is of the form

[TABLE]

where $A,B,C\in\mathbb{R}^{n\times n}$, and $B$ is invertible. We set $n=20$. All the eigenvalues are real positive, and to make the spectrum more clustered in the left end, we multiply the original coeffient matrix $C$ by $n$. The spectrum of the problem is shown in Figure 1.

We first place the initial value $(\mu_{0},x_{0})$ close to the right-most eigenvalue $\lambda\approx 5170$, such that $\mu_{0}=\lambda+5.0$ and $x_{0}=v+a\cdot[1\,\ldots\,1]^{T}$, where $v$ is the exact eigenvector corresponding to $\lambda$ and $a$ is a scalar. We set $c=x_{0}$ and $\sigma=\mu_{0}$. Since we are here mainly concerned with convergence properties, we use computation with high precision arithmetic with a sufficiently high precision, such that round-off errors are not influencing the figures. As shown in Figure 2, the performance of the QN1 is found to be sensitive to the distance of $x_{0}$ from $v$. As expected from Corollary 8, the convergence curves of the Quasi-Newton 2 and the residual inverse iteration are very close to each other. Figure 3 shows the numerically estimated convergence factors and also the a priori convergence factor estimates. For the QN1 the a priori converge factor equals $\rho(A_{1})$, i.e., the spectral radius of the matrix $A_{1}$ given in Theorem 6. For the QN2 and QN3, i.e., the residual inverse iteration, it equals $\rho(B)$, where $B$ is given in

Corollary 8. The estimated convergence factor $\rho_{k}$ at iteration $k$ is computed by

[TABLE]

where $w_{k}=[\begin{smallmatrix}v_{k}\\ \mu_{k}\end{smallmatrix}]$ and $w_{*}$ denotes the exact solution $w_{*}=[\begin{smallmatrix}v\\ \lambda\end{smallmatrix}]$.

Then, we consider the initial value $(\mu_{0},x_{0})$ close to a left-end eigenvalue $\lambda\approx 9.07$, such that again $\mu_{0}=\lambda+5.0$, $x_{0}=v+a\cdot[1\,\ldots\,1]^{T}$ ($a>0$), and $c=x_{0}$ and $\sigma=\mu_{0}$. As can be expected from the bound given in 11, the convergence is now slower since $\lambda$ is in a cluster of eigenvalues. This is depicted by Figure 4. Again, the performance of the QN1 is found to be sensitive to the distance of $x_{0}$ from $v$.

As shown in Figure 4, the convergence is now slower than for the right-most eigenvalue, as is expected from the clustering of the spectrum depicted in Figure 1,

5.2 Quadratic eigenvalue problem

In this section we illustrate the influence of the eigenvalue clustering on the convergence factor. The bounds for the convergence factor illustrating this effect were derived in Section 3 and can be clearly identified from the following example. Consider the quadratic eigenvalue problem

[TABLE]

where $A_{1},A_{2}\in\mathbb{C}^{10\times 10}$ are diagonal matrices. The set eigenvalues of this problem is the uninion of eigenvalues $A_{1}$ and $A_{2}$. We choose the diagonal elements of $A_{1}$ and $A_{2}$ such that eigenvalues of $M(\lambda)$ are $0.1$ and $19$ equally distributed points on the circle of radius $r$, as illustrated in Figure 5 for $r=0.5$.

We first place the initial eigenpair approximation $(\mu_{0},x_{0})$ to the origin, such that $\mu_{0}=0$ and $x_{0}=v+0.1\cdot[1\,\ldots\,1]^{T}$, where $v$ is the exact eigenvector corresponding to the eigevalue 0.1. We set $c=x_{0}$ and $\sigma=\mu_{0}$. Figure 6 shows the convergence of the methods and also the numerically estimated convergence factors for the three first quasi-Newton methods. QN1 is found again to have the slowest convergence, and the a priori convergence factor estimates computed using Theorem 6 and Corollary 8 are again found to be sharp.

Let $\rho_{r}$ be the spectral radius of the iteration matrix $B$ given in Corollary 8, i.e., the convergence factor for QN2 and QN3 (the residual inverse iteration). In Figure 7 we illustrate how $\rho_{r}$ behaves as a function of $r$, the radius of the circle. The value of $\rho_{r}$ is computed for $10$ different values of $r$ varying from $10^{-1/2}$ to $10^{5}$. As expected from Corollary 11 we observe that $\rho_{r}\sim 1/r$.

5.3 A large scale problem

We now consider a large–scale NEP which arises in the study of waves traveling in a periodic medium [16, 36]. More precisely, the waveguide eigenvalue problem, without the Cayley transformation, associated to the waveguide described in [16, Section 5.2] is considered. The problem is formulated as

[TABLE]

We choose the discretization parameters $n_{x}=200$ and $n_{z}=201$, which means that the size of the NEP is $n=n_{x}n_{z}+2n_{z}=40602$. The matrix $C_{2}^{T}$ and the second degree polynomials $Q(\lambda)$ and $C_{1}(\lambda)$ are sparse, and the matrix $P(\lambda)$ is dense and it is defined by nonlinear functions of $\lambda$ involving square roots of polynomials. The matrix–vector product $P(\lambda)w$ is efficiently computed using two Fast Fourier Transforms (FFTs) and a multiplication with a diagonal matrix. The linear systems involving the matrix $M(\sigma)$ can be solved by precomputing a Schur complement. See [16] for a full description of the problem.

We compare QN1, QN2 and QN3 (residual inverse iteration) for approximating a specific eigenpair. The pair $(\lambda,v)$ denotes the accurate approximation of the wanted eigenpair. The shift $\sigma$ is selected close to the wanted eigenvalue $\lambda$, more precisely, $|\sigma-\lambda|\approx 0.42$ (see Figure 9), and the initial guess $x_{0}$ is selected such that $\|x_{0}-v\|\approx 10^{-5}$. The error is computed as the absolute value of the distance between the wanted eigenvalue and the current eigenvalue approximation, namely $|\mu_{k}-\lambda|$. If an initial guess of the eigenvector is provided, all the methods present similar convergence rate and QN2 and QN3 are slightly faster then QN1 (see Figure 9.)

The block $(1,2)$ of the matrix $A_{1}$ defined in (23) has norm approximatively $10^{-3}$. This may suggest that the spectral radius of this matrix is close to the spectral radius of $B$ given in (28). In particular the convergence rate of QN1 is expected to be close to the convergence rate of QN2 and QN3. This is consistent with the numerical simulation. With a random initial guess $x_{0}$ of the eigenvector, QN1 does not converge whereas QN2 and QN3 still converge with the same convergence rate but with they require more iterations since the initial error is larger. See Figure 9. This can be justified observing that, in this case, the block $(1,2)$ of the matrix $A_{1}$ defined in (23) has norm approximatively $10^{-1}$. Therefore, the convergence factor of QN1 is expected to be significantly different from QN2 and QN3 by using the previous reasoning.

6 Conclusions

We have here presented four iterative methods and showed how they can be analyzed with techniques of quasi-Newton methods, and Keldysh’s theorem. We have also illustrated how two well-established methods can be interpreted as quasi-Newton methods.

A secondary conclusion that we wish to stress here is that, in general, methods of this type can in an insightful way be analyzed in the framework of quasi-Newton methods. We have illustrated how linear, as well as higher order convergence, can be approached with quasi-Newton results. There are many results for quasi-Newton methods, which can potentially be applied to these methods, e.g., techniques which improve the convergence basin [9, 11], scaling techniques [6], adaption for non-smooth (or almost non-smooth) problems [24] and other convergence results [44].

The method corresponding to keeping only the (1,1)-block constant, i.e., (7), is to our knowledge new method in the context of NEPs. Moreover, although it has several similarities with residual inverse iteration in terms of convergence and behavior for linear problems, it can be more attractive than residual inverse iteration. Note that unlike QN2, residual inverse iteration requires the solution to a nonlinear scalar problem. For some NEPs, the action of $M(\lambda)$ is only implicitly available, e.g., in the form of a differential equation as in [28]. Therefore, the computation of a solution to the scalar nonlinear equation (19c) is computationally more demanding than computing the product $M^{\prime}(\mu)r$ which is required in (18b).

We do have a negative conclusion in this paper. We have concluded that QN1 does not appear very competitive in practice since better convergence is usually achieved from QN2, although QN1 is based on the most common quasi-Newton approach, i.e., keeping the Jacobian matrix constant. There is another important Newton-like algorithm which we have not considered in this manuscript, the Jacobi-Davidson algorithm ([30, 32]). Although the Jacobi-Davidson algorithm does have interpretations in terms of Newton’s method, as e.g., pointed out in [31, Section 6], our results here are not directly applicable. In [31, Section 6] the authors point out that the correction equation of the Jacobi-Davidson algorithm can be derived from Newton’s method on the nonlinear equation $G(x)=M(p(x))x$ where $p(x)$ is the Rayleigh functional. Note that this nonlinear equation is different from our augmented system (3). Moreover, since the solution is a manifold and the solution (eigenvector) is not isolated in the standard sense which prevents us from directly applying results for quasi-Newton methods.

Finally, we wish to specifically stress that the value of the presented results may be of interested considerably beyond the scope of the presented methods. Several of the methods presented here form the basis of other state-of-the-art such as the subspace accelerated extensions of residual inverse iteration (the nonlinear Arnoldi method [39]), preconditioned versions [8] and inner-outer-iteration constructions as in [43]. Those extensions may possibly also be interpreted in a quasi-Newton setting, although the extension requires attention beyond this manuscript.

Acknowledgements

We thank Wim Michiels (KU Leuven) for valuable discussions regarding partial fraction expansions for time-delay systems.