# Two harmonic Jacobi–Davidson methods for computing a partial generalized singular value decomposition of a large matrix pair

Two harmonic extraction based Jacobi–Davidson (JD) type algorithms are proposed to compute a partial generalized singular value decomposition (GSVD) of a large regular matrix pair. They are called cross product-free (CPF) and inverse-free (IF) harmonic JDGSVD algorithms, abbreviated as CPF-HJDGSVD and IF-HJDGSVD, respectively. Compared with the standard extraction based JDGSVD algorithm, the harmonic extraction based algorithms converge more regularly and suit better for computing GSVD components corresponding to interior generalized singular values. Thick-restart CPF-HJDGSVD and IF-HJDGSVD algorithms with some deflation and purgation techniques are developed to compute more than one GSVD components. Numerical experiments confirm the superiority of CPF-HJDGSVD and IF-HJDGSVD to the standard extraction based JDGSVD algorithm.

There are no comments yet.

## Authors

• 3 publications
• 8 publications
04/29/2020

### A cross-product free Jacobi-Davidson type method for computing a partial generalized singular value decomposition (GSVD) of a large matrix pair

A Cross-Product Free (CPF) Jacobi-Davidson (JD) type method is proposed ...
07/24/2019

### On choices of formulations of computing the generalized singular value decomposition of a matrix pair

For the computation of the generalized singular value decomposition (GSV...
01/10/2020

### The joint bidiagonalization process with partial reorthogonalization

The joint bidiagonalization(JBD) process is a useful algorithm for the c...
11/06/2020

### Efficient Robust Watermarking Based on Quaternion Singular Value Decomposition and Coefficient Pair Selection

Quaternion singular value decomposition (QSVD) is a robust technique of ...
04/05/2020

### New Formulation and Computation for Generalized Singular Values of Grassman Matrix Pair

In this paper, we derive new model formulations for computing generalize...
02/12/2020

### Towards a more robust algorithm for computing the restricted singular value decomposition

A new algorithm to compute the restricted singular value decomposition o...
05/04/2020

### The Multi-Symplectic Lanczos Algorithm and Its Applications to Color Image Processing

Low-rank approximations of original samples are playing more and more an...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

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 [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 , and , , where the superscript denotes the transpose. Then the GSVD of is

 (1.1) {UTAX=ΣA=\diag{C,0l1,q1,Iq2},VTBX=ΣB=\diag{S,Iq1,0l2,q2},

where is nonsingular, and are orthogonal, and the diagonal matrices and satisfy

 0<αi,βi<1andα2i+β2i=1,i=1,…,q

with . Here, and are the zero matrices and identity matrices of order , respectively; see [28]. The GSVD part in (1.1) that corresponds to and can be written as

 (1.2)

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 :

 (1.3) |σ1−τ|≤⋯≤|σℓ−τ|<|σℓ+1−τ|≤⋯≤|σq−τ|.

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 [37] 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 [26] 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 [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 least squares problem with the coefficient matrix 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 with being the machine precision.

Hochstenbach [12] 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 [16]

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 [16]. The results in [16] 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 of

and 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 . 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 [31] 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 [38] exploits the Kronecker canonical form of a regular matrix pair [32] 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

 [uTj,xTj/βj,√σjuTj,√σjvTj]T,[−uTj,−xTj/βj,√σjuTj,√σjvTj]T, [−iuTj,ixTj/βj,√σjuTj,−√σjvTj]T,[iuTj,ixTj/βj,−√σjuTj,−√σjvTj]T.

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 [15] 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 [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 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 [39], where the dimension of the correction equations is no more than half of the dimension 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 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 [15] 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 [30] 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 [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 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 [15] 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

 (2.1) U=AXandV=BX

as the two left searching subspaces, from which approximations to and are computed. It is proved in [15] 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 [15] 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 :

 (2.2) A˜X=˜URAandB˜X=˜VRB,

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 [15] seeks for positive pairs with , normalized vectors , , and vectors that satisfy the Galerkin type conditions:

 (2.3) ⎧⎪⎨⎪⎩A~x−~α~u⊥U,B~x−~β~v⊥V,~βAT~u−~αBT~v⊥X.

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.

It follows from the thin QR factorizations (2.2) of and that and . Write , and . Then (2.3) becomes

 (2.4) RA~d=~α~e,RB~d=~β~f,~βRTA~e=~αRTB~f,

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

 (~α,~β,~u,~v,~x)=(~α,~β,˜U~e,˜V~f,˜X~d)

as an approximation to of . It is straightforward from (2.3) that , and

 (ATA−~θ2 BTB)~x⊥X.

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

 (2.5) r=r(~α,~β,~u,~v,~x)=~βAT~u−αBT~v.

Obviously, is an exact GSVD component of if and only if . The approximate GSVD component is claimed to have converged if

 (2.6) ∥r∥≤(~β∥A∥1+~α∥B∥1)⋅tol,

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

 (2.7) ~y:=(ATA+BTB)~x=~αAT~u+~βBT~v

that satisfies , we first solve the correction equation

 (2.8) (I−~y~xT)(ATA−ρ2BTB)(I−~x~yT)t=−r

with the fixed for until

 (2.9) ∥r∥≤(~β∥A∥1+~α∥B∥1)⋅fixtol

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

 (2.10) ˜Xnew=[˜X, x+]withx+=(I−˜X˜XT)t∥(I−˜X˜XT)t∥,

where is called an expansion vector. We then compute the orthonormal bases and of the expanded left searching subspaces

 Unew=AXnew=span{˜U,Ax+},Vnew=BXnew=span{˜V,Bx+}

by efficiently updating the thin QR factorizations of and , respectively, where

 ˜Unew=[˜U,u+],RA,new=[RArAγA],˜Vnew=[˜V,v+],RB,new=[RBrBγB]

with

 rA=˜UTAx+,γA=∥Ax+−˜UrA∥,u+=Ax+−˜UrAγA, rB=˜VTBx+,γB=∥Bx+−˜VrB∥,v+=Bx+−˜VrBγA.

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 [15], 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 [15] and stop the inner iterations when the inner relative residual norm satisfies

 (2.11) ∥rin∥≤min{2c~ε,0.01},

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

 (3.1) BTB=LLT

is the Cholesky factorization of with being nonsingular and lower triangular, and define the matrix

 (3.2) ˇA=AL−T.

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.

Let be a GSVD component of the regular matrix pair and . Assume that has full column rank and has the Cholesky factorization (3.1), and let be defined by (3.2) and the vector

 (3.3) z∗=1β∗LTx∗.

Then is a singular triplet of :

 (3.4) ˇAz∗=σ∗u∗andˇATu∗=σ∗z∗.
###### Proof.

It follows from the GSVD (1.2) of that with , meaning that . Making use of (3.1), we have

 ∥z∗∥=1β∗∥LTx∗∥=1β∗∥Bx∗∥=1.

By the definitions (3.2) and (3.3) of and , from we obtain

 ˇAz∗=1β∗AL−TLTx∗=1β∗Ax∗=α∗β∗u∗=σ∗u∗,

that is, the first relation in (3.4) holds. From the GSVD (1.2), it is straightforward that . Making use of this relation and (3.1) gives

 ˇATu∗=L−1ATu∗=σ∗β∗L−1BTBx∗=σ∗β∗LTx∗=σ∗z∗,

which proves the second relation in (3.4). ∎

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

 (3.5) [0ˇATˇA0][ˇzˇu]−ϕ[ˇzˇu] ⊥ ([0ˇATˇA0]−τI)R([˜Z˜U]).

This is the harmonic extraction approach for the eigenvalue problem of the augmented matrix

 [0ˇATˇA0]

for the given target [31, 33], where is a harmonic Ritz value and is the harmonic Ritz vector with respect to the searching subspace

 R([˜Z˜U]).

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

 [˜ZT˜UT][−τIˇATˇA−τI][−ϕIˇATˇA−ϕI][˜Z˜U][ˇdˇe]=0.

Decompose , and rearrange the above equation. Then we obtain the generalized eigenvalue problem of a -by- matrix pair:

 (3.6) [˜ZTˇATˇA˜Z+τ2˜ZT˜Z−2τ˜ZTˇAT˜U−2τ˜UTˇA˜Z˜UTˇAˇAT˜U+τ2I][ˇdˇe]=(ϕ−τ)[−τ˜ZT˜Z˜ZTˇAT˜U˜UTˇA˜Z−τI][ˇdˇe].

By (3.2), and the thin QR factorization of in (2.2), we have

 ˇA˜Z=A˜X=˜URA,

showing that

 ˜ZTˇATˇA˜Z=RTARA%and˜ZTˇAT˜U=RTA.

Moreover, exploiting the Cholesky factorization (3.1) of and the thin QR factorization of in (2.2), we obtain

 ˜ZT˜Z = ˜XTLLT˜X=˜XTBTB˜X=RTBRB, ˜UTˇAˇAT˜U = ˜UTA(LLT)−1AT˜U=˜UTA(BTB)−1AT˜U.

Substituting these two relations into (3.6) yields

 (3.7)

For the brevity of presentation, we will denote the symmetric matrices

 HA,B†=˜UTA(BTB)−1AT˜U

and

 (3.8) Gc=[−τRTBRBRTARA−τI],Hc=[RTARA+τ2RTBRB−2τRTA−2τRAHA,B†+τ2I].

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

 (3.9) (ϕ,ˇu,ˇz)=(τ+1μ,˜Uˇe∥ˇe∥,˜Zˇd∥˜Zˇd∥).

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

 (3.10) ˇx=1ˇδ˜Xˇd,

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

 1=1ˇδ2ˇdT˜XT(ATA+BTB)˜Xˇd=1ˇδ2ˇdT(RTARA+RTBRB)ˇd,

from which it follows that

 (3.11) ˇδ=√∥RAˇd∥2+∥RBˇd∥2.

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

 (3.12) ˇu=Aˇx∥Aˇx∥andˇv=Bˇx∥Bˇx∥

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

 Aˇx = 1ˇδA˜Xˇd=1ˇδ˜URAˇd=1ˇδ˜Uˇe, Bˇx = 1ˇδB˜Xˇd=1ˇδ˜VRBˇd=1ˇδ˜Vˇf.

Using them, we can efficiently compute the approximate generalized singular vectors

 (3.13) ˇu=Aˇx∥Aˇx∥=˜Uˇe∥ˇe∥andˇv=Bˇx∥Bˇx∥=˜Vˇf∥ˇf∥<