Shifted Lanczos method for
quadratic forms with Hermitian matrix resolvents

Keiichi Morikuni Faculty of Engineering, Information and Systems, University of Tsukuba, Japan. Email: [email protected]. The work was supported in part by the Japan Society for the Promotion of Science (Grants-in-Aid for Young Scientists (B) 16K17639) and Hattori Hokokai Foundation.

Abstract

Quadratic forms of Hermitian matrix resolvents involve the solutions of shifted linear systems. Efficient iterative solutions use the shift-invariance property of Krylov subspaces. The Hermitian Lanczos method reduces a given vector and matrix to a Jacobi matrix (real symmetric tridiagonal matrix with positive super and sub-diagonal entries) and approximates the quadratic form using the Jacobi matrix. This study develops a shifted Lanczos method that deals directly with the Hermitian matrix resolvent. We derive a matrix representation of a linear operator that approximates the resolvent by solving a Vorobyev moment problem associated with the shifted Lanczos method. We show that an entry of the Jacobi matrix resolvent can approximate the quadratic form, matching the moments. We give a sufficient condition such that the method does not break down, an error bound, and error estimates. Numerical experiments on matrices drawn from real-world applications compare the proposed method with previous methods and show that the proposed method outperforms well-established methods in solving some problems.

1 Introduction

Consider the computation of $m$ quadratic forms

\boldsymbol{v}^{\mathsf{H}}(z_{i}\mathrm{I}-A)^{-1}\boldsymbol{v},\quad i=1,2,\dots,m,

(1.1)

where $\boldsymbol{v}^{\mathsf{H}}$ denotes the complex conjugate transpose of a vector $\boldsymbol{v}\in\mathbb{C}^{n}$ , $A\in\mathbb{C}^{n\times n}$ is a Hermitian matrix that may be indefinite, and $z_{i}\in\mathbb{C}$ . Here, $z_{i}\mathrm{I}-A$ is assumed to be invertible. If $z_{i}$ is not real, then $z_{i}\mathrm{I}-A$ is not Hermitian.

A straightforward approach to the quadratic form (1.1) is to solve the shifted linear systems

\displaystyle(z_{i}\mathrm{I}-A)\boldsymbol{x}^{(i)}=\boldsymbol{v}

(1.2)

and compute $\boldsymbol{v}^{\mathsf{H}}\boldsymbol{x}^{(i)}$ for $i=1$ , $2$ , $\dots$ , $m$ . The development of efficient solutions for shifted linear systems (1.2) with real symmetric $A$ has been on-going for two decades; these solutions include shifted Krylov subspace methods such as variants of the conjugate orthogonal conjugate gradient (COCG) method [45], a version of the quasi-minimal residual (QMR_SYM) method [5], and the conjugate orthogonal conjugate residual (COCR) method [38] proposed in [43], [37], and [39], respectively. These methods use the shift-invariance property of the Krylov subspace $\mathcal{K}_{k}(z_{i}\mathrm{I}-A,\boldsymbol{v})=\mathcal{K}_{k}(A,\boldsymbol{v})=\mathrm{span}\{\boldsymbol{v}_{1},A\boldsymbol{v}_{1},\dots,A^{k-1}\boldsymbol{v}_{1}\}$ and complex symmetry of $z\mathrm{I}-A$ for efficient formulations. A shifted Krylov subspace method forms basis vectors of the Krylov subspace $\mathcal{K}_{k}(A,\boldsymbol{v})$ for a shifted matrix $z\mathrm{I}-A$ with a particular shift $z$ by using a short-term recurrence and determines an iterate for the linear system $(z\mathrm{I}-A)\boldsymbol{x}_{z}=\boldsymbol{v}$ and associated iterates for other shifted linear systems $(z_{i}\mathrm{I}-A)\boldsymbol{x}^{(i)}=\boldsymbol{v}$ simultaneously for different values of $z_{i}$ , $i=1$ , $2$ , $\dots$ , $m$ , without additional matrix-vector products. Typically, shifted Krylov subspace methods take the iterative residual vectors of the shifted linear systems to be collinear to that of the seed linear system. Here, we call a representative linear system $(z\mathrm{I}-A)\boldsymbol{x}_{z}=\boldsymbol{v}$ a seed linear system. These methods may suffer from breakdown, i.e. division by zero, although this rarely occurs in practice. Previous studies on shifted CG and MINRES methods [7, 26] focused on real shifts. The technique in [26] improves the convergence of the methods with preconditioning under the assumption that the factorization or the use of a direct solver for the shifted symmetric matrix is performed efficiently. Previous studies on a shifted Lanczos method [6] for a Hermitian positive definite $A$ with a complex shift $z_{i}$ derived its error bound and estimates. Extensions of the MINRES method [30] to the shifted linear system (1.2) work with the Hermitian $A$ and complex shift $z_{i}$ [11, 36]. See [19] for more shifted Krylov subspace methods. The Padé approximation via the Lanczos process (PVL) [4] has attracted interest for the case $|z|>\rho(A)$ , where $\rho(\cdot)$ is the spectral radius of a square matrix, whereas our interest includes the contrasting case $|z|<\rho(A)$ .

This study focuses on exploiting the advantage of a Krylov subspace method in efficiently approximating the Hermitian matrix resolvent and directly computing the quadratic form (1.1) without solving the shifted linear systems (1.2). The key feature of our approach is the development of a shifted Lanczos method. The Hermitian Lanczos method projects an original model represented by $A$ and $\boldsymbol{v}$ to a lower-order model and matches the lowest-order moments of the original model with those of the reduced model. We show that the shifted Lanczos method retains such a property. The Vorobyev moment problem [44, 2, 40, 24] enables a concise derivation of thisto concisely derive the method. This problem gives a matrix representation of a linear operator that represents the reduced model to approximate the resolvent $(z\mathrm{I}-A)^{-1}$ . We show that the $(1,1)$ entry of the Jacobi matrix resolvent can approximate the quadratic form (1.1) with eight additional operations to the Lanczos method for each shift by using a recursive formula. Moreover, we give a sufficient condition such that the proposed method does not break down. Breakdown may occur in the COCG, COCR, and QMR_SYM methods, although this is rarely seen. Furthermore, we drive an error bound and develop practical error estimates of the shifted Lanczos method.

Practical applications that can benefit from this development include chemistry and physics [24, Section 3.9], eigensolvers using complex moments [35], the computation of Green’s function for a many-particle Hamiltonian [46], the stochastic estimation of the number of eigenvalues [25], samplers for determinantal point processes [23] (see also references therein), the approximation of Markov chains in Bayesian sampling [20], and computational quantum physics [19]. An extension of the single-vector case $\boldsymbol{v}\in\mathbb{R}^{n}$ in (1.1) for real $A\in\mathbb{R}^{n\times n}$ to the multiple-vector case $\boldsymbol{v}_{1}$ , $\boldsymbol{v}_{2}$ , $\dots$ , $\boldsymbol{v}_{\ell}\in\mathbb{R}^{n}$ , namely,

\displaystyle V^{\top}(z_{i}\mathrm{I}-A)^{-1}V,\quad V=[\boldsymbol{v}_{1},\boldsymbol{v}_{2},\dots,\boldsymbol{v}_{\ell}]\in\mathbb{R}^{n\times\ell},\quad i=1,2,\dots,m

(1.3)

can be reduced to the solutions of bilinear forms ${\boldsymbol{v}_{p}}^{\top}(z_{i}\mathrm{I}-A)\boldsymbol{v}_{q}$ , $p$ , $q=1$ , $2$ , $\dots$ , $\ell$ , where $\boldsymbol{v}_{p}$ is the $p$ th column of $V$ . The bilinear form $\boldsymbol{v}_{p}(z_{i}\mathrm{I}-A)^{-1}\boldsymbol{v}_{q}$ can be further reduced to the quadratic form

\displaystyle\boldsymbol{v}_{p}(z_{i}\mathrm{I}-A)^{-1}\boldsymbol{v}_{q}=\frac{1}{4}[\boldsymbol{s}^{\top}(z_{i}\mathrm{I}-A)^{-1}\boldsymbol{s}-\boldsymbol{t}^{\top}(z_{i}\mathrm{I}-A)^{-1}\boldsymbol{t}],

(1.4)

where $\boldsymbol{s}=\boldsymbol{v}_{p}+\boldsymbol{v}_{q}$ and $\boldsymbol{t}=\boldsymbol{v}_{p}-\boldsymbol{v}_{q}$ (cf. [10, p. 114]). This type of problem arises in the analysis of dynamical systems [1].

The rest of this paper is organized as follows. In Section 2, we review the Lanczos method and its moment-matching property. In Section 3, we describe a shifted Lanczos method for computing the quadratic forms (1.1), its moment-matching property, and implementation; discuss related methods; and give a sufficient breakdown-free condition, error bound, and error estimates. In Section 4, we present the results of numerical experiments in which the shifted Lanczos method is compared with previous methods and illustrate the developed error estimates. In Section 5, we conclude the paper.

2 Lanczos method

We review the Lanczos method [22] for Hermitian $A\in\mathbb{C}^{n\times n}$ and $\boldsymbol{v}\in\mathbb{C}^{n}$ . Algorithm 2.1 gives the procedures of the Lanczos method.

Algorithm 2.1 Lanczos method

A\in\mathbb{C}^{n\times n}

\boldsymbol{v}\in\mathbb{C}^{n}

\boldsymbol{v}_{i}\in\mathbb{C}^{n}

\alpha_{i}\in\mathbb{R}

\beta_{i}\in\mathbb{R}

i=1

2

\dots

k

\boldsymbol{v}_{1}=\boldsymbol{v}/\|\boldsymbol{v}\|

\boldsymbol{u}=A\boldsymbol{v}_{1}

\alpha_{1}=\boldsymbol{u}^{\mathsf{H}}\boldsymbol{v}_{1}

2: for

k=1,2,\dots

\boldsymbol{u}=\boldsymbol{u}-\alpha_{k}\boldsymbol{v}_{k}

\beta_{k}=\|\boldsymbol{u}\|

4: if

\beta_{k}=0

then break

\boldsymbol{v}_{k+1}=(\beta_{k})^{-1}\boldsymbol{u}

\boldsymbol{u}=A\boldsymbol{v}_{k+1}-\beta_{k}\boldsymbol{v}_{k}

\alpha_{k+1}=\boldsymbol{u}^{\mathsf{H}}\boldsymbol{v}_{k+1}

6: end for

Here, $\|\cdot\|$ denotes the Euclidean norm. Denote the Lanczos decomposition of $A$ by

AV_{k}=V_{k+1}T_{k+1,k},

(2.1)

where the columns of $V_{k}=[\boldsymbol{v}_{1},\boldsymbol{v}_{2},\dots,\boldsymbol{v}_{k}]\in\mathbb{C}^{n\times k}$ form an orthonormal basis of the Krylov subspace $\mathcal{K}_{k}(A,\boldsymbol{v}_{1})=\mathrm{span}\{\boldsymbol{v}_{1},A\boldsymbol{v}_{1},\dots,A^{k-1}\boldsymbol{v}_{1}\}$ and

	$\displaystyle T_{k+1,k}$	$\displaystyle=\begin{bmatrix}\alpha_{1}&\beta_{1}&&&\smash{\hbox{\bg 0}}\\ \beta_{1}&\alpha_{2}&\beta_{2}\\ &\beta_{2}&\ddots&\ddots\\ &&\ddots&&\beta_{k-1}\\ &&&\beta_{k-1}&\alpha_{k}\\ \smash{\lower 4.30554pt\hbox{\bg 0}}&&&&\beta_{k}\end{bmatrix}$		(2.2)
		$\displaystyle=\begin{bmatrix}T_{k,k}\\ \beta_{k}{\boldsymbol{e}_{k}}^{\top}\end{bmatrix}\in\mathbb{R}^{(k+1)\times k}.$		(2.3)

Here, $\boldsymbol{e}_{k}\in\mathbb{R}^{k}$ is the $k$ th Euclidean basis vector and $T_{k,k}\in\mathbb{R}^{k\times k}$ is the Jacobi matrix (real symmetric tridiagonal matrix with positive super and subdiagonal entries). Then, $V_{k}^{\mathsf{H}}AV_{k}=T_{k,k}$ holds.

2.1 Hamburger moment problem

The Lanczos method projects the original model given by a Hermitian matrix $A\in\mathbb{C}^{n\times n}$ and an initial vector $\boldsymbol{v}\in\mathbb{C}^{n}$ to a lower-order model given by $T_{k,k}\in\mathbb{R}^{k\times k}$ and $\boldsymbol{e}_{1}$ , matching the lowest-order moments of the original model and those of the reduced model, i.e., it approximates a given matrix $A$ via moment matching. A thorough derivation of the well-known moment-matching property (2.12) presented below is given in [24, Chapter 3]. We review this for completeness. Given a sequence of scalars $\xi_{i}$ , $i=0$ , $1$ , $2$ , $\dots$ , $2k-1$ , the problem of finding a nondecreasing real distribution function $w^{(k)}(\lambda)$ , $\lambda\in\mathbb{R}$ , with $k$ points of increase such that the Riemann–Stieltjes integral is equal to the given sequence of scalars

\int_{-\infty}^{\infty}\lambda^{i}\mathrm{d}w^{(k)}(\lambda)=\xi_{i},\quad i=0,1,\dots,2k-1

(2.4)

is called the Hamburger moment problem [12, 13, 14, 15]. The motivation to view the Lanczos method via the Hamburger moment problem instead of the Stieltjes moment problem is to enable working with indefinite Hermitian matrices with a shift. Note that the Stieltjes moment problem goes on the positive real axis and its relation with the Lanczos method on Hermitian positive definite matrices has been well-established [24].

The left-hand side of (2.4) is called the $i$ th moment with respect to the distribution function $w^{(k)}(\lambda)$ . Set another moment

\xi_{i}=\int_{-\infty}^{\infty}\lambda^{i}\mathrm{d}w(\lambda),\quad i=0,1,2,\dots,2k-1,

(2.5)

and the distribution function with $n$ points of increase to

\displaystyle w(\lambda)=\begin{cases}0,\quad&\lambda<\lambda_{1},\\ \sum_{j=1}^{i}w_{j},\quad&\lambda_{i}\leq\lambda<\lambda_{i+1},\quad i=1,2,\dots,n-1,\\ \sum_{j=1}^{n}w_{j}=1,\quad&\lambda_{n}\leq\lambda\end{cases}

(2.6)

associated with weights $w_{j}=(\boldsymbol{v}^{\mathsf{H}}\boldsymbol{u}_{j})^{2}/\|\boldsymbol{v}\|^{2}$ , $j=1$ , $2$ , $\dots$ , $n$ , where $\lambda_{1}<\lambda_{2}<\dots<\lambda_{n}$ are the eigenvalues of $A$ and $\boldsymbol{u}_{i}$ , $i=1$ , $2$ , $\dots$ , $n$ , are the corresponding eigenvectors. Here, for clarity, we assume that the eigenvalues of $A$ are distinct without loss of generality. Thus, the distribution function $w(\lambda)$ is connected with the eigenpairs of $A$ . Then, we can express the moment (2.5) as the Gauss–Christoffel quadrature and the quadratic form

	$\displaystyle\int_{-\infty}^{\infty}\lambda^{i}\mathrm{d}w(\lambda)$	$\displaystyle=\sum_{j=1}^{n}w_{j}\{\lambda_{j}\}^{i}$		(2.7)
		$\displaystyle=\boldsymbol{v}^{\mathsf{H}}A^{i}\boldsymbol{v},\quad i=0,1,2,\dots.$		(2.8)

Thus, the solution of (2.4) is given by

\displaystyle w^{(k)}(\lambda)=\begin{cases}0,\quad&\lambda<\lambda_{1}^{(k)},\\ \sum_{j=1}^{i}w_{j}^{(k)},\quad&\lambda_{i}^{(k)}\leq\lambda<\lambda_{i+1}^{(k)},\quad i=1,2,\dots,k-1,\\ \sum_{j=1}^{k}w_{j}^{(k)}=1,&\lambda_{k}^{(k)}\leq\lambda\end{cases}

(2.9)

associated with weights $w_{j}^{(k)}=({\boldsymbol{e}_{1}}^{\top}\boldsymbol{u}_{j}^{(k)})^{2}$ , $j=1$ , $2$ , $\dots$ , $k$ , where $\lambda_{1}^{(k)}<\lambda_{2}^{(k)}<\dots<\lambda_{k}^{(k)}$ are the eigenvalues of $T_{k,k}$ and $\boldsymbol{u}_{j}^{(k)}$ , $j=1$ , $2$ , $\dots$ , $k$ , are the corresponding eigenvectors. Here, the distribution function $w^{(k)}(\lambda)$ is connected with the eigenpairs of $T_{k,k}$ . Because the Gauss–Christoffel quadrature is exact for polynomials up to degree $2k-1$

	$\displaystyle\int_{-\infty}^{\infty}\lambda^{i}\mathrm{d}w^{(k)}(\lambda)$	$\displaystyle=\sum_{i=1}^{k}w_{j}^{(k)}\{\lambda_{j}^{(k)}\}^{i}$		(2.10)
		$\displaystyle={\boldsymbol{e}_{1}}^{\top}(T_{k,k})^{i}\boldsymbol{e}_{1},\quad i=0,1,\dots,2k-1,$		(2.11)

the first $2k$ moments match

\displaystyle\boldsymbol{v}^{\mathsf{H}}A^{i}\boldsymbol{v}=(\boldsymbol{v}^{\mathsf{H}}\boldsymbol{v}){\boldsymbol{e}_{1}}^{\top}(T_{k,k})^{i}\boldsymbol{e}_{1},\quad i=0,1,\dots,2k-1.

(2.12)

2.2 Model reduction via Vorobyev moment matching

We can state the problem of moment matching in the language of matrices via the Vorobyev moment problem. To derive a linear operator $A_{k}$ that reduces the model order of $A$ , we follow [44] for the derivation (see also [2, 40]). Let

\displaystyle\begin{cases}\boldsymbol{y}_{1}=A\boldsymbol{v},\\ \boldsymbol{y}_{2}=A\boldsymbol{y}_{1}(=A^{2}\boldsymbol{v}),\\ \qquad\vdots\\ \boldsymbol{y}_{k-1}=A\boldsymbol{y}_{k-2}(=A^{k-1}\boldsymbol{v}),\\ \boldsymbol{y}_{k}=A\boldsymbol{y}_{k-1}(=A^{k}\boldsymbol{v}),\end{cases},

(2.13)

for $k=1$ , $2$ , $\dots$ , where $\boldsymbol{v}$ , $\boldsymbol{y}_{1}$ , $\boldsymbol{y}_{2}$ , $\dots$ , $\boldsymbol{y}_{k}$ are assumed to be linearly independent. Then, the Vorobyev moment problem involves determining a sequence of linear operators $A_{k}$ such that

\displaystyle\begin{cases}\boldsymbol{y}_{1}=A_{k}\boldsymbol{v},\\ \boldsymbol{y}_{2}=A_{k}\boldsymbol{y}_{1}(=(A_{k})^{2}\boldsymbol{v}),\\ \qquad\vdots\\ \boldsymbol{y}_{k-1}=A_{k}\boldsymbol{y}_{k-2}(=(A_{k})^{k-1}\boldsymbol{v}),\\ Q_{k}\boldsymbol{y}_{k}=A_{k}\boldsymbol{y}_{k-1}(=(A_{k})^{k}\boldsymbol{v}),\end{cases}

(2.14)

for $k=1$ , $2$ , $\dots$ , where $Q_{k}=V_{k}V_{k}^{\mathsf{H}}$ is the orthogonal projector onto $\mathcal{K}_{k}(A,\boldsymbol{v})$ . A linear operator $A_{k}$ reducing the model order of $A$ is given by

	$\displaystyle A_{k}$	$\displaystyle=Q_{k}AQ_{k}$
		$\displaystyle=V_{k}T_{k,k}{V_{k}}^{\mathsf{H}},$		(2.15)

for $k=1$ , $2$ , $\dots$ , where the sequence $\{A_{k}\}_{k\geq 0}$ is strongly convergent to $A$ [44, Theorem II] (see [2, Section 4.2] for the derivation of (2.15)). Therefore, the first $2k$ moments of the reduced model match those of the original model

	$\displaystyle\boldsymbol{v}^{\mathsf{H}}A^{i}\boldsymbol{v}$	$\displaystyle=\boldsymbol{v}^{\mathsf{H}}(A_{k})^{i}\boldsymbol{v}$
		$\displaystyle=(\boldsymbol{v}^{\mathsf{H}}\boldsymbol{v}){\boldsymbol{e}_{1}}^{\mathsf{H}}(T_{k,k})^{i}\boldsymbol{e}_{1},\quad i=0,1,\dots,2k-1.$		(2.16)

for $k=1$ , $2$ , $\dots$ . We will use this property to derive a shifted Lanczos method in the next section.

3 Shifted Lanczos method

Next, we formulate a shifted Lanczos method to approximate the resolvent $(z_{i}\mathrm{I}-A)^{-1}$ in (1.1). For convenience, we omit subscript $i$ of $z_{i}$ if no confusion can arise. The application of Vorobyev’s method of moments to the shifted matrix $S=z\mathrm{I}-A$ and vector $\boldsymbol{v}$ gives a matrix representation of a linear operator that represents the reduced model to approximate the resolvent $(z\mathrm{I}-A)^{-1}$ . Let

\displaystyle\begin{cases}\boldsymbol{y}_{1}=S\boldsymbol{v},\\ \boldsymbol{y}_{2}=S\boldsymbol{y}_{1}(=S^{2}\boldsymbol{v}),\\ \qquad\vdots\\ \boldsymbol{y}_{k-1}=S\boldsymbol{y}_{k-2}(=S^{k-1}\boldsymbol{v}),\\ \boldsymbol{y}_{k}=S\boldsymbol{y}_{k-1}(=S^{k}\boldsymbol{v}),\end{cases}

(3.1)

\displaystyle\begin{cases}\boldsymbol{y}_{1}=S_{k}\boldsymbol{v},\\ \boldsymbol{y}_{2}=S_{k}\boldsymbol{y}_{1}(=(S_{k})^{2}\boldsymbol{v}),\\ \qquad\vdots\\ \boldsymbol{y}_{k-1}=S_{k}\boldsymbol{y}_{k-2}(=(S_{k})^{k-1}\boldsymbol{v}),\\ Q_{k}\boldsymbol{y}_{k}=S_{k}\boldsymbol{y}_{k-1}(=(S_{k})^{k}\boldsymbol{v}),\end{cases}

(3.2)

for $k=1$ , $2$ , $\dots$ , where $Q_{k}=V_{k}V_{k}^{\mathsf{H}}$ is the orthogonal projector onto $\mathcal{K}_{k}(S,\boldsymbol{v})$ . We first solve the problem for the linear operator $S_{k}$ . Equations (3.2) can be written as

\boldsymbol{y}_{i}=(S_{k})^{i}\boldsymbol{v},\quad i=1,2,\dots,k-1,\quad Q_{k}\boldsymbol{y}_{k}=(S_{k})^{k}\boldsymbol{v}

(3.3)

for $k=1$ , $2$ , $\dots$ . An arbitrary vector $\boldsymbol{u}\in\mathcal{K}_{k}(S,\boldsymbol{v})$ is expanded as

\boldsymbol{u}=\sum_{i=0}^{k-1}a_{i}\boldsymbol{y}_{i},\quad a_{i}\in\mathbb{C},

(3.4)

where $\boldsymbol{y}_{0}=\boldsymbol{v}$ . Multiplying both sides by $S$ gives

S\boldsymbol{u}=\sum_{i=0}^{k-2}a_{i}S^{i+1}\boldsymbol{v}+a_{k-1}\boldsymbol{y}_{k}.

(3.5)

Projecting this onto $\mathcal{K}_{k}(S,\boldsymbol{v})=\mathcal{K}_{k}(A,\boldsymbol{v})$ (shift-invariance property) gives

$\displaystyle Q_{k}S\boldsymbol{u}$	$\displaystyle=\sum_{i=0}^{k-2}a_{i}(S_{k})^{i+1}\boldsymbol{v}+a_{k-1}(S_{k})^{k}\boldsymbol{v}$	(3.6)
	$\displaystyle=\sum_{i=0}^{k-1}a_{i}(S_{k})^{i+1}\boldsymbol{v}$	(3.7)
	$\displaystyle=\sum_{i=0}^{k-1}a_{i}S_{k}\boldsymbol{y}_{i}$	(3.8)
	$\displaystyle=S_{k}\boldsymbol{u}.$	(3.9)

Here, the first equality is due to equations (3.1), (3.2), and

	$\displaystyle Q_{k}(S_{k})^{i+1}\boldsymbol{v}$	$\displaystyle=Q_{k}S^{i+1}\boldsymbol{v}$		(3.10)
		$\displaystyle=S^{i+1}\boldsymbol{v}\in\mathcal{K}_{k}(S,\boldsymbol{v}),\quad i=0,1,\dots,k-2.$		(3.11)

Hence, (3.9) shows that $Q_{k}S=S_{k}$ on $\mathcal{K}_{k}(S,\boldsymbol{v})$ . Because $Q_{k}\boldsymbol{w}\in\mathcal{K}_{k}(S,\boldsymbol{v})$ for any vector $\boldsymbol{w}\in\mathbb{C}^{n}$ , we can obtain the expression

S_{k}=Q_{k}SQ_{k}

(3.12)

by extending the domain to the whole space $\mathbb{C}^{n}$ . Note that the sequence $\{S_{k}\}_{k\geq 0}$ is strongly convergent to $S$ [44, Theorem II]. The expression (3.12) can be obtained from the shifted Lanczos decomposition [6, Lemma 2.1 (i)]

(z\mathrm{I}-A)V_{k}=V_{k+1}\left(z\begin{bmatrix}\mathrm{I}\\ \boldsymbol{0}^{\top}\end{bmatrix}-T_{k+1,k}\right).

(3.13)

Multiplying this by $V_{k}^{\mathsf{H}}$ gives

V_{k}^{\mathsf{H}}SV_{k}=T_{k}^{<}.

(3.14)

This gives the orthogonally projected restriction

	$\displaystyle S_{k}$	$\displaystyle=V_{k}{V_{k}}^{\mathsf{H}}SV_{k}V_{k}^{\mathsf{H}}$		(3.15)
		$\displaystyle=V_{k}T_{k}^{<}V_{k}^{\mathsf{H}}.$		(3.16)

However, this way does not give the insight of strong convergence.

By using the expression $S_{k}$ , we show that the moments of the original model with $S$ and $\boldsymbol{v}$ and those of the reduced model with $T_{k}^{<}=z\mathrm{I}-T_{k,k}$ and $\boldsymbol{e}_{1}$ match. By using the moment-matching property (2.16) and the binomial formula, we have

\displaystyle\boldsymbol{v}^{\mathsf{H}}S^{i}\boldsymbol{v}=\boldsymbol{v}^{\mathsf{H}}(z\mathrm{I}-A_{k})^{i}\boldsymbol{v},\quad i=0,1,\dots,2k-1,

(3.17)

Therefore, the reduced model matches the first $2k$ moments of the original model

	$\displaystyle\boldsymbol{v}^{\mathsf{H}}S^{i}\boldsymbol{v}$	$\displaystyle=\boldsymbol{v}^{\mathsf{H}}(S_{k})^{i}\boldsymbol{v}$		(3.18)
		$\displaystyle=(\boldsymbol{v}^{\mathsf{H}}\boldsymbol{v})\boldsymbol{e}_{1}^{\top}(T_{k}^{<})^{i}\boldsymbol{e}_{1},\quad i=0,1,\dots,2k-1.$		(3.19)

Note that the left-hand side is equal to

\sum_{j=0}^{i}w_{j}(z-\lambda_{j})^{i}=\int_{-\infty}^{\infty}(z-\lambda)^{i}\mathrm{d}w(\lambda)

(3.20)

with distribution function (2.6) and the right-hand side is equal to

\sum_{j=0}^{i}w_{j}(z-\lambda_{j}^{(k)})^{i}=\int_{-\infty}^{\infty}(z-\lambda)^{i}\mathrm{d}w^{(k)}(\lambda)

(3.21)

with distribution function (2.9) (cf. [42, Chapter 15]). Thus, the Hamburger moment problem (Section 2.1) gives these connections that carry over the expressions of moments (2.7) and the exactness of the quadrature (2.10) to the shifted case, whereas Vorobyev’s method of moments translates these connections in the language of matrices.

Consider the approximation of $\boldsymbol{v}^{\mathsf{H}}S^{-1}\boldsymbol{v}$ . Because the matrix representation of the inverse of the reduced-order operator $S_{k}$ restricted onto $\mathcal{K}_{k}(S,\boldsymbol{v})$ [16, p. 79] is given by

{S_{k}}^{-1}=V_{k}(T_{k}^{<})^{-1}{V_{k}}^{\mathsf{H}},

(3.22)

an approximation of $\boldsymbol{v}^{\mathsf{H}}S^{-1}\boldsymbol{v}$ is given by $\boldsymbol{v}^{\mathsf{H}}(S_{k})^{-1}\boldsymbol{v}$ . Therefore, we obtain

\boldsymbol{v}^{\mathsf{H}}(S_{k})^{-1}\boldsymbol{v}=(\boldsymbol{v}^{\mathsf{H}}\!\boldsymbol{v}){\boldsymbol{e}_{1}}^{\top}(T_{k}^{<})^{-1}\boldsymbol{e}_{1}\equiv L_{k}.

(3.23)

3.1 Implementation

The quantity $L_{k}$ in (3.23) is the $(1,1)$ entry of the resolvent $(T_{k}^{<})^{-1}$ of a successively enlarging Jacobi matrix $T_{k,k}$ . An efficient recursive formula for computing such an entry was developed in [10, Section 3.4]. The formula is given by

L_{k+1}=L_{k}+c_{k+1}\pi_{k+1},\quad k=1,2,\dots,

(3.24)

starting with $c_{1}=1$ , $\delta_{1}=z-\alpha_{1}$ , and $\pi_{1}=1/\delta_{1}$ , where

$\displaystyle t_{k}$	$\displaystyle=(\beta_{k})^{2}\pi_{k},$	(3.25)
$\displaystyle\delta_{k+1}$	$\displaystyle=z-\alpha_{k+1}-t_{k},$	(3.26)
$\displaystyle\pi_{k+1}$	$\displaystyle=1/\delta_{k+1},$	(3.27)
$\displaystyle c_{k+1}$	$\displaystyle=c_{k}t\pi_{k}.$	(3.28)

We summarize the procedures for approximating quadratic forms (1.1) in Algorithm 3.1. Here, we denote $L_{k}^{(i)}=(\boldsymbol{v}^{\mathsf{T}}\boldsymbol{v})\boldsymbol{e}_{1}^{\mathsf{T}}(z_{i}\mathrm{I}-T_{k,k})^{-1}\boldsymbol{e}_{1}$ . Quantities denoted with the superscript $(i)$ correspond to the shift $z_{i}$ . The difference from Algorithm 2.1 is the addition of Lines 2 and 7. In particular, when $A$ is a real symmetric matrix and $\boldsymbol{v}$ is a real vector in Algorithm 3.1, only real arithmetic is needed to compute Lines 1, 4, and 6, whereas Lines 2 and 7 require complex arithmetic in general. Note that $\alpha_{k+1}$ is a real number in theory. However, when $A$ is not real but complex Hermitian, due to rounding error in finite precision arithmetic, the imaginary part of $\alpha_{k+1}$ may grow and affect the accuracy. Therefore, it is recommended to explicitly set the real part of the numerically computed $\alpha_{k+1}$ to its value.

Algorithm 3.1 Shifted Lanczos method for quadratic forms

A\in\mathbb{C}^{n\times n}

\boldsymbol{v}\in\mathbb{C}^{n}

z_{i}\in\mathbb{C}

i=1

2

\dots

m

L_{k}^{(i)}\in\mathbb{C}

i=1

2

\dots

m

\boldsymbol{v}_{1}=\boldsymbol{v}/\|\boldsymbol{v}\|

\boldsymbol{s}_{1}=A\boldsymbol{v}_{1}

\alpha_{1}=\boldsymbol{u}^{\mathsf{H}}\boldsymbol{v}_{1}

c_{1}^{(i)}=\boldsymbol{v}^{\top}\boldsymbol{v}

\delta_{1}^{(i)}=z_{i}-\alpha_{1}

\pi_{1}^{(i)}=1/\delta_{1}^{(i)}

L_{1}^{(i)}=c_{1}^{(i)}/(z_{i}-\alpha_{1})

i=1

2

\dots

m

3: for

k=1,2,\dots

, until convergence do

\boldsymbol{t}_{k}=\boldsymbol{s}_{k}-\alpha_{k}\boldsymbol{v}_{k}

\beta_{k}=\|\boldsymbol{t}_{k}\|

5: if

\beta_{k}=0

then break

\boldsymbol{v}_{k+1}=(\beta_{k})^{-1}\boldsymbol{t}_{k}

\boldsymbol{s}_{k+1}=A\boldsymbol{v}_{k+1}-\beta_{k}\boldsymbol{v}_{k}

\alpha_{k+1}=\boldsymbol{s}_{k+1}^{\mathsf{H}}\boldsymbol{v}_{k+1}

t_{k}^{(i)}=(\beta_{k})^{2}\pi_{k}^{(i)}

\delta_{k+1}^{(i)}=z_{i}-\alpha_{k+1}-t_{k}^{(i)}

\pi_{k+1}^{(i)}=1/\delta_{k+1}^{(i)}

c_{k+1}^{(i)}=c_{k}^{(i)}t_{k}^{(i)}\pi_{k}^{(i)}

L_{k+1}^{(i)}=L_{k}^{(i)}+c_{k+1}^{(i)}\pi_{k+1}^{(i)}

i=1

2

\dots

m

8: end for

We compare the shifted Lanczos method with related methods in terms of computational cost. Algorithms 3.2, 3.3, and 3.4 give simple modifications of the shifted COCG, COCR, and MINRES methods, respectively, for computing quadratic forms (1.1). Here, $z_{s}$ is the seed shift. The modifications are given by applying $\boldsymbol{v}^{\top}$ to the $k$ th iterate $\boldsymbol{x}_{k}$ and associated vectors. They produce approximations $G_{k}^{(i)}$ , $R_{k}^{(i)}$ , and $M_{k}^{(i)}$ , respectively, to the quadratic forms (1.1). In the shifted COCG and COCR methods, if a satisfactory seed iterate is obtained, one can use the seed switching technique [46] to choose a different shift as the seed shift $z_{s}$ . Similarly to the shifted MINRES method [36, Section 2.3], the shifted Lanczos method does not need such a seed switching technique.

Remark.

The modifications applied to the shifted COCG, COCR, and MINRES methods for computing quadratic forms to derive Algorithms 3.2, 3.3, and 3.4, respectively, can also be applied to the shifted CG method to compute quadratic forms. However, such a shifted modification of the CG method [19, Section 2.2] is mathematically equivalent to the shifted Lanczos method for Hermitian positive definite $z\mathrm{I}-A$ , whereas it is not mathematically equivalent to the shifted Lanczos method in general. Such a shifted CG method may not work on the case where $z$ is not real, because the shifted coefficient matrix $z\mathrm{I}-A$ is not Hermitian for $z\in\mathbb{C}\backslash\mathbb{R}$ . With the connection between the tridiagonal entries of the Jacobi matrix and the scalar coefficients of vectors in the CG method [33, Section 6.7], we may formulate the shifted CG method for quadratic forms, which is mathematically equivalent to the shifted Lanczos method for quadratic forms. See [27, Section 9.6] for different aspects of the shifted CG method.

For simplicity of comparison, we count basic scalar, vector, and matrix operations. The Lanczos method (Algorithm 2.1) needs one vector scale (scale), one dot product (dot), one vector norm (norm), two scalar–vector additions (axpy) and one matrix–vector product (matvec) per iteration. Table 3.2 gives the number of basic vector and matrix operations of the shifted methods. Table 3.2 gives the number of scalar operations for each shift $z_{i}$ per iteration. These tables show that in terms of the cost per iteration, the shifted Lanczos method is the cheapest of the methods compared.

Algorithm 3.2 Shifted COCG method for quadratic forms

A\in\mathbb{R}^{n\times n}

\boldsymbol{v}\in\mathbb{C}^{n}

z_{s}\in\mathbb{C}

z_{i}\in\mathbb{C}

i=1

2

\dots

m

G_{k}^{(i)}\in\mathbb{C}

i=1

2

\dots

m

\alpha_{-1}=1

\beta_{-1}=0

\boldsymbol{p}_{-1}=\boldsymbol{0}

\boldsymbol{r}_{0}=\boldsymbol{v}

\boldsymbol{p}_{0}=\boldsymbol{r}_{0}

\pi_{-1}^{(i)}=\pi_{0}^{(i)}=1

p_{0}^{(i)}=\boldsymbol{v}^{\mathsf{H}}\boldsymbol{r}_{0}

G_{0}^{(i)}=0

i=1

2

\dots

m

3: for

k=1,2,\dots

, until convergence do

4: if

{\boldsymbol{p}_{k-1}}^{\top}(z_{s}\mathrm{I}-A)\boldsymbol{p}_{k-1}=0

{\boldsymbol{r}_{k-1}}^{\top}\boldsymbol{r}_{k-1}=0

then switch the seed

\alpha_{k-1}=({\boldsymbol{r}_{k-1}}^{\top}\boldsymbol{r}_{k-1})/({\boldsymbol{p}_{k-1}}^{\top}(z_{s}\mathrm{I}-A)\boldsymbol{p}_{k-1})

\boldsymbol{r}_{k}=\boldsymbol{r}_{k-1}-\alpha_{k-1}(z_{s}\mathrm{I}-A)\boldsymbol{p}_{k-1}

\beta_{k-1}=({\boldsymbol{r}_{k}}^{\top}\boldsymbol{r}_{k})/({\boldsymbol{r}_{k-1}}^{\top}\boldsymbol{r}_{k-1})

r_{k}=\boldsymbol{v}^{\mathsf{H}}\boldsymbol{r}_{k}

\boldsymbol{p}_{k}=\boldsymbol{r}_{k}+\beta_{k-1}\boldsymbol{p}_{k-1}

6: for

i=1

2

\dots

m

\pi_{k}^{(i)}=[1+\alpha_{k-1}(z_{i}-z_{s})+(\beta_{k-2}/\alpha_{k-2})\alpha_{k-1}]\pi_{k-1}^{(i)}-(\beta_{k-2}/\alpha_{k-2})\alpha_{k-1}\pi_{k-2}^{(i)}

8: if

\pi_{k}^{(i)}=0

then output

G_{k-1}^{(i)}

\alpha_{k-1}^{(i)}=(\pi_{k-1}^{(i)}/\pi_{k}^{(i)})\alpha_{k-1}

G_{k}^{(i)}=G_{k-1}^{(i)}+\alpha_{k-1}^{(i)}p_{k-1}^{(i)}

\beta_{k-1}^{(i)}=(\pi_{k-1}^{(i)}/\pi_{k}^{(i)})^{2}\beta_{k-1}

p_{k}^{(i)}=r_{k}/\pi_{k}^{(i)}+\beta_{k-1}^{(i)}p_{k-1}^{(i)}

10: end for

11: end for

Algorithm 3.3 Shifted COCR method for quadratic forms

A\in\mathbb{R}^{n\times n}

\boldsymbol{v}\in\mathbb{C}^{n}

z_{s}\in\mathbb{C}

z_{i}\in\mathbb{C}

i=1

2

\dots

m

R_{k}^{(i)}\in\mathbb{C}

i=1

2

\dots

m

\alpha_{-1}=1

\beta_{-1}=0

\boldsymbol{q}_{-1}=\boldsymbol{0}

\boldsymbol{r}_{0}=\boldsymbol{v}

r_{0}=\boldsymbol{v}^{\mathsf{H}}\boldsymbol{r}_{0}

p_{-1}^{(i)}=0

\pi_{-1}^{(i)}=\pi_{0}^{(i)}=1

p_{0}^{(i)}=\boldsymbol{v}^{\mathsf{H}}\boldsymbol{r}_{0}

R_{0}^{(i)}=0

i=1

2

\dots

m

3: for

k=1,2,\dots

, until convergence do

\boldsymbol{q}_{k-1}=(z_{s}\mathrm{I}-A)\boldsymbol{r}_{k-1}+\beta_{k-2}\boldsymbol{q}_{k-2}

5: if

{\boldsymbol{q}_{k-1}}^{\top}\boldsymbol{q}_{k-1}=0

then switch the seed

\alpha_{k-1}=[{\boldsymbol{r}_{k-1}}^{\top}(z_{s}\mathrm{I}-A)\boldsymbol{r}_{k-1}]/({\boldsymbol{q}_{k-1}}^{\top}\boldsymbol{q}_{k-1})

7: for

i=1

2

\dots

m

\pi_{k}^{(i)}=(1+(\beta_{k-2}/\alpha_{k-2})\alpha_{k-1}+\alpha_{k-1}(z_{i}-z_{s}))\pi_{k-1}^{(i)}-(\beta_{k-2}/\alpha_{k-2})\alpha_{k-1}\pi_{k-2}^{(i)}

9: if

\pi_{k}^{(i)}=0

then output

R_{k-1}^{(i)}

10:

\beta_{k-2}^{(i)}=(\pi_{k-2}^{(i)}/\pi_{k-1}^{(i)})^{2}\beta_{k-2}

\alpha_{k-1}^{(i)}=(\pi_{k-1}^{(i)}/\pi_{k}^{(i)})\alpha_{k-1}

p_{k-1}^{(i)}=r_{k-1}/\pi_{k-1}^{(i)}+\beta_{k-2}^{(i)}p_{k-2}^{(i)}

R_{k}^{(i)}=R_{k-1}^{(i)}+\alpha_{k-1}^{(i)}p_{k-1}^{(i)}

11: end for

12:

\boldsymbol{r}_{k}=\boldsymbol{r}_{k-1}-\alpha_{k-1}\boldsymbol{q}_{k-1}

r_{k}=\boldsymbol{v}^{\mathsf{H}}\boldsymbol{r}_{k}

13: if

{\boldsymbol{r}_{k-1}}^{\top}(z_{s}\mathrm{I}-A)\boldsymbol{r}_{k-1}=0

then switch the seed

14:

\beta_{k-1}=[{\boldsymbol{r}_{k}}^{\top}(z_{s}\mathrm{I}-A)\boldsymbol{r}_{k}]/[{\boldsymbol{r}_{k-1}}^{\top}(z_{s}\mathrm{I}-A)\boldsymbol{r}_{k-1}]

15: end for

Algorithm 3.4 Shifted MINRES method for quadratic forms

A\in\mathbb{C}^{n\times n}

\boldsymbol{v}\in\mathbb{C}^{n}

\boldsymbol{x}_{0}\in\mathbb{C}^{n}

z_{s}\in\mathbb{C}

z_{i}\in\mathbb{C}

i=1

2

\dots

m

M_{k}^{(i)}\in\mathbb{C}

i=1

2

\dots

m

\beta_{0}=0

\boldsymbol{q}_{0}=\boldsymbol{0}

\boldsymbol{r}_{0}=\boldsymbol{v}-(z_{s}\mathrm{I}-A)\boldsymbol{x}_{0}

\boldsymbol{q}_{1}={\|\boldsymbol{r}_{0}\|}^{-1}\boldsymbol{r}_{0}

q_{1}=\boldsymbol{v}^{\mathsf{H}}\boldsymbol{q}_{1}

f_{1}^{(i)}=1

p_{-1}^{(i)}=p_{0}^{(i)}=0

M_{0}^{(i)}=0

i=1

2

\dots

m

3: for

k=1,2,\dots

, until convergence do

\boldsymbol{s}_{k}=A\boldsymbol{q}_{k}-\beta_{k-1}\boldsymbol{q}_{k-1}

\alpha_{k}={\boldsymbol{s}_{k}}^{\mathsf{H}}\boldsymbol{q}_{k}

\boldsymbol{t}_{k}=\boldsymbol{s}_{k}-\alpha_{k}\boldsymbol{q}_{k}

\beta_{k}=\|\boldsymbol{t}_{k}\|

q_{k}=\boldsymbol{v}^{\mathsf{H}}\boldsymbol{q}_{k}

5: for

i=1

2

\dots

m

r_{k-2,k}^{(i)}=0

r_{k-1,k}^{(i)}=\beta_{k-1}

r_{k,k}^{(i)}=z_{i}-\alpha_{k}

7: if

k\geq 3

then update

[r_{k-2,k}^{(i)},r_{k-1,k}^{(i)}]^{\top}=G_{k-2}^{(i)}[r_{k-2,k}^{(i)},r_{k-1,k}^{(i)}]^{\top}

8: if

k\geq 2

then update

[r_{k-1,k}^{(i)},r_{k,k}^{(i)}]^{\top}=G_{k-1}^{(i)}[r_{k-1,k}^{(i)},r_{k,k}^{(i)}]^{\top}

9: Compute

G_{k}^{(i)}=\left[\begin{smallmatrix}c_{k}^{(i)}&\bar{s}_{k}^{(i)}\\ -s_{k}^{(i)}&\bar{c}_{k}^{(i)}\end{smallmatrix}\right]

and update

r_{k,k}^{(i)}

such that

[r_{k,k}^{(i)},0]^{\top}=G_{k}^{(i)}[r_{k,k}^{(i)},\beta_{k}]^{\top}

|c_{k}^{(i)}|^{2}+|s_{k}^{(i)}|^{2}=1

c_{k}^{(i)}

s_{k}^{(i)}\in\mathbb{C}

10:

p_{k}^{(i)}=(r_{k,k}^{(i)})^{-1}(q_{k}-r_{k-2,k}^{(i)}p_{k-2}^{(i)}-r_{k-1,k}^{(i)}p_{k-1}^{(i)})

M_{k}^{(i)}=M_{k-1}^{(i)}+\|\boldsymbol{r}_{0}\|c_{k}^{(i)}f_{k}^{(i)}p_{k}^{(i)}

f_{k+1}^{(i)}=-\bar{s}_{k}^{(i)}f_{k}^{(i)}

11: end for

12:

\boldsymbol{q}_{k+1}=(\beta_{k})^{-1}\boldsymbol{t}_{k}

13: end for

Table 3.1: Basic vector and matrix operations

Method	scale	dot	norm	axpy	matvec
Shifted Lanczos	1	1	1	2	1
Shifted COCG	1	1	1	2	1
Shifted COCR	0	3	0	2	1
Shifted MINRES	1	2	1	2	1

Method: name of the method, scale: vector scale, dot: dot product, norm: vector norm, axpy: scalar–vector addition, matvec: matrix–vector product.

Table 3.2: Scalar operations for each shift

z_{i}

per iteration

Method	$+$	$\times$	$/$	$\sqrt{}$	Total
Shifted Lanczos	3	4	1	0	8
Shifted COCG	6	9	3	0	18
Shifted COCR	6	8	3	0	17
Shifted MINRES	11	20	3	1	35

Method: name of the method, $+$ : addition, $\times$ : multiplication, $/$ : division, $\sqrt{}$ : square root, Total: total number of operations.

3.2 Breakdown-free condition

A breakdown resulting from division by zero for $\delta_{k+1}^{(i)}=0$ may occur in Line 7 of Algorithm 3.1 before a solution is obtained. Therefore, we give a sufficient condition such that the shifted Lanczos method does not break down. Let $|T_{0}^{<}|=1$ for convenience, where $|\cdot|$ denotes the determinant of a matrix. For convenience, we omit superscript $(i)$ from the quantities given in Algorithm 3.1 and prepare lemmas on shifted versions of well-known properties of the Jacobi matrix.

Lemma 3.1.

Let $T_{k}^{<}$ and $\delta_{k}$ be defined as above. Assume that $|T_{i}^{<}|\neq 0$ holds for $i=1$ , $2$ , $\dots$ , $k$ , $k\in\mathbb{Z}_{>0}$ . If $\delta_{k+1}=0$ for $k\in\mathbb{Z}_{>0}$ , then we have $|T_{k+1}^{<}|=0$ .

Proof.

The condition $\delta_{k+1}=0$ gives $t_{k}=z-\alpha_{k+1}$ . It follows from [10, Lemma 3.2]

\displaystyle|T_{k+1}^{<}|=(z-\alpha_{k+1})|T_{k}^{<}|-(\beta_{k})^{2}|T_{k-1}^{<}|,\quad k>0

(3.29)

and [10, Section 3.3]

\delta_{k+1}=\frac{|T_{k+1}^{<}|}{|T_{k}^{<}|},\quad k>0

(3.30)

that

$\displaystyle\|T_{k+1}^{<}\|$	$\displaystyle=\|T_{k}^{<}\|\left(z-\alpha_{k+1}-(\beta_{k})^{2}\frac{\|T_{k-1}^{<}\|}{\|T_{k}^{<}\|}\right)$	(3.31)
	$\displaystyle=\|T_{k}^{<}\|\left(z-\alpha_{k+1}-t_{k}\right)$	(3.32)
	$\displaystyle=0.$	(3.33)

∎∎

Then, the breakdown-free condition is described as follows.

Theorem 3.2.

Let $\delta_{k}$ be defined as above. Let $\lambda_{1}$ and $\lambda_{n}$ be the smallest and largest eigenvalues of $A$ , respectively. If $z\in\mathbb{C}$ satisfies $z\not\in[\lambda_{1},\lambda_{n}]$ , then $\delta_{k}\neq 0$ holds for $k\in\mathbb{Z}_{>0}$ .

Proof.

From the interlacing property of eigenvalues [17, Theorem 4.3.17], $T_{k,k}$ does not have an eigenvalue equal to $z$ . Hence, $|T_{k}^{<}|\neq 0$ holds. From Lemma 3.1, the assertion holds. ∎∎

Theorem 3.2 shows that the shifted Lanczos method does not break down whenever each $z_{i}\in\mathbb{C}$ satisfies $z_{i}\not\in[\lambda_{1},\lambda_{n}]$ . The condition $z\not\in[\lambda_{1},\lambda_{n}]$ in Theorem 3.2 is not necessarily equivalent to assuming that $z\mathrm{I}-A$ is positive or negative definite because $z\in\mathbb{C}$ . The shifted MINRES methods [11, 36] do not break down for nonsingular $z\mathrm{I}-A$ , whereas the shifted COCG and COCR methods may break down.

Projection methods for symmetric or Hermitian eigenproblems [34, 31] typically take shift points for quadrature points in a circle and can circumvent taking quadrature points on the real line. Therefore, the shifted Lanczos method can avoid breakdown when applied to quadratic forms in these methods with particular choices of quadrature points.

3.3 Convergence bound

The approximation of the quadratic form $\boldsymbol{v}^{\mathsf{T}}(z\mathrm{I}-A)^{-1}\boldsymbol{v}$ by using the shifted Lanczos method can be viewed as a way of solving the shifted linear system $(z\mathrm{I}-A)\boldsymbol{x}=\boldsymbol{v}$ by using the same method and computing $\boldsymbol{v}^{\mathsf{T}}\boldsymbol{x}_{*}$ , where $\boldsymbol{x}_{*}=S^{-1}\boldsymbol{v}$ . Concretely, to determine the $k$ th iterate $\boldsymbol{x}_{k}$ of the Lanczos method for the linear system $S\boldsymbol{x}=\boldsymbol{v}$ , the method imposes the Galerkin condition

\displaystyle\boldsymbol{v}-S\boldsymbol{x}_{k}\perp\mathcal{K}_{k}(A,\boldsymbol{v}).

(3.34)

The iterate is the same as the CG iterate if $S$ is Hermitian positive definite; otherwise, they may be different. The absolute error of the shifted Lanczos method for quadratic forms for the $k$ th iteration is

\displaystyle\varepsilon_{k}=|L_{k}-\boldsymbol{v}^{\mathsf{H}}S^{-1}\boldsymbol{v}|

(3.35)

and its upper bound is given by using the Cauchy-Schwarz inequality

	$\displaystyle\varepsilon_{k}$	$\displaystyle=\|\boldsymbol{v}^{\mathsf{H}}(\boldsymbol{x}_{k}-S^{-1}\boldsymbol{v})\|$		(3.36)
		$\displaystyle\leq\\|\boldsymbol{v}\\|\\|\boldsymbol{x}_{k}-\boldsymbol{x}_{*}\\|.$		(3.37)

For the Hermitian positive definite linear systems $A\boldsymbol{x}=\boldsymbol{b}$ , $\boldsymbol{b}\in\mathbb{C}^{n}$ , a well-known upper bound [21] of the Lanczos method is related to the $A$ -norm and a recent upper bound depends on the distribution of the eigenvalues [29, Theorem B.1].

The following assertions give an error bound of the shifted Lanczos method for quadratic forms.

Theorem 3.3.

Let $A\in\mathbb{C}^{n\times n}$ be a Hermitian matrix, $z\in\mathbb{C}$ , $S=z\mathrm{I}-A$ , $\boldsymbol{v}\in\mathbb{C}^{n}$ , and $\mathbb{P}_{k}$ be the set of all polynomials with degree less than $k$ . Then, the absolute error of the shifted Lanczos method for the quadratic form $\boldsymbol{v}^{\mathsf{H}}S^{-1}\boldsymbol{v}$ for the $k$ th iteration satisfies

\displaystyle\varepsilon_{k}\leq 2\tau_{k}\|\boldsymbol{v}\|^{2}.

(3.38)

where

\displaystyle\tau_{k}=\min_{p\in\mathbb{P}_{k}}\max_{t\in[\lambda_{1},\lambda_{n}]}|p(z-t)-(z-t)^{-1}|.

(3.39)

Proof.

If $\boldsymbol{x}_{*}$ and $\boldsymbol{x}_{k}$ are defined as above, then the error norm has the expression

	$\displaystyle\\|\boldsymbol{x}_{k}-\boldsymbol{x}_{*}\\|$	$\displaystyle=\\|V_{k}(T_{k,k}^{<})^{-1}(\\|\boldsymbol{v}\\|\boldsymbol{e}_{1})-S^{-1}\boldsymbol{v}\\|$		(3.40)
		$\displaystyle=\\|V_{k}(T_{k,k}^{<})^{-1}\boldsymbol{e}_{1}-S^{-1}\boldsymbol{v}_{1}\\|\\|\boldsymbol{v}\\|.$		(3.41)

For any polynomial $p\in\mathbb{P}_{k}$ , the first factor of the last quantity is bounded as

	$\displaystyle\\|V_{k}(T_{k,k}^{<})^{-1}\boldsymbol{e}_{1}-S^{-1}\boldsymbol{v}_{1}\\|$	(3.42)
$\displaystyle\leq$	$\displaystyle\\|V_{k}[p(T_{k,k}^{<})-(T_{k,k}^{<})^{-1}]\boldsymbol{e}_{1}-[p(S)-S^{-1}]\boldsymbol{v}_{1}\\|+\\|V_{k}p(T_{k,k}^{<})\boldsymbol{e}_{1}-p(S)\boldsymbol{v}_{1}\\|$	(3.43)
$\displaystyle\leq$	$\displaystyle\\|p(T_{k,k}^{<})-(T_{k,k}^{<})^{-1}\\|+\\|p(S)-S^{-1}\\|$	(3.44)

with $V_{k}^{\mathsf{H}}V_{k}=\mathrm{I}$ and $\|\boldsymbol{e}_{1}\|=\|\boldsymbol{v}_{1}\|=1$ . Here, we used the shifted Lanczos decomposition (3.13) and the identity

\displaystyle p(S)\boldsymbol{v}_{1}=V_{k}p(T_{k,k}^{<})\boldsymbol{e}_{1},

(3.45)

cf. [29, Lemma 4.1]. Because of the interlacing property of eigenvalues [17, Theorem 4.3.17]

\displaystyle\|p(S)-(S)^{-1}\|\leq\max_{t\in[\lambda_{1},\lambda_{n}]}|p(z-t)-(z-t)^{-1}|

(3.46)

and

\displaystyle\|p(T_{k,k}^{<})-(T_{k,k}^{<})^{-1}\|\leq\max_{t\in[\lambda_{1},\lambda_{n}]}|p(z-t)-(z-t)^{-1}|.

(3.47)

Therefore, the assertion (3.38) holds. ∎∎

The shifted Lanczos method uses a short-term recurrence and may suffer from rounding errors. Rounding errors in the Lanczos method may quickly lead to loss of orthogonality of the computed basis vectors $\boldsymbol{v}_{i}$ of the Krylov subspace and thus could cause a delay in convergence. We leave an open problem of giving the bound in finite precision, cf. [29, Theorem 6.2].

3.4 Estimation of error

The bound presented in the previous section gives insights on the convergence of the shifted Lanczos method for quadratic forms but is not practically used as a stopping criterion in the iteration of Algorithm 3.1 because it needs to know the condition number of $z\mathrm{I}-A$ . To check if a satisfactory solution is obtained for each iteration in Algorithm 3.1 in practice, we formulate two estimates of the error $\varepsilon_{k}$ .

First, the shifted Lanczos decomposition (3.13) gives

\displaystyle SV_{k}=V_{k}T_{k}^{<}-\beta_{k}\boldsymbol{v}_{k+1}\boldsymbol{e}_{k}.

(3.48)

Multiplying this by $\boldsymbol{v}_{1}^{\mathsf{T}}S^{-1}$ from the left and $(T_{k}^{<})^{-1}\boldsymbol{e}_{1}$ from the right gives

\displaystyle\boldsymbol{e}_{1}(T_{k}^{<})^{-1}\boldsymbol{e}_{1}=\boldsymbol{v}_{1}S^{-1}\boldsymbol{v}_{1}-\beta_{k}(\boldsymbol{v}_{1}S^{-1}\boldsymbol{v}_{k+1})[\boldsymbol{e}_{k}(T_{k}^{<})^{-1}\boldsymbol{e}_{1}].

(3.49)

Together with $\|\boldsymbol{v}\|^{2}$ , we have

\displaystyle\varepsilon_{k}=\beta_{k}\|\boldsymbol{v}\|^{2}|(\boldsymbol{v}_{1}^{\mathsf{H}}S^{-1}\boldsymbol{v}_{k+1})[\boldsymbol{e}_{k}^{\mathsf{T}}(T_{k}^{<})^{-1}\boldsymbol{e}_{1}]|.

(3.50)

With the representation (3.22) and a positive integer $d$ , we approximate $\boldsymbol{v}_{1}^{\mathsf{H}}S^{-1}\boldsymbol{v}_{k+1}$ by

\displaystyle\boldsymbol{v}_{1}^{\mathsf{H}}(S_{k+d})^{-1}\boldsymbol{v}_{k+1}=\boldsymbol{e}_{1}(T_{k+d}^{<})^{-1}\boldsymbol{e}_{k+1}.

(3.51)

Therefore, we obtain the following estimate:

\displaystyle\varepsilon_{k}\simeq\beta_{k}\|\boldsymbol{v}\|^{2}|[\boldsymbol{e}_{1}^{\mathsf{T}}(T_{k}^{<})^{-1}\boldsymbol{e}_{k}][\boldsymbol{e}_{1}(T_{k+d}^{<})^{-1}\boldsymbol{e}_{k+1}]|\equiv\mu_{k,d}.

(3.52)

This estimate $\mu_{k,d}$ requires additional $d$ iterations and needs $\boldsymbol{e}_{1}^{\mathsf{T}}(T_{k}^{<})^{-1}\boldsymbol{e}_{k}$ and $\boldsymbol{e}_{1}(T_{k+d}^{<})^{-1}\boldsymbol{e}_{k+1}$ to compute. The former is the $(1,k)$ entry of $(T_{k}^{<})^{-1}$ and the latter is the $(1,k+1)$ entry of $(T_{k+d}^{<})^{-1}$ . These quantities can be computed recursively [10, Sections 3.2, 3.4] as follows:

	$\displaystyle\boldsymbol{e}_{1}^{\mathsf{T}}(T_{k}^{<})^{-1}\boldsymbol{e}_{k}$	$\displaystyle=(-1)^{k-1}\frac{\beta_{1}\beta_{2}\cdots\beta_{k-1}}{\delta_{1}\delta_{2}\cdots\delta_{k}},$		(3.53)
	$\displaystyle\boldsymbol{e}_{1}^{\mathsf{T}}(T_{k+d})^{-1}\boldsymbol{e}_{k+1}$	$\displaystyle=(-1)^{k}\frac{\beta_{1}\beta_{2}\cdots\beta_{k}}{\delta_{1}\delta_{2}\cdots\delta_{k+d}}\varphi_{2}^{(k+d)}\varphi_{3}^{(k+d)}\cdots\varphi_{d}^{(k+d)},$		(3.54)

where $\delta_{k}$ is defined in (3.26) and

\displaystyle\varphi_{d}^{(k+d)}=z-\alpha_{k+d},\quad\varphi_{j}^{(k+d)}=z-\alpha_{k+j}-\frac{\beta_{k+j}^{2}}{\varphi_{j+1}^{(k+d)}},\quad j=d-1,d-2,\dots,2.

(3.55)

Therefore, we may update $\boldsymbol{e}_{1}^{\mathsf{T}}(T_{k}^{<})^{-1}\boldsymbol{e}_{k}$ by

\displaystyle\boldsymbol{e}_{1}^{\mathsf{T}}(T_{k+1}^{<})^{-1}\boldsymbol{e}_{k+1}=-\frac{\beta_{k}}{\delta_{k+1}}\boldsymbol{e}_{1}^{\mathsf{T}}(T_{k}^{<})^{-1}\boldsymbol{e}_{k},

(3.56)

To update $\boldsymbol{e}_{1}^{\mathsf{T}}(T_{k+1}^{<})^{-1}\boldsymbol{e}_{k+1}$ for each shift per iteration, one multiplication by $\beta_{k-1}$ and one division by $\delta_{k}$ are required. To compute $\varphi_{j}^{(k+d)}$ for $j=2$ , $3$ , $\dots$ , $d$ for each shift per iteration, we require one minus in $z-\alpha_{k+d}$ , $d-2$ minuses and $d-2$ divisions in $(z-\alpha_{k+j})-\beta_{k+j}^{2}/\varphi_{j+1}^{(k+d)}$ if one stores $z-\alpha_{k+j}$ and $\beta_{k+j}^{2}$ for $j=2$ , $3$ , $\dots$ , $d-1$ computed in Line 7 of Algorithm 3.1. To update $\boldsymbol{e}_{1}^{\mathsf{T}}(T_{k+d})^{-1}\boldsymbol{e}_{k+1}$ for each shift per iteration, $d$ multiplications by $\varphi_{j}^{(k+d)}$ , one multiplication by $\beta_{k}$ , and one division by $\delta_{k+d}$ are required.

Second, the triangular equation gives

	$\displaystyle\varepsilon_{k}$	$\displaystyle=\|L_{k}-L_{k+d}+L_{k+d}-\boldsymbol{v}^{\mathsf{H}}S^{-1}\boldsymbol{v}\|$		(3.57)
		$\displaystyle\leq\|L_{k}-L_{k+d}\|+\|L_{k+d}-\boldsymbol{v}^{\mathsf{H}}S^{-1}\boldsymbol{v}\|.$		(3.58)

Under the assumption that the approximation error for the ( $k+d$ )th iteration is significantly smaller than that for the $k$ th iteration, we obtain an estimate for the $k$ th iteration

\displaystyle\varepsilon_{k}\simeq|L_{k}-L_{k+d}|\equiv\nu_{k,d}.

(3.59)

This estimate $\nu_{k,d}$ requires additional $d$ iterations and one minus in $L_{k}-L_{k+d}$ . An analogous estimate can be found in [41, Section 4] for the $A$ -norm error in the CG method under a similar assumption. The shifted coefficient matrix $z\mathrm{I}-A$ may not be Hermitian and may not form the $(z\mathrm{I}-A)$ -norm in the usual sense.

4 Numerical experiments

Numerical experiments were performed to compare the proposed method with the previous methods—shifted COCG method (Algorithm 3.2), shifted COCR method (Algorithm 3.3), shifted MINRES method (Algorithm 3.4), shifted Lanczos method (Algorithm 3.1, and direct solver using the MATLAB function mldivide for solving (1.2)—in terms of the number of iterations and CPU time.

All computations were performed on a computer with an Intel Xeon E5-2670 v2 2.50 GHz CPU, 256 GB of random-access memory (RAM), and CentOS 6.10. All programs were coded and run in MATLAB R2019a in double-precision floating-point arithmetic with unit round-off at $2^{-52}\simeq 2.2\cdot 10^{-16}$ .

Table 4.1 gives information about the test matrices, including the size of each matrix, density of nonzero entries [%], (estimated) condition number, and application from which the matrix arose. The condition number was estimated using the MATLAB function condest. Matrices mhd1280b, apache2, CurlCurl_3 and thermal2 are from [3]; the other matrix, VCNT1000000std, is from [18]. Matrix mhd1280b is complex symmetric and matrix conf5.4-00l8x8-2000 is complex Hermitian and indefinite, whereas the others are real symmetric. The Cholesky factorization [32] proved that matrices mhd1280b and apache2 are positive definite. According to the information in [3], matrix CurlCurl_3 is not positive definite and matrix thermal2 is positive definite.

The condition number for CurlCurl_3 could not be computed because of insufficient computer memory. Vector $\boldsymbol{v}$ was simply set to $\boldsymbol{v}=n^{-1/2}\boldsymbol{e}$ for reproducibility, where $\boldsymbol{e}$ is the all-ones vector. The shifts were set to

z_{i}=\exp\left(-\frac{2i+1}{2m}\pi\mathrm{i}\right),\quad i=1,2,\dots,m,\quad m=16,

(4.1)

where $\mathrm{i}$ is the imaginary unit. Because $\mathrm{Im}(z_{i})\neq 0$ , the conditions in Theorem 3.2 are satisfied and the shifted Lanczos method is breakdown-free. These shifts demonstrate typical choices of quadrature points in the projection method for eigenproblems [35].

Table 4.1: Information on test matrices

Matrix	$n$	Density [%]	condest	Application
mhd1280b	1,280	$1.4\cdot 10^{0\phantom{-}}$	$6.0\cdot 10^{12}$	Magnetohydrodynamics
conf5.4-00l8x8-2000	49,152	$8.2\cdot 10^{-2}$	$3.6\cdot 10^{4\phantom{0}}$	Quantum chromodynamics
apache2	715,176	$9.4\cdot 10^{-4}$	$5.3\cdot 10^{6\phantom{0}}$	Structural analysis
VCNT1000000std	1,000,000	$4.0\cdot 10^{-3}$	$1.2\cdot 10^{7\phantom{0}}$	Quantum mechanics
CurlCurl_3	1,219,574	$9.1\cdot 10^{-4}$	—	Electromagnetic analysis
thermal2	1,228,045	$5.7\cdot 10^{-4}$	$7.5\cdot 10^{6\phantom{0}}$	Steady-state thermal analysis

Matrix: name of the matrix, $n$ : size of the matrix, density: density of nonzero entries, condest: condition number estimated by using the MATLAB function condest, application: application from which the matrix arises.

4.1 Comparisons in terms of CPU time

We compare the methods in terms of CPU time. Figure 4.1 shows the relative error versus the number of iterations for the compared methods on test matrices. The shifted methods terminated when the largest relative error among $i=1$ , $2$ , $\dots$ , $16$ became less than or equal to $10^{-10}$ . Here, we adopted the numerical solution computed by MATLAB’s sparse direct solver as the exact solution. Table 4.2 gives the largest relative residual norm $\max_{i}\|\boldsymbol{v}-(z_{i}\mathrm{I}-A)\tilde{\boldsymbol{x}}^{(i)}\|/\|\boldsymbol{b}\|$ for MATLAB’s sparse direct solver for the test matrices, where $\tilde{\boldsymbol{x}}^{(i)}$ is the numerical solution of the linear system $(z_{i}\mathrm{I}-A)\boldsymbol{x}^{(i)}=\boldsymbol{b}$ computed by the sparse direct solver. The sparse direct solver was accurate on mhd1280b, VCNT1000000std, conf5.4-00l8x8-2000, and thermal2 and not accurate on apache2 and CurlCurl_3. The convergence curve plotted corresponds to the case where the required number of iterations is the largest among shifts $z_{i}$ , $i=1$ , $2$ , $\dots$ , $16$ , because this case determines the total CPU time. Such cases are for $z_{1}$ on mhd1280b, conf5.4-00l8x8-2000, apache2, and thermal2; for $z_{6}$ on CurlCurl_3; and for $z_{1}6$ on VCNT1000000std. The shifted Lanczos method converged faster than the other methods on apache2, CurlCurl_3, and thermal2 and was competitive with the shifted MINRES method on VCNT1000000std. The shifted COCR and MINRES methods were almost identical on CurlCurl_3. The convergence curves of the four methods on VCNT1000000std and CurlCurl_3 decreased monotonically. Although the shifted MINRES method is monotonically non-increasing in terms of the residual norm, Figures 4.1 4.1 and 4.1 4.1 show that the shifted MINRES method does not necessarily converge monotonically in terms of the error. The convergence curves of the four methods are oscillatory on apache2 and thermal2. The shifted Lanczos method did not reach the stopping criterion on CurlCurl_3, whereas the shifted COCR and MINRES methods did. The curves of the shifted MINRES and Lanczos methods for conf5.4-00l8x8-2000, those of the shifted COCG and Lanczos methods for VCNT1000000std, those of the shifted COCG and Lanczos methods for CurlCurl_3, and those of the shifted COCR and MINRES methods for CurlCurl_3 seem to overlap.

Table 4.3 gives the CPU time in seconds taken by the four methods. The symbol * stands for the least CPU time for each matrix. The shifted Lanczos method was competitive with or faster than the MINRES method except on CurlCurl_3, whereas the shifted COCR method took more CPU time. CCOG required fewer iterations than the shifted Lanczos method but took more CPU time on CurlCurl_3 and thermal2. This is because the shifted COCG (and COCR) method used $z_{1}\neq\mathbb{R}$ as the seed, and its complex symmetric Lanczos iterations in Line 5 of Algorithm 3.2 were performed in complex arithmetic. When setting the seed to zero, its complex symmetric Lanczos iterations were performed with real operations; however, the method required more CPU time and iterations than the case with shift $z_{1}$ .

Table 4.2: Maximum relative residual norm for the direct method

Matrix	$\max_{i}\\|\boldsymbol{v}-(z_{i}\mathrm{I}-A)\tilde{\boldsymbol{x}}^{(i)}\\|/\\|\boldsymbol{v}\\|$
mhd1280b	$2.0\cdot 10^{-16}$
conf5.4-00l8x8-2000	$1.0\cdot 10^{-15}$
apache2	$4.6\cdot 10^{-12}$
VCNT1000000std	$9.4\cdot 10^{-16}$
CurlCurl_3	$5.0\cdot 10^{-12}$
thermal2	$6.4\cdot 10^{-16}$

Matrix: name of the matrix, $\max_{i}\|\boldsymbol{v}-(z_{i}\mathrm{I}-A)\tilde{\boldsymbol{x}}^{(i)}\|/\|\boldsymbol{v}\|$ : maximum relative error norm.

Table 4.3: Number of iterations and CPU times [s] on test matrices for the methods compared

Method	mhd1280b		conf5.4-00l8x8-2000		apache2
	iter	time	iter	time	iter	time
mldivide		0.08		6,505		1,978
Shifted COCG	446	0.05	—		6,613	175
Shifted COCR	437	0.04	—		6,461	163
Shifted MINRES	281	0.03	266	1.67	7,283	95.4
Shifted Lanczos	219	*0.02	266	*1.52	5,662	*68.4

Method	VCNT1000000std		CurlCurl_3		thermal2
	iter	time	iter	time	iter	time
mldivide		555		749,696		563
Shifted COCG	16	1.18	4,974	246	214	11.40
Shifted COCR	16	1.16	5,698	382	279	14.57
Shifted MINRES	14	*0.73	5,779	156	216	6.53
Shifted Lanczos	16	0.81	5,139	*135	214	*6.17

4.2 Spectrum and shift

We illustrate the effect of the spectrum of $A$ and shift $z$ on the convergence. The condition number of shifted matrix $z\mathrm{I}-A$ is given by

\displaystyle\kappa=\frac{\max_{i}|z-\lambda_{i}|}{\min_{i}|z-\lambda_{i}|}.

(4.2)

This means that the condition number has a connection between the spectrum and shift. For convenience of the computation of a spectrum, we take a small matrix mhd1280b as an example. The largest and smallest eigenvalues of the matrix are approximately $70.32$ and $1.48\cdot 10^{-11}$ , respectively. Substitute the shift $z=\lambda_{1}+\zeta\mathrm{i}$ for $\zeta=10^{j}$ , $j=-1$ , $-2$ , $-3$ , and $-4$ for controlling the condition number. Then, the corresponding condition numbers of $z\mathrm{I}-A$ are as given in Table 4.4. Table 4.4 gives the required number of iterations and CPU time for the compared methods for each value of $\zeta$ value. The stopping criterion was that the relative error became less than or equal to $10^{-10}$ . The table also shows that the shifted Lanczos method was competitive with or outperformed other methods in terms of the CPU time. It shows that as the condition number increases, the required number of iterations tends to increase.

Table 4.4: Number of iterations and CPU times [s] on test matrices for the methods compared for different values of shift

z=\lambda_{1}+\zeta\mathrm{i}

$\zeta$	$10^{-1}$		$10^{-2}$		$10^{-3}$		$10^{-4}$
condition number	$1.1\cdot 10^{3}$		$1.2\cdot 10^{4}$		$1.0\cdot 10^{5}$		$8.9\cdot 10^{5}$
	iter	time	iter	time	iter	time	iter	time
Shifted COCG	118	*0.01	357	0.03	1,087	0.10	3,013	0.27
Shifted COCR	117	*0.01	331	0.03	1,018	0.10	3,129	0.29
Shifted MINRES	84	*0.01	231	0.03	690	0.07	1,988	0.19
Shifted Lanczos	76	*0.01	226	*0.02	680	*0.06	1,894	*0.17

4.3 Estimation of error

We tested the estimates $\mu_{k,d}$ and $\nu_{k,d}$ developed in Section 3.4 on the matrices used in Section 4.1. Figure 4.3 shows the relative error and its estimates $\mu_{k,d}$ and $\nu_{k,d}$ with $d=5$ versus the number of iterations. Here, the plotted curve corresponds to the case where the required number of iterations was the largest among shifts $z_{i}$ , $i=1$ , $2$ , $\dots$ , $16$ . For presentation, these estimates are normalized by $\boldsymbol{v}^{\mathsf{H}}\tilde{\boldsymbol{x}}^{(i)}$ . Both estimates were accurate on conf5.4-00l8x8-2000, VCNT1000000std, and thermal2. The estimate $\zeta_{k,5}$ was accurate on mdh1280b. Both estimates underestimated the error on apache2 and CurlCurl_3 because the assumption made in Section 3.4 might not hold and the shifted Lanczos method lacks the monotonicity of an error norm such as the one in the CG method (cf. [41]). The oscillations of the estimates are similar to those observed for the errors on apache2 and thermal2. However, we have no theoretical justification of these estimates at this time and when they capture the convergence well, in particular for finite-precision arithmetic, is open.

5 Conclusions

We explored the computation of quadratic forms of Hermitian matrix resolvents. In contrast to previous shifted Krylov subspace methods, our method approximates the matrix resolvent directly. The underlying concept used in approximating the resolvent is to exploit the moment-matching property of a shifted Lanczos method by solving the Vorobyev moment problem. We showed that the shifted Lanczos method matches the first $k$ moments of the original model and those of the reduced model and extended the scope of the problems that the standard Lanczos method can solve. We derived the inverse of a linear operator representing the reduced model and related it to an entry of a Jacobi matrix resolvent. The entry can be efficiently computed by using a recursive formula. Previous shifted Krylov subspace methods work on a real symmetric matrix with a complex shift, whereas the proposed method works on a Hermitian matrix with a complex shift and does not break down, provided that the shift is not in the interior of the extremal eigenvalues of the Hermitian matrix. We gave an error bound and estimates for the shifted Lanczos method for quadratic forms. Numerical experiments on matrices drawn from real-world applications showed that the shifted Lanczos method is competitive with the shifted MINRES method and outperforms it when solving some problems. We illustrated the effect of the spectrum and shift on the convergence and showed that the error estimate is reasonable.

We intend to perform preconditioning for the shifted Lanczos method for quadratic forms in future work. The shifted COCG and COCR methods can use preconditioning if the preconditioned matrix is complex symmetric; however, for the shifted Lanczos method, it is not trivial to incorporate preconditioning. How to monitor the convergence of the proposed method is also not trivial, and we leave the development of a more rigorous and/or sophisticated estimate for future work. With regard to the error estimate, there will be an interesting connection with the Gauss quadratures (cf. [8, 9, 28]).

Acknowledgments

The author would like to thank Editage for English language editing, Professor Zdeněk Strakoš for giving the author a chance to attend his course which motivated the author to study this work, and Professor Ken Hayami for helpful discussions.

References

[1] Bai, Z., Golub, G.: Computation of large-scale quadratic forms and transfer functions using the theory of moments, quadrature and Padé approximation. In: Modern Methods in Scientific Computing and Applications, pp. 1–30. Springer Netherlands, Dordrecht (2002). DOI 10.1007/978-94-010-0510-4_1
[2] Brezinski, C.: The methods of Vorobyev and Lanczos. Linear Algebra Appl. 234, 21–41 (1996). DOI 10.1016/0024-3795(94)00081-6
[3] Davis, T., Hu, Y.: The University of Florida sparse matrix collection. ACM Trans. Math. Software 38(1), 1:1–1:25 (2011). DOI 10.1145/2049662.2049663
[4] Feldmann, P., Freund, R.W.: Efficient linear circuit analysis by Padé approximation via the Lanczos process. IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems 14(5), 639–649 (1995). DOI 10.1109/43.384428
[5] Freund, R.W.: Conjugate gradient-type methods for linear systems with complex symmetric coefficient matrices. SIAM J. Sci. Statist. Comput. 13(1), 425–448 (1992). DOI 10.1137/0913023
[6] Frommer, A., Kahl, K., Lippert, T., Rittich, H.: 2-norm error bounds and estimates for Lanczos approximations to linear systems and rational matrix functions. SIAM J. Matrix Anal. Appl. 34(3), 1046–1065 (2013). DOI 10.1137/110859749
[7] Frommer, A., Maass, P.: Fast CG-based methods for Tikhonov–Phillips regularization. SIAM J. Sci. Comput. 20(5), 1831–1850 (1999). DOI 10.1137/s1064827596313310
[8] Golub, G.H., Meurant, G.: Matrices, moments and quadrature, vol. 303, pp. 105–156. Longman Sci. Tech., Harlow (1993)
[9] Golub, G.H., Meurant, G.: Matrices, moments and quadrature II; How to compute the norm of the error in iterative methods. BIT 37(3), 687–705 (1997). DOI 10.1007/bf02510247
[10] Golub, G.H., Meurant, G.: Matrices, Moments and Quadrature with Applications. Princeton University Press, Princeton, NJ (2010)
[11] Gu, G.D., Liu, G.: On convergence property of MINRES method for solving a complex shifted Hermitian linear system. Int. J. Math. Comput. Sci. 7, 290–294 (2013). DOI 10.5281/ZENODO.1087887
[12] Hamburger, H.: Beiträge zur Konvergenztheorie der Stieltjesschen Kettenbrüche. Math. Z. 4, 186–222 (1919). DOI 10.1007/BF01203012
[13] Hamburger, H.: Über eine Erweiterung des Stieltjesschen Momentenproblems. Math. Ann. 81(2-4), 235–319 (1920). DOI 10.1007/bf01564869
[14] Hamburger, H.: Über eine Erweiterung des Stieltjesschen Momentenproblems. Math. Ann. 82(1-2), 120–164 (1920). DOI 10.1007/bf01457982
[15] Hamburger, H.: Über eine Erweiterung des Stieltjesschen Momentenproblems. Math. Ann. 82(3-4), 168–187 (1921). DOI 10.1007/bf01498663
[16] Hoffman, K., Kunze, R.: Linear Algebra, 2nd edn. Prentice-Hall Inc., Englewood, Cliffs, NJ (1971)
[17] Horn, R.A., Johnson, C.R.: Matrix Analysis, 2nd edn. Cambridge University Press, New York, NY (2013). DOI 10.1017/9781139020411
[18] Hoshi, T., Imachi, H., Kuwata, A., Kakuda, K., Fujita, T., Matsui, H.: Numerical aspect of large-scale electronic state calculation for flexible device material. Jpn. J. Ind. Appl. Math. 36(2), 685–698 (2019). DOI 10.1007/s13160-019-00358-2
[19] Hoshi, T., Kawamura, M., Yoshimi, K., Motoyama, Y., Misawa, T., Yamaji, Y., Todo, S., Kawashima, N., Sogabe, T.: K $\omega$ – Open-source library for the shifted Krylov subspace method of the form $(z{I}-{H})x=b$ . Comput. Phys. Commun. p. 107536 (2020). DOI 10.1016/j.cpc.2020.107536
[20] Johndrow, J.E., Mattingly, J.C., Mukherjee, S., Dunson, D.B.: Optimal approximating Markov chains for Bayesian inference. arXiv 1508.03387 [stat.CO], 1–31 (2017)
[21] Kaniel, S.: Estimates for some computational techniques in linear algebra. Math. Comp. 20, 369–378 (1966). DOI 10.1090/S0025-5718-1966-0234618-4
[22] Lanczos, C.: Solution of systems of linear equations by minimized iterations. J. Res. Nat. Bur. Stand. 49(1), 33 (1952). DOI 10.6028/jres.049.006
[23] Li, C., Sra, S., Jegelka, S.: Gaussian quadrature for matrix inverse forms with applications. In: M.F. Balcan, K.Q. Weinberger (eds.) Proceedings of The 33rd International Conference on Machine Learning, pp. 1766–1775. PMLR, New York, USA (2016)
[24] Liesen, J., Strakoš, Z.: Krylov Subspace Methods, Principles and Analysis. Oxford University Press, Oxford (2013)
[25] Maeda, Y., Futamura, Y., Imakura, A., Sakurai, T.: Filter analysis for the stochastic estimation of eigenvalue counts. JSIAM Lett. 7(0), 53–56 (2015). DOI 10.14495/jsiaml.7.53
[26] Meerbergen, K.: The solution of parametrized symmetric linear systems. SIAM J. Matrix Anal. Appl. 24(4), 1038–1059 (2003). DOI 10.1137/s0895479800380386
[27] Meurant, G.: The Lanczos and Conjugate Gradient Algorithms. SIAM, Philadelphia, PA (2006). DOI 10.1137/1.9780898718140
[28] Meurant, G., Tichý, P.: On computing quadrature-based bounds for the ${A}$ -norm of the error in conjugate gradients. Numer. Algorithms 62(2), 163–191 (2013). DOI 10.1007/s11075-012-9591-9
[29] Musco, C., Musco, C., Sidford, A.: Stability of the Lanczos method for matrix function approximation. In: Proceedings of the Twenty-Ninth Annual ACM-SIAM Symposium on Discrete Algorithms, pp. 1605–1624. Society for Industrial and Applied Mathematics (2018). DOI 10.1137/1.9781611975031.105
[30] Paige, C.C., Saunders, M.A.: Solution of sparse indefinite systems of linear equations. SIAM J. Numer. Anal. 12, 617–629 (1975). DOI 10.1137/0712047
[31] Polizzi, E.: Density-matrix-based algorithm for solving eigenvalue problems. Phys. Rev. B 79(11), 115112–1–6 (2009). DOI 10.1103/physrevb.79.115112
[32] Rump, S.M.: Verification of positive definiteness. BIT 46(2), 433–452 (2006). DOI 10.1007/s10543-006-0056-1
[33] Saad, Y.: Iterative Methods for Sparse Linear Systems, 2nd edn. SIAM, Philadelphia, PA (2003). DOI 10.1137/1.9780898718003
[34] Sakurai, T., Sugiura, H.: A projection method for generalized eigenvalue problems using numerical integration. J. Comput. Appl. Math. 159(1), 119–128 (2003). DOI 10.1016/S0377-0427(03)00565-X
[35] Sakurai, T., Tadano, H.: CIRR: a Rayleigh–Ritz type method with contour integral for generalized eigenvalue problems. Hokkaido Math. J. 36, 745–757 (2007). DOI 10.14492/hokmj/1272848031
[36] Seito, H., Hoshi, T., Yamamoto, Y.: On using the shifted minimal residual method for quantum-mechanical wave packet simulation. JSIAM Lett. 11(0), 13–16 (2019). DOI 10.14495/jsiaml.11.13
[37] Sogabe, T., Hoshi, T., Zhang, S., Fujiwara, T.: On a weighted quasi-residual minimization strategy for solving complex symmetric shifted linear systems. Electron. Trans. Numer. Anal. 31, 126–140 (2008)
[38] Sogabe, T., Zhang, S.L.: A COCR method for solving complex symmetric linear systems. J. Comput. Appl. Math. 199(2), 297–303 (2007). DOI 10.1016/j.cam.2005.07.032
[39] Sogabe, T., Zhang, S.L.: An extension of the COCR method to solving shifted linear systems with complex symmetric matrices. East Asian J. Appl. Math. 1(2), 97–107 (2011). DOI 10.4208/eajam.260410.240510a
[40] Strakoš, Z.: Model reduction using the Vorobyev moment problem. Numer. Algorithms 51(3), 363–379 (2009). DOI 10.1007/s11075-008-9237-0
[41] Strakoš, Z., Tichý, P.: On error estimation in the conjugate gradient method and why it works in finite precision computations. Electron. Trans. Numer. Anal. 13, 56–80 (2002)
[42] Szegö, G.: Orthogonal Polynomials, vol. XXIII. American Mathematical Society, New York (1959). DOI 10.1090/coll/023
[43] Takayama, R., Hoshi, T., Sogabe, T., Zhang, S.L., Fujiwara, T.: Linear algebraic calculation of the Green’s function for large-scale electronic structure theory. Phys. Rev. B 73(165108), 1–9 (2006). DOI 10.1103/physrevb.73.165108
[44] Vorobyev, Y.V.: Method of Moments in Applied Mathematics. Gordon and Breach Science Publishers, New York (1965)
[45] van der Vorst, H.A., Melissen, J.: A Petrov–Galerkin type method for solving ${Ax=b}$ , where ${A}$ is symmetric complex. IEEE Trans. Magn. 26(2), 706–708 (1990). DOI 10.1109/20.106415
[46] Yamamoto, S., Sogabe, T., Hoshi, T., Zhang, S.L., Fujiwara, T.: Shifted conjugate-orthogonal–conjugate-gradient method and its application to double orbital extended Hubbard model. J. Phys. Soc. Jpn. 77(11), 114713 (2008). DOI 10.1143/jpsj.77.114713

	$\displaystyle\\|\boldsymbol{x}_{k}-\boldsymbol{x}_{*}\\|$	$\displaystyle=\\|V_{k}(T_{k,k}^{<})^{-1}(\\|\boldsymbol{v}\\|\boldsymbol{e}_{1})-S^{-1}\boldsymbol{v}\\|$		(3.40)
		$\displaystyle=\\|V_{k}(T_{k,k}^{<})^{-1}\boldsymbol{e}_{1}-S^{-1}\boldsymbol{v}_{1}\\|\\|\boldsymbol{v}\\|.$		(3.41)

Shifted Lanczos method for quadratic forms with Hermitian matrix resolvents

Abstract

1 Introduction

2 Lanczos method

2.1 Hamburger moment problem

2.2 Model reduction via Vorobyev moment matching

3 Shifted Lanczos method

3.1 Implementation

Remark.

3.2 Breakdown-free condition

Lemma 3.1.

Proof.

Theorem 3.2.

Proof.

3.3 Convergence bound

Theorem 3.3.

Proof.

3.4 Estimation of error

4 Numerical experiments

4.1 Comparisons in terms of CPU time

4.2 Spectrum and shift

4.3 Estimation of error

5 Conclusions

Acknowledgments

References

Shifted Lanczos method for
quadratic forms with Hermitian matrix resolvents