For a pair of large and possibly sparse matrices and , the matrix pair is called regular if , i.e., , where and denote the null spaces of and , respectively. The generalized singular value decomposition (GSVD) of was introduced by Van Loan  and developed by Paige and Saunders . Since then, GSVD has become a standard matrix decomposition and has been widely used [2, 3, 4, 9, 10, 25]. Let , and , , where the superscript denotes the transpose. Then the GSVD of is
where is nonsingular, and are orthogonal, and the diagonal matrices and satisfy
where is the th column of and the unit-length vectors and are the th columns of and , respectively. The quintuples , are called nontrivial GSVD components of . Particularly, the numbers or the pairs are called the nontrivial generalized singular values, and and are the corresponding left and right generalized singular vectors, respectively, .
For a given target , we assume that all the nontrivial generalized singular values of are labeled by their distances from :
We are interested in computing the GSVD components corresponding to the nontrivial generalized singular values of closest to . If is inside the nontrivial generalized singular spectrum of , then , are called interior GSVD components of ; 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 is not equal to any generalized singular value of .
Zha  proposes a joint bidiagonalization (JBD) method to compute extreme GSVD components of the large matrix pair . The method is based on a JBD process that successively reduces to a sequence of upper bidiagonal pairs, from which approximate GSVD components are computed. Kilmer, Hansen and Espanol  have adapted the JBD process to the linear discrete ill-posed problem with general-form regularization and developed a JBD process that reduces to lower-upper bidiagonal forms. Jia and Yang  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 least squares problem with the coefficient matrix at each step of the JBD process. Jia and Li  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 with being the machine precision.
Hochstenbach  presents a Jacobi–Davidson (JD) GSVD (JDGSVD) method to compute a number of interior GSVD components of with of full column rank, where, at each step, an -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 -dimensional and upper -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 as the equivalent generalized eigendecomposition of the augmented matrix pair for of full column rank, computes the relevant eigenpairs, and recovers the approximate GSVD components from the converged eigenpairs. The authors 
have shown that the error of the computed eigenvector is bounded by the size of the perturbations times a multiple, where denotes the -norm condition number of with and being the largest and smallest singular values of , respectively. Consequently, with an ill-conditioned , the computed GSVD components may have very poor accuracy, which has been numerically confirmed . The results in  show that if is ill conditioned but has full column rank and is well conditioned then the JDGSVD method can be applied to the matrix pair and computes the corresponding approximate GSVD components with high accuracy. Note that the two formulations require that and
be rectangular or square, respectively. We should also realize that a reliable estimation of the condition numbers ofand may be costly, so that it may be difficult to choose a proper formulation in applications.
Zwaan and Hochstenbach  present a generalized Davidson (GDGSVD) method and a multidirectional (MDGSVD) method to compute an extreme partial GSVD of . These two methods involve no cross product matrices and or matrix-matrix products, and they apply the standard extraction approach, i.e., the Rayleigh–Ritz method  to 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 and , respectively. At iteration of the GDGSVD method, the right searching subspace is spanned by the residuals of the generalized Davidson method [1, Sec. 11.2.4 and Sec. 11.3.6]
applied to the generalized eigenvalue problem of; in the MDGSVD method, an inferior search direction is discarded by a truncation technique, so that the searching subspaces are improved. Zwaan  exploits the Kronecker canonical form of a regular matrix pair  and shows that the GSVD problem of can be formulated as a certain 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 with 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 . 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  have recently proposed a Cross Product-Free JDGSVD method, referred to as the CPF-JDGSVD method, to compute several GSVD components of corresponding to the generalized singular values closest to . The CPF-JDGSVD method is cross products and free when constructing and expanding right and left searching subspaces; it premultiplies the right searching subspace by and to construct two left ones separately, and forms the orthonormal bases of those by computing two thin QR factorizations, as done in . 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 without forming or explicitly. At the subspace expansion stage, an -by- 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 , where the dimension of the correction equations is no more than half of the dimension of those in .
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 . 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 and free, and the other is inversions and 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 and propose our harmonic extraction based JDGSVD type methods. Then by introducing the deflation technique in  into the methods, we present the methods to compute more than one, i.e., , GSVD components. To be practical, combining the thick-restart technique in  and some purgation approach, we develop thick-restart CPF-HJDGSVD and IF-HJDGSVD algorithms to compute the GSVD components associated with the generalized singular values of closest to .
The rest of this paper is organized as follows. In Section 2, we briefly review the CPF-JDGSVD method proposed in . 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 GSVD components of . 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 the column space of a matrix, and by and the - and -norms of a matrix or vector, respectively. As in (1.1), we denote by and the -by- identity and -by- zero matrices, respectively, with the subscripts and dropped whenever they are clear from the context.
2 The standard extraction based JDGSVD method
We review the CPF-JDGSVD method in  for computing the GSVD component of . Assume that a -dimensional right searching subspace is available, from which an approximation to is extracted. Then we construct
as the two left searching subspaces, from which approximations to and are computed. It is proved in  that the distance between and (resp. and ) is as small as that between and , provided that (resp. ) is not very small. In other words, for the extreme and interior GSVD components, and constructed by (2.1) are as good as provided that the desired generalized singular values are neither very small nor very small. It is also proved in  that or is as accurate as for very large or small generalized singular values.
Assume that the columns of form an orthonormal basis of , and compute the thin QR factorizations of and :
where and are orthonormal, and and are upper triangular. Then the columns of and are orthonormal bases of and , respectively. With , , and their orthonormal bases available, we can extract an approximation to the desired GSVD component of with respect to them. The standard extraction approach in  seeks for positive pairs with , normalized vectors , , and vectors that satisfy the Galerkin type conditions:
Among pairs ’s, select closest to , and take as an approximation to . We call or a Ritz value and , and the corresponding left and right Ritz vectors, respectively.
which is precisely the GSVD of the projected matrix pair . Therefore, in the extraction phase, the standard extraction approach computes the GSVD of the -by- matrix pair , picks up the GSVD component with being the generalized singular value of closest to the target , and use
as an approximation to of . It is straightforward from (2.3) that , and
That is, is a standard Ritz pair to the eigenpair of the symmetric definite matrix pair with respect to the subspace . Because of this, we call a standard Ritz approximation in the GSVD context.
Since and , the residual of Ritz approximation is
Obviously, is an exact GSVD component of if and only if . The approximate GSVD component is claimed to have converged if
where is a user prescribed tolerance, and one then stops the iterations.
If has not yet converged, the CFP-JDGSVD method expands the right searching subspace and constructs the corresponding left subspaces and by (2.1). Specifically, the CPF-JDGSVD seeks for an expansion vector in the following way: For the vector
that satisfies , we first solve the correction equation
with the fixed for until
for a user prescribed tolerance , say, , and then solve the modified correction equation with the dynamic for . Note that is an oblique projector onto the orthogonal complement of the subspace .
With the solution of (2.8), we expand to the new -dimensional , whose orthonormal basis matrix is
where is called an expansion vector. We then compute the orthonormal bases and of the expanded left searching subspaces
by efficiently updating the thin QR factorizations of and , respectively, where
CPF-JDGSVD then computes a new approximate GSVD component of with respect to and , 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 , and the outer iterations.
As has been shown in , it suffices to iteratively solve the correction equations approximately with low or modest accuracy and uses an approximate solution to update 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  and stop the inner iterations when the inner relative residual norm satisfies
where is a user prescribed parameter and is a constant depending on and the current approximate generalized singular values.
3 The harmonic extraction based JDGSVD methods
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 with respect to the given left and right searching subspaces , and , and suit better for the computation of interior GSVD components.
3.1 The CPF-harmonic extraction approach
If has full column rank with some special, e.g., banded, structure, from which the inversion can be efficiently applied, we can propose our CPF-harmonic extraction approach to compute a desired approximate GSVD component as follows. For the purpose of derivation, assume that
is the Cholesky factorization of with being nonsingular and lower triangular, and define the matrix
We present the following result, which establishes the relationship between the GSVD of and the SVD of and will be used to propose the CPF-harmonic extraction approach.
Theorem 3.1 motivates us to propose our first harmonic extraction approach to compute the singular triplet of and then recover the desired GSVD component of .
Specifically, take the -dimensional and as the left and right searching subspaces for the left and right singular vectors and of , respectively. Then the columns of form a basis of . Mathematically, we seek for positive and vectors and such that
This is the harmonic extraction approach for the eigenvalue problem of the augmented matrix
We pick up the closest to as the approximation to and take the normalized and as approximations to and , respectively. We will show how to obtain an approximation to afterwards.
Write and with and . Then , and requirement (3.5) amounts to the equation
Decompose , and rearrange the above equation. Then we obtain the generalized eigenvalue problem of a -by- matrix pair:
Substituting these two relations into (3.6) yields
For the brevity of presentation, we will denote the symmetric matrices
In implementations, we compute the generalized eigendecomposition of the symmetric positive definite matrix pair and pick up the largest eigenvalue in magnitude and the corresponding unit-length eigenvector . Then the harmonic Ritz approximation to the desired singular triplet of is
Since is an approximation to the right singular vector of , from (3.3) the vector after some proper normalization is an approximation to the right generalized singular vector of , which we write as
where is a normalizing factor. It is natural to require that the approximate right singular vector be -norm normalized, i.e., , since the exact satisfies by (1.2). With this normalization, from (3.10), we have
from which it follows that
Note that the approximate left generalized singular vector defined by (3.9) is no longer collinear with , as opposed to the collinear and obtained by the standard extraction approach in Section 2. To this end, instead of in (3.9), we take new and defined by
as the harmonic Ritz approximations to and , which are colinear with and , respectively. Correspondingly, define and . Then by (3.11), the parameter in (3.10) becomes . Moreover, by definition (3.10) of and the thin QR factorizations of and in (2.2), we obtain
Using them, we can efficiently compute the approximate generalized singular vectors