Weighted Block Golub-Kahan-Lanczos Algorithms for Linear Response Eigenvalue Problem

Abstract

In order to solve all or some eigenvalues lied in a cluster, we propose a weighted block Golub-Kahan-Lanczos algorithm for the linear response eigenvalue problem. Error bounds of the approximations to an eigenvalue cluster, as well as their corresponding eigenspace, are established and show the advantages. A practical thick-restart strategy is applied to the block algorithm to eliminate the increasing computational and memory costs, and the numerical instability. Numerical examples illustrate the effectiveness of our new algorithms.

1. Introduction

In this paper, we are interested in solving the linear response eigenvalue problem (LREP):

$$\mathbf{H}z:=\left[\begin{array}{cc}0& M\\ K& 0\end{array}\right]\left[\begin{array}{c}u\\ v\end{array}\right]=\lambda \left[\begin{array}{c}u\\ v\end{array}\right]=\lambda z,$$

where K and M are $N\times N$ real symmetric positive definite matrices. Such a problem arises from studying the excitation energy of many particle systems in computational quantum chemistry and physics [1,2,3]. It also known as the Bethe-Salpeter (BS) eigenvalue-problem [4] or the random phase approximation (RPA) eigenvalue problem [5]. There has immense past and recent work in developing efficient numerical algorithms and attractive theories for LREP [6,7,8,9,10,11,12,13,14,15].

Since all the eigenvalues of $\mathbf{H}$ are real nonzero and appear in pairs $\{\lambda ,-\lambda \}$ [6], thus we order the eigenvalues in ascending order, i.e.,

$$-{\lambda }_1\le \cdots \le -{\lambda }_N<{\lambda }_N\le \cdots \le {\lambda }_1.$$

In this paper, we focus on a small portion of the positive eigenvalues for LREP, i.e., ${\lambda }_i$, $i=k,k+1,\cdots ,\ell$ with $1\le k\le \ell \le N$ and $\ell -k+1\ll N$, and their corresponding eigenvectors. We only consider the real case, all the results can be easily applied to the complex case.

The weighted Golub-Kahan-Lanczos method (wGKL) for LREP was introduced in [16]. It produces recursively a much small projection ${\mathbf{B}}_j=\left[\begin{array}{cc}0& B_j\\ B_j^T& 0\end{array}\right]$ of $\mathbf{H}$ at j-th iteration, where $B_j\in {\mathbb{R}}^{j\times j}$ is upper bidiagonal. Afterwards, the eigenpairs of $\mathbf{H}$ can be constructed by the singular value decomposition of $B_j$. The convergence analysis performs that running k iterations of wGKL is equivalently running $2k$ iterations of a weighted Lanczos algorithm for $\mathbf{H}$ [16]. Actually, $B_j$ can be also a lower bidiagonal matrix, and the same discussion can be taken place as in the case of $B_j$ is upper bidiagonal. In the following, we only consider the upper bidiagonal case.

It is well known that the single-vector Lanczos method is widely used for searching a small number of extreme eigenvalues, and it may encounter very slow convergence when the wanted eigenvalues stay in a cluster [17]. Instead, a block Lanczos method with a suitable block size is capable of computing a cluster of eigenvalues including multiple eigenvalues very quickly. Motivated by this idea, we are going to develop a block version of wGKL in [16] in order to find efficiently all or some positive eigenvalues within a cluster for LREP. Based on the standard block Lanczos convergence theory in [17], the error bounds of approximation to an eigenvalue cluster, as well as their corresponding eigenspace are established to illustrate the advantage of our weighted block Golub-Kahan-Lanczos algorithm (wbGKL).

As the increasing size of the Krylov subspace, the storage demands, computational costs, and numerical stability of a simple version of a block Lanczos method may be affected [18]. Several kinds of efficiently restarting strategies to eliminate these effects are developed for the classic Lanczos method, such as, implicitly restart method [19], thick restart method [20]. In order to make our block method more practical, and using the special structure of LREP, we consider the thick restart strategy to our block method.

The rest of this paper is organized as follows. Section 2 gives some necessary preliminaries for our later use. In Section 3, the weighted block Golub-Kahan-Lanczos algorithm (wbGKL) for LREP is presented, and its convergence analysis is discussed. Section 4 proposed the thick restart weighted block Golub-Kahan-Lanczos algorithm (wbGKL-TR). The numerical examples are tested in Section 5 to illustrate the efficiency of our new algorithms. Finally, some conclusions are given in Section 6.

Throughout this paper, ${\mathbb{R}}^{m\times n}$ is the set of all $m\times n$ real matrices, ${\mathbb{R}}^n={\mathbb{R}}^{n\times 1}$, and $\mathbb{R}={\mathbb{R}}^1$. $I_n$ (or simply I if its dimension is clear from the context) is the $n\times n$ identity matrix, and $0_{m\times n}$ is an $m\times n$ matrix of zero. The superscript “${}^T$” denotes transpose. ${\parallel ·\parallel }_F$ denotes the Frobenius norm of a matrix, and ${\parallel ·\parallel }_2$ denotes the 2-norm of a matrix or a vector. For a matrix $X\in {\mathbb{R}}^{m\times n}$, $rank(X)$ denotes the rank of X, and $\mathcal{R}(X)=span(X)$ denotes the column space of X; the submatrices $X_{i:j,:}$ and $X_{:,k:\ell }$ of X composed by the intersections of row i to row j and column k to column ℓ, respectively. For matrices or scalars $X_i$, $diag(X_1,\cdots ,X_k)$ denotes the block diagonal matrix with the i-th diagonal block $X_i$.

2. Preliminaries

For a symmetric positive definite matrix $W\in {\mathbb{R}}^{N\times N}$, the W-inner product is defined as following

$${\langle x,y\rangle }_W:=y^{T}Wx,\forall x,y\in {\mathbb{R}}^N.$$

If ${\langle x,y\rangle }_W=0$, then we denote it by $x{\perp }_{W}y$, and call it with x and y are W-orthogonal. The projector ${\Pi }_W$ is called the W-orthogonal projector onto $\mathcal{Y}$ if for any $y\in {\mathbb{R}}^N$,

$${\Pi }_{W}y\in \mathcal{Y},(I-{\Pi }_W)y{\perp }_W\mathcal{Y}.$$

For two subspaces $\mathcal{X},\mathcal{Y}\subseteq {\mathbb{R}}^N$, and suppose $k=dim(\mathcal{X})\le dim(\mathcal{Y})=\ell$, if $X\in {\mathbb{R}}^{N\times k}$ and $Y\in {\mathbb{R}}^{N\times \ell }$ are W-orthonormal basis of $\mathcal{X}$ and $\mathcal{Y}$, respectively, i.e.,

$$X^{T}WX=I_k,\mathcal{X}=\mathcal{R}(X)\mathrm{and}{Y}^{T}WY=I_{\ell },\mathcal{Y}=\mathcal{R}(Y),$$

and ${\nu }_j$ for $j=1,\cdots ,k$ with ${\nu }_1\le \cdots \le {\nu }_k$ are the singular values of $Y^{T}WX$, then the W-canonical angles ${\theta }_W^{(j)}(\mathcal{X},\mathcal{Y})$ from $\mathcal{X}$ to $\mathcal{Y}$ are defined by

$$0\le {\theta }_W^{(j)}(\mathcal{X},\mathcal{Y})=arccos{\nu }_j\le \pi /2,\mathrm{for}j=1,\cdots ,k.$$

If $k=\ell$, these angles can be said between $\mathcal{X}$ and $\mathcal{Y}$. Obviously, ${\theta }_W^{(1)}(\mathcal{X},\mathcal{Y})\ge \cdots \ge {\theta }_W^{(k)}(\mathcal{X},\mathcal{Y})$. Set

$${\mathsf{\Theta }}_W(\mathcal{X},\mathcal{Y})=diag({\theta }_W^{(1)}(\mathcal{X},\mathcal{Y}),\cdots ,{\theta }_W^{(k)}(\mathcal{X},\mathcal{Y})).$$

Especially, if $k=1$, X is a vector, there is only one W-canonical angle from $\mathcal{X}$ to $\mathcal{Y}$. In the following, we may use a matrix in one or both arguments of ${\mathsf{\Theta }}_W(·,·)$, i.e., ${\mathsf{\Theta }}_W(X,Y)$ with the understanding that it means the subspace spanned by the columns of the matrix argument.

The following two lemmas are important to our later analysis, and for proofs and more details, the reader is referred to [12,16].

Lemma 1.

([12] Lemma 3.2). Let $\mathcal{X}$ and $\mathcal{Y}$ be two subspaces in ${\mathbb{R}}^N$ with equal dimensional $dim(\mathcal{X})=dim(\mathcal{Y})=k$. Suppose ${\theta }_W^{(1)}(\mathcal{X},\mathcal{Y})<\pi /2$. Then, for any set $y_1,y_2,\cdots ,y_{k_1}$ of the basis vectors in $\mathcal{Y}$ where $1\le k_1\le k$, there is a set $x_1,x_2,\cdots ,x_{k_1}$ of linearly independent vectors in $\mathcal{X}$ such that ${\Pi }_W{x}_j=y_j$ for $1\le j\le k_1$, where ${\Pi }_W$ is the W-orthogonal projector onto $\mathcal{Y}$.

Lemma 2.

([16] Proposition 3.1). The matrix $MK$ has N position eigenvalues ${\lambda }_1^2\ge {\lambda }_2^2\ge \cdots \ge {\lambda }_N^2$ with ${\lambda }_j>0$. The corresponding right eigenvectors ${\xi }_1,\cdots ,{\xi }_N$ can be chosen K-orthonormal, and the corresponding left eigenvectors ${\eta }_1,\cdots ,{\eta }_N$ can be chosen M-orthonormal. In particular, for given $\left\{{\xi }_j\right\}$, one can choose ${\eta }_j={\lambda }_j^{-1}K{\xi }_j$, and for given $\left\{{\eta }_j\right\}$, ${\xi }_j={\lambda }_j^{-1}M{\eta }_j$, for $j=1,\cdots ,N$.

3. Weighted Block Golub-Kahan-Lanczos Algorithm

3.1. Weighted Block Golub-Kahan-Lanczos Algorithm

In this section, we plan to introduce the weighted block Golub-Kahan-Lanczos algorithm (wbGKL) for LREP, which is a block version of the weighted Golub-Kahan-Lanczos algorithm [16]. Algorithm 1 gives the process of recursively generating the M-orthonormal matrix ${\mathcal{X}}_n$, the K-orthonormal matrix ${\mathcal{Y}}_n$, and the block bidiagonal matrix ${\mathcal{B}}_n$. Giving $Y_1\in {\mathbb{R}}^{n\times n_b}$ with $Y_1^{T}K{Y}_1=I_{n_b}$, denoting $E_n^T=[0_{n_b\times (n-1)n_b},I_{n_b}]$, and

$${\mathcal{X}}_n=[X_1,\cdots ,X_n],{\mathcal{Y}}_n=[Y_1,\cdots ,Y_n],{\mathcal{B}}_n=\left[\begin{array}{cccc}{A}_1& B_1& & \\ & A_2& \ddots & \\ & & \ddots & B_{n-1}\\ & & & A_n\end{array}\right],$$

then we have the relation from Algorithm 1:

$$K{\mathcal{Y}}_n={\mathcal{X}}_n{\mathcal{B}}_n,M{\mathcal{X}}_n={\mathcal{Y}}_n{\mathcal{B}}_n^T+Y_{n+1}{B}_n^T{E}_n^T,$$

(1)

and

$${\mathcal{X}}_n^{T}M{\mathcal{X}}_n=I_{n{n}_b}={\mathcal{Y}}_n^{T}K{\mathcal{Y}}_n.$$

Remark 1.

In Algorithm 1, we only consider the case that $rank({\tilde{X}}_j)=rank({\tilde{Y}}_{j+1})=n_b$, no further treatment is provided for the cases $rank({\tilde{X}}_j)<n_b$ or $rank({\tilde{Y}}_{j+1})<n_b$. Because K and M are both symmetric positive definite, thus the two W inStep 2are both reversible.

Algorithm 1: wbGKL

1. Choose $Y_1$ satisfying $Y_1^{T}K{Y}_1=I_{n_b}$, and set $W=I_{n_b}$, $B_0=I_{n_b}$, $X_0=0_{n\times n_b}$. Compute     $F=K{Y}_1$.     2. For $j=1,2,\cdots ,n$              ${\tilde{X}}_j=FW-X_{j-1}{B}_{j-1}$              $F=M{\tilde{X}}_j$              Do Cholesky decomposition ${\tilde{X}}_j^{T}F=W^{T}W$              $A_j=W$, $W=inv(W)$, $X_j={\tilde{X}}_{j}W$         %$W=inv(W)$ means $W=W^{-1}$              ${\tilde{Y}}_{j+1}=FW-Y_j{A}_j^T$              $F=K{\tilde{Y}}_{j+1}$              Do Cholesky decomposition ${\tilde{Y}}_{j+1}^{T}F=W^{T}W$              $B_j=W^T$, $W=inv(W)$, $Y_{j+1}={\tilde{Y}}_{j+1}W$          End

Remark 2.

With j increasing inStep 2, the M-orthogonality of ${\mathcal{X}}_j$ and the K-orthogonality of ${\mathcal{Y}}_j$ will slowly lose. Thus, in practice, we can add a re-orthogonalization process in each iteration to eliminate the defect. The same strategy is executed in the following algorithms.

From (1), we have

$$\left[\begin{array}{cc}0& M\\ K& 0\end{array}\right]\left[\begin{array}{cc}{\mathcal{Y}}_n& 0\\ 0& {\mathcal{X}}_n\end{array}\right]=\left[\begin{array}{cc}{\mathcal{Y}}_n& 0\\ 0& {\mathcal{X}}_n\end{array}\right]\left[\begin{array}{cc}0& {\mathcal{B}}_n^T\\ {\mathcal{B}}_n& 0\end{array}\right]+\left[\begin{array}{c}{Y}_{n+1}\\ 0\end{array}\right]B_n^T{E}_{2n}^T$$

with $E_{2n}^T=[0_{n_b\times (2n-1)n_b},I_{n_b}]$. Then, the approximate eigenpairs of $\mathbf{H}$ can be obtained by solving a small eigenvalue problem of $\left[\begin{array}{cc}0& {\mathcal{B}}_n^T\\ {\mathcal{B}}_n& 0\end{array}\right]$. Suppose ${\mathcal{B}}_n$ has an singular value decomposition

$${\mathcal{B}}_n=\mathsf{\Phi }{\Sigma }_n{\mathsf{\Psi }}^T,$$

(2)

where $\mathsf{\Phi }=[{\varphi }_1,{\varphi }_2,\cdots ,{\varphi }_{n{n}_b}]$, $\mathsf{\Psi }=[{\psi }_1,{\psi }_2,\cdots ,{\psi }_{n{n}_b}]$, ${\Sigma }_n=[{\sigma }_1,{\sigma }_2,\cdots ,{\sigma }_{n{n}_b}]$ with ${\sigma }_1\ge \cdots \ge {\sigma }_{n{n}_b}>0$. Thus, we can take $\pm {\sigma }_j(1\le j\le n{n}_b)$ as the Ritz values of $\mathbf{H}$ and

$${\tilde{z}}_j=\frac{1}{\sqrt{2}}\left[\begin{array}{c}{\mathcal{Y}}_n{\psi }_j\\ \pm {\mathcal{X}}_n{\varphi }_j\end{array}\right],1\le j\le n{n}_b,$$

as the corresponding $\mathbf{K}$-orthonormal Ritz vectors, where $\mathbf{K}=\left[\begin{array}{cc}K& 0\\ 0& M\end{array}\right]$.

3.2. Convergence Analysis

In this section, we first consider the convergence analysis when using the first few ${\sigma }_j$ as approximations to the first few ${\lambda }_j$. Then, the similar theories are presented if using the last few ${\sigma }_j$ as approximations to the last few ${\lambda }_j$. Since a block Lanczos method with a suitable block size which is not smaller than the size of an eigenvalue cluster can compute all eigenvalues in the cluster. Now, we are considering the i-th to $(i+n_b-1)$-st eigenpairs of LREP, in which the k-th to ℓ-th eigenvalues form a cluster as in the following figure with $1\le i\le k\le \ell \le i+n_b-1\le n{n}_b$ and $k\le n$.

Here, the squares of the eigenvalues for LREP are listed. Hence, motivated by [12,17], we analyze the convergence of the cluster eigenvalues and their corresponding eigenspace, and give the error bounds of the approximate eigenpairs belonging to eigenvalue cluster together, instead of separately for each individual eigenpair.

We first give some notations and equations, which are critical in our main theorem. Note that from (1), we get

$$MK{\mathcal{Y}}_n={\mathcal{Y}}_n{\mathcal{B}}_n^T{\mathcal{B}}_n+Y_{n+1}{B}_n^T{A}_n{E}_n^T.$$

(3)

Since (2) is the singular value decomposition of ${\mathcal{B}}_n$, thus the eigenvalues of ${\mathcal{B}}_n^T{\mathcal{B}}_n$ are ${\sigma }_j^2$ with the associated eigenvectors ${\psi }_j$ for $1\le j\le n{n}_b$.

From Lemma 2, if we let $\mathsf{\Xi }=[{\xi }_1,\cdots ,{\xi }_N]$, and $\mathsf{\Gamma }=[{\eta }_1,\cdots ,{\eta }_N]$, then $\mathsf{\Gamma }=K\mathsf{\Xi }{\mathsf{\Lambda }}^{-1}$, and

$$MK\mathsf{\Xi }=\mathsf{\Xi }{\mathsf{\Lambda }}^2.$$

(4)

Write $\mathsf{\Xi }$ and ${\mathsf{\Lambda }}^2$ as

Let ${\check{\mathsf{\Xi }}}_2={\mathsf{\Xi }}_{(:,k:\ell )}$ and ${\check{\mathsf{\Lambda }}}_2^2=diag({\lambda }_k^2,\cdots ,{\lambda }_{\ell }^2)$. Denote $C_j$ the first kind Chebyshev polynomial with j-th degree, and $0\le j\le n$.

In the following, we assume

$${\theta }_K^{(1)}(Y_1,{\mathsf{\Xi }}_2)<\pi /2,$$

(5)

i.e., $rank(Y_1^{T}K{\mathsf{\Xi }}_2)=n_b$, then from Lemma 1, we have ∃$Z\in {\mathbb{R}}^{N\times (\ell -k+1)}$ with $\mathcal{R}(Z)\subseteq \mathcal{R}(Y_1)$, s.t.,

$${\mathsf{\Xi }}_2{\mathsf{\Xi }}_2^{T}KZ={\check{\mathsf{\Xi }}}_2.$$

(6)

Theorem 1.

Suppose ${\theta }_K^{(1)}(Y_1,{\mathsf{\Xi }}_2)<\pi /2$, and Z satisfy (6), then we have

$$\parallel diag({\lambda }_k^2-{\sigma }_k^2,\cdots ,{\lambda }_{\ell }^2-{\sigma }_{\ell }^2){\parallel }_F\le ({\lambda }_k^2-{\lambda }_N^2)\frac{{\pi }_{i,k,\ell }^2}{C_{n-k}^2(1+2{\gamma }_{i,\ell })}{\parallel {tan}^2{\mathsf{\Theta }}_K({\check{\mathsf{\Xi }}}_2,Z)\parallel }_F$$

(7)

with

$${\gamma }_{i,\ell }=\frac{{\lambda }_{\ell }^2-{\lambda }_{i+n_b}^2}{{\lambda }_{i+n_b}^2-{\lambda }_N^2},{\pi }_{i,k,\ell }=\frac{\underset{i+n_b\le j\le N}{max}{\displaystyle \prod _{m=1}^{k-1}}|{\sigma }_m^2-{\lambda }_j^2|}{\underset{k\le t\le \ell }{min}{\displaystyle \prod _{m=1}^{k-1}}|{\sigma }_m^2-{\lambda }_t^2|},$$

and

$$\parallel sin{\mathsf{\Theta }}_K({\check{\mathsf{\Xi }}}_2,{\mathcal{Y}}_n{\mathsf{\Psi }}_{(:,k:\ell )}){\parallel }_F\le \frac{{\pi }_{i,k}\sqrt{1+c^2{\parallel A_n^T{B}_n\parallel }_2^2/{\delta }^2}}{C_{n-i}(1+2{\gamma }_{i,\ell })}{\parallel tan{\mathsf{\Theta }}_K({\check{\mathsf{\Xi }}}_2,Z)\parallel }_F$$

(8)

with constant c lies between 1 and $\pi /2$, and $c=1$ if $k=\ell$, and

$$\delta =\underset{\begin{array}{c}k\le j\le \ell \\ p<korp>\ell \end{array}}{min}|{\lambda }_j^2-{\sigma }_p^2|,{\pi }_{i,k}={\displaystyle \prod _{j=1}^{i-1}}\frac{{\lambda }_j^2-{\lambda }_N^2}{{\lambda }_j^2-{\lambda }_k^2}.$$

Particularly if ${\sigma }_{k-1}^2\ge {\lambda }_k^2$, then

$${\pi }_{i,k,\ell }={\displaystyle \prod _{m=1}^{k-1}}\frac{|{\sigma }_m^2-{\lambda }_N^2|}{|{\sigma }_m^2-{\lambda }_k^2|}.$$

Proof. 

Multiplying $L^T$ from left, (4) can be rewritten as $L^{T}ML(L^T\mathsf{\Xi })=(L^T\mathsf{\Xi }){\mathsf{\Lambda }}^2$, so, $({\lambda }_j^2,L^T{\xi }_j)$ is the eigenpair of $L^{T}ML$, for $j=1,\cdots ,N$, and $L^T{\xi }_1,\cdots ,L^T{\xi }_N$ are orthonormal. Do the same process to (3), we have

$$L^{T}ML{\mathcal{V}}_n={\mathcal{V}}_n{\mathcal{B}}_n^T{\mathcal{B}}_n+V_{n+1}{B}_n^T{A}_n{E}_n^T,$$

(9)

where ${\mathcal{V}}_n=L^T{\mathcal{Y}}_n$, $V_{n+1}=L^T{Y}_{n+1}$, and ${\mathcal{V}}_n^T{\mathcal{V}}_n=I_{n{n}_b}$, which can be seen as the relation generalize by using standard Lanczos process to $L^{T}ML$. Thus, ${\sigma }_1^2,\cdots ,{\sigma }_{n{n}_b}^2$ are the Ritz values of $L^{T}ML$, with the corresponding orthonormal Ritz vectors ${\mathcal{V}}_n{\psi }_1,\cdots ,{\mathcal{V}}_n{\psi }_{n{n}_b}$.

Premultiplying $L^T$ to Equation (6), we have $L^T{\mathsf{\Xi }}_2{\mathsf{\Xi }}_2^{T}L(L^{T}Z)=L^T{\check{\mathsf{\Xi }}}_2$. Consequently, the conditions of the block Lanczos convergence Theorem 4.1 and Theorem 5.1 in [17] are satisfied. Thus, using the results Theorem 5.1 in [17], one has

$$\parallel diag({\lambda }_k^2-{\sigma }_k^2,\cdots ,{\lambda }_{\ell }^2-{\sigma }_{\ell }^2){\parallel }_F\le ({\lambda }_k^2-{\lambda }_N^2)\frac{{\pi }_{i,k,\ell }^2}{C_{n-k}^2(1+2{\gamma }_{i,\ell })}{\parallel {tan}^2\mathsf{\Theta }(L^T{\check{\mathsf{\Xi }}}_2,L^{T}Z)\parallel }_F.$$

Then the bound (7) can be easily got by using ([21] Theorem 4.2)

$$\mathsf{\Theta }(L^T{\check{\mathsf{\Xi }}}_2,L^{T}Z)={\mathsf{\Theta }}_K({\check{\mathsf{\Xi }}}_2,Z).$$

(10)

Let ${\Pi }_n={\mathcal{V}}_n{\mathcal{V}}_n^T$, then ${\Pi }_n$ is the orthogonal projection onto ${\mathcal{K}}_n(L^{T}ML,L^{T}Z)$, thus from (9), we have

$$\begin{array}{cc}\hfill \parallel {\Pi }_n{L}^{T}ML(I-{\Pi }_n){\parallel }_2& =\parallel {\mathcal{V}}_n{\mathcal{V}}_n^T{L}^{T}ML(I-{\mathcal{V}}_n{\mathcal{V}}_n^T){\parallel }_2\hfill \\ & =\parallel {\mathcal{V}}_n({\mathcal{B}}_n^T{\mathcal{B}}_n+E_n{A}_n^T{B}_n{V}_{n+1}^T)-{\mathcal{V}}_n{\mathcal{B}}_n^T{\mathcal{B}}_n{\mathcal{V}}_n^T{\parallel }_2\hfill \\ & =\parallel V_n{A}_n^T{B}_n{V}_{n+1}^T{\parallel }_2\hfill \\ & =\parallel A_n^T{B}_n{\parallel }_2.\hfill \end{array}$$

Consequently, applying the results of Theorem 4.1 in [17], we get

$$\begin{array}{cc}\hfill \parallel sin\mathsf{\Theta }(L^T{\check{\mathsf{\Xi }}}_2,{\mathcal{V}}_n{\mathsf{\Psi }}_{(:,k:\ell )}){\parallel }_F& \le \frac{{\pi }_{i,k}\sqrt{1+\parallel {\Pi }_n{L}^{T}ML(I-{\Pi }_n){\parallel }_2^2/{\delta }^2}}{C_{n-i}(1+2{\gamma }_{i,\ell })}{\parallel tan\mathsf{\Theta }(L^T{\check{\mathsf{\Xi }}}_2,L^{T}Z)\parallel }_F\hfill \\ & =\frac{{\pi }_{i,k}\sqrt{1+\parallel A_n^T{B}_n{\parallel }_2^2/{\delta }^2}}{C_{n-i}(1+2{\gamma }_{i,\ell })}{\parallel tan\mathsf{\Theta }(L^T{\check{\mathsf{\Xi }}}_2,L^{T}Z)\parallel }_F.\hfill \end{array}$$

Then the bound (8) can be derived by using $\mathsf{\Theta }(L^T{\check{\mathsf{\Xi }}}_2,{\mathcal{V}}_n{\mathsf{\Psi }}_{(:,k:\ell )})={\mathsf{\Theta }}_K({\check{\mathsf{\Xi }}}_2,{\mathcal{Y}}_n{\mathsf{\Psi }}_{(:,k:\ell )})$ and (10). ☐

Theorem 1 is used to bound the errors of the approximate eigenvalues to an eigenvalue cluster including the multiple eigenvalues. It can be also applied to the single eigenvalue case, the following corollary is derived by setting $k=\ell =i$, except the left equality of (10), which needs to be proved.

Corollary 1.

Suppose ${\theta }_K^{(1)}(Y_1,{\mathsf{\Xi }}_2)<\pi /2$, then for $1\le i\le n{n}_b$, there exits a vector $y\in \mathcal{R}(Y_1)$, s.t., ${\mathsf{\Xi }}_2{\mathsf{\Xi }}_2^{T}y={\xi }_i$, and

$${\lambda }_i^2-{\sigma }_i^2\le ({\lambda }_i^2-{\lambda }_N^2)\frac{{\pi }_{i,j}^2}{C_{n-i}^2(1+2{\gamma }_i)}{tan}^2{\theta }_K({\xi }_i,y)$$

with

$${\gamma }_i=\frac{{\lambda }_i^2-{\lambda }_{i+n_b}^2}{{\lambda }_{i+n_b}^2-{\lambda }_N^2},{\pi }_{i,j}=\underset{i+n_b\le j\le N}{max}{\displaystyle \prod _{m=1}^{i-1}}\frac{|{\sigma }_m^2-{\lambda }_j^2|}{|{\sigma }_m^2-{\lambda }_i^2|},$$

and

$$\begin{array}{cc}& {\left((1-\frac{{\sigma }_i^2}{{\lambda }_i^2})+\frac{{\sigma }_i^2}{{\lambda }_i^2}{sin}^2{\theta }_M({\eta }_i,{\mathcal{X}}_n{\varphi }_i)\right)}^{1/2}\hfill \\ \hfill =& sin{\theta }_K({\xi }_i,{\mathcal{Y}}_n{\psi }_i)\le \frac{{\pi }_i\sqrt{1+\parallel A_n^T{B}_n{\parallel }_2^2/{\delta }^2}}{C_{n-i}(1+2{\gamma }_i)}tan{\theta }_K({\xi }_i,y)\hfill \end{array}$$

(11)

with

$$\delta =\underset{i\ne j}{min}|{\lambda }_j^2-{\sigma }_i^2|,{\pi }_i={\displaystyle \prod _{j=1}^{i-1}}\frac{{\lambda }_j^2-{\lambda }_N^2}{{\lambda }_j^2-{\lambda }_i^2}.$$

Proof. 

We only proof the left equality of (11). From (4) and Lemma 2, we have $\mathsf{\Xi }=MK\mathsf{\Xi }{\mathsf{\Lambda }}^{-2}=M\mathsf{\Gamma }{\mathsf{\Lambda }}^{-1}$. If we let $Z_1={({\mathcal{Y}}_n{\psi }_i)}^{T}K{\xi }_i$, and $Z_2={({\mathcal{X}}_n{\varphi }_i)}^{T}M{\eta }_i$, then we can get $Z_1=\frac{{\sigma }_i}{{\lambda }_i}{Z}_2$ by using $K{\mathcal{Y}}_n\mathsf{\Psi }={\mathcal{X}}_n{B}_n\mathsf{\Psi }={\mathcal{X}}_n\mathsf{\Phi }{\Sigma }_n$. Thus

$$\begin{array}{cc}\hfill {sin}^2{\theta }_K({\xi }_i,{\mathcal{Y}}_n{\psi }_i)& =1-{cos}^2{\theta }_K({\xi }_i,{\mathcal{Y}}_n{\psi }_i)\hfill \\ & =1-Z_1^T{Z}_1\hfill \\ & =1-\frac{{\sigma }_i^2}{{\lambda }_i^2}{Z}_2^T{Z}_2\hfill \\ & =1-\frac{{\sigma }_i^2}{{\lambda }_i^2}{cos}^2{\theta }_M({\eta }_i,{\mathcal{X}}_n{\varphi }_i)\hfill \\ & =1-\frac{{\sigma }_i^2}{{\lambda }_i^2}+\frac{{\sigma }_i^2}{{\lambda }_i^2}{sin}^2{\theta }_M({\eta }_i,{\mathcal{X}}_n{\varphi }_i).\hfill \end{array}$$

Then,

$$sin{\theta }_K({\xi }_i,{\mathcal{Y}}_n{\psi }_i)={\left(1-\frac{{\sigma }_i^2}{{\lambda }_i^2}+\frac{{\sigma }_i^2}{{\lambda }_i^2}{sin}^2{\theta }_M({\eta }_i,{\mathcal{X}}_n{\varphi }_i)\right)}^{1/2}.$$

 ☐

Next, we are going to consider the last few ${\sigma }_j$ to approximate as the last few ${\lambda }_{N-n{n}_b+j}$, $j=k,\cdots ,\ell$, and ${\lambda }_{N-n{n}_b+k}$ to ${\lambda }_{N-n{n}_b+\ell }$ form a cluster in ${\lambda }_{\widehat{i}}$ to ${\lambda }_{\widehat{i}+n_b-1}$, which is described in the following figure, where $N+1-n{n}_b\le \widehat{i}\le \widehat{k}\le \widehat{\ell }\le \widehat{i}+n_b-1\le N$, $n{n}_b-\ell +1\le n$, $\widehat{k}\triangleq N-n{n}_b+k$, and $\widehat{\ell }\triangleq N-n{n}_b+\ell$.

Similar to the above discussion for the first few eigenvalues, we can also obtain the error bounds of the approximate last few eigenpairs belongs to eigenvalue cluster together. We use the same notion, except ${\widehat{\mathsf{\Lambda }}}_2^2=diag({\lambda }_{\widehat{k}}^2,\cdots ,{\lambda }_{\widehat{\ell }}^2)$ and ${\widehat{\mathsf{\Xi }}}_2={\mathsf{\Xi }}_{(:,\widehat{k}:\widehat{\ell })}$. Assuming ${\theta }_K^{(1)}(Y_1,{\mathsf{\Xi }}_2)<\pi /2$, then from Lemma 1, there ∃$\widehat{Z}\in {\mathbb{R}}^{N\times (\ell -k+1)}$ with $\mathcal{R}(\widehat{Z})\subseteq \mathcal{R}(Y_1)$, s.t.,

$${\mathsf{\Xi }}_2{\mathsf{\Xi }}_2^{T}K\widehat{Z}={\widehat{\mathsf{\Xi }}}_2.$$

(12)

Theorem 2.

Suppose ${\theta }_K^{(1)}(Y_1,{\mathsf{\Xi }}_2)<\pi /2$ and $\widehat{Z}$ satisfy (12), then we have

$$\begin{array}{cc}& \parallel diag({\sigma }_k^2-{\lambda }_{\widehat{k}}^2,\cdots ,{\sigma }_{\ell }^2-{\lambda }_{\widehat{\ell }}^2){\parallel }_F\hfill \\ \hfill \le & ({\lambda }_1^2-{\lambda }_{\widehat{\ell }}^2)\frac{{\widehat{\pi }}_{\widehat{i},\widehat{k},\widehat{\ell }}^2}{C_{n-N+\widehat{\ell }-1}^2(1+2{\widehat{\gamma }}_{\widehat{i},\widehat{k}})}{\parallel {tan}^2{\mathsf{\Theta }}_K({\widehat{\mathsf{\Xi }}}_2,\widehat{Z})\parallel }_F\hfill \end{array}$$

(13)

with

$${\widehat{\gamma }}_{\widehat{i},\widehat{\ell }}=\frac{{\lambda }_{\widehat{i}-1}^2-{\lambda }_{\widehat{k}}^2}{{\lambda }_1^2-{\lambda }_{\widehat{i}-1}^2},{\widehat{\pi }}_{\widehat{i},\widehat{k},\widehat{\ell }}=\frac{\underset{1\le j\le \widehat{i}-1}{max}{\displaystyle \prod _{m=\ell +1}^{n{n}_b}}|{\sigma }_m^2-{\lambda }_j^2|}{\underset{\widehat{k}\le t\le \widehat{\ell }}{min}{\displaystyle \prod _{m=\ell +1}^{n{n}_b}}|{\sigma }_m^2-{\lambda }_t^2|},$$

and

$$\parallel sin{\mathsf{\Theta }}_K({\widehat{\mathsf{\Xi }}}_2,{\mathcal{Y}}_n{\mathsf{\Psi }}_{(:,k:\ell )}){\parallel }_F\le \frac{{\widehat{\pi }}_{\widehat{i},\widehat{\ell }}\sqrt{1+{\widehat{c}}^2{\parallel A_n^T{B}_n\parallel }_2^2/{\widehat{\delta }}^2}}{C_{n+\widehat{i}+n_b-N-2}(1+2{\widehat{\gamma }}_{\widehat{i},\widehat{k}})}{\parallel tan{\mathsf{\Theta }}_K({\widehat{\mathsf{\Xi }}}_2,\widehat{Z})\parallel }_F$$

(14)

with constant $\widehat{c}$ lies between 1 and $\pi /2$, and $\widehat{c}=1$ if $k=\ell$, and

$$\widehat{\delta }=\underset{\begin{array}{c}\widehat{k}\le j\le \widehat{\ell }\\ p<korp>\ell \end{array}}{min}|{\lambda }_j^2-{\sigma }_p^2|,{\widehat{\pi }}_{\widehat{i},\widehat{\ell }}={\displaystyle \prod _{j=\widehat{i}+n_b}^N}\frac{{\lambda }_1^2-{\lambda }_j^2}{{\lambda }_{\widehat{\ell }}^2-{\lambda }_j^2}.$$

Remark 3.

Similar to Corollary 1, Theorem 2 can also be applied to the single eigenvalue case, here we omit the detail.

Remark 4.

In Theorem 1 and 2, we use the Frobenius norm to estimate the accuracy of eigenpairs approximations, in fact, any unitary invariant norm can be used to measure.

Remark 5.

Compared with the single-vector type of the weighted Golub-Kahan-Lanczos method in [16], our convergence results show the superiority of the block version. For instance, in Corollary 1, the convergence rate of the approximate eigenvalues ${\sigma }_j$ is proportional to $C_{n-i}^{-2}(1+2{\gamma }_i)$ with ${\gamma }_i=\frac{{\lambda }_i^2-{\lambda }_{i+n_b}^2}{{\lambda }_{i+n_b}^2-{\lambda }_N^2}$, which is obviously better than $C_{n-i}^{-2}(1+2{\tilde{\gamma }}_i)$ with ${\tilde{\gamma }}_i=\frac{{\lambda }_i^2-{\lambda }_{i+1}^2}{{\lambda }_{i+1}^2-{\lambda }_N^2}$ in ([16] Theorem 3.4). While the additional cost caused from the block version can be paid by the improvements generated by ${\gamma }_i$, especially when the desired eigenvalues lie in a well-separated cluster [12].

4. Thick Restart

As the number of iterations increases, Algorithm 1 may encounter the dilemma that the amount of calculation and storage increases sharply and the numerical stability gradually weakens. In this section, we will apply the thick restart strategy [20] to improve the algorithm. After running n iterations, Algorithm 1 derives the following relations for LREP:

$$\left\{\begin{array}{cc}\hfill K{\mathcal{Y}}_n& ={\mathcal{X}}_n{\mathcal{B}}_n,\hfill \\ \hfill M{\mathcal{X}}_n& ={\mathcal{Y}}_n{\mathcal{B}}_n^T+Y_{n+1}{B}_n^T{E}_n^T,\hfill \end{array}\right.$$

(15)

with ${\mathcal{X}}_n^{T}M{\mathcal{X}}_n=I_{n{n}_b}={\mathcal{Y}}_n^{T}K{\mathcal{Y}}_n$.

Recall the SVD (2), let ${\mathsf{\Phi }}_k$ and ${\mathsf{\Psi }}_k$ be the first $k{n}_b$ columns of $\mathsf{\Phi }$ and $\mathsf{\Psi }$, respectively, i.e.,

$${\mathsf{\Phi }}_k=[{\varphi }_1,{\varphi }_2,\cdots ,{\varphi }_{k{n}_b}],{\mathsf{\Psi }}_k=[{\psi }_1,{\psi }_2,\cdots ,{\psi }_{k{n}_b}].$$

Thus it follows that

$${\mathcal{B}}_n{\mathsf{\Psi }}_k={\mathsf{\Phi }}_k{\Sigma }_k\mathrm{and}{\mathcal{B}}_n^T{\mathsf{\Phi }}_k={\mathsf{\Psi }}_k{\Sigma }_k,$$

(16)

where ${\Sigma }_k=diag({\sigma }_1,\cdots ,{\sigma }_{k{n}_b})$.

By using the approximate eigenvectors of $\mathbf{H}$ for thick restart, we post-multiply ${\mathsf{\Psi }}_k$ and ${\mathsf{\Phi }}_k$ to the Equation (15), respectively, and get

$$\left\{\begin{array}{cc}\hfill K{\mathcal{Y}}_n{\mathsf{\Psi }}_k& ={\mathcal{X}}_n{\mathcal{B}}_n{\mathsf{\Psi }}_k,\hfill \\ \hfill M{\mathcal{X}}_n{\mathsf{\Phi }}_k& ={\mathcal{Y}}_n{\mathcal{B}}_n^T{\mathsf{\Phi }}_k+Y_{n+1}{B}_n^T{E}_n^T{\mathsf{\Phi }}_k,\hfill \end{array}\right.$$

(17)

From (16), and let ${\widehat{\mathcal{Y}}}_k={\mathcal{Y}}_n{\mathsf{\Psi }}_k$, ${\widehat{\mathcal{X}}}_k={\mathcal{X}}_n{\mathsf{\Phi }}_k$, ${\widehat{\mathcal{B}}}_k={\Sigma }_k$, ${\widehat{Y}}_{k+1}=Y_{n+1}$, $U^T=E_n^T{\mathsf{\Phi }}_k$, ${\widehat{B}}_k=B_n$, then (17) can be rewritten as

$$\left\{\begin{array}{cc}\hfill K{\widehat{\mathcal{Y}}}_k& ={\widehat{\mathcal{X}}}_k{\widehat{\mathcal{B}}}_k,\hfill \\ \hfill M{\widehat{\mathcal{X}}}_k& ={\widehat{\mathcal{Y}}}_k{\widehat{\mathcal{B}}}_k^T+{\widehat{Y}}_{k+1}{\widehat{B}}_k^T{U}^T,\hfill \end{array}\right.$$

(18)

and ${\widehat{\mathcal{X}}}_k^{T}M{\widehat{\mathcal{X}}}_k=I_{k{n}_b}={\widehat{\mathcal{Y}}}_k^{T}K{\widehat{\mathcal{Y}}}_k$.

Next, ${\widehat{X}}_{k+1}$ and ${\widehat{Y}}_{k+2}$ will be generalized. Firstly, we compute

$$\begin{array}{cc}\hfill {\tilde{X}}_{k+1}& =K{\widehat{Y}}_{k+1}-{\widehat{\mathcal{X}}}_k{\widehat{\mathcal{X}}}_k^{T}MK{\widehat{Y}}_{k+1}\hfill \\ & =K{\widehat{Y}}_{k+1}-{\widehat{\mathcal{X}}}_{k}U{\widehat{B}}_k.\hfill \end{array}$$

From the second equation in (18), we know ${\tilde{X}}_{k+1}^{T}M{\widehat{\mathcal{X}}}_k=0$. Do Cholesky decomposition ${\tilde{X}}_{k+1}^{T}M{\tilde{X}}_{k+1}=W^{T}W$, and set ${\widehat{A}}_{k+1}=W$, $W=inv(W)$. Compute ${\widehat{X}}_{k+1}={\tilde{X}}_{k+1}W$, and let

$${\widehat{\mathcal{X}}}_{k+1}=[{\widehat{\mathcal{X}}}_k,{\widehat{X}}_{k+1}],{\widehat{\mathcal{B}}}_{k+1}=\left[\begin{array}{cc}{\widehat{\mathcal{B}}}_k& U{\widehat{B}}_k\\ 0& {\widehat{A}}_{k+1}\end{array}\right],$$

we have

$$K{\widehat{\mathcal{Y}}}_{k+1}={\widehat{\mathcal{X}}}_{k+1}{\widehat{\mathcal{B}}}_{k+1}\mathrm{with}{\widehat{\mathcal{X}}}_{k+1}^{T}M{\widehat{\mathcal{X}}}_{k+1}=I_{(k+1)n_b}.$$

(19)

Secondly, from the above equation, we can compute

$$\begin{array}{cc}\hfill {\tilde{Y}}_{k+2}& =M{\widehat{X}}_{k+1}-{\widehat{\mathcal{Y}}}_k{\widehat{\mathcal{Y}}}_k^{T}KM{\widehat{X}}_{k+1}-{\widehat{Y}}_{k+1}{\widehat{Y}}_{k+1}^{T}KM{\widehat{X}}_{k+1}\hfill \\ & =M{\widehat{X}}_{k+1}-{\widehat{Y}}_{k+1}{\widehat{A}}_{k+1}^T.\hfill \end{array}$$

Again using (19), it is easily got that ${\tilde{Y}}_{k+2}^{T}K{\widehat{\mathcal{Y}}}_{k+1}=0$. Similarly, do Cholesky decomposition ${\tilde{Y}}_{k+2}^{T}K{\tilde{Y}}_{k+2}=W^{T}W$, and let ${\widehat{B}}_{k+1}=W^T$, $W=inv(W)$. Compute ${\widehat{Y}}_{k+2}={\tilde{Y}}_{k+2}W$, and let ${\widehat{\mathcal{Y}}}_{k+1}=[{\widehat{\mathcal{Y}}}_k,{\widehat{Y}}_{k+1}]$, we get

$$M{\widehat{\mathcal{X}}}_{k+1}={\widehat{\mathcal{Y}}}_{k+1}{\widehat{\mathcal{B}}}_{k+1}^T+{\widehat{Y}}_{k+2}{\widehat{B}}_{k+1}^T{E}_{k+1}^T\mathrm{with}{\widehat{\mathcal{Y}}}_{k+1}^{T}M{\widehat{\mathcal{Y}}}_{k+1}=I_{(k+1)n_b}.$$

Continue the same procedure for ${\widehat{X}}_{k+2},\cdots ,{\widehat{X}}_n$ and ${\widehat{Y}}_{k+3},\cdots ,{\widehat{Y}}_{n+1}$, we can obtain the new M-orthonormal matrix ${\widehat{\mathcal{X}}}_n\in {\mathbb{R}}^{N\times n{n}_b}$, the new K-orthonormal matrix ${\widehat{\mathcal{Y}}}_n\in {\mathbb{R}}^{N\times n{n}_b}$, and the new matrix ${\widehat{\mathcal{B}}}_n\in {\mathbb{R}}^{n{n}_b\times n{n}_b}$, and relations

$$\left\{\begin{array}{cc}\hfill K{\widehat{\mathcal{Y}}}_n& ={\widehat{\mathcal{X}}}_n{\widehat{\mathcal{B}}}_n,\hfill \\ \hfill M{\widehat{\mathcal{X}}}_n& ={\widehat{\mathcal{Y}}}_n{\widehat{\mathcal{B}}}_n^T+{\widehat{Y}}_{n+1}{\widehat{B}}_n^T{E}_n^T,\hfill \end{array}\right.$$

(20)

with ${\widehat{\mathcal{X}}}_n^{T}M{\widehat{\mathcal{X}}}_n=I_{n{n}_b}={\widehat{\mathcal{Y}}}_n^{T}K{\widehat{\mathcal{Y}}}_n$, and

$${\widehat{\mathcal{B}}}_n=\left[\begin{array}{cccc}{\widehat{\mathcal{B}}}_k& U{\widehat{B}}_k& & \\ & {\widehat{A}}_{k+1}& {\widehat{B}}_{k+1}& \\ & & \ddots & {\widehat{B}}_{n-1}\\ & & & {\widehat{A}}_n\end{array}\right].$$

Note that ${\widehat{\mathcal{B}}}_n$ is no longer a block bidiagonal matrix. Algorithm 2 is our thick-restart weighted block Golub-Kahan-Lanczos algorithm for LREP.

Remark 6.

Actually, from the construction of ${\widehat{\mathcal{B}}}_n$, we can know the procedure for getting ${\widehat{X}}_{k+2},\cdots ,{\widehat{X}}_n$ and ${\widehat{Y}}_{k+3},\cdots ,{\widehat{Y}}_{n+1}$ is the same as applying Algorithm 1 to ${\widehat{Y}}_{k+2}$ for $n-k-1$ iterations, thus we use Algorithm 1 directly in restartingStep 2of the following Algorithm 2.

Algorithm 2: wbGKL-TR

1. Given an initial guess $Y_1$ satisfying $Y_1^{T}K{Y}_1=I_{n_b}$, a tolerance $tol$, an integer k that the k blocks approximate eigenvectors we want to add to the solving subspace, an integer n the block dimension of solving subspace, as well as $w_{\ell }$ the desired number of eigenpairs;     2. Apply Algorithm 1 from the current point to generate the rest of ${\mathcal{X}}_n$, ${\mathcal{Y}}_{n+1}$, and ${\mathcal{B}}_n$. If it is the first cycle, the current point is $Y_1$, else $Y_{k+2}$;     3. Compute an SVD of ${\mathcal{B}}_n$ as in (2), select $w_{\ell }(w_{\ell }\le n{n}_b)$ wanted singular values ${\sigma }_j$, and their associated left singular vectors ${\varphi }_j$ and right singular vectors ${\psi }_j$. Form the approximate eigenpairs for $\mathbf{H}$, if the stopping criterion is satisfied, then stop, else continue;     4. Generate new ${\widehat{\mathcal{X}}}_{k+1}$, ${\widehat{\mathcal{Y}}}_{k+2}$ and ${\widehat{\mathcal{B}}}_{k+1}$:     Compute ${\widehat{\mathcal{Y}}}_k={\mathcal{Y}}_n{\mathsf{\Psi }}_k$, ${\widehat{\mathcal{X}}}_k={\mathcal{X}}_n{\mathsf{\Phi }}_k$, ${\widehat{\mathcal{B}}}_k={\Sigma }_k$, ${\widehat{Y}}_{k+1}=Y_{n+1}$, $U^T=E_n^T{\mathsf{\Phi }}_k$, ${\widehat{B}}_k=B_n$;     Compute ${\tilde{X}}_{k+1}=K{\widehat{Y}}_{k+1}-{\widehat{\mathcal{X}}}_{k}U{\widehat{B}}_k$, do Cholesky decomposition ${\tilde{X}}_{k+1}^{T}M{\tilde{X}}_{k+1}=W^{T}W$, set ${\widehat{A}}_{k+1}=W$, $W=inv(W)$, ${\widehat{X}}_{k+1}={\tilde{X}}_{k+1}W$;     Compute ${\tilde{Y}}_{k+2}=M{\widehat{X}}_{k+1}-{\widehat{Y}}_{k+1}{\widehat{A}}_{k+1}^T$, do Cholesky decomposition ${\tilde{Y}}_{k+2}^{T}K{\tilde{Y}}_{k+2}=W^{T}W$, set ${\widehat{B}}_{k+1}=W^T$, $W=inv(W)$, ${\widehat{Y}}_{k+2}={\tilde{Y}}_{k+2}W$;     Let ${\mathcal{X}}_{k+1}={\widehat{\mathcal{X}}}_{k+1}=[{\widehat{\mathcal{X}}}_k,{\widehat{X}}_{k+1}]$, ${\mathcal{B}}_{k+1}={\widehat{\mathcal{B}}}_{k+1}=\left[\begin{array}{cc}{\widehat{\mathcal{B}}}_k& U{\widehat{B}}_k\\ 0& {\widehat{A}}_{k+1}\end{array}\right]$, ${\mathcal{Y}}_{k+2}={\widehat{\mathcal{Y}}}_{k+2}=[{\widehat{\mathcal{Y}}}_k,{\widehat{Y}}_{k+1},{\widehat{Y}}_{k+2}]$, and go to Step 2.

Remark 7.

InStep 3, we compute the harmonic Ritz pairs after n iterations. In practice, we do the computation for each iterations $j=1,\cdots ,n$. When restarting, the information chosen to add to the solving subspaces are the wanted $w_{\ell }$ singular values of ${\mathcal{B}}_n$ with their corresponding left and right singular vectors. Actually, we use MATLAB command “sort” to choose the $w_{\ell }$ smallest ones or the $w_{\ell }$ largest ones, and which singular values to choose depends on the desired eigenvalues of $\mathbf{H}$.

In the end of this section, we list the computational costs in a generic cycle of four algorithms, which are weighted block Golub-Kahan-Lanczos algorithm, thick-restart weighted block Golub-Kahan-Lanczos algorithm, block Lanczos algorithm [12], and thick-restart block Lanczos algorithm [12], and denoted by wbGKL, wbGKL-TR, BLan, and BLan-TR, respectively. The detail pseudocodes of BLan and BLan-TR are be found in [12].

The comparisons are presented in

Table 1. Main computational costs per cycle wbGKL and wbGKL-TR.

	wbGKL	wbGKL-TR	wbGKL-TR
		(1-st Cycle)	(Other Cycle)
mvb	$2n+1$	$2n+1$	$2(n-k)$
dpb	$2n+1$	$2n+1$	$2(n-k)$
saxpyb	$8n$	$8n$	$8(n-k)+2k(2n+1)$
block vector updates	$2n+2$	$2n+2$	$2n+2$
Ep($2n\times 2n$)(with sorting)	0	0	0
Sp($n\times n$)	1	1	1

and

Table 2. Main computational costs per cycle BLan and BLan-TR.

	BLan	BLan-TR	BLan-TR
		(1-st Cycle)	(Other Cycle)
mvb	$2n+1$	$2n+1$	$2(n-k)$
dpb	$2n+1$	$2n+1$	$2(n-k)$
saxpyb	$6n$	$6n$	$6(n-k)+2k(2n+1)$
block vector updates	$2n+2$	$2n+2$	$2n+2$
Ep($2n\times 2n$)(with sorting)	1	1	1
Sp($n\times n$)	0	0	0

. Here, we denote “block vector” a $N\times n_b$ rectangular matrix, denote “mvb” the product number of a $N\times N$ matrix and a block vector. “dpb” denotes the dot product number of two block vectors X and Y, i.e., $X^{T}Y$. “saxpyb” denotes the number of adding two block vectors or multiplying a block vector to a $n_b\times n_b$ small matrix. “Ep($2n\times 2n$)(with sorting)” means the number of $2n\times 2n$ size eigenvalue problem with sorting eigenvalues and their corresponding eigenvectors in one cycle. Similarly, “Sp($n\times n$)” denotes the number of $n\times n$ size singular value decomposition in one cycle. Because wbGKL and BLan are non-restart algorithms, we just count the first n Lanczos iterations.

5. Numerical Examples

In this section, two numerical experiments are carried out by using MATLAB 8.4 (R2014b) on a laptop with an Intel Core i5-6200U CPU 2.3 GHz memory 8 GB under the Windows 10 operating system.

Example 1.

In this example, we check the bounds established in Theorem 1 and 2. For simplicity, we take $N=100$, the number of weighted block Golub-Kahan-Lanczos steps $n=20$, $K=M$ as diagonal matrix $diag({\lambda }_1,{\lambda }_2,\cdots ,{\lambda }_N)$, where

$$\begin{array}{c}{\lambda }_1=11+\rho ,{\lambda }_2=11,{\lambda }_3=11-\rho ,\hfill \\ {\lambda }_{N-2}=1+\rho ,{\lambda }_{N-1}=1,{\lambda }_N=1-\rho ,\hfill \\ {\lambda }_j=5+\frac{5(N-j+1)}{N-3},j=4,\cdots ,N-3,\hfill \end{array}$$

and $i=k=1$, $\ell =3$, $\widehat{i}=\widehat{k}=N-2$, $\widehat{\ell }=N$, $n_b=3$. There are three positive eigenvalue clusters: $\{{\lambda }_1,{\lambda }_2,{\lambda }_3\}$, $\{{\lambda }_4,\cdots ,{\lambda }_{N-3}\}$, or $\{{\lambda }_{N-2},{\lambda }_{N-1},{\lambda }_N\}$. Obviously, $\mathsf{\Xi }=\mathsf{\Gamma }=K^{-\frac{1}{2}}$.

We seek two groups of the approximate eigenpairs, the first is related to the first cluster, the second is related to the last cluster, i.e., $\{{\sigma }_1,{\sigma }_2,{\sigma }_3\}$ approximate $\{{\lambda }_1,{\lambda }_2,{\lambda }_3\}$, and $\{{\sigma }_{n{n}_b-2},{\sigma }_{n{n}_b-1},{\sigma }_{n{n}_b}\}$ approximate $\{{\lambda }_{N-2},{\lambda }_{N-1},{\lambda }_N\}$. In order to see the affect that generated from $\rho$ to the upper bounds of the approximate eigenpairs errors in weighted block Golub-Kahan-Lanczos method for LREP, we change the parameter $\rho >0$ to overmaster the tightness among eigenvalues within $\{{\lambda }_1,{\lambda }_2,{\lambda }_3\}$ and $\{{\lambda }_{N-2},{\lambda }_{N-1},{\lambda }_N\}$. First, we choose the same matrix $Y_0$ as in [12,17], i.e.,

$$Y_0=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\\ \frac{1}{N}& sin1& cos1\\ ⋮& ⋮& ⋮\\ \frac{N-n_b}{N}& sin(N-n_b)& cos(N-n_b)\end{array}\right].$$

Obviously, $rank(Y_0)=n_b$ and $rank(Y_0^{T}K{\mathsf{\Xi }}_{(:,1:3)})=n_b$. Since K symmetric positive definite, thus do Cholesky decomposition $Y_0^{T}K{Y}_0=W^{T}W$, let $Y_1=Y_0{W}^{-1}$, hence, $Y_1$ satisfies (5), i.e., $Y_1^{T}K{\mathsf{\Xi }}_{(:,1:3)}$ is singular. We take $Z=Y_1{({\mathsf{\Xi }}_{(:,1:3)}^{T}K{Y}_1)}^{-1}$, then Z satisfies (6). We execute the weighted block Golub-Kahan-Lanczos method with full re-orthogonalization for LREP in MATLAB, and check the bounds in (7), (8), (13), and (14). Since the approximate eigenvalues are $\{{\sigma }_1,{\sigma }_2,{\sigma }_3\}$ and $\{{\sigma }_{n{n}_b-2},{\sigma }_{n{n}_b-1},{\sigma }_{n{n}_b}\}$, thus ${\pi }_{i,k,\ell }={\pi }_{i,k}={\widehat{\pi }}_{\widehat{i},\widehat{k},\widehat{\ell }}={\widehat{\pi }}_{\widehat{i},\widehat{\ell }}=1$, $c=\widehat{c}=1$, and we measure the following two groups of errors:

$$\begin{array}{c}{\epsilon }_{11}={\parallel diag({\lambda }_1^2-{\sigma }_1^2,{\lambda }_2^2-{\sigma }_2^2,{\lambda }_3^2-{\sigma }_3^2)\parallel }_F,\hfill \\ {\epsilon }_{21}=\frac{{\lambda }_1^2-{\lambda }_N^2}{C_{n-1}^2(1+2{\gamma }_{1,3})}{\parallel {tan}^2{\mathsf{\Theta }}_K({\mathsf{\Xi }}_{(:,1:3)},Z)\parallel }_F,\hfill \\ {\epsilon }_{31}={\parallel sin{\mathsf{\Theta }}_K({\mathsf{\Xi }}_{(:,1:3)},{\mathcal{Y}}_n{\mathsf{\Psi }}_{(:,1:3)})\parallel }_F,\hfill \\ {\epsilon }_{41}=\frac{\sqrt{1+\parallel A_n^T{B}_n{\parallel }_2^2/{\delta }^2}}{C_{n-1}(1+2{\gamma }_{1,3})}{\parallel tan{\mathsf{\Theta }}_K({\mathsf{\Xi }}_{(:,1:3)},Z)\parallel }_F,\hfill \end{array}$$

and

$$\begin{array}{c}{\epsilon }_{12}={\parallel diag({\sigma }_{N-2}^2-{\lambda }_{N-2}^2,{\sigma }_{N-1}^2-{\lambda }_{N-1}^2,{\sigma }_N^2-{\lambda }_N^2)\parallel }_F,\hfill \\ {\epsilon }_{22}=\frac{{\lambda }_1^2-{\lambda }_N^2}{C_{n-1}^2(1+2{\widehat{\gamma }}_{N-2,N})}{\parallel {tan}^2{\mathsf{\Theta }}_K({\mathsf{\Xi }}_{(:,N-2:N)},Z)\parallel }_F,\hfill \\ {\epsilon }_{32}={\parallel sin{\mathsf{\Theta }}_K({\mathsf{\Xi }}_{(:,N-2:N)},{\mathcal{Y}}_n{\mathsf{\Psi }}_{(:,n{n}_b-2:n{n}_b)})\parallel }_F,\hfill \\ {\epsilon }_{42}=\frac{\sqrt{1+\parallel A_n^T{B}_n{\parallel }_2^2/{\widehat{\delta }}^2}}{C_{n-i}(1+2{\widehat{\gamma }}_{N-2,N})}{\parallel tan{\mathsf{\Theta }}_K({\mathsf{\Xi }}_{(:,N-2:N)},Z)\parallel }_F.\hfill \end{array}$$

Actually, ${\epsilon }_{21}$ and ${\epsilon }_{41}$ are upper bounds of ${\epsilon }_{11}$ and ${\epsilon }_{31}$, and ${\epsilon }_{22}$ and ${\epsilon }_{42}$ are upper bounds of ${\epsilon }_{12}$ and ${\epsilon }_{32}$.

Table 3. ${\epsilon }_{11}$, ${\epsilon }_{31}$ together with their upper bounds ${\epsilon }_{21}$, ${\epsilon }_{41}$ of Example 1.

$\mathit{\rho }$	${\mathit{\epsilon }}_{11}$	${\mathit{\epsilon }}_{21}$	${\mathit{\epsilon }}_{31}$	${\mathit{\epsilon }}_{41}$
${10}^{-1}$	$4.0295\times {10}^{-13}$	$2.6773\times {10}^{-10}$	$1.2491\times {10}^{-10}$	$2.6260\times {10}^{-6}$
${10}^{-2}$	$5.1238\times {10}^{-14}$	$5.4555\times {10}^{-11}$	$6.1184\times {10}^{-11}$	$1.1407\times {10}^{-6}$
${10}^{-3}$	$7.1054\times {10}^{-14}$	$4.6711\times {10}^{-11}$	$5.7698\times {10}^{-11}$	$1.0520\times {10}^{-6}$
${10}^{-4}$	$2.4449\times {10}^{-13}$	$4.5993\times {10}^{-11}$	$5.7370\times {10}^{-11}$	$1.0436\times {10}^{-6}$
${10}^{-5}$	$2.1552\times {10}^{-13}$	$4.5922\times {10}^{-11}$	$5.7338\times {10}^{-11}$	$1.0427\times {10}^{-6}$

and

Table 4. ${\epsilon }_{12}$, ${\epsilon }_{32}$ together with their upper bounds ${\epsilon }_{22}$, ${\epsilon }_{42}$ of Example 1.

$\mathit{\rho }$	${\mathit{\epsilon }}_{11}$	${\mathit{\epsilon }}_{21}$	${\mathit{\epsilon }}_{31}$	${\mathit{\epsilon }}_{41}$
${10}^{-1}$	$7.1089\times {10}^{-16}$	$6.0352\times {10}^{-11}$	$1.9393\times {10}^{-10}$	$8.8823\times {10}^{-7}$
${10}^{-2}$	$1.3688\times {10}^{-15}$	$3.5913\times {10}^{-11}$	$1.9562\times {10}^{-10}$	$6.8797\times {10}^{-7}$
${10}^{-3}$	$3.9968\times {10}^{-15}$	$3.4113\times {10}^{-11}$	$1.9580\times {10}^{-10}$	$6.7081\times {10}^{-7}$
${10}^{-4}$	$4.8495\times {10}^{-15}$	$3.3938\times {10}^{-11}$	$1.9582\times {10}^{-10}$	$6.6912\times {10}^{-7}$
${10}^{-5}$	$8.1221\times {10}^{-15}$	$3.3920\times {10}^{-11}$	$1.9582\times {10}^{-10}$	$6.6895\times {10}^{-7}$

report the results of ${\epsilon }_{ij}$, $i=1,\cdots ,4$, $j=1,2$ with the parameter $\rho$ goes to 0. From the two tables, we can see that the bounds for the eigenvalues lie in a cluster and their corresponding eigenvectors are sharp, and they are not sensitive to $\rho$ when $\rho$ goes to 0.

Example 2.

In this example, we are going to test the effectiveness of our weighted block Golub-Kahan-Lanczos algorithms. Four algorithms are tested, i.e.,wbGKL,wbGKL-TR,BLan, andBLan-TR. We choose 3 test problems used in [12,13], which are listed in

Table 5. The matrices K and M in Test 1–3.

Problems	N	K	M
Test 1	1862	$Na2$	$Na2$
Test 2	5660	$SiH4$	$SiH4$
Test 3	9604	$fv1$	$finan512$

. All the matrices K and M in the problems are symmetric positive definite. Specifically, Test 1 and Test 2, which are derived by the turboTDDFT command in QUANTUM ESPRESSO [22], are from the linear response research for Na2 and silane (SiH4) compound, respectively. The matrices K and M in Test 3 are from the University of Florida Sparse Matrix Collection [23], where the order of K is $N=9604$, and M is the leading $N\times N$ principal submatrix of $finan512$.

We aim to compute the smallest 5 positive eigenvalues and the largest 5 eigenvalues, i.e., ${\lambda }_i$ for $i=1,\cdots ,5,N-4,\cdots ,N$, together with their associated eigenvectors. The initial guess is chosen as $V_0=eye(N,n_b)$ with block size $n_b=3$, where $eye$ is the MATLAB command. The same as in Example 1, since K is symmetric positive definite, thus do Cholesky decomposition $Y_0^{T}K{Y}_0=W^{T}W$, let $Y_1=Y_0{W}^{-1}$, hence, $Y_1$ satisfies $Y_1^{T}K{Y}_1=I_{n_b}$. In wbGKL-TR and BLan-TR, we select $n=30$, $k=20$, i.e., the restart will occur once the dimension of the solving subspace is larger than 90, and the information of 60 Ritz vectors are kept. For wbGKL and BLan, because there is no restart, then we compute the approximate eigenpairs when the Lanczos iterations equals to $30+10\times (j-1)$, $j=1,2,\cdots$, hence, the Lanczos iterations are as the same amount as in wbGKL-TR and BLan-TR. The following relative eigenvalue error and relative residual 1-norm for each 10 approximate eigenpairs are calculated:

$$\begin{array}{c}e({\sigma }_j):=\left\{\begin{array}{cc}\frac{|{\lambda }_j-{\sigma }_j|}{{\lambda }_j},\hfill & j=1,\cdots ,5,\hfill \\ \frac{|{\lambda }_{n+j-k}-{\sigma }_j|}{{\lambda }_{n+j-k}},\hfill & j=n{n}_b-4,\cdots ,n{n}_b,\hfill \end{array}\right.\hfill \\ r({\sigma }_j):=\frac{\parallel \mathbf{H}{\tilde{z}}_j-{\sigma }_j{\tilde{z}}_j{\parallel }_1}{{(\parallel \mathbf{H}\parallel }_1+{\sigma }_j)\parallel {\tilde{z}}_j{\parallel }_1},j=1,\cdots ,5,n{n}_b-4,\cdots ,n{n}_b,\hfill \end{array}$$

where the “exact” eigenvalues ${\lambda }_j$ are calculated by the MATLAB code $eig$. The calculated approximate eigenpair $({\sigma }_j,{\tilde{z}}_j)$ is regarded as converged if $r({\sigma }_j)\le tol={10}^{-8}$.

Table 6. Compute 5 smallest positive eigenvalues for Test 1–3.

Algorithms	Test 1		Test 2		Test 3
	$\mathit{CPU}$	$\mathit{iter}$	$\mathit{CPU}$	$\mathit{iter}$	$\mathit{CPU}$	$\mathit{iter}$
wbGKL	1.5070	149	25.7848	319	15.9308	379
wbGKL-TR	1.0746	179	20.3593	359	5.1302	589
BLan	4.6739	149	87.1670	349	43.9506	379
BLan-TR	2.1243	163	39.1306	393	19.9677	592

and

Table 7. Compute 5 largest eigenvalues for Test 1–3.

Algorithms	Test 1		Test 2		Test 3
	$\mathit{CPU}$	$\mathit{iter}$	$\mathit{CPU}$	$\mathit{iter}$	$\mathit{CPU}$	$\mathit{iter}$
wbGKL	0.6387	79	12.4658	179	1.0639	109
wbGKL-TR	0.5284	79	9.9093	179	0.8774	109
BLan	1.4634	79	27.4028	179	6.7574	109
BLan-TR	1.0151	82	18.3415	186	4.1298	113

give the number of the Lanczos iterations (denote by $iter$) and the CPU time in seconds (denote by $CPU$) for the four algorithms, and Table 6 is for the smallest 5 positive eigenvalues, Table 7 is for the largest 5 eigenvalues. From Table 6, one can see that, no matter the smallest or the largest eigenvalues, the iteration number of the four algorithms are competitive, but wbGKL and wbGKL-TR cost significant less time than BLan and BLan-TR, especially, wbGKL-TR consumes the least amount of time. Because BLan and BLan-TR need to compute the eigenvalues of $\left[\begin{array}{cc}0& T_n\\ D_n& 0\end{array}\right]$, which is a nonsymmetric matrix, thus the two algorithms slower than wbGKL and wbGKL-TR. Due to the saving during the orthogonalization procedure and solving a much smaller ${\mathcal{B}}_n$, wbGKL-TR is the faster algorithm.

The accuracy of the last two approximate eigenpairs in Test 1 are shown in Figure 1. From the figure, we can see that, for the last two eigenpairs, wbGKL and BLan require almost the same iterations to obtain the same accuracy, and the case of wbGKL-TR and BLan-TR also need almost the same iterations, which are one or two more restarts than wbGKL and BLan. On one hand, without solving a nonsymmetric eigenproblem, wbGKL and wbGKL-TR can save much more time than BLan and BLan-TR. On the other hand, since the dimension of the solving subspace for wbGKL-TR is bounded by $n{n}_b$, the savings in the process of orthogonalization and a much smaller singular value decomposition problem is sufficient to cover the additional restart steps.

6. Conclusions

In this paper, we present a weighted block Golub-Kahan-Lanczos algorithm to solve the desired small portion of smallest or largest positive eigenvalues which are in a cluster. Convergence analysis is established in Theorems 1 and 2, and bound the errors of the eigenvalue and eigenvector approximations belonging to an eigenvalue cluster. These results also show the advantages of the block algorithm over the single-vector version. To make the new algorithm more practical, we introduced a thick-restart strategy to eliminate the numerical difficulties caused by the block method. Numerical examples are executed to demonstrate the efficiency of our new restart algorithm.