The infinite Lanczos method for symmetric nonlinear eigenvalue problems

TL;DR

This paper introduces an infinite Lanczos method for symmetric nonlinear eigenvalue problems, utilizing an infinite-dimensional linearization and a specialized extraction process to improve computational efficiency and stability.

Contribution

It develops a novel infinite-dimensional linearization and applies the indefinite Lanczos method with a tailored eigenpair extraction for symmetric nonlinear eigenproblems.

Findings

Method achieves finite arithmetic under specific conditions

Exploits problem structure for faster computation

Addresses instability issues in eigenpair extraction

Abstract

A new iterative method for solving large scale symmetric nonlinear eigenvalue problems is presented. We firstly derive an infinite dimensional symmetric linearization of the nonlinear eigenvalue problem, then we apply the indefinite Lanczos method to this specific linearization, resulting in a short-term recurrence. We show how, under specific assumption on the starting vector, this method can be carried out in finite arithmetic and how the exploitation of the problem structure leads to improvements in terms of computation time. The eigenpair approximations are extracted with the nonlinear Rayleigh-Ritz procedure combined with a specific choice of the projection space. We illustrate how this extraction technique resolves the instability issues that may occur due to the loss of orthogonality in many standard Lanczos-type methods.

Tables

Table 1. Table 1: Performance of Algorithm  2 applied to the NEP in Section  5.1 of different sizes. The NEPs are obtained by discretizing ( 30 ) respectively with N = 100 , 300 , 500 𝑁 100 300 500 N=100,300,500 nodes in each direction. The projected problems are solved only at the last iteration with: Beyn contour integral method, IAR  [ 27 ] ( 50 50 50 iterations) and IAR ( 100 100 100 iterations).

		Beyn		IAR: $50$ iter.		IAR: $100$ iter.
prob. size		time	conv. eig.	time	conv. eig.	time	conv. eig.
$10000$		1.997 s	11	1.500 s	9	2.570 s	13
$90000$		14.191 s	11	13.120 s	9	14.295 s	19
$250000$		36.177 s	11	35.496 s	9	36.460 s	17

Equations

$$M(\lambda )x=0,$$

$$M(\lambda )=\sum _{m=1}^p{f}_m(\lambda )A_m,$$

$$M(\lambda {)}^T=M(\lambda )$$

$$A=\left(\begin{array}{cccccc}-M_0& 0& 0& 0& \dots & 0\\ 0& M_2& M_3& M_4& \dots & M_k\\ 0& M_3& M_4& & & 0\\ 0& M_4& & & & 0\\ ⋮& ⋮& & & & ⋮\\ 0& M_k& 0& 0& \dots & 0\end{array}\right),$$

$$\left(\begin{array}{ccccc}-M_0& & & & \\ & I& & & \\ & & I& & \\ & & & \ddots & \\ & & & & I\end{array}\right)\left(\begin{array}{c}x\\ \lambda x\\ {\lambda }^{2}x\\ ⋮\\ {\lambda }^{k-1}x\end{array}\right)=\lambda \left(\begin{array}{ccccc}{M}_1& M_2& \dots & M_{k-1}& M_k\\ I& & & & 0\\ & I& & & 0\\ & & \ddots & & ⋮\\ & & & I& 0\end{array}\right)\left(\begin{array}{c}x\\ \lambda x\\ {\lambda }^{2}x\\ ⋮\\ {\lambda }^{k-1}x\end{array}\right).$$

$$S=\left(\begin{array}{cccccc}I& 0& 0& 0& \dots & 0\\ 0& M_2& M_3& M_4& \dots & M_k\\ 0& M_3& M_4& & & 0\\ 0& M_4& & & & 0\\ ⋮& ⋮& & & & ⋮\\ 0& M_k& 0& 0& \dots & 0\end{array}\right).$$

$$\left(\begin{array}{ccccc}-M_0& & & & \\ & I& & & \\ & & I& & \\ & & & I& \\ & & & & \ddots \end{array}\right)\left(\begin{array}{c}\frac{{\lambda }^0}{0!}x\\ \frac{{\lambda }^1}{1!}x\\ \frac{{\lambda }^2}{2!}x\\ \frac{{\lambda }^3}{3!}x\\ ⋮\end{array}\right)=\lambda \left(\begin{array}{ccccc}{M}_1& \frac{1}{2}{M}_2& \frac{1}{3}{M}_3& \frac{1}{4}{M}_4& \dots \\ \frac{1}{1}I& & & & \\ & \frac{1}{2}I& & & \\ & & \frac{1}{3}I& & \\ & & & \ddots & \end{array}\right)\left(\begin{array}{c}\frac{{\lambda }^0}{0!}x\\ \frac{{\lambda }^1}{1!}x\\ \frac{{\lambda }^2}{2!}x\\ \frac{{\lambda }^3}{3!}x\\ ⋮\end{array}\right).$$

$$\mathbf{A}\mathbf{x}=\lambda \mathbf{B}\mathbf{x},$$

$$\mathbf{S}:=\left[\left(\begin{array}{cc}1& \\ & \mathbf{C}\end{array}\right)\otimes \mathbf{e}{\mathbf{e}}^T\right]\circ \left(\begin{array}{cccccc}I& & & & & \\ & M_2& M_3& M_4& M_5& \dots \\ & M_3& M_4& M_5& & \\ & M_4& M_5& & & \\ & M_5& & & & \\ & ⋮& & & & \end{array}\right)$$

$$\mathbf{S}\mathbf{A}\mathbf{x}=\lambda \mathbf{S}\mathbf{B}\mathbf{x}$$

$$c_{i,1}$$

$$c_{i-1,j}$$

$$(j+1)c_{j+k,j}=$$

$$(j+2)c_{j+k-1,j+1}=$$

$$\cdots =$$

$$(j+s)c_{j+k-s+1,j+s-1}=$$

$$\cdots =$$

$$(j+k)c_{j+1,j+k-1}=$$

$$c_{j+k,j}=\frac{(j+k)(j+k-1)\dots (j+1)}{(j+1)(j+2)\dots (j+k)}{c}_{j,j+k},$$

$$\mathbf{S}\mathbf{A}=\left(\begin{array}{cccccc}-M_0& & & & & \\ & c_{1,1}{M}_2& c_{1,2}{M}_3& c_{1,3}{M}_4& c_{1,4}{M}_5& \dots \\ & c_{2,1}{M}_3& c_{2,2}{M}_4& c_{2,3}{M}_5& & \\ & c_{3,1}{M}_4& c_{3,2}{M}_5& & & \\ & c_{4,1}{M}_5& & & & \\ & ⋮& & & & \end{array}\right)$$

$$\mathbf{S}\mathbf{B}=\left(\begin{array}{cccccc}\frac{1}{1}{M}_1& \frac{1}{2}{M}_2& \frac{1}{3}{M}_3& \frac{1}{4}{M}_4& \frac{1}{5}{M}_5& \dots \\ \frac{c_{1,1}}{1}{M}_2& \frac{c_{1,2}}{2}{M}_3& \frac{c_{1,3}}{3}{M}_4& \frac{c_{1,4}}{4}{M}_5& & \\ \frac{c_{2,1}}{1}{M}_3& \frac{c_{2,2}}{2}{M}_4& \frac{c_{2,3}}{3}{M}_5& & & \\ \frac{c_{3,1}}{1}{M}_4& \frac{c_{3,2}}{2}{M}_5& & & & \\ \frac{c_{4,1}}{1}{M}_5& & & & & \\ ⋮& & & & & \end{array}\right).$$

$$\frac{c_{i-1,j}}{j}=\frac{c_{j-1,i}}{i}=\frac{c_{i,j-1}}{i},$$

$${\mathcal{K}}_k(A^{-1}B,q_1):=\mathrm{span}\left(q_1,A^{-1}B{q}_1,(A^{-1}B{)}^2{q}_1,\dots ,(A^{-1}B{)}^{k-1}{q}_1\right)$$

$$A^{-1}B{Q}_k=Q_{k+1}{T}_{k+1,k},$$

$$Q_{k+1}^{T}B{Q}_{k+1}={\Omega }_{k+1},$$

$$(\lambda ,Q_{k}z),\text{mbox}where{T}_{k}z=\lambda {\Omega }_{k}z,$$

$$[\mathbf{S}\mathbf{A}{]}_{N}x=\lambda [\mathbf{S}\mathbf{B}{]}_{N}x$$

$$W=w_1{e}_1^T+Q_{k}D,$$

$$w_1=M_0^{-1}\sum _{j=1}^k\frac{M_j}{j}{\tilde{q}}_j.$$

$$\left(\begin{array}{c}\mathrm{vec}(W)\\ 0\end{array}\right)$$

$$=([\mathbf{A}{]}_{2N}{)}^{-1}([\mathbf{B}{]}_{2N})\left(\begin{array}{c}\mathrm{vec}([Q_k,0])\\ 0\end{array}\right).$$

$$Z=\sum _{m=1}^p{A}_{m}W(G_{k+1}\circ F_m),$$

$$F_m:=\left(\begin{array}{ccccc}{f}_m^{(1)}(0)& f_m^{(2)}(0)& f_m^{(3)}(0)& \dots & f_m^{(k+1)}(0)\\ f_m^{(2)}(0)& f_m^{(3)}(0)& \dots & & f_m^{(k+2)}(0)\\ f_m^{(3)}(0)& \dots & & & f_m^{(k+3)}(0)\\ ⋮& & & & ⋮\\ f_m^{(k+1)}(0)& f_m^{(k+2)}(0)& f_m^{(k+3)}(0)& \dots & f_m^{(2k+1)}(0)\end{array}\right).$$

$$\mathrm{vec}Z=[\mathbf{S}\mathbf{B}{]}_{k+1}\mathrm{vec}W.$$

Full text

Keywords nonlinear eigenvalue problem  $\cdot$ symmetric  $\cdot$ Lanczos

Mathematics Subject Classification (2000) 35P30  $\cdot$ 65H17  $\cdot$ 65F60  $\cdot$ 15A18  $\cdot$ 65F15

1 Introduction

We consider the nonlinear eigenvalue problem (NEP) which consists of computing $(\lambda,v)\in D\times{\mathbb{C}}^{n}\setminus\left\{0\right\}$ such that

[TABLE]

where

[TABLE]

with $D\subset{\mathbb{C}}$ open disk, $f_{m}:D\rightarrow{\mathbb{C}}$ analytic functions, and $A_{m}\in{\mathbb{C}}^{n\times n}$ for $m=1,\dots,p$. In this work we focus on the symmetric NEP, namely we assume that

[TABLE]

Equivalently, $M(\lambda)$ can be expressed as (2) with symmetric matrices $A_{m}^{T}=A_{m}$ for $m=1,\dots,p$. Notice that the complex matrices $A_{m}$ and $M(\lambda)$ are assumed to be symmetric but not necessarily Hermitian. The NEP arises in many areas such as: stability analysis, control theory, wave propagation, etc, and it has been studied in various settings. See the review papers [19, 37], the PhD theses [14, 50], and the problem collection [6]. Specialized software for NEPs has been recently produced: the package NEP-PACK [24], the library SLEPc [21], and even more open-source software. An approach for solving the NEP consists of constructing a linear eigenvalue problem (linearization) whose eigenvalues approximate, or correspond to, the eigenvalues of the original NEP [2, 30, 3, 33, 18, 13, 45]. When the NEP has specific structures such as being: symmetric, Hermitian, Hamiltonian, palindromic, etc, it is preferable to construct a linearization that preserves these structures. Theoretical and algorithmic aspects of structured linearizations have been extensively analyzed [35, 49, 39, 34, 10, 12, 16]. In particular, it has been shown that methods based on structure preserving linearizations, in certain applications, are more robust than other methods that do not take into account the structure [38, 36]. For the polynomial eigenvalue problem (PEP), i.e., the special case where $f_{m}(\lambda)$ in (1) are polynomials, symmetric linearizations are extensively characterized in [22, 11]. A well established class of methods for solving symmetric eigenvalue problems are Lanczos-like methods. More precisely, the Lanczos method, and its variants, can be applied for solving symmetric generalized eigenvalue problems $Ax=\lambda Bx$ where $A,B\in{\mathbb{C}}^{n\times n}$ are symmetric matrices. The original approach [29] was developed for the case $B=I$, a generalization for $B$ positive definite is presented in [42, Ch.15, Sec.11]. A further extension of this approach, known as indefinite Lanczos method, for the case where $A$ and $B$ are Hermitian or symmetric matrices is discussed in [43] and [4, Section 8.6.1]. Lanczos methods belong to the class of Krylov methods. They exploit the fact that, during the construction of the orthogonal basis of the Krylov space, due to the symmetry of the problem, the orthogonalization can be performed in a more efficient way with a three-term recurrence. The main disadvantage of having a three-term recurrence is that, in floating-point arithmetic and without further specializations, the basis vectors are often affected by loss of orthogonality, resulting in slow convergence of the Ritz pairs and numerical instability [44, 46, 53, 1].

In this work, we present a new symmetric linearization for the symmetric NEP, resulting in a symmetric, linear, and infinite dimensional eigenvalue problem. Symmetric generalized eigenvalue problems can be solved with the indefinite Lanczos method [4, Section 8.6.1]. We present a new method that corresponds to adapting, in an efficient and robust way, the indefinite Lanczos method to the derived linearization. In order to cure the slow convergence, that is due to the loss of orthogonality, we use the nonlinear Rayleigh–Ritz procedure, combined with a proper choice of the projection space, for extracting eigenpair approximations by exploiting the structure of the derived linearization. The proposed method is competitive, in terms of robustness and complexity, with Arnoldi-like methods for NEPs that perform the full orthogonalization.

The paper is organized as follows: in Section 2 we prove that the symmetric NEP is equivalent to a symmetric, linear, and infinite dimensional eigenvalue problem. In Section 3 we derive a method, in finite arithmetic, which consists of applying the indefinite Lanczos method to the derived linearization. In Section 4 we show how the computation time of the resulting method can be considerably reduced by exploiting additional NEP structures. In Section 5 we illustrate the performance of this new approach with numerical simulations by solving large and sparse NEPs. The simulations were carried out in the Julia programming language [9] with NEP-PACK [24], which is an open source Julia package for NEPs.

The method we derive can been seen as an Arnoldi-like method applied to an iteratively expanding linearization, or equivalently to a infinite dimensional linear operator. Other methods that are based on these ideas are: infinite Arnoldi [27] and its tensor variant [26], NLEIGS [20], and CORK [51]. There are also methods based on the bi-orthogonalization procedure, which also lead to a three-term recurrence, presented in [31, 17]. However, these methods and their variations, in the way they are presented and without further research, are not capable of taking advantage of the symmetry of the NEP.

In the rest of this work, vectors and matrices are denoted as $v=[v_{i}]_{i=1}^{n}$ and $A=[a_{i,j}]_{i,j=1}^{n,m}$ respectively, whereas bold letters represent block-vectors and block-matrices with infinite size, namely $\bm{v}=[v_{i}]_{i=1}^{\infty}$ with $v_{i}\in{\mathbb{C}}^{n}$ and $\bm{A}=[A_{i,j}]_{i,j=1}^{\infty}$ with $A_{i,j}\in{\mathbb{C}}^{n\times n}$. The matrix $[\bm{A}]_{k}=[A_{i,j}]_{i,j=1}^{k}\in{\mathbb{C}}^{nk\times nk}$ consists of the main sub-matrix obtained by extracting the first $k$-blocks. The Kronecker product and the Hadamard (element-wise) product are denoted by $\otimes$ and $\circ$ respectively. The vectors $e_{j}$ and $\bm{e}_{j}$ have zeros as elements except one in the $j$-th position whereas $e$ and $\bm{e}$ are the vectors with all ones. Without loss of generality, after a change of variables in (2), we assume that the region of interest $D\subset{\mathbb{C}}$ is a disk centered in the origin. The derivatives of (2) will be denoted as $M_{i}:=M^{(i)}(0)$. We will denote by $A^{T}$ the transpose (not conjugate transpose) of the matrix $A\in{\mathbb{C}}^{n\times n}$.

2 Indefinite Lanczos method in infinite dimensional settings

In order to derive a symmetric linearization for NEPs, we first review a specific linearization technique for PEPs. This technique consists of symmetrizing the companion linearization and it is presented in [39, Theorem 2.1] that we recall below. The approach is previously reported in [28, Ch.4 Sec. 2] for the scalar case.

Theorem 2.1** (Mehrmann and Watkins [39], Lancaster [28]).**

Consider the polynomial eigenvalue problem $M(\lambda)v=0$ where $M(\lambda)=\sum_{j=0}^{k}M_{j}\lambda^{j}$ and $M_{k}$ nonsingular. Then, the pencil $A-\lambda B\in{\mathbb{C}}^{nk\times nk}$, where

[TABLE]

has the same eigenvalues of $M(\lambda)$. If $M_{j}$ are symmetric, i.e., $M(\lambda)$ is symmetric, then $A$ and $B$ are symmetric. If $M(\lambda)v=0$, then $[v^{T},\lambda v^{T},\dots,\lambda^{k-1}v^{T}]^{T}$ is an eigenvector of $A-\lambda B$.

The proof is based on the following argument. A pair $(\lambda,x)$ that fulfill $M(\lambda)x=0$ defines an eigenpair of the companion linearization [18], namely

[TABLE]

To obtain (4) we multiply (5) on the left by the matrix

[TABLE]

The main disadvantage of using this linearization in practice is that the blocks forming the eigenvectors of the pencil, defined by (4), grow or decay exponentially, depending on the value of $|\lambda|$. More precisely, the norm of the $j$-th block is $|\lambda|^{j}\|x\|$ and, if the PEP has high degree, this leads to overflow or underflow when the eigenpairs of (4) are computed numerically. In order to resolve this issue, in this section we consider the scaled companion linearization presented [27, Section 5.1]. We extend the ideas used in Theorem 2.1 to symmetrize the scaled companion linearization. Moreover, we consider the NEP in its general form (1), therefore we derive a linearization that involves vectors and matrices with infinite size.

2.1 An infinite dimensional symmetric linearization

The NEP (1) is equivalent to a linear and infinite dimensional eigenvalue problem, see [27, Section 5.1]. More precisely, if $(\lambda,x)$ is an eigenpair of (1), the following relation between vectors and matrices of infinite size is fulfilled

[TABLE]

The equation (7) defines a linear and infinite dimensional eigenvalue problem

[TABLE]

where $\bm{A},\bm{B},\bm{x}$ are matrices and vector of infinite size defined accordingly. Clearly, the linearization (7) is never symmetric. However, if the NEP is symmetric, i.e., it holds (3), then it is possible to symmetrize (7) with a similar technique as in Theorem 2.1. More precisely, since we consider a scaled and infinite companion linearization, in the following theorem we derive a scaled and infinite version of the matrix (6) that symmetrizes (7).

Theorem 2.2** (Symmetric linearization).**

Assume that the NEP (2) is symmetric, i.e., it holds (3), then there exists a unique matrix $\bm{C}$ such that

[TABLE]

is a symmetrizer for (7), namely

[TABLE]

is a symmetric eigenvalue problem. The vector $\bm{e}$ has infinite length with ones in all the entries. The coefficients of the matrix $\bm{C}$ fulfill the following relations

[TABLE]

Proof.

We start observing that $M_{j}$, for $j\geq 0$, are symmetric matrices as consequence of (3). The relations (11) uniquely define a matrix $\bm{C}$ since the first column is fixed in (11a) and the $j$-th column is computed by the $(j-1)$-th column in (11b). We start by showing that the matrix $\bm{C}$ is symmetric. Let us consider $i>j$, namely $i=j+k$ for some positive integer $k$. By iteratively using (11b) we obtain the relations

[TABLE]

that combined together give

[TABLE]

that is, $c_{i,j}=c_{j,i}$. The case $i<j$ is analogous and we conclude that the matrix $\bm{C}$ is symmetric.

By multiplying (8) on the left by the matrix (9) we get

[TABLE]

and

[TABLE]

The matrix $\bm{S}\bm{A}$ is symmetric because $\bm{C}$ and $M_{j}$, for $j\geq 0$, are symmetric. By using (11a) we get that the first block-row of $\bm{S}\bm{B}$ is equal to its first block column, whereas the equation (11b) and the symmetry of $\bm{C}$ gives the relation

[TABLE]

which directly implies that the $(i,j)$-th and the $(j,i)$-th blocks of $\bm{S}\bm{B}$ are equal. Hence the matrix $\bm{S}\bm{B}$ is symmetric and (10) is a symmetric eigenvalue problem. ∎

Remark 2.3**.**

The eigenvalue problems (7) and (10) have the same eigenpairs if the symmetrizer (9) is nonsingular, namely $\bm{S}\bm{x}=0$ only for $\bm{x}=0$. In the next section we assume that $[\bm{S}]_{2N}$ is invertible for an $N$ large enough. This condition can be phrases in terms of solvability of a specific matrix equation as discussed in Observation 3.4.

The method that we refer to as infinite Lanczos consists of applying the indefinite Lanczos method (Algorithm 1), described in the next section, to the symmetric eigenvalue problem (10).

3 Infinite Lanczos method

3.1 Indefinite Lanczos method

Eigenpair approximations to the generalized eigenvalue problem $Ax=\lambda Bx$, with $A,B\in{\mathbb{C}}^{n\times n}$ symmetric matrices, not necessarily Hermitian, and $A$ nonsingular, can be obtained by using the indefinite Lanczos method [4, Section 8.6.1] that is summarized in Algorithm 1. The method consists of computing an orthogonal basis of the Krylov space

[TABLE]

by using, instead of the (standard) Euclidean scalar product, the indefinite scalar product defined by the matrix $B$, namely $x^{T}By$ is the $B$-product between $x,y\in{\mathbb{C}}^{n}$. The fact that $A^{-1}B$ is self-adjoint, with respect to this indefinite scalar product, leads to the property that the $B$-orthogonal basis of the Krylov space can be computed with a three-term recurrence. In particular, at the $k$-th iteration of Algorithm 1, the following relations are fulfilled

[TABLE]

where the diagonal matrix $\Omega_{k+1}:=\operatorname{diag}(\omega_{1},\dots,\omega_{k+1})$ and the tridiagonal matrix $T_{k+1,k}=[t_{i,j}]_{i,j=1}^{k+1}$ contain the orthogonalization and normalization coefficients. The matrix $Q_{k+1}$ is $B$-orthogonal in the sense of (15b) and its columns, generated with a three-term recurrence, span the Krylov space (14). The Ritz pairs of (15a), defined as follows,

[TABLE]

provide an approximation to the eigenpairs of the original problem. Since the indefinite scalar product defined by $B$ is in general degenerate, there may be cases of break down in Algorithm 1. We refer to [4, Section 8.6.1] and reference therein for a detailed discussion of this issue.

3.2 Infinite Lanczos method in finite arithmetic

We now derive a method that consists of applying the indefinite Lanczos method (Algorithm 1) to the symmetric eigenvalue problem

[TABLE]

obtained by extracting the main block sub-matrices from (10), where $N$ is a non-fixed parameter greater than the number of iterations performed in Algorithm 1. The method we derive is independent on $N$ and, under the assumption that $\bm{S}$ given in (9) is invertibile, corresponds to apply Algorithm 1 directly to the linear infinite dimensional eigenvalue problem (10) with a specific starting vector. This equivalence is formally presented in Theorem 3.5 at the end of this section.

Algorithm 1 can be efficiently applied to (17) by exploiting the structure of the matrices (17). We start with Step 1 that can be performed as stated in the following result.

Theorem 3.1** (Action of $[\bm{S}\bm{A}]_{N}^{-1}[\bm{S}\bm{B}]_{N}$ ).**

Assume that $[\bm{S}]_{2N}$ is invertible, let $q_{k}\in{\mathbb{C}}^{Nn}$ be such that only its first $k$ blocks are nonzero, corresponding to the columns of $Q_{k}:=[\tilde{q}_{1},\dots,\tilde{q}_{k}]$. Then, only the first $k+1$ blocks of $w=[\bm{S}\bm{A}]_{N}^{-1}[\bm{S}\bm{B}]_{N}q_{k}$ are nonzero, corresponding to the columns of $W:=[w_{1},\dots,w_{k+1}]$ given by

[TABLE]

where $D\in{\mathbb{R}}^{k\times(k+1)}$ is a diagonal matrix with coefficients $d_{j,j+1}=1/j$ and

[TABLE]

Proof.

By using the specific structure of the matrices (12) and (13), the nonzero blocks of $w$ fulfill $\operatorname{vec}(W)=([\bm{S}\bm{A}]_{N})^{-1}([\bm{S}\bm{B}]_{N})\operatorname{vec}([Q_{k},\ 0])$. We can then derive the following relations

[TABLE]

Hence $\operatorname{vec}(W)=([\bm{A}]_{N})^{-1}([\bm{B}]_{N})\operatorname{vec}([Q_{k},\ 0])$, this directly implies (18) and (19), c.f., [27, Section 4.2]. ∎

By using the previous result, we conclude that in Algorithm 1, if $q_{1}$ has only the first block which is nonzero, then $q_{k}$ at the $k$-th iteration will have $k$ nonzero blocks. This is due to the fact that, none of the steps, except Step 1, introduce fill-in in the vectors $q_{1},q_{2},\dots,q_{k}$. In Step 1 of Algorithm 1 the products $z^{T}q_{k}$, $z^{T}q_{k-1}$ and $z^{T}w$ are computed. Observe that the vectors multiplied by $z$ have at most $k+1$ nonzero blocks. Therefore, even if $z=[\bm{S}\bm{B}]_{N}w$ is in general a full vector, only the first $k+1$ blocks are required. These blocks can be computed as follows.

Theorem 3.2** (Action of $[\bm{S}\bm{B}]_{N}$).**

Let us consider $z:=[\bm{S}\bm{B}]_{N}w$, where $w$ is given as in Theorem 3.1. Then the first $k+1$ blocks of $z$, corresponding to the columns of $Z:=[z_{1},\dots,z_{k+1}]$, fulfill the relation

[TABLE]

where $G_{k+1}\in{\mathbb{R}}^{(k+1)\times(k+1)}$ has coefficients $g_{j,1}=g_{1,j}=1/j$ for $j=1,\dots,k+1$ and $g_{i,j}=c_{i-1,j}/j$ and

[TABLE]

Proof.

Since $w$ has only the first $k+1$ blocks that are nonzero, we can express

[TABLE]

By using (13) and that $M_{j}=\displaystyle\sum_{m=1}^{p}f_{m}^{(j)}(0)A_{m}$, we can decompose

[TABLE]

Equation (20) follows by combining (22) and (21), and by using the properties of the Kronecker product. ∎

Observation 3.3**.**

The scalar product between the vectorization of two matrices can be carried out directly in matrix form by using the Hadamard product as follows: $(\operatorname{vec}Z)^{T}\operatorname{vec}W=\tilde{e}^{T}(Z\circ W)e$ with $\tilde{e}\in{\mathbb{R}}^{n}$ and $e\in{\mathbb{R}}^{k}$.

Observation 3.4**.**

With the same reasoning as in Theorem 3.2, we can decompose (9) as

[TABLE]

Therefore, we can relate the invertibility of $[\bm{S}]_{2N}$ with the solvability of the following linear matrix equation

[TABLE]

for any $B\in{\mathbb{C}}^{2N\times n}$. Linear matrix equations are extensively studied in recent literature. See the review paper [47] and reference therein. In the numerical examples reported in Section 5 we never encounter a case when $[\bm{S}]_{2N}$ was singular. A case when this matrix is obviously singular is when the NEP is defined by polynomial functions. Although the theory does not cover this case, we have successfully applied the method we are deriving without introducing any breakdown or instability.

Figure 1 illustrates the structure of the matrices and vectors, involved in Algorithm 1, when applied to (17) with a starting vector that has only the first block which is nonzero. At iteration $k$ only the vectors $q_{k-1},q_{k}$ are needed, therefore they are the only vectors that need to be stored.

Algorithm 2 is the combination of the results presented in this section. More precisely, Algorithm 2 is the reformulation of Algorithm 1, applied to (17), where the nonzero blocks of the vectors $q_{k},q_{k-1},w$, and the needed blocks of $z$, are stored as columns of the matrices $Q_{k},Q_{k-1},W$ and $Z$. Moreover, the size of the linearization (17) is implicitly expanded at each iteration. Observe that at iteration $k$ only the first $2k+1$ derivatives $f^{(j)}_{m}(0)$ for $m=1,\dots,p$ and $j=1,\dots,2k+1$ are needed. We now conclude this section by showing the equivalence between Algorithm 1, directly applied to the infinite dimensional problem (10), and Algorithm 2.

Theorem 3.5** (Infinite dimensional equivalence).**

Assume that the matrix $\bm{S}$, given in (9), is invertible and let $\bm{q}_{1}$ be an infinite length vector with only the first block $q_{1}$ nonzero. Then, Algorithm 1, with stating vector $\bm{q}_{1}$, is applicable to (10) and the matrices $Q_{k}$, that have as columns the first $k$ nonzero blocks of $\bm{q}_{k}$, $\underline{T}_{k}$, and $\omega_{k}$ are equal to the homonymous matrices generated by Algorithm 2 with starting matrix $Q_{1}=[q_{1}]$.

Proof.

We denote by $\bm{q}_{1},\bm{q}_{2},\dots,\bm{q}_{k}$ the infinite-length vectors generated by Algorithm 1 and by $Q_{1},Q_{2},\dots,Q_{k}$ the matrices generated by Algorithm 2. The proof is based on induction over the iteration count $k$. The result is trivial for $k=1$. Suppose the results holds for some $k$. In Step 1 of Algorithm 1, by using that $\bm{S}$ is invertible, we have

[TABLE]

By using the induction hypothesis, $\bm{q}_{k}$ has $k$ nonzero blocks, corresponding to the column of the matrix $Q_{k}$ generated at the $(k-1)$-th iteration of Algorithm 2. Because of the structure of the matrices (12) and (13), we get that $\bm{w}$ has only the first $k$ blocks which are nonzero, corresponding to the columns of the matrix $W$ that fulfills (18). Therefore this matrix corresponds to the matrix computed in Step 2 of Algorithm 2. In the Step 1 of Algorithm 1 we compute $\bm{z}=\bm{S}\bm{B}\bm{w}$. This vector is in general full. However, in the Step 1 of Algorithm 1 the products $\bm{z}^{T}\bm{q}_{k}$, $\bm{z}^{T}\bm{q}_{k-1}$ and $\bm{z}^{T}\bm{w}$ are computed. By induction hypothesis, $\bm{q}_{k}$, $\bm{q}_{k-1}$ have respectively $k$ and $k-1$ nonzero blocks, corresponding to the columns of the matrices $Q_{k}$ and $Q_{k-1}$ generated by Algorithm 2. Therefore, only the first $k+1$ blocks of $\bm{z}$ are required and they can be computed with the same reasoning of Theorem 3.2. More precisely, the first $k+1$ blocks of $\bm{z}$ are the columns of the matrix $Z$ that fulfills (20). Therefore this matrix coincides with the matrix generated by Step 2 of Algorithm 2. In order to conclude that $\bm{q}_{k+1}$ has only we first $k+1$ nonzero blocks, corresponding to the columns of the matrix $Q_{k+1}$ generated by Algorithm 2, we only need to use the property $\|M\|_{F}=\|\operatorname{vec}(M)\|_{2}$ for every matrix $M$.

∎

3.3 Robust extraction of eigenpair approximations

We propose to enhance Algorithm 2 as follows. We consider the projected NEP

[TABLE]

where $V\in{\mathbb{C}}^{n\times k}$ is an orthogonal matrix. Under the assumption that $V$ posses good approximation properties, eigenpair approximations to the NEP (1) are given by $(\lambda,Vz)$. This can be seen as the Galerkin projection method that uses the range of $V$ as projection space. This technique for extracting eigenpair approximations, called nonlinear Rayleigh–Ritz procedure or subspace acceleration, is often used to improve properties of more basic algorithms, e.g., the nonlinear Arnoldi method [52], Jacobi–Davidson methods [14, 7], infinite Arnoldi [25], block preconditioned harmonic projection methods [54], and many more.

In our framework, there is a natural choice for the projection space. The matrix $V$ is chosen as the orthogonal matrix whose columns span the subspace of vectors obtained by extracting the first column from $Q_{1},\dots,Q_{k}$ generated by Algorithm 2. The reason this matrix contains good approximation properties is due to the following argument. In the Ritz pairs extraction described in Section 3.1, the eigenvector approximations of (7) (or equivalently of (10)) are given by the first block row of the Ritz vectors (16). Thus, the eigenvector approximations to the NEP (1) are also obtained by the first block of these Ritz vectors and thus by the first block row of the Krylov basis, namely by the columns of the proposed matrix $V$.

The projected problem (23) has size $k$ equal to the number of iterations, which is typically, in the context of Krylov methods, a small number, i.e., $k\ll n$. Therefore, the projected problem (23) has small size and solving (23) is not the computationally dominating part of Algorithm 2. In the numerical experiments we have tested, for solving (23), the following methods: Beyn’s contour integral method [8], NLEIGS [20] and IAR [27]. The choice the method for solving (23) depends on the features of the original problem (1) and there is not a favorite candidate. For example, one may want to exploit that the projected problem (23) is defined by the same nonlinear functions and may inherit several features of the original NEP 1 such as being symmetric or palindromic. Therefore, this problem can be solved in a more robust way with structure preserving methods.

4 Indefinite scalar product computation

Under the assumption that the linear systems with the matrix $M_{0}$ can be efficiently solved, e.g., exploiting the sparsity, the dominating part of Algorithm 2 is the Step 2, namely the computation of (20), which has complexity $\mathcal{O}(k^{2}n)$. In this section we derive efficient methods for computing this quantity.

4.1 Computation of Step 2: General case

The following theorem provides an effective approximation to (20) without any specific assumption on the coefficients of (2).

Theorem 4.1**.**

Let $U,V\in{\mathbb{R}}^{n\times q}$ the factors of the best rank $q$ approximation, with respect to the Euclidean norm, to the matrix $G_{k+1}$. Then

[TABLE]

is such that

[TABLE]

Proof.

The approximation (24) is obtained by replacing $G_{k+1}$ with $UV^{T}$ in (20) and using $u_{j}v_{j}^{T}\circ F_{m}=\operatorname{diag}(u_{j})F_{m}\operatorname{diag}(v_{j})$. The equation (25) follows by the triangular inequality and by the fact that the Frobenius norm is sub-multiplicative with respect to the Hadamard product. ∎

The approximation (24) is effective, with $q$ small, since the matrix $G_{k}$, which is problem independent, has a fast decay in the singular values. In Figure 2 the singular values111The singular values are computed in BigFloat arithmetic using the package GenericSVD. of this matrix are displayed for different sizes $k$. Moreover, the computation of (24), requires less computation time than (20) since the products with the Hankel matrices $F_{j}$ can be efficiently computed with FFTs [32, Section 4]. The complexity for computing (24) is $\mathcal{O}(nk\log k)$.

4.2 Computation of Step 2: delay eigenvalue problem

The stability analysis of delay systems of the form

[TABLE]

is related to solving NEPs, referred to as delay eigenvalue problems, see [23, Ch. 2] and [40, Ch. 1 Sect. 1.2], defined as follows:

[TABLE]

In this case the matrices $F_{m}$ have at most rank one. More precisely, a direct computation leads to $F_{1}=-e_{1}e_{1}^{T}$, $F_{2}=0$ and for $m\geq 3$ we get $F_{m}=-\tau_{m}vv^{T}$ with $v_{j}=(-\tau_{m})^{j-1}$. Therefore, the computation time of (24) is much lower than (20) since it involves only products with low rank matrices. The complexity for computing of (24) is reduced, by exploiting the low-rank structure, to $\mathcal{O}(npk)$.

4.3 Computation of Step 2: polynomial plus low-rank structured NEPs

In certain applications (2) can be written as sum of a polynomial and a low-rank part. See, e.g., [14, Ch. 2 Sec. 4], [50, Ch. 1 Sec. 1.2], [6, gun problem] and [5, Sec. 6.2.2]. More precisely:

[TABLE]

with $U_{m},U_{m}\in{\mathbb{C}}^{n\times r_{m}}$ and $r_{m}\ll n$. In this case we split (20) in the polynomial and low-rank terms, namely $Z=Z_{p}+Z_{lr}$ with

[TABLE]

The term (29a) can be efficiently approximated as in (24) by exploiting the low-rank structure of the matrices $F_{m}$. Moreover, since $U_{m}U_{m}^{T}$ in (29b) have low rank, the computation of $Z$ with (29), respecting the order given by the parentheses in (29b), requires less computation time than (20). The complexity of (29) is $\mathcal{O}((d+r)n)$ where $r=\displaystyle\max_{d+1\leq t\leq p}r_{m}$.

5 Numerical simulations

In the following numerical experiments222All simulations were carried out with Intel octa core i7-4770 CPU 3.40GHz and 24 GB RAM. we use, as error measure, the relative error defined as follows

[TABLE]

An eigenpair approximation is marked as “converged” if $\operatorname{Err}(\lambda,x)<10^{-8}$. The software used in these simulations is implemented in the Julia programming language [9], and publicly available in the Julia package NEP-PACK [24] 333The scripts reproducing several of the presented examples are directly available in the web-page: https://people.kth.se/~gmele/InfLan/.

5.1 Delay eigenvalue problem

We consider the delay eigenvalue problem arising from the spatial discretization of the following partial delay differential equation

[TABLE]

where $a(\xi)=-\xi_{1}\sin(\xi_{1}+\xi_{2})$, resulting in a problem of the form (26) where all the matrices are real and symmetric. The spatial domain is partitioned with a uniform equispaced grid with $N$ points in each direction. The Laplace operator is discretized by the 5-points stencil finite difference approximation, leading to a NEP of size $n=N^{2}$, cf. [5, Section 6.2.1]. More precisely, the NEP is defined as $M(\lambda)=-\lambda I+A_{2}+e^{-\lambda}A_{3}$ where $I\in{\mathbb{R}}^{n\times n}$ and the other matrix coefficients are given by

[TABLE]

with $h:=\pi/(N-1)$ discretization step. We run $50$ iterations of Algorithm 2. In Figure 4 and Figure 4 we illustrate the robustness of the strategy for extracting the eigenpair approximations presented in Section 3.3. More precisely, we compare the standar approach for extracting the eigenpair approximations, which is based on the computation of the Ritz pairs, with the more robust approach consisting of solving the projected NEP. In Figure 4 is displayed the spectrum and the converged eigenvalues, respectively computed with the Ritz and the projected NEP approach. In this example we solve the projected NEP with the Beyn contour integral method444The disk of interest is set with center in the origin and radius $4$ with $N=1000$ discretization points and $\mbox{tol}_{\mbox{\tiny res}}=\mbox{tol}_{\mbox{\tiny rank}}=10^{-8}$. [8]. The error history is presented in Figure 4. As expected, the convergence of the Algorithm 2, with the Ritz pair approximation, appear to be slower with respect to the convergence of the Algorithm 2 with the more robust eigenpair extraction based on solving the projected NEP.

The performance of Algorithm 2 is affected by the method used for solving the projected NEP. In Table 1 we compare the time, and the number of computed eigenvalues, after $50$ iterations of Algorithm 2 combined with three different NEP solves for the projected problem: Beyn contour integral method (with the same settings as before), IAR [27] ($50$ iterations) 555We perform $50$ iteration of IAR for solving the projected NEP and return the converged eigenvalues. and IAR ($100$ iterations).

5.2 A benchmark problem

We now illustrate the performance of Algorithm 2 for solving a NEP that is symmetric but not Hermitian. We consider the gun problem that belong to the problem collection [6]. This NEP has the following form

[TABLE]

where $\sigma_{2}=108.8774$. The matrices $A_{j}\in{\mathbb{R}}^{9956\times 9956}$, for $j=1,\dots,4$, are real and symmetric. The NEP 31 can be written in the form (28) since the matrix coefficients of the nonlinear part have low-rank, namely $\operatorname{rank}(A_{3})=19$ and $\operatorname{rank}(A_{4})=65$. The eigenvalues of interest are located inside the the closed disk centered in $250^{2}$ and with radius $5\cdot 10^{4}$. Before applying Algorithm 2, the problem is shifted and scaled. We set the parameters to $\lambda=\lambda_{0}+\alpha\hat{\lambda}$ where $\lambda_{0}=300^{2}$ and $\alpha=(300-200)^{2}$. This problem has been solved with various methods [27, 20, 51, 31, 17] and, by numerical evidences, there are $21$ eigenvalues in the region of interest. For this problem we use IAR for solving the projected NEP. More precisely we test two variants: IAR (50 iterations) and IAR (200 iterations). As showed in the numerical experiment in Section 5.1, the robustness of the whole Algorithm 2 is effected by the choice of the method used for solving the projected NEP. In Figure 6 we can see that more eigenvalues converge when we solve more accurately the projected problem. The error history is presented in Figure 6. Solving more accurately the projected NEP is necessary to compute the outermost eigenvalues.

5.3 A random symmetrized problem

In conclusion we illustrate how Algorithm 2 can be used for solving a nonsymmetric NEP. We introduce a symmetrization technique, consisting of doubling the problem size, based in the idea presented in [41, Sect. 5]. Namely, we define the symmetrized NEP as

[TABLE]

Observe that if $\left(\lambda,[y^{T},x^{T}]^{T}\right)$ is an eigenpair of (32), then $(\lambda,x)$ is an eigenpair of $M(\lambda)$. We now consider the symmetrization, in the sense of (32), of the following NEP that is artificially constructed:

[TABLE]

where $A_{j}\in{\mathbb{C}}^{500\times 500}$ are defined as follows: $A_{1}$ is the bidiagonal matrix with elements equal to $500$ in the upper and lower diagonal, $A_{2}$ is the identity matrix, $A_{3}=A_{1}/500$ and $A_{4}$ is a diagonal matrix with elements equal to $i$ (complex unit) in the lower diagonal. We perform $50$ iterations of Algorithm 2 and solve the projected NEP with NLEIGS 666We used the static variant with Leja-Bagby points automatically generated from the polygonal target. The shift is in zero and it is kept constant during the iterations. [20] by targeting the eigenvalues contained in the rectangle with opposite vertices in $-1.5-1.5i$ and $0.5+0.5i$. In Figure 8 is illustrated the spectrum and the converged eigenvalues, in the region of interest, after $50$ iterations of Algorithm 2. The error history is illustrated in Figure 8.

6 Conclusions and outlook

We have presented a method for solving symmetric NEPs. We have also illustrated how the problem structure, in particular the structure of the matrices and functions in (2), can be exploited in order to reduce the computation time. However, there are NEPs that cannot be written in the format (2), e.g., the waveguide eigevalue problem [26], the reformulation of the Dirichlet eigenvalue problem with the boundary element method [48, 15], etc. For some of these problems, only a routine for computing $M_{k}x$ is available. We believe that further research can potentially extend the applicability of Infinite Lanczos to such problems.

In the numerical experiment in Section 5.3, we have successfully solved a nonsymmetric NEP in the following way. Firstly we constructed a symmetric NEP (32) whose eigenvalues are also eigenvalues of the original NEP. Then we applied the infinite Lanczos to the symmetrized problem (32). The matrices (32) have clearly a very well defined block structure. We believe that infinite Lanczos can be further specialized for solving nonsymmetric NEPs by exploiting these structures. In conclusion we also believe that similar ideas can be extended to NEPs that are Hermitian, namely, $M(\lambda)^{H}=M(\bar{\lambda})$ where $\bar{\lambda}$ represents the complex conjugate of $\lambda\in{\mathbb{C}}$ and $M(\lambda)^{H}$ the Hermitian, or conjugate transpose, of the matrix $M(\lambda)$.

7 Acknowledgement

The author wishes to thank Sarah Gaaf (Utrecht University) who suggested a first approach to extend Theorem 2.1 to the scaled companion linearization. The author also thanks Elias Jarlebring (KTH Royal Institute of Technology) who has provided feedback and help in the first stage of the project