Verified partial eigenvalue computations using contour integrals for Hermitian generalized eigenproblems

TL;DR

This paper introduces a verified computational approach for partial eigenvalues of Hermitian generalized eigenproblems using contour integrals, providing rigorous error bounds and efficient verification techniques.

Contribution

It develops a method that rigorously bounds errors in contour integral-based eigenvalue computations, improving accuracy and efficiency over standard methods.

Findings

The proposed method provides rigorous error bounds for complex moments.

Numerical experiments demonstrate superior performance compared to standard methods.

The approach is potentially efficient for parallel computations.

Abstract

We propose a verified computation method for partial eigenvalues of a Hermitian generalized eigenproblem. The block Sakurai-Sugiura Hankel method, a contour integral-type eigensolver, can reduce a given eigenproblem into a generalized eigenproblem of block Hankel matrices whose entries consist of complex moments. In this study, we evaluate all errors in computing the complex moments. We derive a truncation error bound of the quadrature. Then, we take numerical errors of the quadrature into account and rigorously enclose the entries of the block Hankel matrices. Each quadrature point gives rise to a linear system, and its structure enables us to develop an efficient technique to verify the approximate solution. Numerical experiments show that the proposed method outperforms a standard method and infer that the proposed method is potentially efficient in parallel.

Figures

Tables

Table 1. Table 1 : Verified eigenvalues for artificially generated problems 1 .

$\ell$	Eigenvalues near $2$
5	${1.715370325}_{576319}^{647292}$,	${1.9048361618}_{37662}^{83239}$,	${2.095163824}_{458061}^{553889}$,	${2.284629679}_{395375}^{449021}$
6	${1.8551304200}_{49649}^{88382}$,	${1.95167256291}_{1146}^{316}$,	${2.048327546}_{241311}^{750658}$,	${2.144869623}_{845262}^{957331}$
7	${1.9269559899}_{20654}^{32997}$,	${1.97564721303}_{2271}^{7844}$,	${2.02435287}_{2964093}^{3629513}$,	${2.073044084}_{849138}^{931074}$
8	${1.963329778}_{797455}^{856981}$,	${1.9877759682}_{78868}^{90228}$,	${2.0122240077}_{57747}^{94568}$,	${2.03667024}_{4513424}^{5464907}$
9	${1.981628386}_{583106}^{639699}$,	${1.993876062}_{151237}^{477787}$,	${2.006123967}_{790464}^{896611}$,	${2.01837164}_{2984799}^{475554}$
10	${1.99080513}_{0338708}^{2328691}$,	${1.9969350}_{32012484}^{40118075}$,	${2.003064970}_{239924}^{538566}$,	${2.00919487}_{492424}^{9374675}$
11	${1.99540031}_{1711072}^{4525073}$,	${1.998466772}_{568391}^{799085}$,	${2.0015332366}_{65236}^{74622}$,	${2.004599697}_{805267}^{937402}$
12	${1.997699590}_{620666}^{95202}$,	${1.99923319}_{1994498}^{2347755}$,	${2.000766798}_{470323}^{533992}$,	${2.002300408}_{831211}^{946773}$
13	${1.998849650}_{636445}^{85915}$,	${1.999616549}_{351151}^{578838}$,	${2.000383446}_{084906}^{393767}$,	${2.0011503}_{36316831}^{4576172}$
14	${1.9994247}_{87629918}^{92481914}$,	${1.999808262}_{201378}^{706644}$,	${2.000191734}_{004361}^{485737}$,	${2.000575205}_{226003}^{599759}$
15	${1.99971238}_{5433325}^{6833255}$,	${1.99990412}_{6923384}^{9815281}$,	${2.0000958}_{6580425}^{74431981}$,	${2.00028761}_{0865148}^{1877939}$
16	${1.999856}_{177687935}^{205050025}$,	${1.9999520}_{36785509}^{89032199}$,	${2.00004}_{7548157389}^{8322762835}$,	${2.00014}_{360912421}^{4007676757}$
17	${1.999928}_{088198068}^{100686499}$,	${1.9999760}_{2728946}^{34322243}$,	${2.00002396}_{5890281}^{8809379}$,	${2.00007190}_{3283288}^{4763985}$
18	${1.9999640}_{4108165}^{5305932}$,	${1.9999880}_{07405001}^{23539268}$,	${2.00001}_{1955135737}^{2012521721}$,	${2.0000359}_{28694293}^{75584551}$
19	${1.99998202}_{1657987}^{5400559}$,	${1.99999}_{399498807}^{4021890356}$,	${2.00000599}_{0718144}^{3112617}$,	${2.00001797}_{5380837}^{6886183}$
20	${1.9999910}_{09237281}^{13961167}$,	${1.99999}_{691993317}^{7090330345}$,	${2.00000}_{2940857038}^{3052063462}$,	${2.0000089}_{79158735}^{96723335}$

Equations

$$A{\mathbf{x}}_i={\lambda }_{i}B{\mathbf{x}}_i,{\mathbf{x}}_i\in {\mathbb{C}}^n\setminus \{0\},i=1,2,\dots ,m$$

$${\mathsf{M}}_p=\frac{1}{2\pi \mathrm{i}}{\oint }_{\Gamma }{\left(z-\gamma \right)}^p{V}^{\mathsf{H}}B(zB-A{)}^{-1}BV\mathrm{d}z,p=0,1,2,\dots ,2M-1$$

$$H_M^{<}$$

$$H_M$$

$$S_p=\frac{1}{2\pi \mathrm{i}}{\oint }_{\Gamma }(z-\gamma {)}^p(zB-A{)}^{-1}BV\mathrm{d}z$$

$${\mathsf{M}}_p$$

$$=\frac{1}{2\pi \mathrm{i}}{\oint }_{\Gamma }{\left(z-\gamma \right)}^p\sum _{k=1}^r\left(\frac{V^{\mathsf{H}}B{\mathbf{x}}_k{\mathbf{x}}_k^{\mathsf{H}}BV}{z-{\lambda }_k}\right)\mathrm{d}z$$

$$=\sum _{k=1}^m({\lambda }_k-\gamma {)}^p{\mathcal{V}}_k,$$

$${\mathsf{M}}_p=V^{\mathsf{H}}B{X}_{\Omega }{\Lambda }_{\Omega }^p{X}_{\Omega }^{\mathsf{H}}BV,X_{\Omega }=[{\mathbf{x}}_1,{\mathbf{x}}_2,\dots ,{\mathbf{x}}_m],{\Lambda }_{\Omega }=\mathrm{diag}({\lambda }_1-\gamma ,{\lambda }_2-\gamma ,\dots ,{\lambda }_m-\gamma ).$$

$$H_M$$

$$=\left[\begin{array}{c}{V}^{\mathsf{H}}B{X}_{\Omega }{X}_{\Omega }^{\mathsf{H}}\\ V^{\mathsf{H}}B{X}_{\Omega }{\Lambda }_{\Omega }{X}_{\Omega }^{\mathsf{H}}\\ ⋮\\ V^{\mathsf{H}}B{X}_{\Omega }{\Lambda }_{\Omega }^{M-1}{X}_{\Omega }^{\mathsf{H}}\end{array}\right]B\left[\begin{array}{cccc}{X}_{\Omega }{X}_{\Omega }^{\mathsf{H}}BV& X_{\Omega }{\Lambda }_{\Omega }{X}_{\Omega }^{\mathsf{H}}BV& \cdots & X_{\Omega }{\Lambda }_{\Omega }^{M-1}{X}_{\Omega }^{\mathsf{H}}BV\end{array}\right]$$

$$=S^{\mathsf{H}}BS,$$

$$S$$

$$S_p$$

$$R(S)=\underset{p=0}{\overset{M-1}{⨁}}R(S_p)\subset \mathrm{span}\left\{{\mathbf{x}}_k:{\mathbf{x}}_k^{\mathsf{H}}BV\ne 0,k=1,2,\dots ,m\right\}.$$

$$\dim \left(R(S)\right)\le \dim \left(\mathrm{span}\left\{{\mathbf{x}}_k:{\mathbf{x}}_k^{\mathsf{H}}BV\ne 0,k=1,2,\dots ,m\right\}\right).$$

$$\Gamma =\{z\in \mathbb{C}|z=\gamma +\rho {\mathrm{e}}^{\mathrm{i}\theta },\theta \in \mathbb{R}\},\gamma =\frac{b+a}{2},\rho =\frac{b-a}{2}$$

$$z_j=\gamma +\rho {\mathrm{e}}^{\mathrm{i}{\theta }_j},{\theta }_j=\frac{2j-1}{N}\pi ,j=1,2,\dots ,N.$$

$${\mathsf{M}}_p^{(N)}=\frac{1}{N}\sum _{j=1}^N{\left(\rho {\mathrm{e}}^{\mathrm{i}{\theta }_j}\right)}^{p+1}\left(\sum _{k=1}^r\frac{{\mathcal{V}}_k}{\rho {\mathrm{e}}^{\mathrm{i}{\theta }_j}-({\lambda }_k-\gamma )}\right)=\sum _{k=1}^r{\mathcal{V}}_k\left(\frac{1}{N}\sum _{j=1}^N{\rho }^p{\mathrm{e}}^{\mathrm{i}p{\theta }_j}\frac{\rho {\mathrm{e}}^{\mathrm{i}{\theta }_j}}{\rho {\mathrm{e}}^{\mathrm{i}{\theta }_j}-({\lambda }_k-\gamma )}\right).$$

$$\frac{1}{N}\sum _{j=1}^N{\rho }^p{\mathrm{e}}^{\mathrm{i}p{\theta }_j}\frac{\rho {\mathrm{e}}^{\mathrm{i}{\theta }_j}}{\rho {\mathrm{e}}^{\mathrm{i}{\theta }_j}-({\lambda }_k-\gamma )}$$

$$=\sum _{\ell =0}^{\infty }{\rho }^{p-\ell }({\lambda }_k-\gamma {)}^{\ell }\left(\frac{1}{N}\sum _{j=1}^N{\mathrm{e}}^{\mathrm{i}(p-\ell ){\theta }_j}\right)$$

$$=\sum _{s=0}^{\infty }{\rho }^{-sN}({\lambda }_k-\gamma {)}^{p+sN}=({\lambda }_k-\gamma {)}^p\left(\frac{1}{1-{\left(\frac{{\lambda }_k-\gamma }{\rho }\right)}^N}\right).$$

$$\frac{1}{N}\sum _{j=1}^N{\mathrm{e}}^{\mathrm{i}h{\theta }_j}=\{\begin{array}{l}1(h\in N\mathbb{Z}),\\ 0(\mathrm{otherwise}).\end{array}$$

$$\frac{1}{N}\sum _{j=1}^N{\rho }^p{\mathrm{e}}^{\mathrm{i}p{\theta }_j}\frac{\rho {\mathrm{e}}^{\mathrm{i}{\theta }_j}}{\rho {\mathrm{e}}^{\mathrm{i}{\theta }_j}-({\lambda }_k-\gamma )}$$

$$=\frac{1}{N}\sum _{j=1}^N{\rho }^p{\mathrm{e}}^{\mathrm{i}p{\theta }_j}\left(-\sum _{\ell =0}^{\infty }{\left(\frac{\rho {\mathrm{e}}^{\mathrm{i}{\theta }_j}}{{\lambda }_k-\gamma }\right)}^{\ell +1}\right)$$

$$=\sum _{\ell =0}^{\infty }-{\rho }^{p+\ell +1}({\lambda }_k-\gamma {)}^{-(\ell +1)}\left(\frac{1}{N}\sum _{j=1}^N{\mathrm{e}}^{\mathrm{i}(p+\ell +1){\theta }_j}\right)$$

$$=\sum _{s=1}^{\infty }-{\rho }^{sN}({\lambda }_k-\gamma {)}^{-(sN-p)}$$

$$=({\lambda }_k-\gamma {)}^p\left(\frac{-{\left(\frac{\rho }{{\lambda }_k-\gamma }\right)}^N}{1-{\left(\frac{\rho }{{\lambda }_k-\gamma }\right)}^N}\right).$$

$${\mathsf{M}}_{p,\mathrm{in}}^{(N)}=\sum _{k=1}^m({\lambda }_k-\gamma {)}^p\left(\frac{1}{1-{\left(\frac{{\lambda }_k-\gamma }{\rho }\right)}^N}\right){\mathcal{V}}_k,{\mathsf{M}}_{p,\mathrm{out}}^{(N)}=\sum _{k=m+1}^r({\lambda }_k-\gamma {)}^p\left(\frac{-{\left(\frac{\rho }{{\lambda }_k-\gamma }\right)}^N}{1-{\left(\frac{\rho }{{\lambda }_k-\gamma }\right)}^N}\right){\mathcal{V}}_k$$

$${\alpha }_k=\frac{1}{1-{\left(\frac{{\lambda }_k-\gamma }{\rho }\right)}^N},k=1,2,\dots ,m,{\beta }_k=\frac{-{\left(\frac{\rho }{{\lambda }_k-\gamma }\right)}^N}{1-{\left(\frac{\rho }{{\lambda }_k-\gamma }\right)}^N},k=m+1,m+2,\dots ,r.$$

$${\mathsf{M}}_p^{(N)}={\mathsf{M}}_{p,\mathrm{in}}^{(N)}+{\mathsf{M}}_{p,\mathrm{out}}^{(N)}=\sum _{k=1}^m({\lambda }_k-\gamma {)}^p{\alpha }_k{\mathcal{V}}_k+\sum _{k=m+1}^r({\lambda }_k-\gamma {)}^p{\beta }_k{\mathcal{V}}_k.$$

$${\mathsf{M}}_p-{\mathsf{M}}_p^{(N)}=\sum _{k=1}^m({\lambda }_k-\gamma {)}^p(1-{\alpha }_k){\mathcal{V}}_k-\sum _{k=m+1}^r({\lambda }_k-\gamma {)}^p{\beta }_k{\mathcal{V}}_k.$$

$$H_{M,\mathrm{in}}^{<,(N)}=\left[\begin{array}{cccc}{\mathsf{M}}_{1,\mathrm{in}}^{(N)}& {\mathsf{M}}_{2,\mathrm{in}}^{(N)}& \cdots & {\mathsf{M}}_{M,\mathrm{in}}^{(N)}\\ {\mathsf{M}}_{2,\mathrm{in}}^{(N)}& {\mathsf{M}}_{3,\mathrm{in}}^{(N)}& \cdots & {\mathsf{M}}_{M+1,\mathrm{in}}^{(N)}\\ ⋮& ⋮& \ddots & ⋮\\ {\mathsf{M}}_{M,\mathrm{in}}^{(N)}& {\mathsf{M}}_{M+1,\mathrm{in}}^{(N)}& \cdots & {\mathsf{M}}_{2M-1,\mathrm{in}}^{(N)}\end{array}\right]\in {\mathbb{C}}^{LM\times LM},$$

$$H_{M,\mathrm{in}}^{(N)}=\left[\begin{array}{cccc}{\mathsf{M}}_{0,\mathrm{in}}^{(N)}& {\mathsf{M}}_{1,\mathrm{in}}^{(N)}& \cdots & {\mathsf{M}}_{M-1,\mathrm{in}}^{(N)}\\ {\mathsf{M}}_{1,\mathrm{in}}^{(N)}& {\mathsf{M}}_{2,\mathrm{in}}^{(N)}& \cdots & {\mathsf{M}}_{M,\mathrm{in}}^{(N)}\\ ⋮& ⋮& \ddots & ⋮\\ {\mathsf{M}}_{M-1,\mathrm{in}}^{(N)}& {\mathsf{M}}_{M,\mathrm{in}}^{(N)}& \cdots & {\mathsf{M}}_{2M-2,\mathrm{in}}^{(N)}\end{array}\right]\in {\mathbb{C}}^{LM\times LM}$$

$${\alpha }_k{\mathcal{V}}_k={\alpha }_k{V}^{\mathsf{H}}B{\mathbf{x}}_k{\mathbf{x}}_k^{\mathsf{H}}BV={\alpha }_k\left(c_k{\mathbf{e}}_k\right){\left(c_k{\mathbf{e}}_k\right)}^{\mathsf{H}}=\left({\alpha }_k^{1/2}{c}_k{\mathbf{e}}_k\right){\left({\alpha }_k^{1/2}{c}_k{\mathbf{e}}_k\right)}^{\mathsf{H}}={V^{\prime }}^{\mathsf{H}}B{\mathbf{x}}_k{\mathbf{x}}_k^{\mathsf{H}}B{V}^{\prime }$$

$$\left|{\mathsf{M}}_{p,\mathrm{in}}^{(N)}-{\tilde{\mathsf{M}}}_p^{(N)}\right|\le \left|{\mathsf{M}}_{p,\mathrm{in}}^{(N)}-{\mathsf{M}}_p^{(N)}\right|+\left|{\mathsf{M}}_p^{(N)}-{\tilde{\mathsf{M}}}_p^{(N)}\right|=\left|{\mathsf{M}}_{p,\mathrm{out}}^{(N)}\right|+\left|{\mathsf{M}}_p^{(N)}-{\tilde{\mathsf{M}}}_p^{(N)}\right|$$

$${\mathsf{M}}_{p,\mathrm{in}}^{(N)}\in ⟨{\mathsf{M}}_p^{(N)},\left|{\mathsf{M}}_{p,\mathrm{out}}^{(N)}\right|⟩\subset ⟨{\tilde{\mathsf{M}}}_p^{(N)},\left|{\mathsf{M}}_{p,\mathrm{out}}^{(N)}\right|+\left|{\mathsf{M}}_p^{(N)}-{\tilde{\mathsf{M}}}_p^{(N)}\right|⟩.$$

Full text

Keywords: partial eigenproblem, contour integral, complex moment, verified numerical computations

MSC 65F15, 65G20, 65G50

1 Introduction

We consider verifying the $m$ eigenvalues $\lambda_{i}$, counting multiplicity, of the Hermitian generalized eigenproblem

[TABLE]

in a prescribed interval $\Omega=[a,b]\subset\mathbb{R}$, where $A=A^{\mathsf{H}}\in\mathbb{C}^{n\times n}$, $B=B^{\mathsf{H}}\in\mathbb{C}^{n\times n}$ is positive semidefinite, and the matrix pencil $zB-A$ ($z\in\mathbb{C}$) is regular111See Appends A for the verification of regularity of a matrix pencil., i.e, $\det(zB-A)$ is not identically equal to zero. We call $\lambda_{i}$ an eigenvalue and $\boldsymbol{x}_{i}$ the corresponding eigenvector of the problem (1) or matrix pencil $zB-A$, $z\in\mathbb{C}$ interchangeably. Throughout, we assume that the number of eigenvalues in the interval $\Omega$ is known to be $m$ and there do not exist eigenvalues of (1) at the end points $a$, $b\in\mathbb{R}$. We also denote the eigenvalues of (1) outside $\Omega$ by $\lambda_{i}$ ($i=m+1,m+2,\dots,r$), where $r=\mathrm{rank}\,B$.

There are plenty of previous works for verification methods of eigenvalue problems (see, e.g., [22] and references therein). These previous works, in particular, for symmetric generalized eigenvalue problems are classified into two kinds: some of them aim at rigorously enclosing specific eigenvalues, and others aim at rigorously enclosing all eigenvalues. For the purposes, different approaches have been taken. Behnke [1] used Temple quotients, their generalizations, and the LDLT decomposition to verify specific eigenvalues. Behnke [2] used the variational principle to verify specific eigenvalues. Watanabe et al. [24] used an approximate diagonalization and generalized Rump’s method, avoiding the Cholesky factorization, to verify the eigenvalue with the maximum magnitude. Yamamoto [25] combined the LDLT decomposition with Sylvester’s law of inertia to verify specific eigenvalues. Maruyama et al. [9] used Geršhgorin’s theorem to verify all eigenpairs. Miyajima et al. [14] used the techniques in [12, 13] and combined it with Rump and Wilkinson’s bounds to verify all eigenpairs. See [11] and references therein for the non-Hermitian case.

In this study, we develop a verification method for partial eigenvalues using the block Hankel-type Sakurai–Sugiura method [7], which receives attentions in recent years by virtue of the scalability in parallel and versatility [8]. We shed light on a new perspective of this method. This method uses contour integrals to form complex moment matrices. Their truncation errors for the trapezoidal rule of numerical quadrature were derived by Miyata et al. [15]. Thanks to their work, we derive a numerically computable enclosure of the complex moment. We point out that our verification method works for multiple eigenvalues in the prescribed region and for semidefinite $B$, whereas the previous methods [1, 2, 25, 9, 12, 13, 14] work only for positive definite $B$. In addition, for each quadrature point, a structured linear system of equations arises to solve. The structure enables us to develop an efficient verification technique in case of $B$ being positive definite. Yamamoto [26] and Rump [23] derived componentwise and normwise bounds, respectively, of the error of the approximate solution. See also [22]. These methods need a numerically computed inverse of the coefficient matrix, whereas the proposed technique does not need such a numerical inverse, and instead needs a lower bound of the smallest eigenvalue of $B$.

In the rest of the paper, we use the following notations: For a real matrix $A=(a_{ij})\in\mathbb{R}^{m\times n}$, a nonnegative matrix consisting of entrywise absolute values is denoted by $|A|=(|a_{ij}|)$. For $B=(b_{ij})\in\mathbb{R}^{m\times n}$ and $\alpha\in\mathbb{R}$, the inequality $A<B$ means $a_{ij}<b_{ij}$ holds entrywise and the inequality $A<\alpha$ means $a_{ij}<\alpha$ holds entrywise.

The rest of this paper is organized as follows: In Section 2, we briefly review the block Sakurai–Sugiura Hankel method and its error analysis derived by Miyata et al. [15]. Thanks to this result, we derive a computable rigorous error bound for complex moment in Section 3. We also put several remarks on the implementation of our method in Section 4. In Section 5, we show two numerical examples illustrating the performance of our method. In Section 6, we conclude the paper for discussing potentials of our method for parallel implementation and future directions.

2 Block Hankel-type Sakurai–Sugiura method

We review the block Sakurai–Sugiura Hankel method [7], which is the basis of the proposed method. The block Sakurai–Sugiura Hankel method has parameters such as the block size $L\in\mathbb{N}_{+}$, the order of moment $M\in\mathbb{N}_{+}$, a random matrix $V\in\mathbb{C}^{n\times L}$ whose column vectors consist of a linear combination of all eigenvectors, the basis vectors of the kernel of $B$, say $\mathrm{Ker}\,B$, and the scaling parameters $(\gamma,\rho)\in\mathbb{R}\times\mathbb{R}$ for the eigenvalues. The $p$ th complex moment matrix is given by

[TABLE]

defined on the closed Jordan curve $\Gamma$ through the end points of the interval $\Omega=[a,b]$, where $\mathrm{i}=\sqrt{-1}$ is the imaginary unit and $\pi$ is the circle ratio. Denote the block Hankel matrices consisting of the moments (2) by

[TABLE]

Then, the following theorem show that the block Sakurai–Sugiura Hankel method can compute eigenvalues in a prescribed domain and their corresponding eigenvectors [7, Theorems 5 and 6].

Theorem 2.1**.**

Let an eigenvalue and the corresponding eigenvector of the regular part of the matrix pencil $zH_{M}-H_{M}^{<}$ be denoted by $\theta_{i}$ and $\boldsymbol{u}_{i}$, respectively. Let

[TABLE]

and $S=[S_{0},S_{1},\dots,S_{M-1}]$. If $\mathrm{rank}(H_{M})=m$ holds, then the eigenvalue in $\Gamma$ and the corresponding eigenvector of $\eqref{eq:gevp}$ are given by $\lambda_{i}=\gamma+\theta_{i}$ and $\boldsymbol{x}_{i}=S\boldsymbol{u}_{i}$ ($i=1,2,\dots,m$), respectively.

We remark that the condition $\mathrm{rank}(H_{M})=m$ implies $LM\geq m$.

Next, we give a relationship between the target eigencomponents in the columns of $V$ and the rank of $H_{M}$. Recall the Weierstrass canonical form of the matrix pencil $zB-A$ [3, Proposition 7.8.3]. There exists a nonsingular matrix $X\in\mathbb{C}^{n\times n}$ such that $X^{\mathsf{H}}(zB-A)X=z\mathrm{I}_{0}-\Lambda$, where $\mathrm{I}_{0}$ is a diagonal matrix whose leading $r$ diagonal entries are one and whose trailing $n-r$ diagonal entries are zeros, and $\Lambda$ is a diagonal matrix whose leading $r$ diagonal entries are the eigenvalues of (1) and whose trailing $n-r$ diagonal entries are one. Note that the columns of $X=[\boldsymbol{x}_{1},\boldsymbol{x}_{2},\dots,\boldsymbol{x}_{n}]$ are the appropriately scaled eigenvectors of matrix pencil $zB-A$, where $\boldsymbol{x}_{1}$, $\boldsymbol{x}_{2}$, …$\boldsymbol{x}_{r}\in\mathbb{C}^{n}$ correspond to the eigenvalues $\lambda_{1}$, $\lambda_{2}$, …, $\lambda_{r}\in\mathbb{R}$, respectively, and $\boldsymbol{x}_{k}$, $\boldsymbol{x}_{k+1}$, …$\boldsymbol{x}_{n}\in\mathbb{C}^{n}$ form a basis of $\mathrm{Ker}\,B$. Then, from $X^{\mathsf{H}}BX=\mathrm{I}_{0}$ and the residue theorem, the complex moment (2) is expressed as

[TABLE]

where $\mathcal{V}_{k}=V^{\mathsf{H}}B\boldsymbol{x}_{k}\boldsymbol{x}_{k}^{\mathsf{H}}BV\in\mathbb{C}^{L\times L}$. This is represented by

[TABLE]

Using this form, we have

[TABLE]

where

[TABLE]

Meanwhile, it follows that the range of $S$ satisfies

[TABLE]

This implies

[TABLE]

If we set $V$ such that $\boldsymbol{x}_{k}^{\mathsf{H}}BV=0$ for some $k=1,2,\dots,m$, then $\dim\left(R(S)\right)<m$ and the Hankel matrix $H_{M}$ becomes singular.

In practice, the method uses the $N$-point trapezoidal rule to approximate the complex moment (2) multiplied by $\rho^{-(p+1)}$. We take a domain of integration $\Gamma$ in (2) as the circle

[TABLE]

and approximate the complex moment (2) with the following equi-distributed quadrature points:

[TABLE]

We review the error analysis in [15] to derive a rigorous error bound of the complex moment (2) in Section 3. The trapezoidal rule with the equi-distributed quadrature points (7) approximates the complex moment (2) as

[TABLE]

Since the number of eigenvalues inside $\Gamma$ is $m$, $\left|(\lambda_{k}-\gamma)/\rho\right|<1$ holds for $k=1$, $2$, $\dots$, $m$. Noting the sum of geometric series, the quantity in the parentheses in (8) for $k=1$, $2$, $\dots$, $m$ is written as

[TABLE]

Here, we set $p-\ell=-sN$ ($s=0$, $1$, $2$, $\dots$), due to the property

[TABLE]

The other $r-m$ eigenvalues $\lambda_{k}$ ($k=m+1$, $m+2$, $\dots$, $r$) outside the domain $\Gamma$ satisfy the inequalities $\left|\rho/(\lambda_{k}-\gamma)\right|<1$. Noting the sum of geometric series, the quantity in the parentheses in (8) for $k=m+1,m+2,\dots,r$ is written as

[TABLE]

Here, we set $p+\ell+1=sN$ ($s=1,2,\dots$). It follows from (8), (9), and (10) that the approximated complex moment is split into two parts $\mathsf{M}_{p}^{(N)}=\mathsf{M}_{p,\mathrm{in}}^{(N)}+\mathsf{M}_{p,\mathrm{out}}^{(N)}$, where

[TABLE]

are regarding the inside and outside of $\Gamma$, respectively. Together with (3), we have the truncation error analysis of the $N$-point trapezoidal rule $\mathsf{M}_{p}^{(N)}-\mathsf{M}_{p}$.

3 Error bound of the complex moment

Based on the error analysis in the previous section, we derive a rigorous error bound for each complex moment $\mathsf{M}_{p}$. Let

[TABLE]

Then, the rightmost sides of (9) and (10) become $(\lambda_{k}-\gamma)^{p}\alpha_{k}$ ($k=1,2,\dots,m$) and $(\lambda_{k}-\gamma)^{p}\beta_{k}$ ($k=m+1,\dots,n$), respectively. Then, we simplify the expressions of the approximated complex moment (11)

[TABLE]

The truncation error is given by

[TABLE]

We note that the following identities of the eigenvalues of a Hankel matrix pencil are useful for our verification methods.

Lemma 3.1**.**

Assume that $\mathrm{rank}(H_{M})=m$ holds. Then, the Hankel matrix pencil $zH_{M}-H_{M}^{<}$ consisting of $\mathsf{M}_{p}$ and the Hankel matrix pencil $zH_{M,\mathrm{in}}^{(N)}-H_{M,\mathrm{in}}^{<,(N)}$ with

[TABLE]

consisting of $\mathsf{M}_{p,\mathrm{in}}^{(N)}$ have the same eigenvalues.

Proof.

Let $V=[\boldsymbol{v}_{1},\boldsymbol{v}_{2},\dots,\boldsymbol{v}_{L}]$, $\boldsymbol{v}_{i}=\sum_{j=1}^{n}c_{j}\boldsymbol{x}_{j}$ and $V^{\prime}=[\boldsymbol{v}^{\prime}_{1},\boldsymbol{v}^{\prime}_{2},\dots,\boldsymbol{v}^{\prime}_{L}]$, $\boldsymbol{v}^{\prime}_{i}=\sum_{j=1}^{n}\alpha_{j}^{1/2}c_{j}\boldsymbol{x}_{j}$. Then, we have the equalities

[TABLE]

for $k=1$, $2$, …, $m$. Since Theorem 2.1 holds irrespective of the scaling regarding the eigenvectors in the columns of $V$, the lemma holds. ∎

Hence, we derive an enclosure of $\mathsf{M}_{p,\mathrm{in}}^{(N)}$ instead of an enclosure of $\mathsf{M}_{p}$. We can enclose $\mathsf{M}_{p,\mathrm{in}}^{(N)}$ by using the quantity $|\mathsf{M}_{p,\mathrm{out}}^{(N)}|$ and computing the truncated complex moment $\mathsf{M}_{p}^{(N)}$ with interval arithmetic. Let us denote a numerical approximation of $\mathsf{M}_{p}^{(N)}$ by $\tilde{\mathsf{M}}_{p}^{(N)}$. Hereafter, we denote a numerically computed quantity that may suffer from rounding errors with a tilde. Then, it follows from $\mathsf{M}_{p}^{(N)}-\mathsf{M}_{p,\mathrm{in}}^{(N)}=\mathsf{M}_{p,\mathrm{out}}^{(N)}$ that the inequality

[TABLE]

holds. Let us denote the interval matrix with radius $r\in\mathbb{R}_{+}^{L\times L}$ centered at $c\in\mathbb{C}^{L\times L}$ by $\langle c,r\rangle$. To sum up the above discussion, we have the following theorem:

Theorem 3.2**.**

The computable rigorous enclosure of $\mathsf{M}_{p,\mathrm{in}}^{(N)}$ is given by

[TABLE]

The proof is already completed by the above discussions. We can enclose $\left|\mathsf{M}_{p}^{(N)}-\tilde{\mathsf{M}}_{p}^{(N)}\right|$ using standard verification methods using interval arithmetic, whereas the complex moment $\mathsf{M}_{p,\mathrm{out}}^{(N)}$ regarding the outside of $\Gamma$ is bounded as follows:

Theorem 3.3**.**

Let $V\in\mathbb{C}^{n\times L}$ be an arbitrary matrix such that

[TABLE]

where the columns of $X_{0}=\left[\boldsymbol{x}_{r+1},\boldsymbol{x}_{r+2},\dots,\boldsymbol{x}_{n}\right]$ form a basis of $\mathrm{Ker}\,B$ ($r=\mathrm{rank}(B)$), the columns of $X_{1}=\left[\boldsymbol{x}_{1},\boldsymbol{x}_{2},\dots,\boldsymbol{x}_{r}\right]$ form a basis of $\mathrm{Ker}\,B^{\perp}$, $C_{0}\in\mathbb{C}^{(n-r)\times L}$, and $C_{1}\in\mathbb{C}^{r\times L}$ has at least one nonzero entry in each column and each row. Suppose $N>2M-1\geq p$ and that $\hat{\lambda}$ satisfies $\left|\hat{\lambda}-\gamma\right|=\min_{k=m+1,m+2,\dots,r}\left|\lambda_{k}-\gamma\right|$. Then, the complex moment (11) is bounded above as

[TABLE]

for $p=0$, $1$, …, $2M-1$, where $\|\cdot\|_{\mathsf{F}}$ denotes the Frobenius norm.

Proof.

Regarding the fraction factor in (11) as the geometric series, we have

[TABLE]

Note that the property $BV=BX_{1}C_{1}+BX_{0}C_{0}=BX_{1}C_{1}$ gives

[TABLE]

Hence, we have

[TABLE]

where $\boldsymbol{e}_{k}$ is the $k$ th standard basis vector of $\mathbb{R}^{n}$, i.e., the $k$ th entry is one and the remaining entries are zero. Therefore, we obtain (13). ∎

4 Implementation

In this section, we present an implementation of the block Sakurai–Sugiura Hankel-based method for numerically verifying the partial eigenvalues $\lambda_{i}\in\Omega$, $i=1$, $2$, …, $m$. Suppose that the number of the eigenvalues in $\Gamma$ is $m$. We set $L$ and $M$ such that $m=LM$. Note that if $m$ is a prime number, either $L$ or $M$ must be one and the other must be $m$. To rigorously enclose the eigenvalues, we verify each block $\mathsf{M}_{p,\mathrm{in}}^{(N)}$ of the block Hankel matrices by using Theorem 3.2, and then apply the verified eigenvalue computation methods [10, 19] to the small eigenproblem of regular Hankel matrix pencil consisting of $\mathsf{M}_{p,\mathrm{in}}^{(N)}$. The matrix $\left|\mathsf{M}_{p,\mathrm{out}}^{(N)}\right|$ in (12) can be bounded by using (13). The number of quadrature points can be automatically determined from the error bound (13) by

[TABLE]

where $\delta$ denotes the tolerance of quadrature error. Hence, there is a trade-off between the accuracy for the quadrature and the central processing unit (CPU) time.

The matrix $\left|\mathsf{M}_{p}^{(N)}-\tilde{\mathsf{M}}_{p}^{(N)}\right|$ in (12) can be also bounded by evaluating the numerical error. To rigorously bound the numerical error, we need verification of a numerical solution of the linear system with multiple right-hand sides, that is $(z_{j}B-A)Y_{j}=BV$, which comes from

[TABLE]

The enclosure of $Y_{j}^{*}$ can be obtained by standard verification methods, e.g., [23], whereas we consider efficiently enclosing the solution $Y_{j}^{*}$ for positive definite $B$.

Theorem 4.1**.**

Let $A$ be a Hermitian matrix and $B$ a Hermitian positive definite matrix. The quadrature points $z_{j}$, $j=1$, $2$, …, $N$ are defined as in (7). Denote the $i$ th entries of the solution $\boldsymbol{y}^{\ast}=(z_{j}B-A)^{-1}\boldsymbol{b}$ and an approximate solution $\tilde{\boldsymbol{y}}$ of $(z_{j}B-A)\boldsymbol{y}=\boldsymbol{b}$ by $\tilde{{y}}_{i}$ and $y_{i}^{\ast}$, respectively. If we denote the residual by $\tilde{\boldsymbol{r}}=\boldsymbol{b}-(z_{j}B-A)\tilde{\boldsymbol{y}}$, then the error $\tilde{\boldsymbol{y}}-\boldsymbol{y}^{\ast}$ satisfies

[TABLE]

for all $i=1$, $2$, …, $n$, where $\lambda_{\min}(\cdot)$ is the smallest eigenvalue of a matrix and $\|\cdot\|_{2}$ denotes the Euclidean norm.

Proof.

Denote the square root of $B$ by $B^{1/2}$. Then, for all $i=1$, $2$, …, $n$ we have

[TABLE]

The bound $\|(z_{j}\mathrm{I}-B^{-1/2}AB^{-1/2})^{-1}\|_{2}\leq(\operatorname{Im}z_{j})^{-1}$ can be geometrically interpreted as in Figure 1. Namely, the distance from the quadrature point $z_{j}$ to the nearest eigenvalue of $B^{-1/2}AB^{-1/2}$ is bounded below by the absolute value of the imaginary part of $z_{j}$.

∎

Note that $z_{j}B-A$ is nonsingular for $j=1$, $2$, …, $N$, since $z_{j}$ is not in the real axis (7). Hence, we do not need to verify the regularity of the coefficient matrix $z_{j}B-A$ such as in [23]. In addition, the bound (15) can be efficiently evaluated for sparse $A$ and $B$. On the other hand, the bound (15) shows that, if $\lambda_{\mathrm{min}}(B)$ is very small, the verification of $\tilde{\boldsymbol{y}}_{j}$ will be loose and the subsequent verification may fail. This indicates that Theorem 4.1 works well for well-conditioned $B$. For ill-conditioned $B$, applying iterative refinements with multi-precision arithmetics [18] to the linear system will potentially remedy the bound (15). Furthermore, if each entry of $z_{j}B-A$ and $\boldsymbol{b}$ is not wide interval, one can use a staggered correction [17, Section 4.3]. That is,

[TABLE]

where $\tilde{\boldsymbol{d}}$ solves $(z_{j}B-A)\tilde{\boldsymbol{d}}\approx\tilde{\boldsymbol{r}}$ in a numerical (non-rigorous) sense and $\tilde{d}_{i}$ denotes the $i$ th entry of $\tilde{\boldsymbol{d}}$. This technique is expected to give sharper error bounds than (15) in Theorem 4.1.

We summarize the above procedures in Algorithm 1. In this implementation, we scale the target interval $\Omega$ into $[-1,1]$ by $A^{\prime}=\frac{1}{\rho}\left(A-\gamma B\right)$ and compute the eigenvalues of $A^{\prime}\boldsymbol{x}=\lambda^{\prime}B\boldsymbol{x}$ for simplicity. Here, we denote interval quantities with squares brackets.

The verification in line 4 of Algorithm 1 can be done by, e.g., the following steps:

Compute a numerical approximation $\tilde{\lambda}$ of $\hat{\lambda}$ (defined in Section 3) using MATLAB function eigs. 2. 2.

Set $c\in(0,1)$ such that $1<c|\tilde{\lambda}|$. 3. 3.

Verify regularity of the interval matrix $[A]-[1,c|\tilde{\lambda}|]B$ by using INTLAB function isregular. 4. 4.

Adopt $c|\tilde{\lambda}|$ as the lower bound of $|\hat{\lambda}|$.

5 Numerical examples

To illustrate effectiveness of the proposed method, we show three numerical examples (two artificially generated eigenproblems and one practical eigenproblem). In first and third examples, we compared the proposed method with INTLAB’s function verifyeig in terms of the CPU time. The second example was set for illustrating the performance of the proposed method under the case that the matrix $B$ is positive semidefinite or ill-conditioned. All computations were carried out on Ubuntu 16.04, Intel(R) Xeon(R) Gold 6128 CPU @ 3.40 gigahertz (GHz) with 12 cores, 256 gigabytes (GB) random-access memory (RAM). All programs were coded and run in MATLAB R2018a for double precision floating operation arithmetic with unit roundoff $2^{-53}\simeq 1.1\cdot 10^{-16}$ and with INTLAB version 10.2 [20]. The matrix $V\in\mathbb{R}^{n\times L}$ was generated by using built-in MATLAB function randn. The tolerance of quadrature error was $\delta=10^{-15}$. We determined the smallest $N$ that satisfies (14). Note again that the number of eigenvalues in the interval is given in advance.

In this example, numerically computed solutions of linear systems $(z_{j}B-A)Y_{j}=BV$ were obtained by using MATLAB function mldivide. The eigenvalues of $H_{M}^{<,\text{in}}\boldsymbol{y}=\lambda^{\prime}H_{M}^{\text{in}}$ in line 9 of Algorithm 1 were verified by using INTLAB function verifyeig.

Artificially generated eigenproblems 1

The test matrix pencil $zB-A$ used was given by

[TABLE]

where $\mathrm{tridiag}(\cdot,\cdot,\cdot)$ denotes the tridiagonal Toeplitz matrix consisting of a triplet and the value of $b_{i}$ normally distributes with mean $1$ and variance $10^{-7}$. The generalized eigenproblem of matrix pencil (16) models harmonic oscillators consisting of mass points and springs. In particular, the matrix pencil (16) arises from an equation of motion of mass points in one dimension. Let $u_{i}(t)$ be the displacement of the $i$ th point from the equilibrium of spring $i$ at time $t$ with mass $b_{i}$ and connected with two springs with stiffnesses $k_{i}=k_{i+1}=1$. Then, we have the equation $i=1$, $2$, …, $n$

[TABLE]

Suppose that the mass point has a simple harmonic oscillation $u_{i}(t)=x_{i}\sin(wt+\phi)$, where $w$ is the angular rate, $\phi$ is the phase, and the homogeneous Dirichlet boundary condition $u_{0}(t)=u_{n+1}(t)=0$ is imposed. Then, we have the eigenproblem $A\boldsymbol{x}=\omega^{2}B\boldsymbol{x}$, where $\boldsymbol{x}=[-x_{1},-x_{2},\dots,-x_{n}]^{\mathsf{T}}$.

The verification targets were four eigenvalues near $2$ for $n=2^{\ell}$, $\ell=5$, $6$, …, $20$ of matrix pencil (16). We set the parameters $L=2$ and $M=2$. It is well-known that the eigenvalue of $A$ is given by $\lambda_{i}(A)=2-2\cos(i\pi/(n+1))$ for $i=1,2,\dots,n$. Perturbation theory of Hermitian generalized eigenproblems [16, Theorem 8.3] gives the following bound between $\lambda_{i}$ and $\lambda_{i}(A)$:

[TABLE]

where $\Delta B=I-B$. Then, we derived the lower bound of $|\hat{\lambda}|$ using the eigenvalue $\lambda_{i}(A)$ with its bound (17).

Figure 2 shows the CPU time of the proposed method (Algorithm 1) and a standard method for verifying specific eigenvalues in MATLAB (build-in MATLAB function eigs for the solution of the eigenproblem and INTLAB function verifyeig for eigenvalue verification).

As shown in Figure 2, the efficient verification technique based on Theorem 4.1 achieved a substantial improvement of the proposed method in the CPU time, and the proposed method using the technique based on Theorem 4.1 was faster than the standard method when the size of matrix $n$ is larger than $2^{10}$. Furthermore, due to the limit of RAM, the standard method did not run for $\ell>16$. The proposed method tended to be more effective, as the size of the matrix $n$ becomes large and sparse. On the other hand, the proposed method diminished more than verifyeig in terms of the error bounds. Table 1 gives the verified eigenvalues for the proposed method for each $\ell$. For each $\ell$, the digits in single lines are the same as those of the exact eigenvalues, whereas the digits in double lines denote the supremum and infimum of the exact eigenvalues. Table 1 shows that the proposed method succeeded in verifying the eigenvalues at least 5 digits up to $\ell=20$. For example, for $n=2^{10}$, verifyeig displayed correct 13 digits of the target eigenvalues

[TABLE]

This is mainly due to an overestimation of the error $\tilde{Y}_{j}-Y_{j}$ and in particular $\|(z_{j}B-A)^{-1}\|_{2}$ (see Theorem 4.1). In addition, we remark that this example (16) is very ideal to show the effectiveness of the proposed method, thanks to the sparsity of $A$ and $B$ and the simple structure of $B$.

Artificially generated eigenproblems 2

Another test matrix pencil $zB-A$ was considered for second numerical example, which is defined by

[TABLE]

where “$\mathrm{pentadiag}$” denotes the pentadiagonal Toeplitz matrix. We changed $b_{100}$ as [math],$10^{-16}$, $10^{-15}$, …, $10^{0}$ for illustrating the performance of our method under the case that $B$ is positive semidefinite or ill-conditioned.

We considered six ($m=6$) eigenvalues in $\Omega=[0.95,1.05]$. We set the parameters $L=3$, $M=2$. For the scaled eigenproblem, we verified $|\hat{\lambda}|>1.36$ by using INTLAB’s function isregular.

Figure 3 shows a transition of verified partial eigenvalues with respect to $b_{100}$ entry. The six target eigenvalues were plotted in Figure 3 (a). Changing $b_{100}$ entry, these values slightly moves between $b_{100}=1$ and $10^{-2}$. Our proposed method succeeded in including these eigenvalues with the radius up to $10^{-9}$ as shown in Figure 3 (b). This result implies that our proposed algorithm works well in the case of the martix $B$ being semidefinite or ill-conditioned. Finally, we remark that Theorem 4.1 cannot work in this case because $\lambda_{\mathrm{min}}(B)^{-1}$ becomes very large or infinity. One should use INTLAB’s function verifylss or another verification methods for linear systems.

Practical eigenproblems

Finally, we consider a practical eigenproblem in quantum mechanics. The verification targets are 52 eigenvalues in the interval $[-0.530,-0.425]$ of the Hermitian generalized eigenvalue problem for VCNT900 [4, 5, 6], which is associated with a vibrating carbon nanotube within a supercell with spd orbitals. Both matrices $A$ and $B$ have nonzero density 42.8% and are not sparse. Figure 4 shows the distribution of the 52 eigenvalues and the outer eigenvalues nearest to $[a,b]$. To verify the lower bound of $\lambda_{\mathrm{min}}(B)$, we used Rump’s method [21] using the INTLAB function isspd. That is, we firstly computed an approximate smallest eigenvalue of $B$ (e.g., by built-in MATLAB function eigs), say $\tilde{\lambda}_{\mathrm{min}}(B)$. We secondly checked the positive definiteness of $B-c\tilde{\lambda}_{\mathrm{min}}(B)I$ using isspd for a certain $c\in(0,1)$. If the matrix is positive definite, then we adopt $c\tilde{\lambda}_{\mathrm{min}}(B)$ as the desired lower bound of $\lambda_{\mathrm{min}}(B)$. Furthermore, for the scaled eigenproblem, we verified $|\hat{\lambda}|>1.19$ by using INTLAB’s function isregular. Execution time of this part is about 120 seconds because of our naive implementation. Indeed, there is a room to improve this part. For example, we can use an efficient technique given in [25], which is based on Sylvester’s law of inertia, to verify non-existence of the eigenvalues in the prescribed interval.

The proposed method based on Theorem 4.1 successfully verified 37 of 52 eigenvalues in 7.9 seconds and failed to obtain the inclusion of the rest 15 eigenvalues. This is due to an overestimation of the entries of $\tilde{Y}_{j}-Y_{j}$. When using a verification method in INTALB (so-called backslash ‘\’) for the linear systems, the proposed method successfully verified all 52 eigenvalues in 36.0 seconds. The standard method (eigs+verifyeig) also succeeded in verifying all 52 eigenvalues in 5.2 seconds, since the sizes of the matrices are not so large. Although the most expensive part in Algorithm 1 is the verification of $\tilde{Y}_{j}$, we note that this can be done in parallel for all $j=1$, $2$, …, $N$.

6 Conclusions

We proposed a verified computation method for partial eigenvalues of a Hermitian generalized eigenproblem. A contour integral-type eigensolver, the block Sakurai–Sugiura Hankel method, reduces a given eigenproblem into a generalized eigenproblem of block Hankel matrices consisting of complex moments. The error of the complex moment can split into the error of numerical quadrature and the rounding error of numerical computations, which should be controlled rigorously. We derived a truncation error bound of the quadrature and developed an efficient technique to verify the rounding error in the numerical solution of a linear system arising from each quadrature point. Numerical experiments showed that, as the sizes of matrices become large and sparse, the proposed method outperforms a standard method on artificially generated eigenproblems. It is also shown that proposed methods is applicable for practical eigenproblems. We left the issue of how to verify the number of the eigenvalues in the prescribed interval. Finally, we remark that the proposed method will be potentially efficient in parallel. This is one of future directions for this research.

Acknowledgements

We would like to thank Prof. Yusaku Yamamoto for letting us know the work [15]. We also would like to thank Prof. Katsuhisa Ozaki for helpful discussions of parallel implementations. This work was supported in part by the Faculty of Engineering, Information and Systems, University of Tsukuba. The work of the first author was supported in part by JST/ACT-I (No. JPMJPR16U6) and JSPS KAKENHI Grant Numbers 17K12690 and 18H03250. The work of the second author was supported in part by JSPS KAKENHI Grant Number 16K17639 and Hattori Hokokai Foundation. The work of the third author was supported in part by JSPS KAKENHI Grant Number 18K13453.

Appendix A Regularity of a matrix pencil

Consider verifying the regularity of matrix pencil $zB-A$ for Hermitian $A$ and Hermitian positive semidefinite $B$. Recall that a matrix pencil $zB-A$ is said to be singular for square matrices $A$ and $B$ if $\det(zB-A)$ is identically equal to zero; regular otherwise. The matrix pencil $zB-A$ is regular if and only if $\mathrm{Ker}(\left[\begin{smallmatrix}A\\ B\end{smallmatrix}\right])=\{\boldsymbol{0}\}$ [3, Proposition 7.8.4]. Hence, we can guarantee the regularity of matrix pencil $zB-A$ by proving positive definiteness of $[B,A][B,A]^{\mathsf{H}}$ by using the INTLAB function isspd in [21].