A Chebyshev--Jackson series based block SS--RR algorithm for computing partial eigenpairs of real symmetric matrices

TL;DR

This paper introduces a Chebyshev--Jackson series-based block SS--RR algorithm for efficiently computing partial eigenpairs of large symmetric matrices within a specific interval, avoiding expensive linear system solutions.

Contribution

It develops a novel CJ series expansion approach for moments in the block SS--RR method, improving efficiency over traditional contour integral techniques.

Findings

The CJ series converges pointwise to moments with quantifiable rates.

The new algorithm outperforms the trapezoidal rule-based method in numerical tests.

Theoretical convergence results extend to higher-order moments.

Abstract

This paper considers eigenpair computations of large symmetric matrices with the desired eigenvalues lying in a given interval using the contour integral-based block SS--RR method, a Rayleigh--Ritz projection onto a certain subspace generated by moment matrices. Instead of using a numerical quadrature to approximately compute the moments by solving a number of large shifted complex linear systems at each iteration, we make use of the Chebyshev--Jackson (CJ) series expansion to approximate the moments, which only involves matrix-vector products and avoids expensive solutions of the linear systems. We prove that the CJ series expansions pointwise converge to the moments as the series degree increases, but at different convergence rates depending on point positions and moment orders. These extend the available convergence results on the zeroth moment of CJ series expansions to higher order ones. Based on the results established, we develop a CJ--SS--RR algorithm. Numerical experiments illustrate that the new algorithm is more efficient than the contour integral-based block SS--RR algorithm with the trapezoidal rule.

Tables

Table 1. Table 1: The computational cost of one iteration of Algorithm   2

Steps	MVs	Other flops
3	$d\ell$	$2dn\ell +2dnM\ell$
4		$2n{\left(M\ell \right)}^2$
5	$M\ell$	$4n{\left(M\ell \right)}^2+8⁤\frac{2}{3}{\left(M\ell \right)}^3$
6		$2n{\left(M\ell \right)}^2$
Total cost	$\left(\frac{d}{M}+1\right)M\ell$	$2\left(d+\frac{d}{M}\right)nM\ell +8n{\left(M\ell \right)}^2+8⁤\frac{2}{3}{\left(M\ell \right)}^3$

Table 2. Table 2: Properties of the test matrices, where λ min \lambda_{\min} , λ max \lambda_{\max} and n e ​ v n_{ev} are computed by eig .

Matrix A	Size	$nnz\left(A\right)$	${\lambda }_{\min }$	${\lambda }_{\max }$	$[a,b]$	$n_{ev}$
SiH4	5041	171903	-0.9956	36.79	$[3.7,4.7]$	115
SiNa	5743	198787	-0.7042	25.62	$[2,2.8]$	134
delaunay_n13	8192	49094	-3.564	6.457	$[2.4,2.8]$	214
benzene	8219	242669	-0.7297	58.40	$[4,5.1]$	115
stokes64	12546	140034	-1.696e-3	3.998	$[1.9,2.1]$	232
Pres_Poisson	14822	715804	1.278e-5	26.03	$[10,12]$	25
Si10H16	17077	875923	-1.149	36.93	$[3.7,4.7]$	371
brainpc2	27607	179395	-2000	4465	$[-10,10]$	88
rgg_n_2_15_s0	32768	320480	-5.119	17.36	[10, 11]	221
SiO	33401	1317655	-1.675	84.32	$[2,4]$	332

Table 3. Table 3: The condition numbers of the moment matrices generated by different polynomials p k ​ ( t ) p_{k}(t) .

$M$	$p_k\left(t\right)$	The matrices
		SiH4	SiNa	delaunay_n13	benzene	stokes64
4	$t^k$	2.5e+10	2.1e+10	3.6e+8	1.9e+11	5.6e+7
	${\left(\tilde{l}\left(t\right)\right)}^k$	1.8e+5	2.0e+5	2.6e+4	3.4e+5	9.5e+3
	$T_k\left(\tilde{l}\left(t\right)\right)$	5.7e+4	5.3e+4	7.8e+3	1.2e+5	3.4e+3
8	$t^k$	7.0e+15	6.0e+15	5.7e+15	4.4e+15	2.1e+14
	${\left(\tilde{l}\left(t\right)\right)}^k$	7.8e+6	6.6e+6	9.5e+5	2.0e+7	3.6e+5
	$T_k\left(\tilde{l}\left(t\right)\right)$	1.2e+5	1.1e+5	1.2e+4	3.3e+5	5.3e+3
16	$t^k$	6.6e+16	5.1e+16	8.1e+22	2.0e+16	7.0e+28
	${\left(\tilde{l}\left(t\right)\right)}^k$	7.7e+9	8.6e+9	6.7e+9	2.6e+10	1.6e+10
	$T_k\left(\tilde{l}\left(t\right)\right)$	1.8e+5	2.0e+5	1.8e+5	5.4e+5	3.8e+5

Table 4. Table 4: The numerical ranks of the moment matrices generated by different polynomials p k ​ ( t ) p_{k}(t) , where the boldface numbers denote numerical rank deficiency.

$M$	$p_k\left(t\right)$	The matrices
		SiH4	SiNa	delaunay_n13	benzene	stokes64
4	$t^k$	176	208	336	176	352
	${\left(\tilde{l}\left(t\right)\right)}^k$	176	208	336	176	352
	$T_k\left(\tilde{l}\left(t\right)\right)$	176	208	336	176	352
8	$t^k$	134	163	296	132	344
	${\left(\tilde{l}\left(t\right)\right)}^k$	176	208	336	176	352
	$T_k\left(\tilde{l}\left(t\right)\right)$	176	208	336	176	352
16	$t^k$	88	104	168	79	176
	${\left(\tilde{l}\left(t\right)\right)}^k$	176	208	336	176	352
	$T_k\left(\tilde{l}\left(t\right)\right)$	176	208	336	176	352
$M\ell$		176	208	336	176	352
$n_{ev}$		115	134	214	115	232

Table 5. Table 5: The performance of Algorithm   2 on the test matrix stokes64.

$M$	Using the constant $d$					Using the changing $d$
	$d$	Iter	Residual	Total MVs		$d$	Iter	Residual	Total MVs
1	310	5	1.34e-12	642345.0		211	8	3.00e-11	700433.8
2	310	5	1.49e-12	345798.8		212	8	1.71e-11	379261.0
4	310	5	9.70e-13	197525.7		220	8	1.88e-11	225104.7
8	310	5	2.84e-12	123389.1		259	7	3.30e-12	144730.8
16	310	12	6.92e-11	207169.9		433	5	3.73e-12	119872.4

Table 6. Table 6: The results of block SS–RR and Algorithm   2 with t ​ o ​ l = 10 − 12 tol=10^{-12} . The stopping tolerance of minres in block SS–RR is 10 − 12 10^{-12} . Speedup ratios are the ratios of MVs used by Algorithm   2 and those used by minres .

Matrix	The total MVs						Speedup ratios
	SS–RR			CJ–SS–RR
	$M=4$	$M=8$		$M=4$	$M=8$		$M=4$	$M=8$
SiH4	1.1e+6	7.6e+5		2.1e+5	9.9e+4		5.1	7.6
SiNa	3.2e+6	1.7e+6		2.0e+5	1.3e+5		15.9	12.7
delaunay_n13	1.3e+7	6.7e+6		5.0e+5	3.3e+5		26.1	20.3
benzene	2.5e+6	1.8e+6		3.4e+5	2.1e+5		7.3	8.3
stokes64	8.0e+6	4.1e+6		2.8e+5	2.1e+5		28.5	19.7
Pres_Poisson	7.7e+4	3.9e+4		1.4e+4	7.3e+3		5.5	5.3
Si10H16	2.2e+7	1.1e+7		6.4e+5	3.9e+5		33.8	27.5
brainpc2	1.6e+7	1.1e+7		3.1e+6	2.2e+6		5.0	5.0
rgg_n_2_15_s0	9.0e+6	5.8e+6		3.2e+5	2.4e+5		28.1	24.6
SiO	9.8e+6	4.9e+6		6.9e+5	4.3e+5		14.2	11.4

Equations

$$\hat{S}=({\hat{S}}_0,{\hat{S}}_1,\dots ,{\hat{S}}_{M-1})$$

$${\hat{S}}_k=\frac{1}{2\pi \mathrm{i}}{\oint }_{\partial \Omega }{z}^k(zI-A{)}^{-1}V\mathrm{d}z,k=0,1,\dots ,M-1,$$

$${\hat{S}}_k\approx S_k:=\sum _{j=1}^q{w}_j{z}_j^k(z_{j}I-A{)}^{-1}V,k=0,1,\dots ,M-1,$$

$$S=(S_0,S_1,\dots ,S_{M-1})$$

$$l(t)=\frac{2t-{\lambda }_{\max }-{\lambda }_{\min }}{{\lambda }_{\max }-{\lambda }_{\min }}$$

$$h(t)=\{\begin{array}{ll}1,& t\in (a,b),\\ \frac{1}{2},& t\in \{a,b\},\\ 0,& t\in [-1,1]\setminus [a,b],\end{array}$$

$$\frac{1}{2}\left(h(a+0)+h(a-0)\right)=\frac{1}{2}\left(h(b+0)+h(b-0)\right)=h(a)=h(b)=\frac{1}{2}.$$

$$\frac{1}{2\pi \mathrm{i}}{\oint }_{\partial \Omega }{z}^k(z-t{)}^{-1}\mathrm{d}z=\{\begin{array}{ll}{t}^k,& t\in (a,b),\\ 0,& t\in [-1,1]\setminus [a,b].\end{array}$$

$${\hat{S}}_k=A^k{\hat{S}}_0=A^{k}h(A)V,\mathcal{R}(\hat{S})={\mathcal{K}}_M(A,h(A)V),$$

$$\hat{S}=K_M(A,h(A)V)=X_{n_{ev}}{K}_M({\Lambda }_{n_{ev}},h({\Lambda }_{n_{ev}})X_{n_{ev}}^{T}V),$$

$$\mathcal{R}(\hat{S})\subseteq {\mathcal{X}}_{n_{ev}},$$

$$\mathrm{rank}\left(K_M({\Lambda }_{n_{ev}},h({\Lambda }_{n_{ev}})X_{n_{ev}}^{T}V)\right)=n_{ev},$$

$$\mathcal{R}(\hat{S})={\mathcal{K}}_M(A,h(A)v)=\{p(A)h(A)v\mid p\in {\mathcal{P}}_{M-1}\},$$

$$\tilde{h}(t)=p(t)h(t)\approx F_d(p)(t):=\frac{c_0}{2}+\sum _{j=1}^d{\rho }_{j,d}{c}_j{T}_j(t),$$

$$c_j=\frac{2}{\pi }{\int }_{-1}^1\frac{\tilde{h}(t)T_j(t)}{\sqrt{1-t^2}}\mathrm{d}t=\frac{2}{\pi }{\int }_a^b\frac{p(t)T_j(t)}{\sqrt{1-t^2}}\mathrm{d}t,$$

$${\rho }_{j,d}=\frac{\sin ((j+1){\alpha }_d)}{(d+2)\sin ({\alpha }_d)}+\left(1-\frac{j+1}{d+2}\right)\cos (j{\alpha }_d),\text{where }{\alpha }_d=\frac{\pi }{d+2}.$$

$$F_d({\beta }_1{p}_1+{\beta }_2{p}_2)={\beta }_1{F}_d(p_1)+{\beta }_2{F}_d(p_2).$$

$$g(\theta )=\tilde{h}(\cos \theta )=\{\begin{array}{ll}p(\cos \theta ),& \theta \in (-\alpha ,-\beta )\cup (\beta ,\alpha ),\\ \frac{1}{2}p(\cos \theta ),& \theta \in \{\pm \alpha ,\pm \beta \},\\ 0,& \theta \in [-\pi ,-\alpha )\cup (-\beta ,\beta )\cup (\alpha ,\pi ],\end{array}$$

$$\alpha =\arccos (a),\beta =\arccos (b)$$

$$q_d(\theta )=F_d(p)(\cos \theta )=\frac{c_0}{2}+\sum _{j=1}^d{\rho }_{j,d}{c}_j\cos (j\theta ).$$

$$q_d(\theta )=\frac{1}{\pi }{\int }_{-\pi }^{\pi }g(\theta +\varphi )u_d(\varphi )\mathrm{d}\varphi ,$$

$$u_d(\varphi )=\frac{1}{2}+\sum _{j=1}^d{\rho }_{j,d}\cos (j\varphi ).$$

$$\frac{1}{\pi }{\int }_{-\pi }^{\pi }g(\varphi )\cos (j\varphi )\mathrm{d}\varphi =c_j,\frac{1}{\pi }{\int }_{-\pi }^{\pi }g(\varphi )\sin (j\varphi )\mathrm{d}\varphi =0.$$

$$q_d(\theta )$$

$$=\frac{1}{\pi }{\int }_{-\pi }^{\pi }g(\varphi )u_d(\varphi -\theta )\mathrm{d}\varphi =\frac{1}{\pi }{\int }_{-\pi -\theta }^{\pi -\theta }g(\theta +\varphi )u_d(\varphi )\mathrm{d}\varphi .$$

$$u_d(\varphi )\ge 0,$$

$$\frac{1}{\pi }{\int }_{-\pi }^{\pi }{u}_d(\varphi )\mathrm{d}\varphi =1,$$

$$\frac{1}{\pi }{\int }_{-\pi }^{\pi }|\varphi |u_d(\varphi )\mathrm{d}\varphi \le \frac{{\pi }^2}{2(d+2)},$$

$$\frac{1}{\pi }{\int }_{-\pi }^{\pi }{\varphi }^2{u}_d(\varphi )\mathrm{d}\varphi \le \frac{{\pi }^4}{4(d+2{)}^2},$$

$$\frac{1}{\pi }{\int }_{-\pi }^{\pi }{\varphi }^4{u}_d(\varphi )\mathrm{d}\varphi \le \left(\frac{1}{4}+\frac{(d+1){\pi }^2}{16(d+2{)}^2}\right)\frac{{\pi }^6}{(d+2{)}^3}\le \frac{{\pi }^6}{2(d+2{)}^3}.$$

$$u_d(\varphi )=\left(\sum _{k=0}^d{t}_k{\mathrm{e}}^{\mathrm{i}k\varphi }\right)\left(\sum _{k=0}^d{t}_k{\mathrm{e}}^{-\mathrm{i}k\varphi }\right)\ge 0,\text{where }{t}_k=\frac{\sin \frac{k+1}{d+2}\pi }{\sqrt{2{\sum }_{k=0}^d{\sin }^2\frac{k+1}{d+2}\pi }},$$

$$\frac{1}{\pi }{\int }_{-\pi }^{\pi }|\varphi |u_d(\varphi )\mathrm{d}\varphi \le {\left(\frac{1}{\pi }{\int }_{-\pi }^{\pi }{\varphi }^2{u}_d(\varphi )\mathrm{d}\varphi \right)}^{1/2},$$

$$\frac{1}{\pi }{\int }_{-\pi }^{\pi }{\varphi }^4{u}_d(\varphi )\mathrm{d}\varphi \lesssim \frac{{\pi }^6}{4(d+2{)}^3}.$$

Full text

{AMS}

65F15, 15A18, 65F10, 41A10

1 Introduction

The problem of computing a number of specific eigenvalues of a large matrix is common in vast scientific and engineering applications. In this paper, we consider the numerical solution of the following symmetric eigenvalue problem $Ax=\lambda x$ with $A\in\mathbb{R}^{n\times n}$: Given a real interval $[a,b]$ contained in the spectrum interval of $A$, determine all the $n_{ev}$ eigenpairs $(\lambda,x)$ with the eigenvalues $\lambda\in[a,b]$ counting multiplicities. Such kind of eigenvalue problems arises from numerous applications, e.g., dynamic analysis of structures [10], linear response eigenvalue problems in quantum chemistry [3], and the density functional theory based electronic structure calculations [31], to mention only a few.

Over the past two decades, starting with the initial the Sakurai–Sugiura (SS) method [32], called the SS–Hankel method, and the FEAST eigensolver [26], a new class of numerical methods that are based on contour integration and rational filtering has emerged for computing the eigenvalues of a large matrix in a given region and the associated eigenvectors or invariant subspaces. The SS–Hankel method uses certain complex moment matrices constructed by a contour integral with the contour containing the region of interest that the eigenvalues lie in, but it and its block variant [15] are numerically unstable because of the ill-conditioning of Hankel matrices [4]. To improve numerical stability and compute eigenpairs more accurately and reliably, based on the Rayleigh–Ritz projection, the SS–RR method [33] and its block variant [14] have been proposed; the block SS–Arnoldi method based on the block Arnoldi process has also been presented in [16], and it is mathematically equivalent to the SS-RR method when the moments are exactly computed. The FEAST eigensolver, initially introduced by Polizzi [26] in 2009 for Hermitian eigenvalue problems, only computes an approximation to the zeroth moment, which is the exact spectral projector corresponding to the eigenvalues inside the interval of interest. It works on subspaces of a fix dimension, and uses subspace iteration [9, 24, 31, 34] on an approximate spectral projector to generate a sequence of subspaces [25], onto which the Rayleigh–Ritz projection is realized. The FEAST eigensolver has been generalized to non-Hermitian matrices [23, 36], and it can be regarded as a restarted variant of the block SS–RR method, in which only the zeroth moment is required and the subspace dimension is no smaller than the number $n_{ev}$ of the desired eigenvalues.

In the SS eigensolver [8] and FEAST eigensolver [27], one needs to use some suitably chosen numerical quadratures to approximate the zeroth to higher order moment matrices and the zeroth order moment matrix, respectively, which requires the solutions of a number of large shifted complex linear systems at each iteration, where the shifts are the quadrature nodes. Besides, several researchers have proposed different rational approximation or filtering methods [1, 11, 35], which also require to solve a number of large shifted linear systems. These linear systems are generally solved by iterative methods or solved by a direct solver, e.g., the sparse LU factorizations [9, 30] if they are computationally feasible. For truly large matrices, if they are structured, e.g., banded, then the sparse LU factorizations are preferable and can solve the shifted linear systems very efficiently. But if they are generally sparse or dense, the storage requirement and computational cost of the LU factorizations is prohibitive, and we have no way but use an iterative method, e.g., a Krylov subspace iterative method such as BiCGStab method [30], to solve the linear systems. We suppose that the shifted linear systems are solved iteratively throughout this paper. However, notice that these shifted linear systems are highly indefinite when the region of the interest is inside the spectrum of $A$. It is well known that, for highly indefinite linear systems, Krylov subspace solvers are generally inefficient and can be too slow. Furthermore, it lacks a general effective preconditioning technique for highly indefinite linear systems.

Noticing that the contour integral corresponds to a specific step or piecewise continuous function, Jia and Zhang have proposed a new FEAST SVDsolver based on Chebyshev–Jackson (CJ) series expansion [21, 22]; instead of numerical quadrature or any rational filtering, they make use of the polynomial CJ-series expansion to construct an approximation to the step function and obtain an approximation to the underlying spectral projector without solving any linear system. The resulting CJ-FEAST solver is directly adaptable to real symmetric eigenvalue problems. They have established the pointwise convergence, error estimates of the CJ series, and accuracy estimates for the approximate spectral projector. The results show that the convergence of CJ polynomial series expansion is uniquely determined by the distances between the interval ends and the closest eigenvalues inside and outside the interval and is independent of the interval location in the spectrum of $A$. By contrast, iterative solutions of the shifted linear systems in numerical quadratures and other rational approximations critically depends on both the interval location and the distances between the ends and the closest eigenvalues inside and outside, which determine the degrees of indefiniteness and ill conditioning of the shifted linear systems, respectively, and affect the convergence of Krylov solvers.

In this paper, following CJ series expansion approximations [21, 22] to the step function and the zeroth order moment, we nontrivially extend them to those higher order moments in the block SS–RR methods, and propose a new block SS–RR method, called the CJ–SS–RR method. We prove that the CJ series expansion pointwise converges to the moments as the series degree increases, but convergence rates depend on point positions and moment orders. Given the property that the subspace generated in the block SS–RR method is a block Krylov subspace [15, 18] when the moments are exactly computed, we extend the results of the block Krylov subspace in [20, 34] to the CJ–SS–RR method, and establish the convergence results that are different from those of the SS-type methods in [15, 17, 18]. These results show that, just as the CJ-FEAST SVDsolver, the approximation accuracy of the CJ series expansion does not depend on the interval location in the spectrum of $A$. We consider the practical choice of series degree and the issue of more reliable computation of orthonormal basis of a given search subspace by replacing its possibly ill conditioned basis matrix formed by the monomial moments with a better conditioned basis matrix formed by the better conditioned modified Chebyshev moments, such that the subspace generated by the computed orthonormal basis is a more accurate approximation to the true search subspace, thereby retaining more effective information and improving the computational accuracy. Finally, we will report numerical experiments to demonstrate the effectiveness and efficiency of the CJ–SS–RR method; in the meantime, we illustrate that the new method is much more efficient than the contour integral-based block SS–RR method using the numerical quadrature with the trapezoidal rule.

The paper is organized as follows. In Section 2, we briefly describe the block SS–RR method and its theoretical background. In Section 3, we establish compact quantitative pointwise convergence results on the CJ series expansion. Then we develop a CJ–SS–RR method in Section 4 and analyze its convergence. In Section 5, we present a complete CJ–SS–RR method with several practical issues addressed. In Section 6, we report numerical experiments to illustrate the performance of the CJ–SS–RR method and to compare it with the contour integral-based block SS–RR method using the trapezoidal rule. Finally, we conclude the paper in Section 7.

Throughout this paper, denote by $\|\cdot\|_{2}$ the 2-norm of a vector or matrix, by $\mathcal{R}(X)$ the column space of a matrix $X$, by $\lambda(A)$ the set of all the eigenvalues of $A$, and by $\lambda_{\min}$ and $\lambda_{\max}$ the minimum and maximum eigenvalues of $A$, respectively. We denote by $(\lambda_{i},x_{i}),i=1,\ldots,n$ the eigenpairs of $A$ with $\|x_{i}\|_{2}=1$, label the eigenvalues in the interval $[a,b]$ as $b\geq\lambda_{1}\geq\cdots\geq\lambda_{n_{ev}}\geq a$, and write $X_{n_{ev}}=(x_{1},\ldots,x_{n_{ev}})$, $\Lambda_{n_{ev}}=\operatorname{diag}\{\lambda_{1},\ldots,\lambda_{n_{ev}}\}$ and $\mathcal{X}_{n_{ev}}=\mathcal{R}(X_{n_{ev}})$.

2 The block SS–RR method

Define $\Omega=\{|z-(a+b)/2|\leq(b-a)/2\}\subset\mathbb{C}$ to be the circle that contains $[a,b]$ and goes through $a$ and $b$. Given two positive integers $M,\ell$ and a matrix $V\in\mathbb{R}^{n\times\ell}$ of full column rank, define the moment matrix

[TABLE]

with the monomial moments

[TABLE]

where $\partial\Omega$ is the positively-oriented boundary of $\Omega$. The block SS–RR method [14] is the Rayleigh–Ritz projection onto the subspace $\mathcal{R}(\widehat{S})$ and computes the Ritz pairs with respect to the subspace. Since the exact $\widehat{S}$ is not available, the block SS–RR method uses a numerical quadrature to compute the integrals approximately:

[TABLE]

where the $w_{j}$ and $z_{j}$ are the quadrature weights and nodes, respectively. Therefore,

[TABLE]

is an approximation to $\widehat{S}$ in (1).

Next we review some basic results and provide some preliminaries. For brevity, suppose $[\lambda_{\min},\lambda_{\max}]=[-1,1]$, and let $[a,b]\subset[\lambda_{\min},\lambda_{\max}]$; otherwise, we use the linear transformation

[TABLE]

to map $[\lambda_{\min},\lambda_{\max}]$ to $[-1,1]$. Define the step function

[TABLE]

where $a$ and $b$ are the discontinuity points of $h(t)$, and the values of $h(t)$ at $a$ and $b$ equal the means of respective right and left limits:

[TABLE]

The Cauchy integral formula tells us

[TABLE]

Thus, according to [15, Theorem 4], we have

[TABLE]

where $\mathcal{K}_{M}(A,h(A)V)$ is the block Krylov subspace generated by $A$ and $h(A)V$.

Since $h(A)=X_{n_{ev}}h(\Lambda_{n_{ev}})X_{n_{ev}}^{T}$, we can write the matrix $\widehat{S}$ defined in Eq. 1 as

[TABLE]

where $K_{M}(A,V)=(V,AV,\ldots,A^{M-1}V)$ denotes the block Krylov matrix generated by $A$ and the starting block $V$. The following theorem (cf. [17, Theorem 1]) shows that the block SS–RR method onto $\mathcal{R}(\widehat{S})$ solves the eigenvalue problem under consideration exactly.

Theorem 2.1**.**

Let $V\in\mathbb{R}^{n\times\ell}$. Then

[TABLE]

*and the equality is attained if and only if $\operatorname{rank}(\widehat{S})=n_{ev}$. *

The condition $\operatorname{rank}(\widehat{S})=n_{ev}$ is equivalent to

[TABLE]

which indicates that $M\ell\geq M\cdot\operatorname{rank}(X_{n_{ev}}^{T}V)\geq n_{ev}$ and requires that no multiplicity of $\lambda_{i}\in[a,b]$ be greater than $\ell$, a necessary condition for a block Krylov to determine the multiplicities of desired eigenvalues [20].

Particularly, for $V=v\in\mathbb{R}^{n}$, suppose that $\widehat{S}$ is of column full rank. Then

[TABLE]

where $\mathcal{P}_{M-1}$ is the set of real polynomials of degree not exceeding $M-1$.

3 The CJ series expansion

Let $p\in\mathcal{P}_{M-1}$ be a real polynomial, and define the operator $F_{d}$ as a mapping from $p$ to an approximate expansion of $\tilde{h}(t)=p(t)h(t)$ obtained by the $d$-degree CJ polynomial series expansion:

[TABLE]

where $T_{j}$ is the $j$-degree Chebyshev polynomial of the first kind [28], the coefficients

[TABLE]

and the Jackson damping factors (cf. [7, 19])

[TABLE]

Since the coefficients $c_{j}$ are linear with respect to the polynomial $p$, the mapping $F_{d}$ is linear in $p$ too: for all real numbers $\beta_{1},\beta_{2}$ and real polynomials $p_{1},p_{2}$, it holds that

[TABLE]

Define the step function $g$ with the period $2\pi$ to be

[TABLE]

where

[TABLE]

and $g(\theta)$ at $\pm\alpha$ and $\pm\beta$ are the means of respective right and left limits, and the trigonometric polynomial

[TABLE]

Analogously to Lemma 1.4 of [28, Section 1.1.2], we have the following result.

Lemma 3.1**.**

Let $g$ and $q_{d}$ be defined as Eqs. 12 and 14, respectively. Then

[TABLE]

where

[TABLE]

Proof 3.2**.**

Since $g$ is even, from Eq. 9 we obtain

[TABLE]

Substituting the above relations into Eq. 14 yields

[TABLE]

*Because of the periodicity of $g$ and $u_{d}$, Eq. 15 follows directly from the above relation. *

For $u_{d}$, there are the following properties.

Proposition 3.3**.**

Let $u_{d}$ be defined as Eq. 16 and $d\geq 2$. Then

[TABLE]

Proof 3.4**.**

It is known from [28, Section 1.1.2] that

[TABLE]

which yields Eq. 17. Relation Eq. 18 is from the proof of [21, Theorem 3.1]. By (18) and the Cauchy–Schwarz inequality, relation (1.1.24) of [28, Section 1.1.2] shows that

[TABLE]

*and (1.1.25) in this book establishes an expression of its bound, whose estimates Eqs. 19 and 20 are from the proof of [21, Theorem 3.2]. Finally, Eq. 21 is straightforward from the proof of [21, Theorem 3.2]. *

The intermediate bound in Eq. 21 can be asymptotically simplified. For a modestly sized $d$, e.g., say 50, ignoring the higher order term in it gives

[TABLE]

We next establish quantitative pointwise convergence results of the approximation $q_{d}$ to $g$.

Theorem 3.5**.**

Let $g$ and $q_{d}$ be as defined in Eqs. 12 and 14, respectively, and $\alpha$ and $\beta$ be as defined in Eq. 13. Define

[TABLE]

and the $L_{\infty}$-norm of $g$ on an interval $(\phi_{1},\phi_{2})$ as

[TABLE]

Then for $\theta\in[0,\pi]\setminus[\beta,\alpha]$ and $d\geq 2$, we have

[TABLE]

Proof 3.6**.**

Due to Eq. 15, Eq. 17 and $g(\theta)=0$ for $\theta\in[0,\pi]\setminus[\beta,\alpha]$, we have

[TABLE]

According to Eq. 12, $g(\varphi)=0$ for $\varphi\in(-\beta,\beta)\cup(\alpha,2\pi-\alpha)$. For $|\phi|<\Delta_{\theta}$, since $\theta\in[0,\pi]\setminus[\beta,\alpha]$, we have $\theta+\phi\in(-\beta,\beta)\cup(\alpha,2\pi-\alpha)$ and hence $g(\theta+\phi)=0$. Therefore,

[TABLE]

For $\Delta_{\theta}\leq|\phi|\leq\pi$, we have

[TABLE]

Then combining the above three relations and Eq. 21, we obtain

[TABLE]

*which proves (23). *

Remark 3.7**.**

*As is seen from the proof, the difference $g(\theta+\phi)-g(\theta)$ on $(-\Delta_{\theta},\Delta_{\theta})$ is zero for $\theta\in[0,\pi]\setminus[\beta,\alpha]$ while, if $p$ in Eq. 12 is not constant, it is nonzero for $\theta\in[\beta,\alpha]$, which implies that $|g(\theta+\phi)-g(\theta)|$ may be larger when $\theta\in[\beta,\alpha]$ than it is when $\theta\in[0,\pi]\setminus[\beta,\alpha]$. This is indeed the case, as the following two theorems indicate. *

Theorem 3.8**.**

With the same notations as in Theorem 3.5, for $\theta\in(\beta,\alpha)$ and $d\geq 2$, we have

[TABLE]

Proof 3.9**.**

From Eqs. 15 and 18, we obtain

[TABLE]

Below we analyze the above two terms in $q_{d}(\theta)-g(\theta)$ separately. From the nonnegativity of $u_{d}$ shown in Eq. 17, we have

[TABLE]

For $\Delta_{\theta}\leq|\phi|\leq\pi$, it holds that

[TABLE]

Thus from the above relations and Eq. 21, we obtain

[TABLE]

On the other hand, since $u_{d}$ defined in Eq. 16 is even, we have

[TABLE]

Let $G(\phi)=g(\theta+\phi)+g(\theta-\phi)-2g(\theta)$. For $|\phi|<\Delta_{\theta}$, we have $\theta\pm\phi\in(\beta,\alpha)$, meaning that $g(\theta\pm\phi)=p(\cos(\theta\pm\phi))$ and $G(\phi)$ is thus smooth. Since $G(0)=G^{\prime}(0)=0$, there is a $\mu\in(0,1)$ such that

[TABLE]

which shows that

[TABLE]

Therefore,

[TABLE]

Making use of Eq. 20 and observing that the integration region is half, we obtain

[TABLE]

*which, together with Eq. 25, establishes Eq. 24. *

Theorem 3.10**.**

With the same notations as in Theorem 3.5, let

[TABLE]

Then for $d\geq 2$, we have

[TABLE]

Proof 3.11**.**

Since $g(\alpha)=\frac{1}{2}(g(\alpha+0)+g(\alpha-0))$, from Eqs. 15 and 18, we obtain

[TABLE]

We first consider the integral over $[0,\pi]$. For $0<\phi<\Delta_{\alpha}$, we have $\alpha+\phi\in(\alpha,2\pi-\alpha)$. Therefore, it follows from Eq. 12 that $g(\alpha+\phi)=g(\alpha+0)=0$. For $\Delta_{\alpha}<\phi<\pi$, we have

[TABLE]

Then by Eq. 21, we obtain

[TABLE]

For $-\Delta_{\alpha}<\phi<0$, we have $\alpha+\phi\in(\beta,\alpha)$. Since $p$ is smooth, there is a $\mu\in(0,1)$ such that

[TABLE]

Then combining the above relation with Eq. 19 yields

[TABLE]

Similarly to (25), we can establish the following bound:

[TABLE]

*Thus Eq. 26 follows from the above two relations and Eq. 28. The proof of Eq. 27 is analogous. *

Remark 3.12**.**

*If $M=1$, then $p\in\mathcal{P}_{M-1}$ is a constant and the second terms in (24), (26) and (27) vanish. As a result, the preceding two theorems degenerate to [21, Theorems 3.2 and 3.3], where the error bounds for $\theta\in[0,\pi]$ tend to zero as fast as $(d+2)^{-3}$. However, as long as $p$ is non-constant, the $L_{\infty}$-norms of $g^{\prime}$ and $g^{\prime\prime}$ are nonzero, and the second terms in (24), (26) and (27) will dominate the error bounds for sufficiently large $d$, and the convergence becomes slow. *

We present an example to illustrate that the convergence rates in above theorems are indeed order optimal. Take $[a,b]=[-0.2,0.4]\subset[-1,1]$ and three points $t=-0.6,-0.2,0.1$, of which $-0.6$ and $0.1$ are outside and inside $[a,b]$, respectively. Polynomials $p(t)$ are the Chebyshev polynomials of the first kind on $[a,b]$. Note that $q_{d}(\theta)=F_{d}(p)(\cos\theta)$, $g(\theta)=p(\cos\theta)h(\cos\theta)$, $\alpha=\arccos(-0.2)$, $\beta=\arccos(0.4)$ and $\theta=\arccos(t)$. We approximately replace the $L_{\infty}$-norms of $g^{\prime}$ and $g^{\prime\prime}$ by $\|p^{\prime}\|_{L_{\infty}(a,b)}$ and $\|p^{\prime}\|_{L_{\infty}(a,b)}+\|p^{\prime\prime}\|_{L_{\infty}(a,b)}$, respectively. For each of $t$ and $p$, we plot the true errors $|F_{d}(p)(t)-p(t)h(t)|$ and the corresponding bounds for $d=2,3,\ldots,10^{4}$ in Fig. 1.

Clearly, our bounds are order optimal since they accurately reflect the convergence rates of the true errors as the degree $d$ increases. Precisely, for $t\not\in[a,b]$ or $\deg p=0$, the convergence is as fast as $(d+2)^{-3}$; for $t\in[a,b]$ and $\deg p\geq 1$, the convergence rates are the same as the second terms in the right-hand sides of Eqs. 24 and 26.

4 The polynomial approximation in block SS–RR and its convergence

Although Theorem 2.1 shows that the search subspace $\mathcal{R}(\widehat{S})$ is the desired eigenspace if $\operatorname{rank}(\widehat{S})=n_{ev}$, the exact computation of integrals $\widehat{S}_{k}$ in (6) is impractical, so is $p(A)h(A)$. Therefore, we have to consider their approximations and the resulting approximation to $p(A)h(A)$ for $p\in\mathcal{P}_{M-1}$, and study the convergence of the block SS–RR method.

Given a matrix $V\in\mathbb{R}^{n\times\ell}$, rather than using the rational approximations $S_{k}$ in (3) to construct an approximation $S$ to $\widehat{S}$ (cf. (4) and Eq. 1), based on the CJ series expansion (8), we construct a polynomial approximation $S$ to $\widehat{S}$, where

[TABLE]

Relationships (6) and (29) between $\mathcal{R}(S)$ and the block Krylov subspace $\mathcal{R}(\widehat{S})$ enable us to establish the following theorem, which is an analogue of the results on the block Krylov subspace [20, Lemma 1 and Theorem 2], where the assumption on the matrix $B$ and the definition of the index set $N_{i}$ are exactly the same as those in Lemma 1 of [20] when the underlying matrix $A$ is symmetric.

Theorem 4.1**.**

Let $S$ be defined as Eq. 29 and $N_{i}=\{1,2,\ldots,n\}\setminus\{i,\ldots,i+\ell-1\}$, and suppose that $B=(x_{i},\ldots,x_{i+\ell-1})^{T}V$ is nonsingular, where $x_{1},\ldots,x_{n}$ are the orthonormal eigenvectors of $A$ corresponding to the eigenvalues $\lambda_{1},\ldots,\lambda_{n}$ labeled in any desired order. Then there exists a unique $v_{i}=VB^{-1}e_{i}\in\operatorname{\mathcal{R}}(V)$ such that

[TABLE]

and

[TABLE]

Proof 4.2**.**

Let $p(t)=\beta_{0}+\beta_{1}t+\cdots+\beta_{M-1}t^{M-1}\in\mathcal{P}_{M-1}$. Then by Eq. 11, we have

[TABLE]

Therefore,

[TABLE]

Relation (30) is the restatement of [20, Lemma 1 and Theorem 2] for $A$ symmetric. Therefore, we can express $v_{i}=x_{i}+\sum_{j\in N_{i}}\xi_{j}x_{j}$, where the two terms are orthogonal and the coefficients $\xi_{j}=x_{j}^{T}v_{i}$.

Since $F_{d}(p)(A)x_{j}=F_{d}(p)(\lambda_{j})x_{j}$, we have

[TABLE]

where the two terms are orthogonal. As a result,

[TABLE]

*Rhelation Eq. 31 follows from the above relation and Eq. 32. *

To analyze the bound in Eq. 31, we next exploit Theorems 3.5, 3.8, and 3.10 to establish error bounds for the difference between $F_{d}(p)$ and $p\cdot h$ at the eigenvalues of $A$. Note that these three theorems concern the error bounds for $|q_{d}(\theta)-g(\theta)|$, while our ultimate focus is error bounds for $|F_{d}(p)(\lambda)-p(\lambda)h(\lambda)|$. To this end, we first present the following two lemmas to estimate the $L_{\infty}$-norms involved.

Lemma 4.3**.**

Suppose $a<b$ and $M\geq 2$. Then for $p\in\mathcal{P}_{M-1}$, it holds that

[TABLE]

where

[TABLE]

Proof 4.4**.**

The value of $T_{M-1}^{(k)}(1)$ sees relation (1.97) of [29]. For $a=-1$ and $b=1$, Eq. 33 is the Markov’s theorem (cf. [29, p. 123]):

[TABLE]

Notice that the linear transformation $l(t)=(2t-a-b)/(b-a)$ maps $[a,b]$ to $[-1,1]$. Then for a given $p\in\mathcal{P}_{M-1}$, we have the one-one correspondence $p(t)=\tilde{p}(l(t))$ for some $\tilde{p}\in\mathcal{P}_{M-1}$. Therefore,

[TABLE]

which leads to

[TABLE]

Applying the Markov’s theorem to $\tilde{p}$, we obtain

[TABLE]

*Thus Eq. 33 follows from the above two relations. *

Lemma 4.5**.**

Let $-1\leq a<b\leq 1$, $M\geq 1$ and $p\in\mathcal{P}_{M-1}$. Then

[TABLE]

Proof 4.6**.**

For $M=1$, we have $p^{\prime}(t)\equiv p^{\prime\prime}(t)\equiv 0$, and thus relation Eq. 34 is trivial. For $M\geq 2$, we apply Eq. 33 to $p$ and get

[TABLE]

Since $b-a\leq 2$ and $M-1\geq 1$, we have

[TABLE]

*From the above two relations, we obtain Eq. 34. *

With the preceding two lemmas, we have the following theorem.

Theorem 4.7**.**

Let $\lambda_{ia},\lambda_{ib}$ and $\lambda_{oa},\lambda_{ob}$ be the eigenvalues of $A$ closest to the ends $a$ and $b$ inside and outside $[a,b]$, respectively. Define

[TABLE]

where $\alpha=\arccos(a),\beta=\arccos(b)$. Then for $\lambda\in\lambda(A)$ and $d\geq 2$, the error

[TABLE]

Proof 4.8**.**

Let $\lambda$ be an eigenvalue of $A$, and denote $\theta=\arccos(\lambda)$. We obtain

[TABLE]

It follows from Theorems 3.5, 3.8, and 3.10 that

[TABLE]

where $g(\theta)=p(\cos(\theta))h(\cos(\theta))=p(\lambda)h(\lambda)$.

Next, we investigate sizes of the $L_{\infty}$-norms of the two derivatives in (36). By definition Eq. 12, we have $\|g^{\prime}\|_{L_{\infty}(\beta,\alpha)}=\|g^{\prime\prime}\|_{L_{\infty}(\beta,\alpha)}=0$ for $M=1$ and

[TABLE]

and these immediately yield Eq. 35 when $M=1$. Consequently, we only need to consider the case $M\geq 2$. For $\phi\in(\beta,\alpha)$, by definition (12), we have

[TABLE]

which, together with Eqs. 33 and 34, yield

[TABLE]

*Substituting Eqs. 37 and 38 into Eq. 36 proves Eq. 35. *

Since the ratio $|F_{d}(p)(\lambda_{j})/F_{d}(p)(\lambda_{i})|$ is independent of the scaling of $p$, i.e.,

[TABLE]

for any real number $c\not=0$, it is insightful to obtain more from Theorem 4.7: Define

[TABLE]

where $\Delta_{\min}$ is defined as in Theorem 4.7. Then

[TABLE]

where $\lambda_{i}\in[a,b]$. We see from Eq. 39 that the parameters $\gamma_{i}(p)$, $\delta_{i}(p)$ and $\eta_{i}(p)$ tend to zero as $d$ increases. Thus for sufficiently large $d$, the ratio $|F_{d}(p)(\lambda_{j})/F_{d}(p)(\lambda_{i})|$ is governed by $h(\lambda_{j})|p(\lambda_{j})/p(\lambda_{i})|$, which vanishes for $\lambda_{j}\not\in[a,b]$. In this case, we can present the following result, which is analogous to [20, Theorem 3] but is stated in a slightly different form.

Theorem 4.9**.**

Assume that $b\geq\lambda_{1}\geq\cdots\geq\lambda_{n_{ev}}\geq a$ and the multiplicities of $\lambda_{i},\,i=1,\ldots,n_{ev}$ do not exceed $\ell$. Given $1\leq i\leq n_{ev}$, let $\Xi_{i}$ be the set of the distinct ones among $\lambda_{1},\ldots,\lambda_{i-1}$, and suppose that $\lambda_{i-1}\neq\lambda_{i}$ for $i>1$ and $|\Xi_{i}|\leq M-1$, where $|\Xi_{i}|$ denotes the cardinality of $\Xi_{i}$. Define

[TABLE]

Then

[TABLE]

Proof 4.10**.**

The case $i=1,\ldots,n_{ev}-\ell$ has been proved in [20, Theorem 3]. Here we only need to show the case $i=n_{ev}-\ell+1,\ldots,n_{ev}$. For $i=n_{ev}-\ell+1,\ldots,n_{ev}$, there is $j\in\{1,\ldots,n_{ev}\}\setminus\{i,\ldots,i+\ell-1\}=\{1,\ldots,i-1\}$ so that $\lambda_{j}\in\Xi_{i}$. In this case, taking $p(t)=m_{i}(t)$ gives

[TABLE]

*Therefore, (41) holds for $1\leq i\leq n_{ev}$. *

With the preceding theorems, we can establish the following estimate for the bound in Eq. 31. From Eq. 40, we can see that for $\lambda_{j}\not\in[a,b]$ the upper bound of $|F_{d}(p)(\lambda_{j})/p(\lambda_{i})|$ is $\gamma_{i}(p)$, independently of the position of $\lambda_{j}$. Thus for the concerning interval eigenvalue problem of $A$, we label the eigenvalues of $A$ inside $[a,b]$ in the same order as in Theorem 4.9 and those outside $[a,b]$ at will.

Theorem 4.11**.**

Label the eigenvalues of $A$ in $[a,b]\subset[\lambda_{\min},\lambda_{\max}]$ as $\lambda_{1},\ldots,\lambda_{n_{ev}}$ with $b\geq\lambda_{1}\geq\cdots\geq\lambda_{n_{ev}}\geq a$, suppose that the multiplicity of each of them does not exceed $\ell$, and write eigenvalues outside $[a,b]$ as $\lambda_{n_{ev}+1},\ldots,\lambda_{n}$ at will. Given $1\leq i\leq n_{ev}$, assume that $\lambda_{i-1}\neq\lambda_{i}$ for $i>1$ and the cardinality $|\Xi_{i}|\leq M-1$, where $\Xi_{i}$ is the set of the distinct eigenvalues among $\lambda_{1},\ldots,\lambda_{i-1}$. Define $\sigma_{i}$, $m_{i}(t)$, $\kappa_{i}$ and $\varepsilon_{i}$ as in Theorem 4.9, and

[TABLE]

where $\Delta_{\min}$ is defined as in Theorem 4.7. Then for sufficiently large CJ series degree $d$ such that $\nu_{i}$ defined in (44) is uniformly positive, we have

[TABLE]

with $N_{i}=\{1,2,\ldots,n\}\setminus\{i,\ldots,i+\ell-1\}$ and

[TABLE]

Proof 4.12**.**

Case I*: $i=1,\ldots,n_{ev}-\ell$. For $j\in N_{i}$, there are three possibilities: $1\leq j\leq i-1$, $i+\ell\leq j\leq n_{ev}$ and $j>n_{ev}$, which correspond to the cases $\lambda_{j}\in\Xi_{i}$, $\lambda_{j}\in[a,\lambda_{i+\ell}]$ and $\lambda_{j}\not\in[a,b]$, respectively. Let $\hat{l}(t)=1+\frac{2(t-\lambda_{i+\ell})}{\lambda_{i+\ell}-a}$ be a linear transformation mapping $[a,\lambda_{i+\ell}]$ to $[-1,1]$ and*

[TABLE]

Then $|T_{M-|\Xi_{i}|-1}(\hat{l}(t))|\leq 1$ for $t\in[a,\lambda_{i+\ell}]$, $|T_{M-|\Xi_{i}|-1}(\hat{l}(t))|\leq T_{M-|\Xi_{i}|-1}(\hat{l}(b))$ for $t\in[a,b]$, and $\hat{l}(\lambda_{i})=\sigma_{i}$. Therefore, for $\lambda_{j}\in[a,\lambda_{i+\ell}]$, we obtain

[TABLE]

Since $p(\lambda_{j})=0$ for $\lambda_{j}\in\Xi_{i}$ and $h(\lambda_{j})=0$ for $\lambda_{j}\not\in[a,b]$, we have

[TABLE]

Notice that the only difference between Eq. 39 and Eq. 42 is the substitution of $\frac{\|p\|_{L_{\infty}(a,b)}}{|p(\lambda_{i})|}$ with $\tau_{i}$. Therefore,

[TABLE]

indicating that

[TABLE]

Substituting Eqs. 45 and 46 into Eq. 40 yields

[TABLE]

that is,

[TABLE]

Case II*: $i=n_{ev}-\ell+1,\ldots,n_{ev}$. In this case, for $j\in N_{i}$, there is $j\in\{1,\ldots,i-1\}$ or $j\in\{i+\ell,\ldots,n\}$ such that $\lambda_{j}\in\Xi_{i}$ or $\lambda_{j}\not\in[a,b]$. Take $p(t)=m_{i}(t)$. Then for $j\in N_{i}$, we have*

[TABLE]

since $p(\lambda_{j})=0$ for $\lambda_{j}\in\Xi_{i}$ and $h(\lambda_{j})=0$ for $\lambda_{j}\not\in[a,b]$, and

[TABLE]

*which yields $\gamma_{i}(p)=\hat{\gamma}_{i}$, $\delta_{i}(p)=\hat{\delta}_{i}$ and $\eta_{i}(p)=\hat{\eta}_{i}$. Substitute these relations into Eq. 40. We obtain Eq. 47 again. *

Remark 4.13**.**

*The quantity $\varepsilon_{i}$ is an upper bound of the ratio $|p(\lambda_{j})/p(\lambda_{i})|$ and $\tau_{i}$, which, together with the degree $d$, forms $\hat{\gamma}_{i}$, $\hat{\delta}_{i}$ and $\hat{\eta}_{i}$ and gives a bound for $\|p\|_{L_{\infty}(a,b)}/|p(\lambda_{i})|$. Moreover, the parameters $\hat{\gamma}_{i}$, $\hat{\delta}_{i}$ and $\hat{\eta}_{i}$ decide the upper bounds of $|F_{d}(p)(t)-p(t)h(t)|/|p(\lambda_{i})|$ for $t$ outside, inside $[a,b]$ and at the ends $a$ and $b$, respectively. *

Remark 4.14**.**

*This theorem considers only the case $M\geq 2$. For $M=1$, we set $N_{i}=\{n_{ev}+1,\ldots,n\}$; Eq. 43 still holds for $\mu_{i}=\hat{\gamma}_{i}$ and $\nu_{i}=h(\lambda_{i})-\hat{\gamma}_{i}$. In this manner, relation Eq. 43 follows from the convergence of the FEAST solver based on subspace iteration [21]. *

Remark 4.15**.**

We can write the ratio $\mu_{i}/\nu_{i}$ as

[TABLE]

*Since $\lim_{d\to\infty}\mu_{i}=\varepsilon_{i}$ and $\lim_{d\to\infty}\nu_{i}=h(\lambda_{i})$, the first term in the above relation reduces to Theorem 4.9 as $d\to\infty$, and it tends to zero as $M$ increases. The second and third terms tend to zero as $d$ increases, and are the deviations caused by the CJ series expansion approximation. Therefore, the ratio $\mu_{i}/\nu_{i}$ decreases as $d$ increases, meaning faster convergence. However, unlike the series degree $d$, the parameter $M$ is not so free to change and must be of small to modest size, and actually one only requires that $M\ell\geq n_{ev}$ slightly in the block SS–RR method. *

5 The CJ–SS–RR algorithm

Recall that the linear transformation (5) maps the spectral interval $[\lambda_{\min},\lambda_{\max}]$ to $[-1,1]$. When using this transformation in computations, it is necessary to give reasonably good estimates for $\lambda_{\max}$ and $\lambda_{\min}$. As commonly suggested in a large number of literature, e.g., [7, 37], one can run the symmetric Lanczos method on $A$ for several steps, say $30\sim 50$, to get estimates of $\lambda_{\max}$ and $\lambda_{\min}$. For instance, $\lambda_{\max}\approx\theta_{1}+\|r_{1}\|_{2}$ and $\lambda_{\min}\approx\theta_{s}-\|r_{s}\|_{2}$, where $\theta_{1}$ and $\theta_{s}$ are the largest and smallest Ritz values and $r_{1},r_{s}$ are the corresponding residuals of the $s$-step Lanczos method. In what follows, we always assume that the spectrum interval of $A$ has been transformed into $[-1,1]$.

5.1 The construction of an orthonormal basis of the search subspace

Given the search subspace $\mathcal{R}(S)$ with $S$ defined in Eq. 29, to obtain approximate eigenpairs, let the columns of $U$ form an orthonormal basis of $\mathcal{R}(S)$. We use the standard Rayleigh–Ritz projection

[TABLE]

to obtain the approximate eigenpairs $(\tilde{\lambda},\tilde{x})$, called the Ritz pairs of $A$ with respect to $\mathcal{R}(S)$. Theoretically, a commonly used approach is to compute a rank-revealing QR decomposition or the singular value decomposition (SVD) of $S$ to get an orthonormal basis of $\mathcal{R}(S)$. In finite precision arithmetic, however, notice that

[TABLE]

The matrix $\widehat{S}$ is typically ill conditioned even for $M$ small and its condition number increases at exponential rate with $M$ increasing [4], so does that of $S$. As a consequence, the subspace generated by the computed orthonormal basis vectors of $\mathcal{R}(S)$ may be far from the true $\mathcal{R}(S)$ [12, pp. 373-4]. Realize the ill conditioning of $S$. It is common in the SS–RR type methods [14, 17, 18] to use the truncated SVD and take those dominant right singular vectors associated with the large singular values as a replacement of $\mathcal{R}(S)$, whose dimension is generally at most a couple of dozens for $\ell$ small when $\mathcal{R}(S)$ is a good approximation to $\mathcal{R}(\widehat{S})$ and may be much lower than the dimension $M\ell$ of $\mathcal{R}(S)$. Consequently, such approach possibly discards too much effective information, and impairs the strength of the SS–RR method artificially and severely. Even worse, it is well possible that the dimension of effective subspace is smaller than the number $n_{ev}$ of desired eigenpairs when $n_{ev}\gg 1$, thereby leading to the failure of the block SS–RR method when computing the $n_{ev}$ desired eigenpairs.

In order to fix this prevailing severe implementation deficiency, computationally, it is very attempting to take new well-conditioned basis vectors of $\mathcal{R}(S)$ rather than the ill-conditioned ones generated by the monomials $1,t,\ldots,t^{M-1}$ and then orthonormalize them to obtain a new orthonormal basis vectors that retain all the information of $\mathcal{R}(S)$. In finite precision arithmetic, such orthonormal basis vectors generate $\mathcal{R}(S)$ to the working precision when the condition number of the new basis matrix is of size $O(1)$. The Chebyshev basis vectors of $\mathcal{R}(S)$ are better choices, and they form a better conditioned matrix. Precisely, we replace $t^{k},k=0,1,\ldots,M-1$ by the $k$-degree Chebyshev polynomials $p_{k}=T_{k}((2t-a-b)/(b-a))$ in (29) and obtain the modified moments. Numerical experiments in Section 6.1 will confirm that the new basis vectors are always numerically linearly independent and thus retain all the information of $\mathcal{R}(S)$.

Algorithm 1 describes the construction of the new basis matrix

[TABLE]

where we need to compute the Chebyshev coefficients of $p_{k},k=0,\ldots,M-1$:

[TABLE]

In computations, we compute them by the MATLAB built-in function integral in the MATLAB computing environment.

5.2 The estimate for the number of desired eigenvalues

By Theorem 2.1, the condition $M\ell\geq n_{ev}$ is indispensable for computing the $n_{ev}$ desired eigenpairs. But the value of $n_{ev}$ is generally unknown in practice. Therefore, we have to estimate $n_{ev}$ in advance.

Note that

[TABLE]

Combine the above relation with $F_{d}(1)(A)\approx h(A)$, and suppose that the eigenvalues of $A$ equal to $a$ or $b$ are simple. Then

[TABLE]

We estimate the trace by Monte-Carlo simulation [2, 5, 13]. Given $s$ random vectors $v_{i}$ with the components $v_{ij}$ being i.i.d. Rademacher random variables, i.e., $v_{ij}=\pm 1$ with probability $1/2$, we estimate $n_{ev}$ using

[TABLE]

More details of the estimation can be found in [21].

5.3 A practical choice of degree $d$ of the CJ series expansion

We now investigate how to select the degree $d$ of the CJ series expansion to ensure that $\mu_{i}/\nu_{i}<\zeta$ for a given $\zeta\in(0,1)$. The following derivations apply to $i=1,\ldots,n_{ev}$. Therefore, we only need to consider $i=1$, and drop the subscripts for brevity. For $M\geq 2$, (44) means

[TABLE]

From (44), suppose that $d$ is sufficiently large so that $\frac{1}{2}>\hat{\gamma}+\hat{\eta}>\hat{\eta}$, which ensures $\nu>0$. Then the expressions of $\hat{\delta}$ and $\hat{\eta}$ in (42) mean

[TABLE]

The sufficient conditions for (50) are

[TABLE]

For $M=1$, according to Remark 4.14, we have $2\hat{\gamma}<\frac{\zeta}{1+\zeta}$.

Substituting Eq. 42 into the above relations yields

[TABLE]

where, by Eq. 22, we rewrite $\hat{\gamma}$ as $\hat{\gamma}=\frac{\pi^{6}}{\Delta_{\min}^{4}(d+2)^{3}}\cdot\frac{\tau}{2}$, and $\tau=1$ when $M=1$ or $\lambda_{1}=b$.

Without a priori knowledge of the eigenvalue distribution, suppose that the desired $\lambda_{i}$ are uniformly distributed in $[a,b]$. Then

[TABLE]

Since $M\ell\geq n_{ev}$, for $M\geq 2$ we have

[TABLE]

which shows

[TABLE]

Thus, setting the maximum $\zeta=1$, we obtain a sufficient condition for Eq. 51:

[TABLE]

However, as illustrated in Fig. 1, the above bound, especially its first term, is overestimated; that is, the estimated degree $d$ may be too large to be feasible in practice. We will propose a practical estimation of $d$ as follows: Choose

[TABLE]

where $D\in[1,8]$ and $K\in[1,10]$.

5.4 A complete algorithm

With the algorithmic implementation details described previously, we can run the CJ–SS–RR algorithm for given $M$ and $\ell$. If it converges, we stop the algorithm and output the converged approximate eigenpairs. Otherwise, for $M\ell$ fixed, it is necessary to restart the algorithm to compute better approximations until convergence. In this paper, we use the restart approach proposed in [17, 18]: At the $k$th restart, we take $V^{(k)}=F_{d}(1)(A)V^{(k-1)}$, or equivalently, the first $\ell$ columns of $\widetilde{S}^{(k-1)}$. We use the stopping criteria in [36] to test convergence, and describe the complete algorithm as Algorithm 2. The computational cost of one iteration of Algorithm 2 is listed in Table 1, where MVs are the number of matrix-vector products and other flops are the flops spent on the other computations.

Let $nnz(A)$ be the number of nonzero entries in $A$ and suppose $d>M\ell$; in fact, generally $d\gg M\ell$ in computations. Notice that $n\gg M\ell$. Then it is seen from Table 1 that one iteration costs $O(nnz(A)d\ell)+O(dnM\ell)+O(n(M\ell)^{2})$ flops, where the first two terms consumed by Step 3 dominate the overall cost. As a consequence, we use the cost of Step 3 to represent the overall cost of one iteration. For ease of counting and comparison, we measure the computational cost in terms of MVs. Since one matrix-vector product of $A$ costs $2nnz(A)$ flops and

[TABLE]

the cost of the other flops in Step 3 is approximately $\frac{(M+1)n}{nnz(A)}d\ell$ MVs, which are not necessarily integers. We call them the equivalent MVs. As a result, we measure the overall efficiency of Algorithm 2 by the total MVs, i.e., the actual MVs of one iteration plus the equivalent MVs of Step 3.

6 Numerical experiments

We report numerical experiments to confirm our theoretical results and illustrate the performance of Algorithm 2. The test matrices are from the SuiteSparse Matrix Collection [6], and we list some of their basic properties and the intervals of interest in Table 2. All the numerical experiments were performed on an Intel Core i7-9700, CPU 3.0GHz, 8GB RAM using MATLAB R2024b with the machine precision $\epsilon_{\rm mach}=2.22\times 10^{-16}$ under the Microsoft Windows 10 64-bit system. To make a fair comparison, for each test problem and given $\ell$, we used the same starting $n\times\ell$ matrix $V^{(0)}$ for the algorithms tested, which was generated randomly in a normal distribution. For reference, we use the built-in function eig in MATLAB to compute $\lambda_{\min}$, $\lambda_{\max}$ and the exact $n_{ev}$ and obtain $\|A\|_{2}=\max\{|\lambda_{\min}|,|\lambda_{\max}|\}$.

Based on Section 5, in the experiments, we apply the linear transformation $l(t)$ in LABEL:eq:_lineartransformation to $A$ such that the spectrum interval of $l(A)$ is $[-1,1]$ and $[a,b]$ is mapped to $[\tilde{a},\tilde{b}]\subset[-1,1]$, where $\tilde{a}=l(a)$ and $\tilde{b}=l(b)$. An approximate eigenpair $(\tilde{\lambda},\tilde{x})$ is claimed to have converged if its relative residual norm satisfies

[TABLE]

6.1 Condition numbers of different bases and effective dimensions of the search subspaces

We first investigate the condition numbers of different kinds of matrices $S$ and effective dimensions of $\mathcal{S}$ constructed by the CJ series approximations to the moments corresponding to different kinds of polynomials $p_{k}(t),k=0,1,\ldots,M-1$, where the numerical ranks, i.e., the effective dimensions of $\mathcal{R}(S)$’s, are determined by the truncated SVDs with the tolerance $10\|S\|\epsilon_{\rm mach}$. Tables 4 and 3 display the condition numbers and effective subspace dimensions obtained by using the monomial moments, the shifted-and-scaled $(\tilde{l}(t))^{k}$ with $\tilde{l}(t)=(2t-\tilde{a}-\tilde{b})/(\tilde{b}-\tilde{a})$ and the shifted-and-scaled Chebyshev moments for the first five matrices in Table 2 with $M=4,8,16$, where, for each test matrix, we set the true dimension $M\ell=16\lceil 1.5n_{ev}/16\rceil$ of $\mathcal{R}(S)$, and determined $\ell$ by the equations $M\ell=16\lceil 1.5n_{ev}/16\rceil$. We took the series degree $d$ in Eq. 53 as $d=\lceil\pi^{2}(b-a)^{-4/3}+\pi^{2}(b-a)^{-1}\rceil$.

As is seen from the tables, the condition numbers of $S$’s using the monomial moments are already bigger than $O(1/\sqrt{\epsilon_{\rm mach}})$ even for $M=4$, and for $M=8$ these $S$ are already numerically rank deficient. For $M=16$, the effective search subspaces obtained by the truncated SVD’s of $S$ loses too much information on $\mathcal{R}(S)$ and effective dimensions of $\mathcal{R}(S)$’s are already smaller than $n_{ev}$’s, causing that Algorithm 2 cannot compute the $n_{ev}$ desired eigenpairs of $A$.

In contrast, the situation changes very substantially when using the shifted-and-scaled moments, and the condition numbers are reduced greatly. However, it is still far from good enough to solve the concerning eigenvalue problems because the condition numbers are around $10^{10}$ for $M=16$ and the errors of the computed subspaces by the rank-revealing QR decompositions and the true subspaces $\mathcal{R}(S)$ are the condition numbers times $\epsilon_{\rm mach}$ and thus around $10^{-5}$.

Strikingly, the shifted-and-scaled Chebyshev moments is very superior to the shifted-and-scaled moments, and the condition numbers are further significantly reduced and are around $10^{5}$, approximately the square roots of the latter ones for $M=16$. As a result, the computed subspaces are much more accurate with the errors $10^{-10}$, so that we can compute the desired eigenpairs much more accurately. For better stability and accuracy, we always use the Chebyshev polynomials $T_{k}(\tilde{l}(t))$ to generate $S$ in the subsequent experiments.

6.2 The performance of the CJ–SS–RR algorithm

We test the matrix stokes64 for $M=1,2,4,8,16$ by fixing the subspace dimension $M\ell=16\lceil 1.5n_{ev}/16\rceil$ and using the stopping tolerance $tol=10^{-10}$ in (55). In the experiments, two options are available for the series degree $d$. The first is to use the constant degree $d=\lceil\pi^{2}(b-a)^{-4/3}+\pi^{2}(b-a)^{-1}\rceil$ for all the moment matrices, and the second is to use the changing series degree $d$ determined by Eq. 53 with $D=1$ and $K=10$ so that the moment matrices are generated using the CJ series expansions of the changing degree $d$. The number of restarts denoted by Iter, the maximum relative residuals and the total MVs are listed in Table 5.

Relation Eq. 53 indicates that the degree $d$ has to increase quadratically with $M$ so as to ensure a similar convergence. Consequently, the fixed constant degree $d$ may slow down convergence when $M>1$ considerably. The results in Table 5 confirm the considerable slowdown of convergence of the constant $d$ for $M=16$. The results also demonstrate that for a given $M$, a higher degree $d$ leads to a faster convergence, and vice versa, which is consistent with Remark 4.15.

Table 5 illustrate that the total MVs used by Algorithm 2 are notably fewer for $M=2$ than for $M=1$. On the contrary, the total MVs for $M=8$ and $M=16$ are comparable. This can be explained by the computational cost list in Table 1, where, for relatively large $M$, the reduction in total MVs due to the increase in $M$ is offset by increasing the degree $d$. Since $d$ increases quadratically with $M$, the total MVs of an excessively large $M$ may instead be bigger. We thus suggest that a modest value of $M$, e.g., $4\sim 8$, be used.

6.3 A comparison of the CJ–SS–RR algorithm and the SS–RR algorithm

In this subsection we compare Algorithm 2 with the block SS–RR method [8, 14]. Here, for a given interval $[a,b]$ of interest, the contour is a circle with the center $(a+b)/2$ and radius $(b-a)/2$. The numerical quadrature rule used is the trapezoidal rule with sixteen nodes, which are default parameters in [8]. Notice that $A$ is symmetric and the quadrature nodes are symmetric with respect to the real axis, the block SS–RR method manages to solve only $M\ell\cdot\frac{q}{2}$ linear systems at each iteration, where $q$ is the number of quadrature nodes. The shifted linear systems are solved by the MATLAB built-in function minres with the stopping tolerance $10^{-12}$.

We tested the two algorithms for $M=4$ and $8$ by taking $M\ell=8\lceil 1.5n_{ev}/8\rceil$ to ensure that $M\ell\geq 1.5n_{ev}$. The convergence tolerance in (55) is $tol=10^{-12}$. In Algorithm 2, we chose the degree $d$ using Eq. 53 with $D=2$ and $K=7$. Since the generally dominating cost of the contour-based block SS–RR method is iterative solutions of shifted linear systems, we only record the total MVs with $A$ used by minres, which means that we ignore other costs of the contour-based block SS–RR method when comparing the overall efficiency.

We record the total MVs of both algorithms and the speedup ratios of Algorithm 2 over the block SS–RR method in Table 6. As we see from Table 6, all the speedup ratios are no less than 5, indicating that the CJ–SS–RR algorithm is significantly more efficient than the block SS–RR method. For Pres_Poisson, only 47 of the total 14822 eigenvalues are on the right of the interval of interest, which means that the desired eigenvalues are near to the largest ones and the shifted linear systems are slightly indefinite, so that MINRES has relatively quick convergence. In this case, the speedup ratios are relatively not so substantial, though already quite considerable. However, for stokes64, Si10H16 and rgg_n_2_15_s0, where the intervals of interest correspond to the interior eigenvalues, the speedup ratios of the total MVs are significantly big. These indicate that when the interval of interest corresponds to interior eigenvalues, the high indefiniteness of the coefficient matrices of the shifted linear systems in the block SS–RR method leads to much more slow convergence of MINRES.

7 Conclusion

We have studied the problem of approximating the piecewise continuous function $p(t)h(t)$ by the CJ series expansion and proved its pointwise convergence to the function $p(t)h(t)$, showing that the convergence rates critically depend on the point positions and the degree of $p(t)$. Making use of these results, we have extended the theoretical results of the block Krylov subspace to the subspace constructed by the CJ series expansion. Based on them, we have developed a practical restart CJ–SS–RR algorithm, where we have replaced ill-conditioned basis vectors of the search subspace generated by the monomial moments by the much better conditioned Chebyshev basis vectors, which enable us to make full of the search subspace and improve the robustness of the SS–RR method and computational accuracy very greatly. In the meantime, we have proposed a practical selection strategy of the series degree $d$.

We have numerically tested our algorithm on several problems to confirm the theoretical results and analysis, and shown that the CJ–SS–RR algorithm is effective and efficient. The experiments have demonstrated that the CJ–SS–RR algorithm outperforms the block SS–RR method and often improve the overall efficiency dozens of times in terms of the total MVs, especially when the interval of interest corresponds to interior eigenvalues. In summary, the developed CJ–SS–RR algorithm is a much more efficient and robust alternative of the block SS–RR method when computing the eigenpairs of large sparse real symmetric matrices.

Declarations

The authors declare that they have no financial interests, and they read and approved the final manuscript. The algorithmic Matlab code is available upon reasonable request from the corresponding author.