Two harmonic Jacobi–Davidson methods for computing a partial generalized singular value decomposition of a large matrix pair^††thanks: Supported by the National Natural Science Foundation of China (No. 12171273).

For a pair of large and possibly sparse matrices $A\in\mathbb{R}^{m\times n}$ and $B\in\mathbb{R}^{p\times n}$, the matrix pair $(A,B)$ is called regular if $\mathcal{N}(A)\cap\mathcal{N}(B)=\{\bm{0}\}$, i.e., $\mathrm{rank}\left(\begin{bmatrix}\begin{smallmatrix}A\\ B\end{smallmatrix}\end{bmatrix}\right)=n$, where $\mathcal{N}(A)$ and $\mathcal{N}(B)$ denote the null spaces of $A$ and $B$, respectively. The generalized singular value decomposition (GSVD) of $(A,B)$ was introduced by Van Loan [34] and developed by Paige and Saunders [28]. Since then, GSVD has become a standard matrix decomposition and has been widely used [2, 3, 4, 9, 10, 25]. Let $q_{1}=\dim(\mathcal{N}(A))$, $q_{2}=\dim(\mathcal{N}(B))$ and $l_{1}=\dim(\mathcal{N}(A^{T}))$, $l_{2}=\dim(\mathcal{N}(B^{T}))$, where the superscript $T$ denotes the transpose. Then the GSVD of $(A,B)$ is

where $X=[X_{q},X_{q_{1}},X_{q_{2}}]$ is nonsingular, $U=[U_{q},U_{l_{1}},U_{q_{2}}]$ and $V=[V_{q},V_{q_{1}},V_{l_{2}}]$ are orthogonal, and the diagonal matrices $C=\mathop{\operator@font diag}\nolimits\{\alpha_{1},\dots,\alpha_{q}\}$ and $S=\mathop{\operator@font diag}\nolimits\{\beta_{1},\dots,\beta_{q}\}$ satisfy

with $q=n-q_{1}-q_{2}$. Here, $\mathbf{0}_{l_{i},q_{i}}$ and $I_{q_{i}},i=1,2,$ are the $l_{i}\times q_{i}$ zero matrices and identity matrices of order $q_{i}$, respectively; see [28]. The GSVD part in (1.1) that corresponds to $\alpha_{i}$ and $\beta_{i}$ can be written as

where $x_{i}$ is the $i$th column of $X_{q}$ and the unit-length vectors $u_{i}$ and $v_{i}$ are the $i$th columns of $U_{q}$ and $V_{q}$, respectively. The quintuples $(\alpha_{i},\beta_{i},u_{i},v_{i},x_{i})$, $i=1,\dots,q$ are called nontrivial GSVD components of $(A,B)$. Particularly, the numbers $\sigma_{i}=\frac{\alpha_{i}}{\beta_{i}}$ or the pairs $(\alpha_{i},\beta_{i})$ are called the nontrivial generalized singular values, and $u_{i},v_{i}$ and $x_{i}$ are the corresponding left and right generalized singular vectors, respectively, $i=1,\dots,q$.

For a given target $\tau>0$, we assume that all the nontrivial generalized singular values of $(A,B)$ are labeled by their distances from $\tau$:

We are interested in computing the GSVD components $(\alpha_{i},\beta_{i},u_{i},v_{i},x_{i})$ corresponding to the $\ell$ nontrivial generalized singular values $\sigma_{i}$ of $(A,B)$ closest to $\tau$. If $\tau$ is inside the nontrivial generalized singular spectrum of $(A,B)$, then $(\alpha_{i},\beta_{i},u_{i},v_{i},x_{i})$, $i=1,\dots,\ell$ are called interior GSVD components of $(A,B)$; otherwise, they are called the extreme, i.e., largest or smallest, ones. A large number of GSVD components, some of which are interior ones [5, 6, 7], are required in a variety of applications. Throughout this paper, we assume that $\tau$ is not equal to any generalized singular value of $(A,B)$.

Zha [37] proposes a joint bidiagonalization (JBD) method to compute extreme GSVD components of the large matrix pair $(A,B)$. The method is based on a JBD process that successively reduces $(A,B)$ to a sequence of upper bidiagonal pairs, from which approximate GSVD components are computed. Kilmer, Hansen and Espanol [26] have adapted the JBD process to the linear discrete ill-posed problem with general-form regularization and developed a JBD process that reduces $(A,B)$ to lower-upper bidiagonal forms. Jia and Yang [24] have developed a new JBD process based iterative algorithm for the ill-posed problem and considered the convergence of extreme generalized singular values. In the GSVD computation and the solution of discrete ill-posed problem, one needs to solve an $(m+p)\times n$ least squares problem with the coefficient matrix $[A^{T},B^{T}]^{T}$ at each step of the JBD process. Jia and Li [22] have recently considered the JBD process in finite precision and proposed a partial reorthogonalization strategy to maintain numerical semi-orthogonality among the generated basis vectors so as to avoid ghost approximate GSVD components, where the semi-orthogonality means that two unit-length vectors are numerically orthogonal to the level of $\epsilon_{\rm mach}^{1/2}$ with $\epsilon_{\rm mach}$ being the machine precision.

Hochstenbach [12] presents a Jacobi–Davidson (JD) GSVD (JDGSVD) method to compute a number of interior GSVD components of $(A,B)$ with $B$ of full column rank, where, at each step, an $(m+n)$-dimensional linear system, i.e., the correction equation, needs to be solved iteratively with low or modest accuracy; see [14, 15, 20, 21]. The lower $n$-dimensional and upper $m$-dimensional parts of the approximate solution are used to expand the right and one of the left searching subspaces, respectively. The JDGSVD method formulates the GSVD of $(A,B)$ as the equivalent generalized eigendecomposition of the augmented matrix pair $\left(\begin{bmatrix}\begin{smallmatrix}&A\\ A^{T}&\end{smallmatrix}\end{bmatrix},\begin{bmatrix}\begin{smallmatrix}I&\\ &B^{T}B\end{smallmatrix}\end{bmatrix}\right)$ for $B$ of full column rank, computes the relevant eigenpairs, and recovers the approximate GSVD components from the converged eigenpairs. The authors [16] have shown that the error of the computed eigenvector is bounded by the size of the perturbations times a multiple $\kappa(B^{T}B)=\kappa^{2}(B)$, where $\kappa(B)=\sigma_{\max}(B)/\sigma_{\min}(B)$ denotes the $2$-norm condition number of $B$ with $\sigma_{\max}(B)$ and $\sigma_{\min}(B)$ being the largest and smallest singular values of $B$, respectively. Consequently, with an ill-conditioned $B$, the computed GSVD components may have very poor accuracy, which has been numerically confirmed [16]. The results in [16] show that if $B$ is ill conditioned but $A$ has full column rank and is well conditioned then the JDGSVD method can be applied to the matrix pair $(\begin{bmatrix}\begin{smallmatrix}&B\\ B^{T}&\end{smallmatrix}\end{bmatrix},\begin{bmatrix}\begin{smallmatrix}I&\\ &A^{T}A\end{smallmatrix}\end{bmatrix})$ and computes the corresponding approximate GSVD components with high accuracy. Note that the two formulations require that $B$ and $A$ be rectangular or square, respectively. We should also realize that a reliable estimation of the condition numbers of $A$ and $B$ may be costly, so that it may be difficult to choose a proper formulation in applications.

Zwaan and Hochstenbach [39] present a generalized Davidson (GDGSVD) method and a multidirectional (MDGSVD) method to compute an extreme partial GSVD of $(A,B)$. These two methods involve no cross product matrices $A^{T}A$ and $B^{T}B$ or matrix-matrix products, and they apply the standard extraction approach, i.e., the Rayleigh–Ritz method [31] to $(A,B)$ directly and compute approximate GSVD components with respect to the given left and right searching subspaces, where the two left subspaces are formed by premultiplying the right one with $A$ and $B$, respectively. At iteration $k$ of the GDGSVD method, the right searching subspace is spanned by the $k$ residuals of the generalized Davidson method [1, Sec. 11.2.4 and Sec. 11.3.6] applied to the generalized eigenvalue problem of $(A^{T}A,B^{T}B)$; in the MDGSVD method, an inferior search direction is discarded by a truncation technique, so that the searching subspaces are improved. Zwaan [38] exploits the Kronecker canonical form of a regular matrix pair [32] and shows that the GSVD problem of $(A,B)$ can be formulated as a certain $(2m+p+n)\times(2m+p+n)$ generalized eigenvalue problem without involving any cross product or any other matrix-matrix product. Such formulation currently is mainly of theoretical value since the nontrivial eigenvalues and eigenvectors of the structured generalized eigenvalue problem are always complex: the generalized eigenvalues are the conjugate quaternions $(\sqrt{\sigma_{j}},-\sqrt{\sigma_{j}},\mathrm{i}\sqrt{\sigma_{j}},-\mathrm{i}\sqrt{\sigma_{j}})$ with $\mathrm{i}$ the imaginary unit, and the corresponding right generalized eigenvectors are

Clearly, the size of the generalized eigenvalue problem is much bigger than that of the GSVD of $(A,B)$. The conditioning of eigenvalues and eigenvectors of this problem is also unclear. In the meantime, no structure-preserving algorithm has been found for such kind of complicated structured generalized eigenvalue problem. Definitely, it will be extremely difficult and highly challenging to seek for a numerically stable structure-preserving algorithm for this problem.

The authors [15] have recently proposed a Cross Product-Free JDGSVD method, referred to as the CPF-JDGSVD method, to compute several GSVD components of $(A,B)$ corresponding to the generalized singular values closest to $\tau$. The CPF-JDGSVD method is cross products $A^{T}A$ and $B^{T}B$ free when constructing and expanding right and left searching subspaces; it premultiplies the right searching subspace by $A$ and $B$ to construct two left ones separately, and forms the orthonormal bases of those by computing two thin QR factorizations, as done in [39]. The resulting projected problem is the GSVD of a small matrix pair without involving any cross product or matrix-matrix product. Mathematically, the method implicitly deals with the equivalent generalized eigenvalue problem of $(A^{T}A,B^{T}B)$ without forming $A^{T}A$ or $B^{T}B$ explicitly. At the subspace expansion stage, an $n$-by-$n$ correction equation is approximately solved iteratively with low or modest accuracy, and the approximate solution is used to expand the searching subspaces. Therefore, the subspace expansion is fundamentally different from that used in [39], where the dimension $n$ of the correction equations is no more than half of the dimension $m+n$ of those in [12].

Just like the standard Rayleigh–Ritz method for the matrix eigenvalue problem and the singular value decomposition (SVD) problem, the CPF-JDGSVD method suits better for the computation of some extreme GSVD components, but it may encounter some serious difficulties for the computation of interior GSVD components. Remarkably, adapted from the standard extraction approach for the eigenvalue problem and SVD problem to the GSVD computation, an intrinsic shortcoming of a standard extraction based method is that it may be hard to pick up good approximate generalized singular values correctly even if the searching subspaces are sufficiently good. This potential disadvantage may make the resulting algorithm expand the subspaces along wrong directions and converge irregularly, as has been numerically observed in [15]. To this end, inspired by the harmonic extraction based methods that suit better for computing interior eigenpairs and SVD components [11, 13, 14, 17, 23, 21, 27], we will propose two harmonic extraction based JDGSVD methods that are particularly suitable for the computation of interior GSVD components. One method is cross products $A^{T}A$ and $B^{T}B$ free, and the other is inversions $(A^{T}A)^{-1}$ and $(B^{T}B)^{-1}$ free. As will be seen, the derivations of the two harmonic extraction methods are nontrivial, and they are subtle adaptations of the harmonic extraction for matrix eigenvalue and SVD problems. In the sequel, we will abbreviate Cross Product-Free and Inverse-Free Harmonic JDGSVD methods as CPF-HJDGSVD and IF-HJDGSVD, respectively.

We first focus on the case $\ell=1$ and propose our harmonic extraction based JDGSVD type methods. Then by introducing the deflation technique in [15] into the methods, we present the methods to compute more than one, i.e., $\ell>1$, GSVD components. To be practical, combining the thick-restart technique in [30] and some purgation approach, we develop thick-restart CPF-HJDGSVD and IF-HJDGSVD algorithms to compute the $\ell$ GSVD components associated with the generalized singular values of $(A,B)$ closest to $\tau$.

The rest of this paper is organized as follows. In Section 2, we briefly review the CPF-JDGSVD method proposed in [15]. In Section 3, we propose the CPF-HJDGSVD and IF-HJDGSVD methods. In Section 4, we develop thick-restart CPF-HJDGSVD and IF-HJDGSVD with deflation and purgation to compute $\ell$ GSVD components of $(A,B)$. In Section 5, we report numerical experiments to illustrate the performance of CPF-HJDGSVD and IF-HJDGSVD, make a comparison of them and CPF-JDGSVD, and show the superiority of the former two to the latter one. Finally, we conclude the paper in Section 6.

Throughout this paper, we denote by $\mathcal{R}(\cdot)$ the column space of a matrix, and by $\|\cdot\|$ and $\|\cdot\|_{1}$ the $2$- and $1$-norms of a matrix or vector, respectively. As in (1.1), we denote by $I_{i}$ and $\bm{0}_{i,j}$ the $i$-by-$i$ identity and $i$-by-$j$ zero matrices, respectively, with the subscripts $i$ and $j$ dropped whenever they are clear from the context.

We review the CPF-JDGSVD method in [15] for computing the GSVD component $(\alpha_{*},\beta_{*},u_{*},v_{*},x_{*}):=(\alpha_{1},\beta_{1},u_{1},v_{1},x_{1})$ of $(A,B)$. Assume that a $k$-dimensional right searching subspace $\mathcal{X}\subset\mathbb{R}^{n}$ is available, from which an approximation to $x_{*}$ is extracted. Then we construct

as the two left searching subspaces, from which approximations to $u_{*}$ and $v_{*}$ are computed. It is proved in [15] that the distance between $u_{*}$ and $\mathcal{U}$ (resp. $v_{*}$ and $\mathcal{V}$) is as small as that between $x_{*}$ and $\mathcal{X}$, provided that $\alpha_{*}$ (resp. $\beta_{*}$) is not very small. In other words, for the extreme and interior GSVD components, $\mathcal{U}$ and $\mathcal{V}$ constructed by (2.1) are as good as $\mathcal{X}$ provided that the desired generalized singular values $\sigma_{*}=\frac{\alpha_{*}}{\beta_{*}}$ are neither very small nor very small. It is also proved in [15] that $\mathcal{U}$ or $\mathcal{V}$ is as accurate as $\mathcal{X}$ for very large or small generalized singular values.

Assume that the columns of $\widetilde{X}\in\mathbb{R}^{n\times k}$ form an orthonormal basis of $\mathcal{X}$, and compute the thin QR factorizations of $A\widetilde{X}$ and $B\widetilde{X}$:

where $\widetilde{U}\in\mathbb{R}^{m\times k}$ and $\widetilde{V}\in\mathbb{R}^{p\times k}$ are orthonormal, and $R_{A}\in\mathbb{R}^{k\times k}$ and $R_{B}\in\mathbb{R}^{k\times k}$ are upper triangular. Then the columns of $\widetilde{U}$ and $\widetilde{V}$ are orthonormal bases of $\mathcal{U}$ and $\mathcal{V}$, respectively. With $\mathcal{X}$, $\mathcal{U}$, $\mathcal{V}$ and their orthonormal bases available, we can extract an approximation to the desired GSVD component $(\alpha_{*},\beta_{*},u_{*},v_{*},x_{*})$ of $(A,B)$ with respect to them. The standard extraction approach in [15] seeks for positive pairs $(\tilde{\alpha},\tilde{\beta})$ with $\tilde{\alpha}^{2}+\tilde{\beta}^{2}=1$, normalized vectors $\tilde{u}\in\mathcal{U}$, $\tilde{v}\in\mathcal{V}$, and vectors $\tilde{x}\in\mathcal{X}$ that satisfy the Galerkin type conditions:

Among $k$ pairs $(\tilde{\alpha},\tilde{\beta})$’s, select $\tilde{\theta}=\tilde{\alpha}/\tilde{\beta}$ closest to $\tau$, and take $(\tilde{\alpha},\tilde{\beta},\tilde{u},\tilde{v},\tilde{x})$ as an approximation to $(\alpha_{*},\beta_{*},u_{*},v_{*},x_{*})$. We call $(\tilde{\alpha},\tilde{\beta})$ or $\tilde{\theta}=\frac{\tilde{\alpha}}{\tilde{\beta}}$ a Ritz value and $\tilde{u}$, $\tilde{v}$ and $\tilde{x}$ the corresponding left and right Ritz vectors, respectively.

It follows from the thin QR factorizations (2.2) of $A\widetilde{X}$ and $B\widetilde{X}$ that $R_{A}=\widetilde{U}^{T}A\widetilde{X}$ and $R_{B}=\widetilde{V}^{T}A\widetilde{X}$. Write $\tilde{u}=\widetilde{U}\tilde{e}$, $\tilde{v}=\widetilde{V}\tilde{f}$ and $\tilde{x}=\widetilde{X}\tilde{d}$. Then (2.3) becomes

which is precisely the GSVD of the projected matrix pair $(R_{A},R_{B})$. Therefore, in the extraction phase, the standard extraction approach computes the GSVD of the $k$-by-$k$ matrix pair $(R_{A},R_{B})$, picks up the GSVD component $(\tilde{\alpha},\tilde{\beta},\tilde{e},\tilde{f},\tilde{d})$ with $\tilde{\theta}=\frac{\tilde{\alpha}}{\tilde{\beta}}$ being the generalized singular value of $(R_{A},R_{B})$ closest to the target $\tau$, and use

as an approximation to $(\alpha_{*},\beta_{*},u_{*},v_{*},x_{*})$ of $(A,B)$. It is straightforward from (2.3) that $A\tilde{x}=\tilde{\alpha}\tilde{u}$, $B\tilde{x}=\tilde{\beta}\tilde{v}$ and

That is, $(\tilde{\theta}^{2},\tilde{x})$ is a standard Ritz pair to the eigenpair $(\sigma_{*}^{2},x_{*})$ of the symmetric definite matrix pair $(A^{T}A,B^{T}B)$ with respect to the subspace $\mathcal{X}$. Because of this, we call $(\tilde{\alpha},\tilde{\beta},\tilde{u},\tilde{v},\tilde{x})$ a standard Ritz approximation in the GSVD context.

Since $A\tilde{x}=\tilde{\alpha}\tilde{u}$ and $B\tilde{x}=\tilde{\beta}\tilde{v}$, the residual of Ritz approximation $(\tilde{\alpha},\tilde{\beta},\tilde{u},\tilde{v},\tilde{x})$ is

Obviously, $(\tilde{\alpha},\tilde{\beta},\tilde{u},\tilde{v},\tilde{x})$ is an exact GSVD component of $(A,B)$ if and only if $\|r\|=0$. The approximate GSVD component $(\tilde{\alpha},\tilde{\beta},\tilde{u},\tilde{v},\tilde{x})$ is claimed to have converged if

where $tol>0$ is a user prescribed tolerance, and one then stops the iterations.

If $(\tilde{\alpha},\tilde{\beta},\tilde{u},\tilde{v},\tilde{x})$ has not yet converged, the CFP-JDGSVD method expands the right searching subspace $\mathcal{X}$ and constructs the corresponding left subspaces $\mathcal{U}$ and $\mathcal{V}$ by (2.1). Specifically, the CPF-JDGSVD seeks for an expansion vector $t$ in the following way: For the vector

that satisfies $\tilde{y}^{T}\tilde{x}=1$, we first solve the correction equation

with the fixed $\rho=\tau$ for $t\perp\tilde{y}$ until

for a user prescribed tolerance $fixtol>0$, say, $fixtol=10^{-4}$, and then solve the modified correction equation with the dynamic $\rho=\tilde{\alpha}/\tilde{\beta}$ for $t\perp\tilde{y}$. Note that $I-\tilde{y}\tilde{x}^{T}$ is an oblique projector onto the orthogonal complement of the subspace ${\rm span}\{x\}$.

With the solution $t$ of (2.8), we expand $\mathcal{X}$ to the new $(k+1)$-dimensional $\mathcal{X}_{\rm new}={\rm span}\{\widetilde{X},t\}$, whose orthonormal basis matrix is

where $x_{+}$ is called an expansion vector. We then compute the orthonormal bases $\widetilde{U}_{\mathrm{new}}$ and $\widetilde{V}_{\mathrm{new}}$ of the expanded left searching subspaces

by efficiently updating the thin QR factorizations of $A\widetilde{X}_{\mathrm{new}}=\widetilde{U}_{\mathrm{new}}R_{A,\mathrm{new}}$ and $B\widetilde{X}_{\mathrm{new}}=\widetilde{V}_{\mathrm{new}}R_{B,\mathrm{new}}$, respectively, where

with

CPF-JDGSVD then computes a new approximate GSVD component of $(A,B)$ with respect to $\mathcal{U}_{\rm new},\mathcal{V}_{\rm new}$ and $\mathcal{X}_{\rm new}$, and repeat the above process until the convergence criterion (2.6) is achieved. We call iterative solutions of (2.8) the inner iterations and the extractions of the approximate GSVD components with respect to $\mathcal{U}$, $\mathcal{V}$ and $\mathcal{X}$ the outer iterations.

As has been shown in [15], it suffices to iteratively solve the correction equations approximately with low or modest accuracy and uses an approximate solution to update $\mathcal{X}$ in the above way, in order that the resulting inexact CPF-JDGSVD method and its exact counterpart with the correction equations solved accurately behave similarly. Precisely, for the correction equation (2.8), we adopt the inner stopping criteria in [15] and stop the inner iterations when the inner relative residual norm $\|r_{in}\|$ satisfies

where $\tilde{\varepsilon}\in[10^{-4},10^{-3}]$ is a user prescribed parameter and $c$ is a constant depending on $\rho$ and the current approximate generalized singular values.

We shall make use of the principle of the harmonic extraction [31, 33] to propose the CPF-harmonic and IF-harmonic extraction based JDGSVD methods in Section 3.1 and Section 3.2, respectively. They compute new approximate GSVD components of $(A,B)$ with respect to the given left and right searching subspaces $\mathcal{U}$, $\mathcal{V}$ and $\mathcal{X}$, and suit better for the computation of interior GSVD components.