Theoretical and Computable Optimal Subspace Expansions for Matrix Eigenvalue Problems^††thanks: Supported by the National Natural Science Foundation of China (Nos. 11771249 and 12171273).

Consider the large scale matrix eigenproblem $Ax=\lambda x$ with $A\in\mathcal{C}^{n\times n}$ and $\|x\|=1$, where $\|\cdot\|$ is the 2-norm of a vector or matrix. We are interested in an exterior eigenvalue $\lambda$ and the associated eigenvector $x$. A number of numerical methods have been available for solving this kind of problem [1, 14, 15, 21, 19]. Of them, Lanczos or Arnoldi type methods [1, 15, 19] are typical, and they are projection methods on the sequence of nested $k$-dimensional Krylov subspaces

and compute approximations to $\lambda$ and $x$. For $A$ Hermitian, Arnoldi type methods reduce to Lanczos type methods [1, 14]. Suppose that $V_{k}=(v_{1},\ldots,v_{k})$ is column orthonormal and is generated by the Lanczos process in the Hermitian case or the Arnoldi process in the non-Hermitian case, and take the expansion vector $w_{k}=v_{k}\in\mathcal{V}_{k}$. Then the next basis vector $v_{k+1}$ is obtained by orthonormalizing $Aw_{k}$ against $V_{k}$, the projection matrix $V_{k}^{H}AV_{k}$ of $A$ onto $\mathcal{V}_{k}$ is Hermitian tridiagonal for $A$ Hermitian or upper Hessenberg for $A$ non-Hermitian, and the columns of $(V_{k},v_{k+1})$ form an orthonormal basis of the $(k+1)$-dimensional Krylov subspace $\mathcal{V}_{k+1}$, where the superscript $H$ denotes the conjugate transpose of a vector or matrix. If the expansion vector $v_{k}$ in $Av_{k}$ is replaced by any $w_{k}\in\mathcal{V}_{k}$ that is not deficient in $v_{k}$, we deduce from Proposition 1 of Ye [25] that the same $\mathcal{V}_{k+1}$ and $V_{k+1}$ are generated.

However, if we start with a non-Krylov subspace $\mathcal{V}_{k}$, then the expanded subspace $\mathcal{V}_{k+1}$ critically depends on the choice of expansion vector $w_{k}\in\mathcal{V}_{k}$. It is well known that the success of a general projection method requires that an underlying subspace contain sufficiently accurate approximations to $x$ but a sufficiently accurate subspace cannot guarantee the convergence of the Ritz and harmonic Ritz vectors obtained by the standard and harmonic Rayleigh–Ritz methods with respect to the subspace [4, 10]. To fix this deficiency, the refined and refined harmonic Rayleigh–Ritz methods have been proposed, which are shown to ensure the convergence of refined and refined harmonic Ritz vectors computed by the refined and refined harmonic Rayleigh–Ritz methods when the subspace is good enough [5, 7, 10, 19, 23].

The accuracy of a projection subspace $\mathcal{V}$ for eigenvector approximation can be measured by its angle with the desired unit length eigenvector $x$ [14, 15, 19]:

where $\angle(z,x)$ denotes the acute angle between $x$ and the nonzero vector $z$. For a general $k$-dimensional non-Krylov subspace $\mathcal{V}$ with the subscript $k$ dropped, in this paper we will consider the following optimal subspace expansion problem: Suppose $\cos\angle(\mathcal{V},x)\not=0$, for any nonzero $w\in\mathcal{V}$, multiply it by $A$, and define the $(k+1)$-dimensional expanded subspace

Then which vector $w_{opt}\in\mathcal{V}$ is an optimal expansion vector that, up to scaling, solves

This problem was first considered by Ye [25] and later paid attention by Wu and Zhang [24]. The choice of $w\in\mathcal{V}$ is essential to subspace expansion, and different $w$ may affect the quality of $\mathcal{V}_{w}$ substantially. At expansion iteration $k>1$, the Lanczos or Arnoldi type expansion takes $w=v_{k}$, the last column of $V_{k}$. The Ritz expansion method [24, 25], a variant of the Arnoldi method, uses the approximating Ritz vector as $w$ to expand the subspace. It is mathematically equivalent to the Residual Arnoldi (RA) method proposed in [11, 12].

The Shift-Invert Residual Arnoldi (SIRA) method is an alternative of the RA method for computing the eigenvalue $\lambda$ closest to a given target $\tau$. More variants have been developed in [8, 9], and they fall into the categories of the harmonic Rayleigh–Ritz, refined Rayleigh–Ritz, and refined harmonic Rayleigh–Ritz methods [1, 14, 15, 19, 21]. The SIRA type methods are mathematically equivalent to the counterparts of the Jacobi–Davidson (JD) type methods; see Theorem 2.1 of [8]. Just as the shift-invert Arnoldi (SIA) type methods applied to $B=(A-\tau I)^{-1}$, the SIRA type and JD type methods [18] both use $B$ to construct nested subspaces but, unlike the SIA type methods, they project the original $A$ rather than $B$ onto the subspaces when finding approximations of $(\lambda,x)$ [8]. At iteration $k$, these methods and the SIA type methods are mathematically common in solving a linear system

i.e., computing $u=Bw$, where $w=v_{k}$ in the SIA type methods and is mathematically equal to the current approximate eigenvector in SIRA and JD type methods; see [8, Theorem 2.1 and its proof] for details.

Ye [25] finds an explicit expression of $w_{opt}$ for a general $A$, as stated below.

Ye’s result [25, Theorem 1]: Let the columns of $V$ form an orthonormal basis of $\mathcal{V}$, define the residual

and assume that the rank $\mathop{\operator@font rank}\nolimits(R)\geq 2$. Then

for any $c\in\mathcal{N}(R)$, the null space of $R$, where the superscript $\dagger$ denotes the Moore-Penrose generalized inverse of a matrix.

Unfortunately, Ye shows that (6) cannot be exploited to obtain a computable (nearly) optimal replacement of $w_{opt}$ as it involves the a priori $x$: (i) $Vc$ makes no contribution to the expansion of $\mathcal{V}$ as $AVc=0$; (ii) since $R^{H}V=0$, we must have $R^{\dagger}z=0$ when replacing $x$ by any $z\in\mathcal{V}$ and taking $c=0$, leading to $w=VR^{\dagger}z=0$. To this end, he gives up (6) and instead turns to deriving a maximization characterization of $\cos\angle(\mathcal{V}_{w},x)$ for a given $w\in\mathcal{V}$ under the assumption that $A$ is Hermitian; see Theorem 2 of [25] and the first result of Theorem 2.1 to be presented in this paper. But he does not find a solution to the maximization characterization problem for a given $w\in\mathcal{V}$. Notice that one cannot obtain $w_{opt}$ from the maximization characterization. Without any other information, e.g., $w_{opt}$, Ye has made an approximate analysis on the maximization characterization of $\cos\angle(\mathcal{V}_{w},x)$ and argued that the Ritz vector from $\mathcal{V}$ might be a good approximate maximizer of (3) and thus might be a nearly optimal expansion vector. Unfortunately, as we shall see, Ye’s analysis and arguments have evident defects; see a detailed explanation after Remark 3 in this paper. We stress that Ye’s proof of the maximization characterization of $\cos\angle(\mathcal{V}_{w},x)$ holds only for $A$ Hermitian.

Wu and Zhang [24] do not consider the optimal expansion problem (3) itself. Instead, they focus on showing that the refined Ritz vector from $\mathcal{V}$ can be a better expansion vector than the Ritz vector from $\mathcal{V}$ for $A$ general. Based on the analysis, they have developed a refined variant of the RA method for a general $A$, which has been numerically confirmed to be more efficient than the RA method.

In this paper, we shall revisit problem (3) for $A$ general. Our theoretical contributions consists of two parts. The first part includes the generalization of Ye’s major result in [25], i.e., Theorem 2, to the non-Hermitian case. We first prove that $\cos\angle(\mathcal{V}_{w},x)$ for a given $w\in\mathcal{V}$ can be formulated as a maximization characterization problem, which extends Theorem 2 of [25] to the general case. Then we find its maximizer. This result is new even for $A$ Hermitian and can be exploited to explain the defect of Ye’s approximate analysis. These results are secondary, and our major theoretical contributions are in the second part. We establish informative expressions on $w_{opt}$, $(I-P_{V})Aw_{opt}$ and the optimally expanded $\mathcal{V}_{w_{opt}}$, where $P_{V}=VV^{H}$ is the orthogonal projector onto $\mathcal{V}$. The results show that (i) generally $w_{opt}\not=P_{V}x$ where $P_{V}x$ is the orthogonal projection of $x$ onto $\mathcal{V}$ and $P_{V}x/\|P_{V}x\|$ is the best approximation to $x$ from $\mathcal{V}$, (ii) $(I-P_{V})Aw_{opt}=RR^{\dagger}x$, which is the orthogonal projection of $x$ onto the subspace ${\rm span}\{R\}$ and whose normalization $RR^{\dagger}x/\|RR^{\dagger}x\|$ is the best approximation to $x$ from ${\rm span}\{R\}$, (iii) $\mathcal{V}_{w_{opt}}=\mathcal{V}\oplus{\rm span}\{RR^{\dagger}x\}$ with $\oplus$ the orthogonal direct sum, and (iv) the orthogonal projection of $x$ onto $\mathcal{V}_{w_{opt}}$ is $(P_{V}+RR^{\dagger})x$.

For the purpose of practical computations, we make a clear and rigorous analysis on the theoretical results and obtain a number of computable optimally expanded subspaces $\mathcal{V}_{\widetilde{w}_{opt}}$, which depend on chosen projection methods. As has already seen from (6) and the comments followed, it is hard to interpret $w_{opt}$, let alone a computable optimal replacement of $w_{opt}$. Fortunately, we are able to take a completely new approach to consider a computable optimal subspace expansion. Our key observation is: when expanding $\mathcal{V}$ to a computable optimal subspace, it follows from the fundamental identity

that seeking a computable optimal replacement of $w_{opt}$ is unnecessary. This identity provides us a new perspective and motivates us to find out a computable optimal replacement of $(I-P_{V})Aw_{opt}$ as a whole rather than $w_{opt}$ itself by a certain computable optimal one. As it will be clear, such an optimal replacement can be defined precisely and is deterministic within the framework of each of the standard, harmonic, refined, and refined harmonic Rayleigh–Ritz methods, respectively. As it turns out, computable optimal replacements are the Ritz vector, refined Ritz vector, harmonic Ritz vector and refined harmonic Ritz vector of $A$ from the subspace ${\rm span}\{R\}$ rather than $\mathcal{V}$ for these four kinds of projection methods. Taking the standard Rayleigh–Ritz method as an example, we will describe how to expand $\mathcal{V}$ to the computable optimal subspace $\mathcal{V}_{\widetilde{w}_{opt}}$. We will also show that, based on our results, there is some other novel optimal expansion that is independent of the desired $x$ in practical computations and is promising and more robust.

The paper is organized as follows. In Section 2, we consider the solution of the optimal subspace expansion problem (3) for a general $A$, present our theoretical results, and make an analysis on them. In Section 3, we show how to obtain computable optimal replacements of $(I-P_{V})Aw_{opt}$ and construct optimally expanded subspaces $\mathcal{V}_{\widetilde{w}_{opt}}$. We also present some other robust optimal expansion approach. In Section 4, we report numerical experiments to demonstrate the effectiveness of theoretical and computable optimal subspace expansion approaches, and compare them with the Lanczos or Arnoldi type expansion with $w=v_{k}$ and the the RA method, i.e., the Ritz expansion [24, 25] with $w$ being the Ritz vector from $\mathcal{V}$. In Section 5, we conclude the paper with some problems and issues that deserve future considerations.

Throughout the paper, denote by $I$ the identity matrix with the order clear from the context, by $\mathcal{C}^{k}$ the complex space of dimension $k$, and by $\mathcal{C}^{n\times k}$ or $\mathcal{R}^{n\times k}$ the set of $n\times k$ complex or real matrices.

For a general $A$, we first establish two results on $\cos\angle(\mathcal{V}_{w},x)$ for a given $w$. The first characterizes it as a maximization problem and generalizes Ye’s Theorem 2 in [25] and the second gives a maximizer of this problem. We remind the reader that these results are secondary. After them, we will present our major results.

From the definition of $r_{w}$, we have

By assumptions, we have $r_{w}\not=0$. Notice that any nonzero $a\in\mathcal{V}_{w}$ but $a\not\in\mathcal{V}$ can be uniquely written as

where $b\in\mathcal{V}$ and $\beta$ is a nonzero scalar. Since $Aw=r_{w}+\phi w$ and $x^{H}r_{w}=0$, we obtain

which establishes the maximization characterization in (8).

We now seek a maximizer of the maximization problem in (8). Since $x^{H}r_{w}=0$, by the definition of $Q_{w}$, there exists a vector $z_{w}\in\mathcal{C}^{n-1}$ such that the unit length eigenvector $x=Q_{w}z_{w}$ with $\|z_{w}\|=1$. As a result, we obtain $Q_{w}^{H}x=z_{w}$, $Q_{w}Q_{w}^{H}x=Q_{w}z_{w}=x$, and

Since $\frac{\|Q_{w}^{H}b\|}{\|b\|}=\sin\angle(b,r_{w})$, it follows from the above that

Therefore, from (10), (11), the orthonormality of $Q_{w}$ and $x=Q_{w}z_{w}$, writing a nonzero $b=Vy$, we obtain

since $(Q_{w}Q_{w}^{H}V)(Q_{w}Q_{w}^{H}V)^{\dagger}x$ is the orthogonal projection of $x$ onto the subspace ${\rm span}\{Q_{w}Q_{w}^{H}V\}$. As a result, $y_{w}=(Q_{w}Q_{w}^{H}V)^{\dagger}x$ solves $\max\limits_{y\in\mathcal{C}^{m},y\neq 0}\cos\angle(Q_{w}Q_{w}^{H}Vy,x)$, and $b_{w}=Vy_{w}$ is the first relation defined in (9), which proves the second relation in the right-hand side of (8).

We next prove that $Q_{w}Q_{w}^{H}V$ has full column rank. This amounts to showing that the cross-product matrix $V^{H}Q_{w}Q_{w}^{H}V$ is positive definite by noting that $(Q_{w}Q_{w}^{H})^{2}=Q_{w}Q_{w}^{H}$. To this end, it suffices to prove that the solution of the homogenous linear system $Q_{w}Q_{w}^{H}Vz=0$ is zero. Since $Q_{w}Q_{w}^{H}=I-\frac{r_{w}r_{w}^{H}}{\|r_{w}\|^{2}}$, the system becomes

which yields

Since

taking norms in the two sides of (12) gives

which holds only if $z=0$ or $\cos\angle(\mathcal{V},r_{w})=1$. The latter means that $r_{w}=Aw-\phi w\in\mathcal{V}$, i.e, $Aw\in\mathcal{V}$, a contradiction to our assumption. Hence we must have $z=0$, and $Q_{w}Q_{w}^{H}V$ has full column rank. Exploiting $Q_{w}Q_{w}^{H}x=x$, we obtain

which proves the second relation in (9).

As a matter of fact, Ye’s analysis in the Hermitian case does not seem theoretically sound and has some arbitrariness, causing that his claim may be problematic, as will be shown below.

Since

for any $w\in\mathcal{V}$, Ye attempts to seek a good approximation to the maximizer $w_{opt}$ of $\max_{w\in\mathcal{V}}\cos\angle(\mathcal{V}_{w},x)$. To this end, in the first relation of the right-hand side in (8), Ye takes $b\in\mathcal{V}$ to be the Ritz vector $z_{1}$ of $A$ from $\mathcal{V}$ that is used to approximate the desired $x$. Then (8) shows that

where $b_{w}$ is defined by (9). Notice that $b_{w}$ is a function of $w\in\mathcal{V}$. However, since $z_{1}$ is a constant vector and independent of $w\in\mathcal{V}$, there is no reason that $z_{1}$ is close to $b_{w}$ unless $b_{w}\approx z_{1}$ for some specific $w$ and $\mathcal{V}$, as we argue below.

Notice that $Q_{w}Q_{w}=I-\frac{r_{w}r_{w}^{H}}{\|r_{w}\|^{2}}$ and $r_{w}^{H}x=0$. Applying the Sherman–Morrison formula to $(V^{H}Q_{w}Q^{H}V)^{-1}$ and making use of $r_{w}^{H}P_{V}x=-r_{w}^{H}(I-P_{V})x$, by some manipulation we can justify that

where

with $\phi=\frac{x^{H}Aw}{x^{H}w}$. Observe that $\|\alpha_{w}w\|$ and $\|\beta_{w}Aw\|$ do not depend on the size of $\|w\|$ and are generally not small. Suppose that $w$ is normalized. In this case, $|\alpha_{w}|$ and $|\beta_{w}|$ are generally not small. Recall that $P_{V}x/\|P_{V}x\|$ is the best approximation to $x$ from $\mathcal{V}$ and $z_{1}$ is the computable approximation to $x$ obtained by the standard Rayleigh–Ritz method. Therefore, it is seen from (14) that the constant vector $z_{1}$ is generally not a good approximation to $b_{w}$ in direction for all $w\in\mathcal{V}$ unless $w\approx x$, which can occur only when $P_{V}x/\|P_{V}x\|\approx x$, that is, $\cos\angle(\mathcal{V},x)$ is already sufficiently small. These mean that the functions

generally have no similarity and their difference is not close to zero unless $z_{1}\approx x$. As a result, the maximizer of $\max_{w\in\mathcal{V}}\frac{\cos\angle(z_{1},x)}{\sin\angle(z_{1},r_{w})}$ generally bears no relation to the maximizer $w_{opt}$ of

The above analysis indicates that $b=z_{1}$ in the first relation of (8) is generally not a good approximation to $b_{w}$ and the resulting $\max_{w\in\mathcal{V}}\frac{\cos\angle(z_{1},x)}{\sin\angle(z_{1},r_{w})}$ is generally not a good replacement of $\max_{w\in\mathcal{V}}\cos\angle(\mathcal{V}_{w},x)$. In the Hermitian case, by using a few approximate and heuristic arguments, Ye [25] argues that $w=z_{1}$ may be a good approximate maximizer of $\max_{w\in\mathcal{V}}\frac{\cos\angle(z_{1},x)}{\sin\angle(z_{1},r_{w})}$ and then uses $z_{1}$ as a replacement of $w_{opt}$ that solves (3).

In what follows we give up any possible further analysis on Theorem 2.1 and instead consider (3) from new perspectives. We will establish a number of important and insightful results. We derive a new expression of $w_{opt}$, which is essentially the same as but formally different from (6) with $c=0$. The new form of $w_{opt}$ will play a critical role in establishing explicit expressions of the a priori $(I-P_{V})w_{opt}$, the optimally expanded subspace $\mathcal{V}_{w_{opt}}$ and some other important quantities.

Note that $P_{V}=VV^{H}$ is the orthogonal projector onto $\mathcal{V}$. Assume that $x\not\in\mathcal{V}$, $w\in\mathcal{V}$ and $Aw\not\in\mathcal{V}$. Lemma 1 of [25] states that

where $(I-P_{V})Aw=r=Ry$ with $R$ defined by (32) and $w=Vy$. We should remind the reader that Lemma 1 of [25] uses the form $r=Ry$, but we use the different form $(I-P_{V})Aw$ here by writing $R$ in the form of $(I-P_{V})AV$. (7) motivates us this change in form and enables us to establish more informative theoretical results on the optimal subspace expansion problem (3). Relation (15) indicates that

Let the matrix $(V,V_{\perp})$ be unitary. Then $I-P_{V}=V_{\perp}V_{\perp}^{H}$ and $(I-P_{V})AV=V_{\perp}V_{\perp}^{H}AV$. Define the vector

which is the orthogonal projection of $x$ onto the orthogonal complement of $\mathcal{V}$ with respect to $\mathcal{C}^{n}$ and is nonzero for $x\not\in\mathcal{V}$. Then the matrix pair

restricted to the orthogonal complement $\mathcal{N}^{\perp}(V_{\perp}^{H}AV)$ of $\mathcal{N}(V_{\perp}^{H}AV)$ is Hermitian definite, that is, the range restricted $V^{H}A^{H}(I-P_{V})AV=V^{H}A^{H}V_{\perp}V_{\perp}^{H}AV$ is Hermitian positive definite. Note that

the column space of $V^{H}A^{H}V_{\perp}$. We write the range restricted matrix pair (18) as

Write $w=Vy$. Then by assumption, $\mathcal{V}$ is not an invariant subspace of $A$, so that $(I-P_{V})AV=V_{\perp}V_{\perp}^{H}AV\not=0$. From (17) and $(I-P_{V})^{2}=I-P_{V}$, we obtain

which is the Rayleigh quotient of the range restricted Hermitian definite pair (19) with respect to the nonzero vector $y$. Keep in mind that $w_{opt}$ satisfies (16). Therefore, by the min-max characterization of the eigenvalues of a Hermitian definite pair (cf. e.g., [20, pp. 281-2]), the solution to the maximization problem

is an eigenvector $y_{opt}$ of the range restricted definite pair (19) associated with its largest eigenvalue $\mu_{opt}$, and the optimal expansion vector $w_{opt}=Vy_{opt}$ up to scaling.

Notice that $V^{H}A^{H}x_{\perp}x_{\perp}^{H}AV\big{|}_{\mathcal{R}(V^{H}A^{H}V_{\perp})}$ is a rank-one Hermitian positive semi-definite matrix, and denote

Therefore, the range restricted Hermitian definite pair (19) has exactly one positive eigenvalue, and the other eigenvalues are zeros. By (28), the eigenvalues $\mu$ of the matrix pair (19) are identical to those of the rank-one Hermitian matrix

and the eigenvectors $y$ of the pair (19) are related to the eigenvectors $\hat{y}$ of the above matrix by

It is known that the unique positive eigenvalue $\mu_{opt}$ of the matrix in (29) is

and its associated (unnormalized) eigenvector is

Substituting this into (30) and $w_{opt}=Vy_{opt}$ gives the first relation in the right-hand side of (20). Exploiting $P_{V}=VV^{H}$ and $I-P_{V}=V_{\perp}V_{\perp}^{H}$, by the definition of Moore–Penrose generalized inverse, we obtain the second relation in the right-hand side of (20).

By (5), the residual

which shows that $(I-P_{V})AP_{V}=RV^{H}$. As a result, we have

which, together with the second relation in (20), proves

i.e., the third relation in (20). Therefore, from (32) and (34), we obtain

which establishes (21).