Backward Stability of Explicit External Deflation for the Symmetric Eigenvalue Problem

Let $A$ be an $n\times n$ real symmetric matrix with ordered eigenvalues $\lambda_{1}\leq\lambda_{2}\leq\cdots\leq\lambda_{n}$, and the corresponding eigenvectors $v_{1},v_{2},\ldots,v_{n}$. Let $\mathcal{I}=[\lambda_{\rm low},\lambda_{\rm upper}]$ be an interval containing the eigenvalues $\lambda_{1},\lambda_{2},\ldots,\lambda_{n_{e}}$ at the lower end of the spectrum of $A$. We are interested in computing the partial eigenvalue decomposition

$$AV_{n_{e}}=V_{n_{e}}\Lambda_{n_{e}},$$

(1.1)

where $\Lambda_{n_{e}}=\mbox{diag}(\lambda_{1},\lambda_{2},\ldots,\lambda_{n_{e}})$, $V_{n_{e}}=[v_{1},v_{2},\ldots,v_{n_{e}}]$ and $V^{\mathrm{T}}_{n_{e}}V_{n_{e}}=I_{n_{e}}$. We are particularly interested in the case where the interval is of the width $|\mathcal{I}|=\lambda_{\rm upper}-\lambda_{\rm low}\leq\frac{1}{2}\|A\|_{2}$, and may contain a large number of eigenvalues of $A$.

A majority of subspace projection methods for computing the partial eigenvalue decomposition (1.1), such as the Lanczos methods [13, 1], typically converge first to the eigenvalues on the periphery of the spectrum. A projection subspace of a large dimension is necessary to compute the eigenvalues located deep inside the spectrum. Maintaining the orthogonality of a large projection subspace is computationally expensive. To reduce the dimensions of projection subspaces and accelerate the convergence, a number of techniques such as the restarting [17, 8, 20] and the filtering [11] have been developed.

This paper revisits a classical deflation technique known as Hotelling’s deflation [6, 19, 13]. Hotelling’s deflation computes the partial eigenvalue decomposition (1.1) through a sequence of low-rank updates. In the simplest case, supposing that the lowest eigenpair $(\lambda_{1},v_{1})$ of $A$ is computed by an eigensolver, by choosing a real shift $\sigma_{1}$ and defining the rank-one update matrix

$\displaystyle A_{1}\equiv A+\sigma_{1}v_{1}v_{1}^{\mathrm{T}},$

the eigenpair $(\lambda_{1},v_{1})$ of $A$ is displaced to an eigenpair $(\lambda_{1}+\sigma_{1},v_{1})$ of $A_{1}$, while all the other eigenpairs of $A$ and $A_{1}$ are the same. If the shift $\sigma_{1}>\lambda_{\rm upper}-\lambda_{1}$, then the eigenpair $(\lambda_{1},v_{1})$ is shifted outside the interval $\mathcal{I}$ and the eigenpair $(\lambda_{2},v_{2})$ of $A$ becomes the lowest eigenpair of $A_{1}$. Subsequently, we can use the eigensolver again to compute the lowest eigenpair $(\lambda_{2},v_{2})$ of $A_{1}$. The procedure is repeated until all the eigenvalues in the interval $\mathcal{I}$ are found. Since the low-rank update of the original $A$ can be done outside of the eigensolver, we refer to this strategy as an explicit external deflation (EED) procedure, to distinguish from the deflation techniques used inside an eigensolver, such as TRLan [20] and ARPACK [9]. Recently the EED has been combined with Krylov subspace methods to progressively compute the eigenpairs toward the interior of the spectrum of symmetric eigenvalue problems [12, 21], and extended to structured eigenvalue problems [2, 3].

In the presence of finite precision arithmetic, the EED procedure is susceptible to numerical instability due to the accumulation of numerical errors from previous computations; see, e.g., [19, p. 585] and [15, p. 96]. In [13, Chap. 5.1], Parlett showed that the change to the smallest eigenvalue in magnitude by deflating out the largest computed eigenvalue would be at the same order of the round-off error incurred. In [14], for nonsymmetric matrices, Saad derived an upper bound on the backward error norm and showed that the stability is determined by the angle between the computed eigenvectors and the deflated subspace. Saad concluded that “if the angle between the computed eigenvector and the previous invariant subspace is small at every step, the process may quickly become unstable. On the other hand, if this is not the case then the process is quite safe, for small number of requested eigenvalues.”

However, none of the existing studies had examined the effect of the choice of shifts on the numerical stability. In this paper, we present a thorough backward stability analysis of the EED procedure. We first derive governing equations of the EED procedure in the presence of finite precision arithmetic, and then derive upper bounds of the loss of the orthogonality of the computed eigenvectors and the symmetric backward error norm of the computed eigenpairs. From these upper bounds, under mild assumptions, we conclude that the EED procedure is backward stable if the shifts $\sigma_{j}$ are dynamically chosen such that (i) the spectral gaps, defined as the separation between the computed eigenvalues and shifted ones, are maintained at the same order of the norm of the matrix $A$, and (ii) the shift-gap ratios, defined as the ratios between the largest shift in magnitude and the spectral gaps, are at the order of one. Our analysis implies that the re-orthogonalization of computed eigenvectors and the large angles between the computed eigenvectors and the previous invariant subspaces are unnecessary. Based on theoretical analysis, we provide a practical shift selection scheme to guarantee the backward stability. For a set of test matrices from applications, we demonstrate that the EED in combination with an established eigensolver can compute a large number of eigenvalues with backward stability.

The rest of the paper is organized as follows. In Section 2, we present the governing equations of the EED in finite precision arithmetic. In Section 3, we derive the upper bounds on the loss of orthogonality of computed eigenvectors and the symmetric backward error norm of computed eigenpairs. From these upper bounds, we derive conditions for the backward stability of the EED procedure. In Section 4, we discuss how to properly choose the shifts to satisfy the stability conditions and outline an algorithm for combining the EED procedure with an eigensolver. Numerical results to verify the sharpness of the upper bounds and to demonstrate the backward stability of the EED for symmetric eigenvalue problems from applications are presented in Section 5. Concluding remarks are given in section 6.

Following the convention of matrix computations, we use the upper case letters for matrices and the lower case letters for vectors. In particular, we use $I_{n}$ for the identity matrix of dimension $n$. If not specified, the dimensions of matrices and vectors conform to the dimensions used in the context. We use $\cdot^{\mathrm{T}}$ for transpose, and $\|\cdot\|_{2}$ and $\|\cdot\|_{\rm F}$ for $2$-norm and the Frobenius norm, respectively. The range of a matrix $A$ is denoted by $\mathcal{R}(A)$. We also use $\sigma_{\min}(\cdot)$ to denote the minimal singular values of a matrix. Other notations will be explained as used.

Suppose that we use an established eigensolver called EIGSOL, such as TRLan [20] and ARPACK [9], to successfully compute the lowest eigenpair $(\widehat{\lambda},\widehat{v})$ of $A$ with

$$A\widehat{v}=\widehat{\lambda}\widehat{v}+\eta,$$

where $\|\widehat{v}\|_{2}=1$ and the residual vector $\eta$ satisfies

$$\|\eta\|_{2}\leq tol\cdot\|A\|_{2}$$

(2.1)

for a prescribed convergence tolerance $tol$. The EED procedure starts with computing the lowest eigenpair $(\widehat{\lambda}_{1},\widehat{v}_{1})$ of $A$ by EIGSOL satisfying

$\displaystyle A\widehat{v}_{1}=\widehat{\lambda}_{1}\widehat{v}_{1}+\eta_{1},$

where $\widehat{\lambda}_{1}\in\mathcal{I}=[\lambda_{\rm low},\lambda_{\rm upper}]$, and the residual vector $\eta_{1}$ satisfies (2.1). At the first EED step, we choose a shift $\sigma_{1}$ and define

$\displaystyle\widehat{A}_{1}\equiv A+\sigma_{1}\widehat{v}_{1}\widehat{v}_{1}^{\mathrm{T}}.$

By choosing the shift $\sigma_{1}>\lambda_{\rm upper}-\widehat{\lambda}_{1}$, the lowest eigenpair of $\widehat{A}_{1}$ is an approximation of the second eigenpair $({\lambda}_{2},{v}_{2})$ of $A$. Subsequently, we use EIGSOL to compute the lowest eigenpair $(\widehat{\lambda}_{2},\widehat{v}_{2})$ of $\widehat{A}_{1}$ satisfying

$\displaystyle\widehat{A}_{1}\widehat{v}_{2}=\widehat{\lambda}_{2}\widehat{v}_{2}+\eta_{2},$

where the residual vector $\eta_{2}$ satisfies (2.1). Meanwhile, expressing the computed eigenpair $(\widehat{\lambda}_{1},\widehat{v}_{1})$ in terms of $\widehat{A}_{1}$, we have

$$\widehat{A}_{1}\widehat{v}_{1}=(\widehat{\lambda}_{1}+\sigma_{1})\widehat{v}_{1}+\eta_{1}.$$

Proceeding to the second EED step, we choose a shift $\sigma_{2}$ and define

$$\widehat{A}_{2}\equiv\widehat{A}_{1}+\sigma_{2}\widehat{v}_{2}\widehat{v}_{2}^{\mathrm{T}}=A+\widehat{V}_{2}\Sigma_{2}\widehat{V}_{2}^{\mathrm{T}},$$

where $\widehat{V}_{2}=[\widehat{v}_{1},\widehat{v}_{2}]$ and $\Sigma_{2}=\operatorname{diag}(\sigma_{1},\sigma_{2})$. By choosing the shift $\sigma_{2}>\lambda_{\rm upper}-\widehat{\lambda}_{2}$, the lowest eigenpair of $\widehat{A}_{2}$ is an approximation of the third eigenpair $({\lambda}_{3},{v}_{3})$ of $A$. Then we use EIGSOL again to compute the lowest eigenpair $(\widehat{\lambda}_{3},\widehat{v}_{3})$ of $\widehat{A}_{2}$ satisfying

$\displaystyle\widehat{A}_{2}\widehat{v}_{3}=\widehat{\lambda}_{3}\widehat{v}_{3}+\eta_{3},$

where the residual vector $\eta_{3}$ satisfies (2.1). Meanwhile, expressing the computed eigenpairs $(\widehat{\lambda}_{1},\widehat{v}_{1})$ and $(\widehat{\lambda}_{2},\widehat{v}_{2})$ in terms of $\widehat{A}_{2}$, we have

$\displaystyle\widehat{A}_{2}\widehat{V}_{2}=\widehat{V}_{2}(\widehat{\Lambda}_{2}+\Sigma_{2})+\widehat{V}_{2}\Sigma_{2}\Phi_{2}+E_{2},$

where $\widehat{\Lambda}_{2}=\operatorname{diag}(\widehat{\lambda}_{1},\widehat{\lambda}_{2})$, $E_{2}=[\eta_{1},\eta_{2}]$, and $\Phi_{2}\in{\mathbb{R}}^{2\times 2}$ is the strictly lower triangular part of the matrix $\widehat{V}_{2}^{\mathrm{T}}\widehat{V}_{2}-I_{2}$, i.e., $\Phi_{2}+\Phi^{\mathrm{T}}_{2}=\widehat{V}_{2}^{\mathrm{T}}\widehat{V}_{2}-I_{2}$.

In general, at the $j$-th EED step, we choose a shift $\sigma_{j}$ and define

$\displaystyle\widehat{A}_{j}\equiv\widehat{A}_{j-1}+\sigma_{j}\widehat{v}_{j}\widehat{v}_{j}^{\mathrm{T}}=A+\widehat{V}_{j}\Sigma_{j}\widehat{V}_{j}^{\mathrm{T}},$

(2.2)

where $\widehat{V}_{j}\equiv[\widehat{v}_{1},\ldots,\widehat{v}_{j}]$ and $\Sigma_{j}=\operatorname{diag}(\sigma_{1},\ldots,\sigma_{j})$ with $\widehat{A}_{0}\equiv A$. Then by choosing the shift $\sigma_{j}>\lambda_{\rm upper}-\widehat{\lambda}_{j}$, the lowest eigenpair of $\widehat{A}_{j}$ is an approximation of the $(j+1)$-th eigenpair $({\lambda}_{j+1},{v}_{j+1})$ of $A$. We use EIGSOL to compute the lowest eigenpair $(\widehat{\lambda}_{j+1},\widehat{v}_{j+1})$ of $\widehat{A}_{j}$ satisfying

$\displaystyle\widehat{A}_{j}\widehat{v}_{j+1}=\widehat{\lambda}_{j+1}\widehat{v}_{j+1}+\eta_{j+1},$

(2.3)

where $\|\widehat{v}_{j+1}\|_{2}=1$ and the residual vector $\eta_{j+1}$ satisfies (2.1), i.e.,

$\displaystyle\|\eta_{j+1}\|_{2}\leq tol\cdot\|A\|_{2}.$

(2.4)

Meanwhile, for the computed eigenpairs $(\widehat{\lambda}_{j},\widehat{v}_{j})$ in terms of $\widehat{A}_{j}$, we have

$\displaystyle\widehat{A}_{j}\widehat{v}_{j}=\widehat{A}_{j-1}\widehat{v}_{j}+\sigma_{j}\widehat{v}_{j}=(\lambda_{j}+\sigma_{j})\widehat{v}_{j}+\eta_{j},$

(2.5)

and for the computed eigenpairs $(\widehat{\lambda}_{i},\widehat{v}_{i})$ with $1\leq i\leq j-1$ in terms of $\widehat{A}_{j}$, we have

	$\displaystyle\widehat{A}_{j}\widehat{v}_{i}$	$\displaystyle=\left(\widehat{A}_{j-1}+\sigma_{j}\widehat{v}_{j}\widehat{v}_{j}^{\mathrm{T}}\right)\widehat{v}_{i}$
		$\displaystyle=\left(\widehat{A}_{i-1}+\sigma_{i}\widehat{v}_{i}\widehat{v}_{i}^{\mathrm{T}}+\widehat{V}_{i+1:j}\Sigma_{i+1:j}\widehat{V}_{i+1:j}^{\mathrm{T}}\right)\widehat{v}_{i}$
		$\displaystyle=\left(\widehat{A}_{i-1}+\sigma_{i}\widehat{v}_{i}\widehat{v}_{i}^{\mathrm{T}}\right)\widehat{v}_{i}+\widehat{V}_{i+1:j}\Sigma_{i+1:j}\widehat{V}_{i+1:j}^{\mathrm{T}}\widehat{v}_{i}$
		$\displaystyle=\widehat{\lambda}_{i}\widehat{v}_{i}+\eta_{i}+\sigma_{i}\widehat{v}_{i}+\widehat{V}_{i+1:j}\Sigma_{i+1:j}\widehat{V}_{i+1:j}^{\mathrm{T}}\widehat{v}_{i}$
		$\displaystyle=(\widehat{\lambda}_{i}+\sigma_{i})\widehat{v}_{i}+\widehat{V}_{i+1:j}\Sigma_{i+1:j}\widehat{V}_{i+1:j}^{\mathrm{T}}\widehat{v}_{i}+\eta_{i}$
		$\displaystyle=(\widehat{\lambda}_{i}+\sigma_{i})\widehat{v}_{i}+\widehat{V}_{j}\Sigma_{j}\left[\begin{array}[]{c}0\\ \widehat{V}_{i+1:j}^{\mathrm{T}}\widehat{v}_{i}\end{array}\right]+\eta_{i},$		(2.8)

where $\widehat{V}_{i+1:j}\equiv[\widehat{v}_{i+1},\ldots,\widehat{v}_{j}]$ and $\Sigma_{i+1:j}\equiv\operatorname{diag}(\sigma_{i+1},\ldots,\sigma_{j})$.

Combining (2.5) and (2.8), we have

$\displaystyle\widehat{A}_{j}\widehat{V}_{j}=\widehat{V}_{j}(\widehat{\Lambda}_{j}+\Sigma_{j})+\widehat{V}_{j}\Sigma_{j}\Phi_{j}+E_{j},$

(2.9)

where $\widehat{\Lambda}_{j}=\operatorname{diag}(\widehat{\lambda}_{1},\ldots,\widehat{\lambda}_{j})$, $E_{j}=[\eta_{1},\ldots,\eta_{j}]$, $\Phi_{j}$ is the strictly lower triangular part of the matrix $\widehat{V}_{j}^{\mathrm{T}}\widehat{V}_{j}-I_{j}$ and $\Phi_{j}+\Phi^{\mathrm{T}}_{j}=\widehat{V}_{j}^{\mathrm{T}}\widehat{V}_{j}-I_{j}$. Eqs. (2.3) and (2.9) are referred to as the governing equations of $j$ steps of (inexact) EED procedure.

For the backward stability analysis of the EED procedure, we introduce the following two quantities associated with the shifts $\sigma_{1},\dots,\sigma_{j}$ for a $j$-step EED:

• •

  the spectral gap of $\widehat{A}_{j}$, defined as the separation between the computed eigenvalues and the shifted ones:

  $$\gamma_{j}\equiv\min_{\lambda\in\mathcal{I}_{j+1},\theta\in\mathcal{J}_{j}}|\lambda-\theta|>0,$$

  (2.10)

  where $\mathcal{I}_{j+1}\equiv\{\widehat{\lambda}_{1},\ldots,\widehat{\lambda}_{j},\widehat{\lambda}_{j+1}\}$, the set of computed eigenvalues, and $\mathcal{J}_{j}\equiv\{\widehat{\lambda}_{1}+\sigma_{1},\ldots,\widehat{\lambda}_{j}+\sigma_{j}\}$, the set of computed eigenvalues after shifting (see Figure 1 for an illustration);

• •

  the shift-gap ratio, defined as the ratio of the largest shift and the spectral gap $\gamma_{j}$:

  $$\tau_{j}\equiv\frac{1}{\gamma_{j}}\cdot\max_{1\leq i\leq j}|\sigma_{i}|.$$

  (2.11)

We will see that $\gamma_{j}$ and $\tau_{j}$ are crucial quantities to characterize the backward stability of the EED procedure.

In this section, we first review the notion of backward stability for the computed eigenpairs. Then we derive upper bounds on the loss of orthogonality of computed eigenvectors and the backward error norm of computed eigenpairs by the EED procedure, and reveal conditions for the backward stability of the procedure.