Nonlinearizing two-parameter eigenvalue problems

TL;DR

This paper introduces a method to convert linear two-parameter eigenvalue problems into nonlinear eigenvalue problems by eliminating one equation, enabling the use of NEP techniques and simplifying certain structured problems.

Contribution

The paper presents a novel transformation technique that converts two-parameter eigenvalue problems into NEPs, with theoretical characterization and computational strategies, including special cases and error analysis.

Findings

Transformation is effective for problems with smaller eliminated equations.

The method generalizes companion linearization for polynomial eigenvalue problems.

Error analysis shows the approach is stable with backward stable solvers.

Abstract

We investigate a technique to transform a linear two-parameter eigenvalue problem, into a nonlinear eigenvalue problem (NEP). The transformation stems from an elimination of one of the equations in the two-parameter eigenvalue problem, by considering it as a (standard) generalized eigenvalue problem. We characterize the equivalence between the original and the nonlinearized problem theoretically and show how to use the transformation computationally. Special cases of the transformation can be interpreted as a reversed companion linearization for polynomial eigenvalue problems, as well as a reversed (less known) linearization technique for certain algebraic eigenvalue problems with square-root terms. Moreover, by exploiting the structure of the NEP we present algorithm specializations for NEP methods, although the technique also allows general solution methods for NEPs to be directly applied. The nonlinearization is illustrated in examples and simulations, with focus on problems where the eliminated equation is of much smaller size than the other two-parameter eigenvalue equation. This situation arises naturally in domain decomposition techniques. A general error analysis is also carried out under the assumption that a backward stable eigenvalue solver method is used to solve the eliminated problem, leading to the conclusion that the error is benign in this situation.

Equations

$$0$$

$$1$$

$$M(\lambda )x=(A_1+\lambda A_2+g_i(\lambda )A_3)x=0.$$

$${\Delta }_0$$

$${\Delta }_1$$

$${\Delta }_2$$

$${\Delta }_{1}z$$

$${\Delta }_{2}z$$

$$C_0:=\left[\begin{array}{cc}{v}^H{A}_{2}x& v^H{A}_{3}x\\ w^H{B}_{2}y& w^H{B}_{3}y\end{array}\right],$$

$$-(B_1+\lambda B_2)y=\mu B_{3}y,$$

$$J(\lambda ,\mu ,y):=\left[\begin{array}{cc}B(\lambda ,\mu )& B_{3}y\\ c^T& 0\end{array}\right]$$

$$f(\lambda ,\mu ,y):=\left[\begin{array}{c}B(\lambda ,\mu )y\\ c^{T}y-1\end{array}\right].$$

$$B(\lambda ,{\mu }_i)u+B_3{y}_i=0,$$

$$B(\lambda ,{\mu }_i)z+B_3{y}_i\alpha =0,$$

$$c^{T}z=0.$$

$$r=\min\big{\{}|\lambda-s|:\textrm{ such that }g_{i}(s)\textrm{ is a double eigenvalue of the GEP }\eqref{eq:B_eig}\big{\}}.$$

$$0$$

$$0=A_{1}x+\lambda A_{2}x+{\lambda }^2{A}_{3}x.$$

$$\left[\begin{array}{cc}-A_1& 0\\ 0& A_3\end{array}\right]z=\lambda \left[\begin{array}{cc}{A}_2& A_3\\ A_3& 0\end{array}\right]z,$$

$$0$$

$$0=(\mu +a)(\mu +d)-(c+\lambda f)(b+\lambda e).$$

$$\mu =g_{\pm }(\lambda )=-\frac{a+d}{2}\pm \sqrt{\frac{(a+d{)}^2}{4}-ad+(b+\lambda e)(c+\lambda f)}.$$

$$(A_1+\lambda A_2+\sqrt{p(\lambda )}{A}_3)x,$$

$$y:=\left[\begin{array}{c}\sqrt{b+\lambda e}\\ -\sqrt{c+\lambda f}\end{array}\right].$$

$$\left[\begin{array}{cc}{A}_1& -b{A}_3\\ -c{A}_3& A_1\end{array}\right]z=\lambda \left[\begin{array}{cc}-A_2& e{A}_3\\ f{A}_3& -A_2\end{array}\right]z.$$

$$x_{k+1}$$

$${\alpha }_k^{-1}$$

$$\left[\begin{array}{c}{y}_i^{(k)}(\lambda )\\ g_i^{(k)}(\lambda )\end{array}\right]=J(\lambda ,{\mu }_i,y_i{)}^{-1}\left[\begin{array}{c}-b_k\\ 0\end{array}\right],$$

$$b_k=B_2{y}_i^{(k-1)}(\lambda )+\sum _{j=1}^{k-1}\left(\genfrac{}{}{0ex}{}{k}{j}\right)g_i^{(k-j)}(\lambda )B_3{y}_i^{(j)}(\lambda ).$$

$$0=\frac{d^k}{d{\lambda }^k}\left[\begin{array}{c}{B}_1{y}_i(\lambda )\\ c^T{y}_i(\lambda )-1\end{array}\right]+\frac{d^k}{d{\lambda }^k}\left[\begin{array}{c}\lambda B_2{y}_i(\lambda )\\ 0\end{array}\right]+\frac{d^k}{d{\lambda }^k}\left[\begin{array}{c}{g}_i(\lambda )B_3{y}_i(\lambda )\\ 0\end{array}\right].$$

Full text

1 Introduction

This paper concerns the two-parameter eigenvalue problem: Determine non-trivial quadruplets $(\lambda,x,\mu,y)\in\mathbb{C}\times\mathbb{C}^{n}\times\mathbb{C}\times\mathbb{C}^{m}$ such that

[TABLE]

where $A_{1},A_{2},A_{3}\in\mathbb{C}^{n\times n}$, and $B_{1},B_{2},B_{3}\in\mathbb{C}^{m\times m}$. We denote the corresponding functions $A(\lambda,\mu):=A_{1}+\lambda A_{2}+\mu A_{3}$ and $B(\lambda,\mu):=B_{1}+\lambda B_{2}+\mu B_{3}$. This problem has been extensively studied in the literature, see, e.g., the fundamental work of Atkinson [2], and the summary of recent developments below. We assume that $m\ll n$ and that $A_{1}$, $A_{2}$ and $A_{3}$ are large and sparse matrices, although several theoretical contributions of this paper are valid without this assumption.

The main idea of our approach can be described as follows. We view (1b) as a parameterized generalized linear eigenvalue problem, where $\lambda$ is the parameter. Due to perturbation theory for eigenvalue problems, there is a family of continuous functions $\{g_{i}(\lambda)\}$, defined by the eigenvalues of (1b), where $\mu$ is the eigenvalue, of a generalized eigenvalue problem (GEP). More formally, for a fixed value of $\lambda$ the functions $g_{i}(\lambda)$ and $y_{i}(\lambda)$ can be defined, as the solution to

[TABLE]

for a given vector $c\in\mathbb{C}^{m}$. We explicitly introduced the normalization condition (2b), to uniquely define a corresponding eigenvector. The condition (2b) is not a restriction of generality except for the rare situation that the eigenvector is orthogonal to $c$. We prefer this condition over the standard Euclidean normalization, since the right-hand side of (2b) is an analytic function.

By insertion of $\mu=g_{i}(\lambda)$ into (1a), we see that a solution to (1) will satisfy

[TABLE]

Note that we have now eliminated $\mu$ and (1b), at the cost of the introduction of a nonlinear function into the eigenvalue problem. The problem $M(\lambda)x=0$ is called a nonlinear eigenvalue problem (NEP). In our setting it is rather a family of NEPs, since we have a different nonlinearity for each function $g_{1},\ldots,g_{m}$. The study of NEPs is a mature field within numerical linear algebra, and there are considerable theoretical results, as well as algorithms and software for NEPs.

The main contributions of this paper consist of a theoretical characterization of the elimination procedure (Section 2), analysis of structured perturbations corresponding to the elimination (Section 4), as well as new algorithms for (1) based on NEP-algorithms (Section 3). We provide software for the simulations, both for MATLAB and for Julia [5]. The Julia software is implemented using the data structures of the NEP-PACK software package [18], including adaption of theory for how to compute derivatives and projections. This provides new ways to solve (1), using the large number of NEP-solvers available in NEP-PACK. Some contributions are also converse, i.e., we provide insight to NEPs based on the equivalence with two-parameter eigenvalue problems. For instance, in Sections 2.2–2.3 we show how to transform certain NEPs with square-root nonlinearities to two-parameter eigenvalue problems. This in turn (using the operator determinants described below) allows us to transform the problem to a standard generalized eigenvalue problem, similar to companion linearization techniques for polynomial and rational eigenvalue problems.

We now summarize the NEP-results relevant for our approach. For a broad overview see the summary papers [34, 26, 46, 9], as well as the benchmark collection [3] and software packages with NEP-solvers [33, 11, 12, 18]). There is considerable theoretical works available for the NEP, in particular for polynomial eigenvalue problems. Techniques to transform polynomial NEPs to standard eigenvalue problems (known as linearization) have been completely characterized in a number of works, e.g., [24, 25] and [29]. We relate our approach to this type of linearization in Section 2.2. In our derivation, we make explicit use of the implicit function theorem applied to the NEP. This has been done in the context of sensitivity analysis, leading to eigenvector free formulas for conditioning [1]. There are a number of algorithms available for NEPs, of which many seem to be applicable to (3). More specifically, we characterize the specialization of residual inverse iteration [30], which forms the basis of more recent methods such as the nonlinear Arnoldi method [45]. We also show how the infinite Arnoldi method [20] can be adapted to (3).

In Section 5.2 we illustrate how two-parameter eigenvalue problems of this type can arise by the separation of domains of a partial-differential equation (PDE). The domains are decoupling in a way that the discretization leads to a two-parameter eigenvalue problem. In this context, the elimination corresponds to an elimination of one of the domains. The elimination of an outer domain, in a way that directly leads to NEPs, by introduction of artificial boundary conditions is the origin of several standard NEPs in the literature, e.g., [40] and the electromagnetic cavity model in [44].

Relevant results for two-parameter eigenvalues can be summarized as follows. Many results for two-parameter eigenvalue problems are phrased in the more general setting of multiparameter eigenvalue problems. There are a number of recent efficient algorithms available, e.g., based on the Jacobi–Davidson approach [15, 16], including subspace methods in [14]. A number of generalizations of inverse iteration are derived in [32]. Our approach is based on an eigenvalue parameterization viewpoint. Eigenvalue parameterization and continuation techniques (but with an additional parameter) have been studied, e.g., in [31].

One of the most fundamental properties of two-parameter eigenvalue problems is the fact that solutions are given by the solution to a larger linear (generalized) eigenvalue problem. This is also often used in the numerical algorithms mentioned above, and to our knowledge first proposed as a numerical method in [36]. More precisely, we associate with (1) the operator determinants

[TABLE]

where $\otimes$ denotes the Kronecker product. The solutions to (1) are (under certain assumptions) equivalent to the solutions to the two generalized eigenvalue problems

[TABLE]

where $z=y\otimes x$. In practice, the application of a general purpose eigenvalue solver on one of the GEPs in (7) yields an accurate solution for small systems. The equivalence between (7) and (1) holds under non-singularity assumption; in particular the problem is singular if $A_{3}$ and $B_{3}$ both are singular; or $A_{2}$ and $B_{2}$ both are singular. See [2] for a precise characterization, and [21, 17] for more recent formulations.

The following matrix is often used in theory for eigenvalue multiplicity and eigenvalue conditioning, and will be needed throughout the paper. We denote

[TABLE]

where $v$ and $w$ are left eigenvectors associated with (1a) and (1b) respectively. In particular, for an (algebraically) simple eigenvalue of the two-parameter eigenvalue problem (1), the matrix $C_{0}$ is nonsingular; see [21, Lemma 3], [15, Lemma 1.1], and [17, Lemma 1]. For a simple eigenvalue, the normwise condition number for the two-parameter eigenvalue problem is expressed as a special induced matrix norm of $C_{0}^{-1}$, see [17, Section 4].

2 Nonlinearization

2.1 Existence and equivalence

The elimination of the $B$-equation in the two-parameter eigenvalue problem can be explicitly characterized as we describe next. We show how the existence of a nonlinearization can be explicitly related to the Jordan structure of the (parametrized) GEP

[TABLE]

which we can also use in practice for the computation of $\mu=g_{i}(\lambda)$ for a given $\lambda$. The existence of analytic functions is formalized in the following lemma. The invertibility assumption in the lemma will be further characterized in Theorem 3.

Lemma 1** (Existence of implicit functions).**

Let $\lambda\in\mathbb{C}$ be given and assume that $(\mu,y)$ is such that (1b) is satisfied with $y$ normalized as $c^{T}y=1$. Moreover, assume that

[TABLE]

is nonsingular. Then, there exist functions $g_{i}:\mathbb{C}\rightarrow\mathbb{C}$ and $y_{i}:\mathbb{C}\rightarrow\mathbb{C}^{m}$ such that

• •

$g_{i}$* and $y_{i}$ are analytic in $\lambda$,*

• •

$g_{i}$* and $y_{i}$ satisfy (2) in a neighborhood of $\lambda$, and*

• •

$\mu=g_{i}(\lambda)$* and $y=y_{i}(\lambda)$.*

Proof.

The proof is based on the complex implicit function theorem. Consider the analytic function $f:\mathbb{C}^{m+2}\rightarrow\mathbb{C}^{m+1}$ given by

[TABLE]

Then $J$ as in (10) is the partial Jacobian of $f$ with respect to the variables $y$ and $\mu$, i.e., $J=\partial f/\partial(y,\mu)$. Since $f(\lambda,\mu_{i},y_{i})=0$ and $J(\lambda,\mu_{i},y_{i})$ is nonsingular, the existence of the desired functions, analytic in the point $\lambda$, follows from the complex implicit function theorem [8, Theorem I.7.6]. ∎

Under the same conditions that the implicit functions exist we have the following equivalence between the solutions to the NEP (3) and the solutions to the two-parameter eigenvalue problem (1).

Theorem 2** (Equivalence).**

Suppose the quadruplet $(\lambda,x,\mu,y)\in\mathbb{C}\times\mathbb{C}^{n}\times\mathbb{C}\times\mathbb{C}^{m}$ is such that $c^{T}y=1$ and $J(\lambda,\mu,y)$ defined in (10) is nonsingular. Then, $(\lambda,x,\mu,y)$ is a solution to (1) if and only if $(\lambda,x)$ is a solution to the NEP (3) for one pair of functions $(g(\lambda),y(\lambda))=(\mu,y)$ which satisfies (2).

Proof.

To prove the forward implication direction suppose $(\lambda,x,\mu,y)$ is a solution to (1). From Lemma 1, there are functions $g(\lambda)$ and $y(\lambda)$ such that $g(\lambda)=\mu$ and $y(\lambda)=y$. Therefore, (3) is satisfied for that pair $(g(\lambda),y(\lambda))$.

To prove the backward implication direction suppose $(\lambda,x)$ is a solution to (3) for a given pair $(g(\lambda),y(\lambda))$. Then $(\lambda,x,\mu,y)=(\lambda,x,g(\lambda),y(\lambda))$ is a solution to (1). More precisely, (1a) is satisfied since (3) and (1b) is satisfied due to (2). ∎

The situation that the Jacobian matrix in (10) is singular is a non-generic situation. It turns out, as we show in the following theorem, that it is singular (essentially) if and only if the GEP (9) has a Jordan chain of length two or more. Therefore, our technique in general works if a solution to the two parameter eigenvalue problem (1) corresponds to a simple eigenvalue $\mu$ of the GEP (9).

Theorem 3** (Singularity and Jordan structure).**

Let $\lambda\in\mathbb{C}$ be given and assume that $c$ is not orthogonal to any eigenvector of the GEP (9). Moreover, assume that $(\mu_{i},y_{i})$ is a solution to the GEP with $y_{i}$ normalized as $c^{T}y_{i}=1$. Then $J(\lambda,\mu_{i},y_{i})$ is singular if and only if there exists a vector $u\in\mathbb{C}^{m}$ such that

[TABLE]

i.e., there exists a Jordan chain of length at least two corresponding to the GEP (9) and eigenpair $(\mu_{i},y_{i})$.

Proof.

We start by proving that singularity implies the existence of a Jordan chain. Assume that $J(\lambda,\mu_{i},y_{i})$ is singular. Then there exists a non-trivial vector $\begin{bmatrix}z&\alpha\end{bmatrix}^{T}\in\mathbb{C}^{m+1}$ such that $J(\lambda,\mu_{i},y_{i})\begin{bmatrix}z&\alpha\end{bmatrix}^{T}=0$. The first row gives

[TABLE]

and the second row gives

[TABLE]

The cases $\alpha=0$ and $\alpha\neq 0$ are investigated separately. Assume that $\alpha=0$, then $z\neq 0$ and thus (12) implies that $z$ is an eigenvector to the GEP. However, (13) gives a contradiction. If $\alpha\neq 0$ then (12) gives that there exists a Jordan chain of length at least two, with $u=z/\alpha$. Hence, singularity implies the existence of a Jordan chain.

To prove the converse we assume that there exists a Jordan chain of length at least two. Let $z:=u-(c^{T}u)y_{i}$. Note that from construction of $z$ (13) holds. Moreover, from the Jordan chain we know that (12) holds for the constructed $z$ with $\alpha=1$. Hence, the vector $\begin{bmatrix}z&1\end{bmatrix}^{T}$ is a non-trivial vector in the null-space of $J(\lambda,\mu_{i},y_{i})$. Thus the existence of a Jordan chain implies singularity. ∎

From a practical point of view we know that if we compute simple eigenvalues of the GEP (9) such that $c$ is not orthogonal to the corresponding eigenvector, then the Jacobian is nonsingular. Hence, the nonlinearization exists. The result is formalized in the following corollary to Theorem 3.

Corollary 4**.**

Let $\lambda\in\mathbb{C}$ be given. Assume that $\mu_{i}\in\mathbb{C}$ is a simple eigenvalue of the GEP (9), and that $y_{i}$ is a corresponding right eigenvector normalized as $c^{T}y_{i}=1$. Then $J(\lambda,\mu_{i},y_{i})$ is nonsingular, where $J$ is defined in (10).

Many algorithms for NEPs depend on analyticity in a target domain. From the above theory we directly conclude the following result for the convergence radius of the implicitly constructed analytic function.

Corollary 5** (Convergence radius).**

Let $\lambda\in\mathbb{C}$ be given and assume that $(\mu_{i},y_{i})$ is a solution to the GEP (9) with $y_{i}$ normalized as $c^{T}y_{i}=1$. Moreover, assume that $J(\lambda,\mu_{i},y_{i})$ is invertible, where $J$ is defined in (10). The functions $g_{i}$ and $y_{i}$ in Theorem 1 can be chosen such that they are analytic in the open disk with radius $r$ centered in $\lambda$, where $r$ is defined by

[TABLE]

Proof.

By the definition of $r$, Theorem 3 ensures that $J(\lambda,g_{i}(\lambda),y_{i})$ is nonsingular in all points in the open disk. Application of Theorem 1 to all those points establishes the result. ∎

As discussed above, the choice of solution to the GEP (9) corresponds to the choice of function $g_{i}(\lambda)$. We note that from Corollary 4 it is clear that the existence of the nonlinearization only relies on that the chosen eigenvalue, of the GEP, is simple. Similarly from Corollary 5 it is clear that the convergence radius is dependent only on the specific function $g_{i}$ in consideration. Hence, the existence, and the convergence radius, of the NEP (3) only depends on the behavior of $g_{i}$ and not the complete eigenstructure of the GEP.

2.2 Nonlinearizations leading to quadratic eigenvalue problems

We first illustrate the theory in the previous section with an implicitly defined function which can be derived explicitly. Consider the two-parameter eigenvalue problem

[TABLE]

for general matrices $A_{1}$, $A_{2}$ and $A_{3}$. The second row in (14b) implies that the elements in the vector $y^{T}=\begin{bmatrix}y_{1}&y_{2}\end{bmatrix}$ are related by $y_{2}=\lambda y_{1}$. The first row in (14b) becomes $\lambda^{2}y_{1}-\mu y_{1}=0$. Hence, since $y_{1}\neq 0$, we have $\mu=\lambda^{2}$ and (14a) becomes

[TABLE]

This problem is commonly known as the quadratic eigenvalue problem, which has been extensively studied in the literature [41]. The example shows that the two-parameter eigenvalue problem (14) can be nonlinearized to a quadratic eigenvalue problem. Moreover, the determinant operator equation (7a) leads to the equation

[TABLE]

which is a particular companion linearization of (15). (It is in fact a symmetry preserving linearization [41, Section 3.4].) Many of the linearizations of polynomial eigenvalue problems given in [25] can be obtained in a similar fashion. Since, the second equation (1b) can be expressed as $\det(B(\lambda,\mu))=0$, which is a bi-variate polynomial, this example is consistent with the bivariate viewpoint of companion linearizations in [29]. Some higher-degree polynomials can be constructed analogously to above, e.g., the polynomial eigenvalue problem $A_{1}+\lambda A_{2}+\lambda^{m}A_{3}$. However, the general higher-degree polynomial eigenvalue problem does not seem to fit into the class of two-parameter eigenvalue problems.

2.3 Nonlinearization leading to algebraic functions

The previous example can be modified in a way that it leads to algebraic functions, which is also the generic situation. Nontrivial solutions to (1b) satisfy $\det(B(\lambda,\mu))=0$, which is a bivariate polynomial. Therefore, the functions $g_{i}(\lambda)$ are roots of a polynomial, where the coefficients are polynomials in $\lambda$, i.e., $g_{i}$ are algebraic functions. The generic situation can be seen from the case where $m=2$:

[TABLE]

We obtain that $\mu$ is the root of a polynomial, where the coefficients depend on $\lambda$, i.e.,

[TABLE]

The explicit solutions to this quadratic equation are given by

[TABLE]

We see by insertion of $\mu=g_{\pm}$ into (16a) that the nonlinearization of (16) is a NEP with an algebraic nonlinearity. The function $g_{+}$ is illustrated in Figure 1.

Several general conclusions can be made from this example. Note that the variables $a,b,c,d,e,f$ can be used for fitting of any function $\sqrt{p(\lambda)}$ where $p$ is a polynomial of degree two. Therefore, we can now reverse the nonlinearization, and for the trivial case $a=d=0$ we directly obtain the following characterization.

Lemma 6** (Two-parameterization of an algebraic NEP).**

Suppose $p(\lambda)=(b+\lambda e)(c+\lambda f)$ is given, and let $a=d=0$. If $(\lambda,x)$ is a solution to the NEP

[TABLE]

then $(\lambda,x,\mu,y)$ satisfies the two-parameter eigenvalue problem equation (16) with $\mu:=\sqrt{p(\lambda)}$ and

[TABLE]

A further consequence of the lemma is that problems of the type (17) can be linearized to a standard GEP using the determinant operators (7). More precisely, the combination of Lemma 6 and (7) shows that (17) can be solved by computing solutions to

[TABLE]

The fact that algebraic NEPs can be linearized was already pointed out in the conference presentation [35], for a specific case using techniques not involving two-parameter eigenvalue problems.

Also note that the functions $g_{i}(\lambda)$ have branch-point singularities. This is the generic situation and we can therefore never expect that the nonlinearizations are entire functions in general. The singularities restrict the performance of many methods, as we will see in the simulations. The implications of singularities in practice is well-known in quantum chemistry, where parameterized eigenvalue problems is a fundamental tool and the singularities are referred to as intruder states [10, Chapter 14]. In that context, methods for computing the closest singularity (which limits the performance of the method) are given in [19, 22].

3 Algorithm specializations

3.1 Derivative based algorithms

Many NEP-algorithms are based on derivatives of $M$. We will now illustrate how to efficiently and reliably access the derivatives of the NEP stemming from a nonlinearization of a two-parameter eigenvalue problem. As a representative first situation we consider the augmented Newton method, see [34, 43]. It can be derived by an elimination of the correction equation in Newtons method, and leads to separate eigenvalue and eigenvector update formulas expressed as

[TABLE]

and $\lambda_{k+1}=\lambda_{k}-\alpha_{k}$, where $d\in\mathbb{C}^{n}$ is a normalization vector. In an implementation, one takes advantage of the fact that the same linear system appears twice, and only needs to be computed once. The iteration has appeared in many variations with different names, e.g., inverse iteration [37] and Newtons method [42].

In order to apply (18) we clearly need the derivative of $M$ defined in (3), which can be obtained directly if we can compute the derivative of the implicitly defined function $g_{i}$. Note that the functions $g_{i}(\lambda)$ (as well as the auxiliary vector $y_{i}(\lambda)$) can be evaluated by solving the GEP (9), and normalizing according to $c^{T}y_{i}=1$. Since the functions are analytic in general, their respective derivatives exist. They can be computed according to the following result, which gives a recursion that can compute the $k$th derivative by solving $k$ linear systems of dimension $(m+1)\times(m+1)$. The adaption of the theorem and (18) into an algorithm results in Algorithm 1.

Theorem 7** (Explicit recursive form for derivatives).**

Let $\lambda\in\mathbb{C}$ be given and assume that $(\mu_{i},y_{i})$ is a solution to the GEP (9) with $y_{i}$ normalized as $c^{T}y_{i}=1$. Moreover, assume that $J(\lambda,\mu_{i},y_{i})$ is invertible, where $J$ is defined in (10). Let $g_{i}$ and $y_{i}$ be the functions defined in Lemma 1, then the $k$th derivative, $k=1,2,\dots$, of $g_{i}$ and $y_{i}$ are given by

[TABLE]

where

[TABLE]

Proof.

We again consider the analytic function $f$ given by (11). By Lemma 1 we know that $g_{i}$ and $y_{i}$ are analytic around $\lambda$, and that $f(\lambda,g_{i}(\lambda),y_{i}(\lambda))=0$ in a neighborhood of $\lambda$. Taking the $k$th implicit derivative with respect to $\lambda$ gives

[TABLE]

The two first terms are found directly as

[TABLE]

and

[TABLE]

The third term can be calculated, by using Leibniz derivation rule for products, to be

[TABLE]

We emphasize the recursion: All derivatives up to order $k-1$ can be considered known since these do not depend on the higher derivatives. Collecting the known terms in the right-hand-side gives the result. ∎

Remark 8**.**

As a special case of Theorem 7, for $k=1$, we find that $g_{i}^{\prime}(\lambda)=-\frac{w_{i}^{H}B_{2}y_{i}}{w_{i}^{H}B_{3}y_{i}}$ where $w_{i}$ is the corresponding left eigenvector to the eigenpair $(\mu_{i},y_{i})$. It follows from multiplying the first block-row of equation system (19) from the left with $w_{i}^{H}$. The result is a special case of well known perturbation analysis for generalized eigenvalue problems [13, Theorem 2.5]. In our case $g^{\prime}(\lambda)$ is the perturbation of the eigenvalue $\mu$ with respect to $\lambda$ in the GEP (9). More precisely, a perturbation of the matrix $-(B_{1}+\lambda B_{2})$ with the structured perturbation $\varepsilon B_{2}$.

Specifically, the closed form of $g_{i}^{\prime}(\lambda)$ means that the derivative of the NEP (3) can be written in closed form, as

[TABLE]

For methods only requiring the first derivative of $M(\lambda)$, the above expression can be used instead of (19). However, that requires the computations of the left eigenvector of the GEP. We will need the expression for theoretical purposes in Section 4.

The family of methods in [20, 6, 28] (flavors of the infinite Arnoldi method) also requires derivative information. These methods require computation of quantities such as

[TABLE]

where $x_{1},\dots,x_{p}$ are given vectors. The computation requires higher derivatives of $g_{i}$. However, $\sigma$ is unchanged throughout the iteration and therefore the matrix in the linear system for derivative computation (19) is unchanged. Hence, all needed derivatives can be computed by solving an additional linear system. If $m\ll n$, this will in general not be computationally demanding. We also note that these fixed-shift methods choose a branch $g_{i}$ in the initial solution of the GEP (9), and then stay on that branch. Hence, convergence properties will depend on the convergence radius of that function $g_{i}$, as mentioned in the end of Section 2.1.

3.2 Projection methods

Many NEP-algorithms require the computation of a projected problem

[TABLE]

where $V,W\in\mathbb{C}^{n\times p}$ are orthogonal matrices. The problem (20) is again a NEP, but of smaller size. This can be viewed as a Petrov–Galerkin projection of the spaces spanned by the columns of $V$ and $W$. The projection is sometimes called subspace acceleration (or the nonlinear Rayleigh–Ritz procedure), since it is often used to improve properties of a more basic algorithm, e.g., the nonlinear Arnoldi method [45], Jacobi–Davidson methods [7, 4], block preconditioned harmonic projection methods [47], the infinite Lanczos method [27], and many more.

In order to give access to these methods, we need to provide a way to solve (20) for our nonlinearized problem. Fortunately, the projected problem stemming from the nonlinearized two-parameter eigenvalue problem, i.e.,

[TABLE]

has a structure which suggests straightforward methods for the projected problem. This is because the projected NEP has the same structure as the nonlinearized two-parameter eigenvalue problem, and can therefore be lifted back to a two-parameter eigenvalue problem, but now of much smaller size. We can then use general methods for two-parameter eigenvalue problems. This is directly observed from the fact that (21) is the nonlinearization of a two-parameter eigenvalue problem with projected $A$-matrices. It is made more precise in the following result.

Corollary 9** (Projected nonlinearized problem).**

Suppose the quadruplet $(\lambda,z,\mu,y)\in\mathbb{C}\times\mathbb{C}^{p}\times\times\mathbb{C}\times\mathbb{C}^{m}$ is such that $c^{T}y=1$ and suppose $J(\lambda,\mu,y)$ with $J$ as defined in (10) is nonsingular. Then, $(\lambda,z,\mu,y)$ is a solution to the the two-parameter eigenvalue problem

[TABLE]

if and only if $(\lambda,z)$ is a solution to (21) for one pair of functions $(g(\lambda),y(\lambda))=(\mu,y)$ which satisfies (2).

Proof.

This follows directly from the application of Theorem 2 on the projected problem (22) and the NEP (21). ∎

If the projection space is small $p\ll n$, and $m\ll n$, we may even solve the two-parameter eigenvalue problem using the operator determinant eigenvalue equations (7) or [15, Algorithm 2.3].

The situation $p=1$ implies that the projected problem is a scalar problem, and reduces to the so-called Rayleigh functional. There are several methods based on the Rayleigh functional, e.g., residual inverse iteration [30], and variational principle based approaches such as [39] and references therein. The fact that the projected problem is scalar and linear allows us to eliminate $\mu$, and we find that $\lambda$ is a solution to the generalized eigenvalue problem. The following corollary specifies the formulas more precisely, and the adaption of the result into the residual inverse iteration is given in Algorithm 2.

Corollary 10**.**

The solution to the projected NEP (21) with $p=1$ is given by $\lambda$, $\mu\in\mathbb{C}$ and $y\in\mathbb{C}^{m}$, where $(\lambda,y)$ is a solution to the GEP

[TABLE]

and $\mu$ is given by

[TABLE]

Proof.

This is derived from a special case of Corollary 9 where $p=1$. The relation (22a) with $W=w$ and $V=v$ can be solved for $\mu$ resulting in the relation (24). By inserting this relation into (22b) we obtain the GEP (23). ∎

4 Conditioning and accuracy

In order to characterize when the elimination procedure works well, we now analyze how the technique behaves subject to perturbations. As a consequence of this we can directly conclude how backward stable computation of $g$ influences the accuracy (Section 4.2).111For notational convenience the $i$ index on $g_{i}$ is dropped in this section.

4.1 Conditioning as a nonlinear eigenvalue problem

Standard results for the condition number of NEPs can be used to analyze perturbations with respect to the $A$-matrices. More precisely, if we define

[TABLE]

where $\alpha_{j}$ are scalars for $j=1,2,3$, and $\Delta\lambda$ is such that

[TABLE]

then we know (see, e.g., [1]) that

[TABLE]

where $v,x$ are the corresponding left and right eigenvectors. In the following we will establish how this formula is modified when we also consider perturbations in the $B$-matrices. Note that this implies that the function $g$ is also perturbed and we cannot directly use the standard result. We therefore define the condition number

[TABLE]

where $\beta_{j}$ are scalars for $j=1,2,3$, and $\Delta\lambda$ fulfills (25) but with a perturbed $g$, i.e., $\mu+\Delta\mu=g(\lambda+\Delta\lambda)$, such that

[TABLE]

The definitions can be used both for absolute and relative condition numbers by setting $\alpha_{j}=\beta_{j}=1$ or $\alpha_{j}=\|A_{j}\|$, $\beta_{j}=\|B_{j}\|$ for $j=1,2,3$ respectively.

As an intermediate step we first consider the perturbation of $\mu$ subject to perturbations in the $B$-matrices and fixed perturbations in $\lambda$ by analyzing

[TABLE]

where $\gamma$ is a scalar, and $\Delta\mu$ satisfies (27) for a given $\lambda$ . The following result shows that $\kappa_{g}$ can be expressed as a sum of perturbations associated with the $B$-matrices and perturbations associated with $\lambda$.

Lemma 11**.**

Let $\lambda\in\mathbb{C}$ be given and suppose $\mu=g(\lambda)$ is a simple eigenvalue of the GEP (9) with $w$ and $y$ being corresponding left and right eigenvectors. Then,

[TABLE]

where

[TABLE]

Proof.

Since $\mu$ is a simple eigenvalue of the GEP (9), the eigenvalue and eigenvector are analytic, and therefore $\Delta y=\mathcal{O}(\varepsilon)$ when all the perturbations are $\mathcal{O}(\varepsilon)$. By collecting all the higher order terms the perturbed GEP (27a) can thus be written as

[TABLE]

Multiplying with $w^{H}$ from the left, solving for $\Delta\mu$, and dividing with $\varepsilon$ gives that

[TABLE]

An upper bound is thus found as

[TABLE]

It remains to show that the bound can be attained. This follows from considering $\hat{B}=wy^{H}/\|w\|\|y\|$, and inserting

[TABLE]

into (28). ∎

Using the intermediate result we can now show that the condition number $\kappa(\lambda)$ is the sum of the standard condition number of NEPs and a term representing perturbations in $g$ generated by perturbations in the $B$-matrices, i.e., $\kappa_{g,B}(\lambda)$.

Theorem 12**.**

Let $\lambda$ be a simple eigenvalue of the NEP (3) with $v$ and $x$ being corresponding left and right eigenvectors. Moreover, for this $\lambda$, suppose $\mu=g(\lambda)$ is a simple eigenvalue of the GEP (9) with $w$ and $y$ being corresponding left and right eigenvectors. Then,

[TABLE]

where $\kappa_{A}(\lambda)$ is given by (26).

Proof.

Recall the assumptions that the NEP (3), i.e., $M$, is analytic, that $\lambda$ is a simple eigenvalue of the NEP, and that $\mu$ is a simple eigenvalue of the GEP (9). Hence, the eigenvalues and eigenvectors are analytic, and therefore $\Delta x=\mathcal{O}(\varepsilon)$ when all the perturbations are $\mathcal{O}(\varepsilon)$. By using that $g(\lambda+\Delta\lambda)=g(\lambda)+\Delta\mu$ and collecting all the higher order terms, the perturbed NEP (25) can therefore be written as

[TABLE]

Multiplying with $v^{H}$ from the left, expanding $\Delta\mu$ according to (28), solving for $\Delta\lambda$, and dividing with $\varepsilon$, gives that

[TABLE]

where $\theta_{g,B}(\lambda):=-(w^{H}\Delta B_{1}y+\lambda w^{H}\Delta B_{2}y+g(\lambda)w^{H}\Delta B_{3}y)/(w^{H}B_{3}y)$. Based on Remark 8 we observe that the denominator of (29) is equal to $\varepsilon v^{H}M^{\prime}(\lambda)x$. An upper bound is is therefore

[TABLE]

It remains to show that the bound can be attained. Similar to the proof of Lemma 11, this follows from considering $\hat{B}=wy^{H}/\|w\|\|y\|$ and $\hat{A}=vx^{H}/\|v\|\|x\|$, and inserting

[TABLE]

into (29). ∎

4.2 Backward stable computation of $g$

The nonlinearization is based on solving a GEP to evaluate the function $g(\lambda)$. We analyze the effects on the accuracy in the computed $\lambda$ when the GEP is solved numerically with a backward stable method. The analysis assumes the two triplets $(\lambda,x,v)\in\mathbb{C}\times\mathbb{C}^{n}\times\mathbb{C}^{n}$ and $(\mu,y,w)\in\mathbb{C}\times\mathbb{C}^{m}\times\mathbb{C}^{m}$ are such that $\lambda$ is a simple eigenvalue of the NEP (3), $\mu$ is a simple eigenvalue of the GEP (9), and $v,w$ and $x,y$ are corresponding left and right eigenvectors respectively.

From the assumption that the GEP (9) is solved by a backward stable method we know that $\mu$ can be characterized as the exact solution to a nearby problem. More precisely, $\mu$ solves

[TABLE]

where $C_{1}=-(B_{1}+\lambda B_{2})$, $C_{2}=B_{3}$, with perturbations, $\Delta C_{1}$ and $\Delta C_{2}$, that are proportional to the errors in our GEP solver. Specifically, there are non-negative $\beta_{1},\beta_{3}\in\mathbb{R}$ such that $\|\Delta C_{1}\|=\beta_{1}\varepsilon$ and $\|\Delta C_{2}\|=\beta_{3}\varepsilon$. Thus, the perturbation in $g$ is precisely captured by $\kappa_{g,B}(\lambda)$ from Lemma 11, with $\beta_{2}=0$ and $\beta_{1}$ and $\beta_{3}$ given above, i.e., by the specific choice of GEP solver. Hence, by application of Theorem 12 with $\alpha_{j}=0$ for $j=1,2,3$ we can conclude that the forward error in $\lambda$, induced by the inexact but backward stable computation of $g(\lambda)$ is bounded by

[TABLE]

Without loss of generality we now assume that $\|x\|=\|v\|=\|y\|=\|w\|=1$.

The upper bound (30) is related to the condition number for multiparameter eigenvalue problems as follows. As mentioned in the introduction, the condition number for the two-parameter eigenvalue problem can be directly expressed with the inverse of $C_{0}$ defined in (8). First note that our assumptions imply that $C_{0}$ is invertible.

Lemma 13**.**

Under the conditions of Theorem 12 the matrix $C_{0}$ is nonsingular, where $C_{0}$ is defined in (8).

Proof.

By using the expression for $M^{\prime}(\lambda)$ from Remark 8 we thus have

[TABLE]

Since the eigenvalues $\lambda$ and $\mu$ are simple we know that $w^{H}B_{3}y\neq 0$, and that $v^{H}M^{\prime}(\lambda)x\neq 0$. Hence, $\det(C_{0})\neq 0$. ∎

From (31) we can conclude that the bound (30) on $|\Delta\lambda|$ can be written as

[TABLE]

Moreover, for a nonsingular $C_{0}$ it is shown in [17, Theorem 6] that the condition number of the two-parameter eigenvalue is

[TABLE]

where the $\theta$-norm, i.e., $\|\cdot\|_{\theta}$, is an induced norm defined in [17, Equation (5)]. In our case we can explicitly bound the condition number by using bounds following directly from the definition of the $\theta$-norm:

[TABLE]

The parameter $\theta_{2}$ is the second component of the $\theta$-vector used in the definition of the $\theta$-norm. Hence, the bound in (32) can be further bounded by

[TABLE]

The typical choices of $\theta$ corresponding to the absolute respectively relative condition number of the two-parameter eigenvalue problem are $|\theta_{2}|=1+|\lambda|+|g(\lambda)|$ and $|\theta_{2}|=\|B_{1}\|+|\lambda|\|B_{2}\|+|g(\lambda)|\|B_{3}\|$. From the bounds in (33) we therefore conclude: The error generated by a backward stable method is benign for well conditioned two-parameter eigenvalue problems.

5 Simulations

5.1 Random example

We generate an example similar to the example in [15], but with $m\ll n$. More precisely, we let

[TABLE]

where $n=5000$ and $m=20$. The matrices $V_{A_{i}}$, $U_{A_{i}}$, $V_{B_{i}}$, $U_{B_{i}}$ have randomly normal distributed elements and $F_{i}$, $G_{i}$ are diagonal matrices with randomly normal distributed diagonal elements. The scalars $\alpha_{i}$ and $\beta_{i}$ were selected such that the eigenvalues closest to the origin were of order of magnitude one in modulus ($\alpha_{1}=\beta_{1}=1$, $\alpha_{2}=\beta_{2}=1/500$, $\alpha_{3}=\beta_{3}=1/50$). The simulations were carried out using the Julia language [5] (version 1.1.0), but implementations of the algorithms are available online for both Julia and MATLAB.222The matrices and the simulations are provided online for reproducibility: http://www.math.kth.se/~eliasj/src/nonlinearization

Since $m=20$, we in general obtain $20$ different functions $g_{1},\ldots,g_{20}$, which we order by magnitude in the origin, each corresponding to a different NEP. Some of the nonlinear functions $g_{i}$ are visualized in Figure 2. The solutions closest to the origin, for the NEPs corresponding to the functions $g_{1},g_{2},g_{3}$, are given in Figure 3.

Note that this problem is of such size that the naive approach with operator determinants (7) is not feasible, since we cannot even store them in memory on the computers we use for the simulations333The simulations were carried out using Ubuntu Linux, Intel(R) Core(TM) i7-8550U CPU @ 1.80GHz, 16 GB of RAM.

We illustrate our algorithms and compare with several other single-vector state-of-the-art algorithms in [32]. As starting values we use $\lambda_{0}=0.15+0.1i$ and $\mu_{0}=35+0.25i$, and a starting vector with an elementwise absolute error (from a nearby solution) less than $0.05$. The iteration history of Algorithm 1 is given in Figure 4. We observe an asymptotic fast convergence for Algorithm 1, which is expected since the solution point is analytic and simple. The error is measured at Step 3 in Algorithm 1 which implies that by construction, the error in the $B$-equation is (numerically) zero. This is a property of the elimination in our approach. We compare (with the same starting values) with the inverse iteration Newton approach proposed in [32]. Note that this method is designed for more general problems, and not specifically our situation where $m\ll n$ and also multiparameter nonlinear problems. Since [32, InvIter] requires several linear solves per iteration, our algorithm is faster in this case. The comparison between the two algorithms as a function of iteration is inconclusive, as can be seen in Figure 4a.

The convergence of our adaption of residual inverse iteration (Algorithm 2) initiated in the same way (except the starting vector is chosen as a vector of ones) is illustrated in Figure 5. We clearly see the expected linear convergence, since it is equivalent to residual inverse iteration for NEPs and the convergence theory in [30, Section 3–4] is directly applicable. We compare with a proposed generalization of residual inverse iteration [32, InvIter], again noting that it has a much wider applicability domain than our approach. In this case, our method has a smaller convergence factor, intuitively motivated by the fact that we solve the $B$-equation exactly.

The problem can also be solved with the tensor infinite Arnoldi method [6]. More specifically, we use the implementation of the method available in the Julia package NEP-PACK [18] (version 0.2.7). By directly using Theorem 7 we can compute the 60 first derivatives. The convergence of the first ten eigenvalues are visualized in Figure 6, for two branches. The solutions are visualized in Figure 3.

5.2 Domain decomposition example

We consider a PDE-eigenvalue problem, which we separate into two domains in a way that it leads to a two-parameter eigenvalue problem. Although domain decomposition (and coupling via boundary conditions) is not new for standard eigenvalue problems, the fact that this type of domain decomposition can be phrased as a two-parameter eigenvalue problem has, to our knowledge, not previously been observed. Although the technique seems applicable in several physical dimensions, we consider a one-dimensional problem for reproducibility.

Consider the Helmholtz eigenvalue problem defined in the domain $x\in[x_{0},x_{2}]$,

[TABLE]

with a wavenumber $\kappa$ which is discontinuous in one part of the domain and smooth in another, as in Figure 7. We take a point $x_{1}$ such that $\kappa$ is smooth for $x>x_{1}$, assume that the solution is non-zero in the interface point $x=x_{1}$, and define

[TABLE]

This means we have two separate PDEs for the two domains:

[TABLE]

and

[TABLE]

These are standard linear PDEs (with robin boundary conditions) and the uniqueness of these PDEs implies an equivalence with the original PDE (34). See [23] and references therein for domain decomposition methods for PDE eigenvalue problems.

The wavenumber is given as in Figure 7, i.e., it is discontinuous at several points in $[x_{0},x_{1}]$ and with a high frequency decaying sine-curve in $[x_{1},x_{2}]$, representing a inhomogeneous periodic medium. We invoke different discretizations in the two domains, for the following reasons. Since $\kappa$ is discontinuous in $[x_{0},x_{1}]$ spectral discretization in that domain will not be considerably faster than a finite difference approximation. We therefore use a uniform second order finite difference for (35) to obtain sparse matrices and one sided second order finite different scheme for the boundary condition. A spectral discretization is used for $[x_{1},x_{2}]$ where the wavenumber is smooth. Since $\mu$ appears linearly in the boundary condition, the discretization leads to a two-parameter eigenvalue problem of the type (1). In our setting $A_{1},A_{2},A_{3}$ are large and sparse, and $B_{1},B_{2},B_{3}$ are full matrices of smaller size. We use the discretization parameters such that $n=10^{6}$ and $m=30$, and $x_{0}=0$, $x_{1}=4$ and $x_{2}=5$. In order to make the measurement of error easier, we use left diagonal scaling of the problem such that the diagonal elements of $A(1.0,1.0)$ and $B(1.0,1.0)$ are equal to one.

The eigenvalues and some corresponding eigenfunctions are plotted in Figure 8 and Figure 9. The nonlinear function $g_{1}$ of this problem is given in Figure 10. Clearly the function has singularities for some real $\lambda$-values. The convergence of Algorithm 1 and Algorithm 2 are again compared to [32] in Figure 12. We again conclude that both our approaches are competitive, although not always faster in terms of iterations, but our approach is generally faster in terms of CPU-time. The algorithms are initiated with approximate rounded eigenvectors and eigenvalues close to the eigenvalue $\lambda\approx 18$. We note that our methods do not require a starting value for $\mu$ (in contrast to [32]) which is an attractive feature from an application point of view, since the value $\mu=u^{\prime}(x_{1})/u(x_{1})$ is artificially introduced parameter and may not be easy to estimate.

We apply the tensor infinite Arnoldi method also for this problem. Since this family of methods is based on a power series expansion, one can only expect to be able to compute eigenvalues on the same side of the singularities as the shift. We therefore run the algorithm several times for different shifts, and select the shifts far away from the singularities, as described in Figure 8. The convergence of the two runs are illustrated in Figure 11. Note that the convergence corresponding to one eigenvalue for the shift $\sigma=12$ stagnates. This is because the eigenvalue is close to the singularity, and therefore difficult to compute, as can be seen in Figure 8.

6 Conclusions and outlook

We have presented a general framework to approach two-parameter eigenvalue problems, by nonlinearization to NEPs. Several steps in this technique seem to be generalizable (but beyond the scope of the paper), e.g., to general multiparameter eigenvalue problems essentially by successive application of the elimination. One such elimination leads to a nonlinear two-parameter eigenvalue problem as considered, e.g., in [32].

Our paper uses the assumption $m\ll n$ and that $A_{1}$, $A_{2}$ and $A_{3}$ are large and sparse. We made this assumption mostly for convenience, since this allows us to apply a general purpose method for the parameterized eigenvalue problem (9). If, on the other hand, we wish to solve two-parameter eigenvalue problems where these assumptions are not satisfied, the ideas may still be useful. The GEP (9) may for instance be approached with structured algorithms (exploiting sparsity, low-rank properties and symmetry), or iterative solution methods for the GEP, where early termination is coupled with the NEP-solver.

The generated nonlinear functions $g_{i}$ are algebraic functions, and can therefore contain singularities (e.g. branch point singularities as characterized in Section 2). These can be problematic in the numerical method, and therefore it would useful with transformations that remove singularities. Linearization which do not lead to singularities have been established for rational eigenvalue problems [38].

The problem in Section 5.2 is such that we obtain one large and sparse parameterized matrix $A(\lambda,\mu)$ which is coupled with a small and dense system. The setting matches the assumptions of the paper and is a representative example of cases where the behavior is different in the two physical domains. The example may be generalizable, to other coupled physical systems where the modeling in one domain leads to a much smaller matrix, e.g., using domain decomposition with more physical dimensions. Note however that the presented methods seem mostly computationally attractive if the discretization of one domain is much smaller. If we apply the same technique to domains of equal size, other generic two-parameter eigenvalue methods (such as those in [32]) may be more effective.

Acknowledgment

We thank Olof Runborg, KTH Royal Institute of Technology, for discussions regarding the Helmholtz eigenvalue problem.