A skew-symmetric Lanczos bidiagonalization method for computing several largest eigenpairs of a large skew-symmetric matrix

Let $A\in\mathbb{R}^{n\times n}$ be a large scale and possibly sparse skew-symmetric matrix, i.e., $A^{T}=-A$, where the superscript $T$ denotes the transpose of a matrix or vector. We consider the eigenvalue problem

where $\lambda\in\mathbb{C}$ is an eigenvalue of $A$ and $x\in\mathbb{C}^{n}$ with the 2-norm $\|x\|=1$ is an corresponding eigenvector. It is well known that the nonzero eigenvalues $\lambda$ of $A$ are purely imaginary and come in conjugate pairs. The real skew-symmetric eigenvalue problem arises in a variety of applications, such as matrix function computations [4, 7], the solution of linear quadratic optimal control problems [19, 39], model reductions [22, 37], wave propagation solution for repetitive structures [12, 38], crack following in anisotropic materials [1, 21], and some others [20, 26, 36].

For small to medium sized skew-symmetric matrices, several numerical algorithms have been developed to compute their spectral decompositions in real arithmetic [8, 23, 26, 35]. For a large scale skew-symmetric $A$, however, it still lacks a high-performance numerical algorithm that can make use of its skew-symmetric structure to compute several eigenvalues and/or the corresponding eigenvectors only in real arithmetic. In this paper, we close this gap by proposing a Krylov subspace type method to compute several extreme conjugate eigenvalues in magnitude and the corresponding eigenvectors of $A$ that only uses real arithmetic.

Our work originates from a basic fact that, as will be shown in Section 2, the spectral decomposition of the skew-symmetric $A$ has a close relationship with its specifically structured singular value decomposition (SVD), where each conjugate pair of purely imaginary eigenvalues of $A$ corresponds to a double singular value of $A$ and the mutually orthogonal real and imaginary parts of associated eigenvectors are the corresponding left and right singular vectors of $A$. Therefore, the solution of eigenvalue problem of the skew-symmetric $A$ is equivalent to the computation of its structured SVD, and the desired eigenpairs can be recovered from the converged singular triplets by exploiting this equivalence.

It has been well realized in [35] that the spectral decomposition of the skew-symmetric $A$ can be reduced to the SVD of a certain bidiagonal matrix $B$ whose size is $n/2$ for $n$ even. Exploiting such property, Ward and Gray [35] have designed an eigensolver for dense skew-symmetric matrices. For $A$ large scale, however, the computation of $B$ using orthogonal similarity transformations is prohibitive due to the requirement of excessive storage and computations. Nevertheless, with the help of Lanczos bidiagonalization (LBD) [24], whose complete bidiagonal reduction is originally due to [9], we are able to compute a sequence of leading principal matrices of $B$, based on which, approximations to the extreme singular values of $A$ and/or the corresponding singular vectors can be computed. Along this line, a number of LBD based methods have been proposed and intensively studied, and several practical explicitly and implicitly restarted LBD type algorithms have been well developed [2, 3, 5, 11, 14, 15, 17, 18, 34].

For the skew-symmetric matrix $A$, LBD owns some unique properties in both mathematics and numerics. However, direct applications of the aforementioned algorithms do not make use of these special properties, and will thus encounter some severe numerical difficulties. We will fully exploit these properties to propose a skew-symmetric LBD (SSLBD) method for computing several extreme singular triplets of $A$, from which its extreme conjugate eigenvalues in magnitude and the associate eigenvectors are recovered. Attractively, our algorithm is performed in real arithmetic when computing complex eigenpairs of $A$. Particularly, one conjugate pair of approximations to two conjugate pairwise eigenpairs can be recovered from one converged approximate singular triplet. We estimate the distance between the subspace generated by the real and imaginary parts of a pair of desired conjugate eigenvectors and the Krylov subspaces generated by SSLBD, and prove how it tends to zero as the subspace dimension increases. With these estimates, we establish a priori error bounds for the approximate eigenvalues and approximate eigenspaces obtained by the SSLBD method. One of them is the error bound for the two dimensional approximate eigenspace, which extends [29, pp.103, Theorem 4.6], a classic error bound for the Ritz vector in the real symmetric (complex Hermitian) case. These results show how fast our method converges when computing several extreme singular triplets.

As the subspace dimension increases, the SSLBD method will eventually become impractical due to the excessive memory and computational cost. For a practical purpose, it is generally necessary to restart the method by selecting an increasingly better starting vector when the maximum iterations allowed are attained, until the restarted algorithm converges ultimately. We will develop an implicitly restarted SSLBD algorithm. Initially, implicit restart was proposed by Sorensen [31] for large eigenvalue problems, and then it has been developed by Larsen [18] and Jia and Niu [14, 15] for large SVD computations. Our implicit restart is based on them but specialized to the skew-symmetric $A$. We will show that two sets of left and right Lanczos vectors generated by SSLBD are meanwhile biorthogonal for the skew-symmetric $A$. However, in finite precision arithmetic, this numerical biorthogonality loses gradually even if they themselves are kept numerically orthonormal. As a matter of fact, we have found that the numerical orthogonality of the computed left and right Lanczos vectors in an available state-of-the-art Lanczos bidiagonalization type algorithm does not suffice for the SVD problem of the skew-symmetric $A$, and unfortunate ghosts appear when the numerical biorthogonality loses severely; that is, a singular value of $A$ may be recomputed a few times. Therefore, certain numerical biorthogonality is crucial to make the method work regularly.

As has been proved by Simon [30] for the symmetric Lanczos method on the eigenvalue problem and extended by Larsen [17] to the LBD method for the SVD problem, as far as the accurate computation of singular values is concerned, the numerical semi-orthogonality suffices for a general $A$, and the numerical orthogonality to the level of machine precision is not necessary. Here semi-orthogonality means that two vectors are numerically orthogonal to the level of the square root of machine precision. To maintain the numerical stability and make the SSLBD method behave like it in exact arithmetic, we shall prove that the semi-orthogonality of each set of Lanczos vectors and the semi-biorthogonality of left and right Lanczos vectors suffice for the accurate computation of singular values of the skew-symmetric $A$. We will seek for an effective and efficient reorthogonalization strategy for this purpose. Specifically, besides the commonly used partial reorthogonalization to maintain the desired semi-orthogonality of each set of Lanczos vectors, we will propose a partial reorthogonalization strategy to maintain the numerical semi-biorthogonality of two sets of Lanczos vectors in order to avoid the ghost phenomena and make the algorithm converge regularly. We will introduce such kind of reorthogonalization strategy into the implicitly restarted SSLBD algorithm and design a simple scheme to recover the desired conjugate eigenpairs of $A$ pairwise from the converged singular triplets.

The rest of this paper is organized as follows. In Section 2, we describe some properties of the spectral decomposition of $A$, and show its mathematical equivalence with the SVD of $A$. Then we establish a close relationship between the SVD of $A$ and that of an upper bidiagonal matrix whose order is half of that of $A$ when the order of $A$ is even. In Section 3, we present some properties of the SSLBD process, and propose the SSLBD method for computing several extreme singular values and the corresponding singular vectors of $A$, from which conjugate approximate eigenpairs of $A$ are pairwise reconstructed. In Section 4, we establish convergence results on the SSLBD method. In Section 5, we design a partial reorthogonalization scheme for the SSLBD process, and develop an implicitly restarted SSLBD algorithm to compute several extreme singular triplets of $A$. Numerical experiments are reported in Section 6. We conclude the paper in Section 7.

Throughout this paper, we denote by $\mathrm{i}$ the imaginary unit, by $X^{H}$ and $X^{{\dagger}}$ the conjugate transpose and pseudoinverse of $X$, by $\mathcal{R}(X)$ the range space of $X$, and by $\|X\|$ and $\|X\|_{F}$ the $2$- and Frobenius norms of $X$, respectively. We denote by $\mathbb{R}^{k}$ and $\mathbb{C}^{k}$ the real and complex spaces of dimension $k$, by $I_{k}$ the identity matrix of order $k$, and by $\bm{0}_{k}$ and $\bm{0}_{k,\ell}$ the zero matrices of orders $k\times k$ and $k\times\ell$, respectively. The subscripts of identity and zero matrices are omitted when they are clear from the context.

For the skew-symmetric matrix $A\in\mathbb{R}^{n\times n}$, the following five properties are well known and easily justified: Property (i): $A$ is diagonalizable by a unitary similarity transformation; Property (ii): the eigenvalues $\lambda$ of $A$ are either purely imaginary or zero; Property (iii): for an eigenpair $(\lambda,z)$ of $A$, the conjugate $(\bar{\lambda},\bar{z})$ is also an eigenpair of $A$; Property (iv): the real and imaginary parts of the eigenvector $z$ corresponding to a purely imaginary eigenvalue of $A$ have the same length, and are mutually orthogonal; Property (v): $A$ of odd order must be singular, and has at least one zero eigenvalue. It is easily deduced that a nonsingular $A$ must have even order, i.e., $n=2\ell$ for some integer $\ell$, and its eigenvalues $\lambda$ are purely imaginary and come in conjugate, i.e., plus-minus, pairs.

We formulate the following basic results as a theorem, which establishes close structure relationships between the spectral decomposition of $A$ and its SVD and forms the basis of our method and algorithm to be proposed. For completeness and our frequent use in this paper, we will give a rigorous proof by exploiting the above five properties.

By Properties (i) and (ii) and the assumption that $A$ is nonsingular, we denote by $\Lambda_{+}\in\mathbb{C}^{\ell\times\ell}$ the diagonal matrix consisting of all the eigenvalues of $A$ with positive imaginary parts, and by $X_{+}\in\mathbb{C}^{n\times\ell}$ the corresponding orthonormal eigenvector matrix. By Property (iv), we can write $(\Lambda_{+},X_{+})$ as

where the columns of $U$ and $V$ have unit-length, and $AX_{+}=X_{+}\Lambda_{+}$ and $X_{+}^{H}X_{+}=I$. By these and Property (iii), the conjugate of $(\Lambda_{+},X_{+})$ is

and $AX_{-}=X_{-}\Lambda_{-}$ and $X_{-}^{H}X_{-}=I$. As a result, $X$ and $\Lambda$ defined by (2.2) are the eigenvector and eigenvalue matrices of $A$, respectively, and

Premultiplying the two sides of $AX_{+}=X_{+}\Lambda_{+}$ by $X_{+}^{T}$ delivers $X_{+}^{T}AX_{+}=X_{+}^{T}X_{+}\Lambda_{+}$, which, by the skew-symmetry of $A$, shows that

Since $\Lambda_{+}=\mathrm{i}\Sigma$ and $\Sigma$ is a positive diagonal matrix, it follows from the above equation that

By definition (2.2) of $X$ and the fact that $X_{+}$ and $X_{-}$ are column orthonormal, $X$ is unitary. Postmultiplying (2.7) by $X^{H}$ yields (2.1), from which it follows that (2.4) holds.

Inserting $X_{+}=\frac{1}{\sqrt{2}}(U+\mathrm{i}V)$ into $X_{+}^{H}X_{+}=I$ and $X_{+}^{T}X_{+}=\bm{0}$ and solving the resulting equations for $U^{T}U,V^{T}V$ and $U^{T}V$, we obtain (2.3). This shows that both $[U,V]$ and $[V,-U]$ are orthogonal. Therefore, (2.4) is the SVD of $A$.

For $U$ and $V$ in (2.4), write

On the basis of Theorem 2.1, in the sequel, we denote by

The SVD (2.4) indicates that $(\sigma_{j},u_{j},v_{j})$ and $(\sigma_{j},v_{j},-u_{j})$ are two singular triplets of $A$ corresponding to the multiple singular value $\sigma_{j}$. Therefore, in order to obtain the conjugate eigenpairs $(\lambda_{\pm j},x_{\pm j})$ of $A$, one can compute only the singular triplet $(\sigma_{j},u_{j},v_{j})$ or $(\sigma_{j},v_{j},-u_{j})$, and then recovers the desired eigenpairs from the singular triplet by (2.9).

Next we present a bidiagonal decomposition of the nonsingular skew-symmetric $A$ with order $n=2\ell$, which essentially appears in [35]. For completeness, we include a proof.

It is known from, e.g., [35] that $A$ is orthogonally similar to a skew-symmetric tridiagonal matrix, that is, there exists an orthogonal matrix $F\in\mathbb{R}^{n\times n}$ such that

Denote $\Pi_{1}=[e_{2,n},e_{4,n},\dots,e_{2\ell,n}]$ and $\Pi_{2}=[e_{1,n},e_{3,n},\dots,e_{2\ell-1,n}]$ with $e_{j,n}$ being the $j$th column of $I_{n}$, $j=1,\dots,n$. Then it is easy to verify that

Combining these two relations with (2.11), we obtain

which proves (2.10), where we have denoted $P=F\Pi_{1}$ and $Q=F\Pi_{2}$. Note that $F$ is orthogonal and $\Pi_{1}$, $\Pi_{2}$ are not only orthonormal but also biorthogonal: $\Pi_{1}^{T}\Pi_{1}=\Pi_{2}^{T}\Pi_{2}=I$ and $\Pi_{1}^{T}\Pi_{2}=\mathbf{0}$. Therefore, $P^{T}P=Q^{T}Q=I$ and $P^{T}Q=\bm{0}$.

Combining Remark 2 with Remark 1 after Theorem 2.1, for ease of presentation, in the sequel we will always assume that $A$ has even order and is nonsingular. However, our method and algorithm to be proposed apply to a singular skew-symmetric $A$ with only some slight modifications.

Decomposition (2.10) is a bidiagonal decomposition of $A$ that reduces $A$ to $\mathop{\operator@font diag}\nolimits\{B,B^{T}\}$ using the structured orthogonal matrices $\begin{bmatrix}P&Q\end{bmatrix}^{T}$ and $\begin{bmatrix}Q&-P\end{bmatrix}$ from the left and right, respectively. The singular values of $A$ are the union of those of $B$ and $B^{T}$ but the SVD of $B$ alone delivers the whole SVD of $A$: Let $B=U_{B}\Sigma_{B}V_{B}^{T}$ be the SVD of $B$. Then Theorem 2.2 indicates that $(\Sigma,U,V)=(\Sigma_{B},PU_{B},QV_{B})$ is an $\ell$-dimensional partial SVD of $A$ with $U$ and $V$ being orthonormal and $U^{T}V=\mathbf{0}$. By Theorem 2.1, the spectral decomposition (2.1) of $A$ can be restored directly from $(\Sigma,U,V)$. This is exactly the working mechanism of the eigensolver proposed in [35] for small to medium sized skew-symmetric matrices.

For $A$ large, using orthogonal transformations to construct $P$, $Q$ and $B$ in (2.10) and computing the SVD of $B$ are unaffordable. Fortunately, based on decomposition (2.10), we can perform the LBD process on the skew-symmetric $A$, i.e., the SSLBD process, to successively generate the columns of $P$ and $Q$ and the leading principal matrices of $B$, so that the Rayleigh–Ritz projection comes into the computation of a partial SVD of $A$.

For ease of presentation, from now on, we always assume that the eigenvalues $\lambda_{\pm j}=\pm\mathrm{i}\sigma_{j}$ of $A$ in (2.9) are simple, and that $\sigma_{1},\dots,\sigma_{\ell}$ are labeled in decreasing order:

Our goal is to compute the $k$ pairs of extreme, e.g., largest conjugate eigenvalues $\lambda_{\pm 1},\dots,\lambda_{\pm k}$ and/or the corresponding eigenvectors $x_{\pm 1},\dots,x_{\pm k}$ defined by (2.9). By Theorem 2.1, this amounts to computing the following partial SVD of $A$:

For later use, for each $j=1,\dots,\ell$, we denote

Then $AX_{j}=X_{j}\Lambda_{j}$, and $\mathcal{X}_{j}:=\mathcal{R}(X_{j})=\mathrm{span}\{x_{\pm j}\}$ is the two dimensional eigenspace of $A$ associated with the conjugate eigenvalues $\lambda_{\pm j}$, $j=1,\dots,\ell$. We call the pair $(\Lambda_{j},X_{j})$ a block eigenpair of $A$, $j=1,\dots,\ell$. From (2.3), $\widehat{X}=[X_{1},\dots,X_{\ell}]$ is orthogonal, and (2.1) can be written as

To make a convergence analysis on the SSLBD method, we review some preliminaries. Let the columns of $Z$, $W$ and $W_{\perp}$ form orthonormal base of $\mathcal{Z}$, $\mathcal{W}$ and the orthogonal complement of $\mathcal{W}$, respectively. Denote by $\angle(\mathcal{W},\mathcal{Z})$ the canonical angles between $\mathcal{W}$ and $\mathcal{Z}$ [10, pp.329–330]. The 2-norm distance between $\mathcal{Z}$ and $\mathcal{W}$ is defined by

which is the sine of the largest canonical angle between $\mathcal{W}$ and $\mathcal{Z}$ (cf. Jia and Stewart [16]). It is worthwhile to point out that if the dimensions of $\mathcal{W}$ and $\mathcal{Z}$ are not equal then this measure is not symmetric in its arguments. In fact, if the dimension of $\mathcal{W}$ is greater than that of $\mathcal{Z}$, then $\|\sin\angle(\mathcal{Z},\mathcal{W})\|=1$, although generally $\|\sin\angle(\mathcal{W},\mathcal{Z})\|<1$. In this paper, we will use the $F$-norm distance

which equals the square root of the squares sum of sines of all the canonical angles between the two subspaces. Correspondingly, we define $\|\tan\angle(\mathcal{W},\mathcal{Z})\|_{F}$ to be the square root of the squares sum of the tangents of all the canonical angles between the subspaces. These tangents are the generalized singular values of the matrix pair $\{W_{\perp}^{H}Z,W^{H}Z\}$.

Notice that for $\mathcal{W}$ and $\mathcal{Z}$ with equal dimensions, if $W^{H}Z$ is nonsingular then the tangents of canonical angles are the singular values of $W_{\perp}^{H}Z(W^{H}Z)^{-1}$, so that

which holds when the $F$-norm is replaced by the 2-norm. We remark that the inverse in the above cannot be replaced by the pseudoinverse $\dagger$ when $W^{H}Z$ is singular. In fact, by definition, $\|\tan\angle(\mathcal{W},\mathcal{Z})\|_{F}$ is infinite when $W^{H}Z$ is singular, and the generalized singular values of $\{W_{\perp}^{H}Z,W^{H}Z\}$ are not the singular values of $W_{\perp}^{H}Z(W^{H}Z)^{{\dagger}}$ in this case.

Note that the real and imaginary parts of the eigenvectors $x_{\pm j}$ of $A$ can be interchanged since $\mathrm{i}x_{\pm j}$ are the eigenvectors of $A$ and its left and right singular vectors $u_{j}$ and $v_{j}$ are the right and left ones, too (cf. (2.4)). As we have pointed out, any pair of approximate left and right singular vectors extracted from the biorthogonal $\mathcal{U}_{m}$ and $\mathcal{V}_{m}$ are mutually orthogonal, using which we construct the real and imaginary parts of an approximate eigenvector of $A$. As a consequence, in the SVD context of the skew-symmetric $A$, because of these properties, when analyzing the convergence of the SSLBD method, we consider the orthogonal direct sum $\mathcal{U}_{m}\oplus\mathcal{V}_{m}$ as a whole, and estimate the distance between a desired two dimensional eigenspace $\mathcal{X}_{j}$ of $A$ and the $2m$-dimensional subspace $\mathcal{U}_{m}\oplus\mathcal{V}_{m}$ for $j$ small. We remind that their dimensions are unequal for $m>1$. Using these distances and their estimates, we can establish a priori error bounds for approximate eigenvalues and approximate eigenspaces obtained by the SSLBD method, showing how fast they converge when computing largest eigenvalues in magnitudes and the associated eigenvectors.

In terms of the definition of $\|\tan\angle(\mathcal{W},\mathcal{Z})\|_{F}$, we present the following estimate for $\|\tan\angle(\mathcal{U}_{m}\oplus\mathcal{V}_{m},\mathcal{X}_{j})\|_{F}$.

For a fixed $1\leq j<m$, denote by $\widehat{\Lambda}_{j}=\mathop{\operator@font diag}\nolimits\{\Lambda_{1},\dots,\Lambda_{j-1},\Lambda_{j+1},\dots,\Lambda_{\ell}\}$ and $\widehat{X}_{j}=\![X_{1},\dots,X_{j-1},X_{j+1},\dots,X_{\ell}]$ which delete $\Lambda_{j}$ and $X_{j}$ from $\widehat{X}$ and $\widehat{\Lambda}$ in (4.2), respectively. Then relation (4.2) can be written as

Notice that $Y=[p_{1},q_{1}]$ is orthonormal by Theorem 3.1, and $\mathcal{Y}=\mathcal{R}(Y)=\mathcal{U}_{1}\oplus\mathcal{V}_{1}$. Since $\widehat{X}$ in (4.2) is orthogonal, there exists an orthonormal matrix $G=[G_{1}^{T},\dots,G_{\ell}^{T}]^{T}$ with $G_{1},\dots,G_{\ell}\in\mathbb{R}^{2\times 2}$ such that

For an arbitrary polynomial $\rho(\cdot)$ in $\mathcal{P}^{m-1}$, the set of polynomials of degree no more than $m-1$ that satisfies $\rho(-\sigma_{j}^{2})\neq 0$, write $Y_{\rho}=\rho(A^{2})Y$ and $\mathcal{Y}_{\rho}=\mathcal{R}(Y_{\rho})$. Then

From (4.1), it is known that $\Lambda_{i}^{2}=-\sigma_{i}^{2}I_{2}$, $i=1,\dots,\ell$.

Since $G_{j}=X_{j}^{H}Y$ is nonsingular, $\rho(-\sigma_{j}^{2})\neq 0$, and the columns of $Y_{\rho}(Y_{\rho}^{T}Y_{\rho})^{-1/2}$ form an orthonormal basis of $\mathcal{Y}_{\rho}$, by definition (4.4) and relation (4.9), we obtain