# Characterizing GSVD by singular value expansion of linear operators and its
computation

Haibo Li School of Mathematics and Statistics, The University of Melbourne,
Parkville, VIC 3010, Australia<EMAIL_ADDRESS>
###### Abstract

The generalized singular value decomposition (GSVD) of a matrix pair
$\\{A,L\\}$ with $A\in\mathbb{R}^{m\times n}$ and $L\in\mathbb{R}^{p\times n}$
generalizes the singular value decomposition (SVD) of a single matrix. In this
paper, we provide a new understanding of GSVD from the viewpoint of SVD, based
on which we propose a new iterative method for computing nontrivial GSVD
components of a large-scale matrix pair. By introducing two linear operators
$\mathcal{A}$ and $\mathcal{L}$ induced by $\\{A,L\\}$ between two finite-
dimensional Hilbert spaces and applying the theory of singular value expansion
(SVE) for linear compact operators, we show that the GSVD of $\\{A,L\\}$ is
nothing but the SVEs of $\mathcal{A}$ and $\mathcal{L}$. This result
characterizes completely the structure of GSVD for any matrix pair with the
same number of columns. As a direct application of this result, we generalize
the standard Golub-Kahan bidiagonalization (GKB) that is a basic routine for
large-scale SVD computation such that the resulting generalized GKB (gGKB)
process can be used to approximate nontrivial extreme GSVD components of
$\\{A,L\\}$, which is named the gGKB_GSVD algorithm. We use the GSVD of
$\\{A,L\\}$ to study several basic properties of gGKB and also provide
preliminary results about convergence and accuracy of gGKB_GSVD for GSVD
computation. Numerical experiments are presented to demonstrate the
effectiveness of this method.

#### Keywords

GSVD, linear operator, singular value expansion, generalized Golub-Kahan
bidiagonalization, Krylov subspace

#### MSC codes

15A22, 47A05, 65F99

##  1 Introduction

The generalized singular value decomposition (GSVD) of a matrix pair is an
extension of the singular value decomposition (SVD) of a single matrix. First
introduced by Van Loan [52] and further developed by many others [41, 49, 53],
now the GSVD has been a standard matrix decomposition [8, 19]. The GSVD
provides an important mathematical tool for analyzing relationships between
two sets of variables or matrices, which is particularly useful in various
applications, including signal processing [37, 48], statistics [43, 39],
computational biology [1], and many others [7, 17, 27, 29, 22].

Let $I_{k}$ denote the identity matrix of order $k$ and $\mathbf{0}$ denote
the zero matrix or vector with dimensions clarified by the context. For any
two matrices with the same number of columns, the general-form GSVD is stated
as follows [41]:

######  Theorem 1.1 (GSVD)

Let $A\in\mathbb{R}^{m\times n}$ and $L\in\mathbb{R}^{p\times n}$ with
$\mathrm{rank}((A^{\top},L^{\top})^{\top})=r$. Then the GSVD of $\\{A,L\\}$ is

$A=P_{A}C_{A}X^{-1},\ \ \ L=P_{L}S_{L}X^{-1},$ (1.1a) with
$C_{A}=\vtop{\hbox{$\left(\vbox{\vbox{\halign{ \hfil$#$\hfil\kern
2\p@\kern\@leftbwd& \thinspace\hfil$#$\hfil&&
\quad\hfil$#$\hfil\cr\Sigma_{A}$\hfil\kern 2.0pt\kern
6.66669pt&\mathbf{0}&m\crcr\kern 2.0pt\cr r$\hfil\kern 2.0pt\kern
6.66669pt&n-r\cr} }\kern-27.45218pt}\kern 0.0pt\right)\kern-6.66669pt\kern
0.0pt$}\kern 27.45218pt}\ ,\ \ \ S_{L}=\vtop{\hbox{$\left(\vbox{\vbox{\halign{
\hfil$#$\hfil\kern 2\p@\kern\@leftbwd& \thinspace\hfil$#$\hfil&&
\quad\hfil$#$\hfil\cr\Sigma_{L}$\hfil\kern 2.0pt\kern
6.66669pt&\mathbf{0}&p\crcr\kern 2.0pt\cr r$\hfil\kern 2.0pt\kern
6.66669pt&n-r\cr} }\kern-27.48329pt}\kern 0.0pt\right)\kern-6.66669pt\kern
0.0pt$}\kern 27.48329pt}$ (1.1b) and
$\Sigma_{A}=\vtop{\hbox{$\left(\vbox{\vbox{\halign{ \hfil$#$\hfil\kern
2\p@\kern\@leftbwd& \thinspace\hfil$#$\hfil&& \quad\hfil$#$\hfil\cr
I_{q_{1}}$\hfil\kern 2.0pt\kern 6.66669pt&&&q_{1}\crcr\kern
2.0pt\cr$\hfil\kern 2.0pt\kern 6.66669pt&C_{q_{2}}&&q_{2}\cr$\hfil\kern
2.0pt\kern 6.66669pt&&\mathbf{0}&m-q_{1}-q_{2}\cr q_{1}$\hfil\kern 2.0pt\kern
6.66669pt&q_{2}&q_{3}\crcr} }\kern-47.29436pt}\kern
0.0pt\right)\kern-6.66669pt\kern 0.0pt$}\kern 47.29436pt},\ \ \
\Sigma_{L}=\vtop{\hbox{$\left(\vbox{\vbox{\halign{ \hfil$#$\hfil\kern
2\p@\kern\@leftbwd& \thinspace\hfil$#$\hfil&&
\quad\hfil$#$\hfil\cr\mathbf{0}$\hfil\kern 2.0pt\kern
6.66669pt&&&p-r+q_{1}\crcr\kern 2.0pt\cr$\hfil\kern 2.0pt\kern
6.66669pt&S_{q_{2}}&&q_{2}\cr$\hfil\kern 2.0pt\kern
6.66669pt&&I_{q_{3}}&q_{3}\cr q_{1}$\hfil\kern 2.0pt\kern
6.66669pt&q_{2}&q_{3}\crcr} }\kern-47.29436pt}\kern
0.0pt\right)\kern-6.66669pt\kern 0.0pt$}\kern 47.29436pt},$ (1.1c)

where $q_{1}+q_{2}+q_{3}=r$, and $P_{A}\in\mathbb{R}^{m\times m}$,
$P_{L}\in\mathbb{R}^{p\times p}$ are orthogonal, $X\in\mathbb{R}^{n\times n}$
is invertible, and $\Sigma_{A}^{T}\Sigma_{A}+\Sigma_{L}^{T}\Sigma_{L}=I_{r}$.
The values of $q_{1}$, $q_{2}$ and $q_{3}$ are defined internally by the
matrices $A$ and $L$.

If $r=n$, then $\\{A,L\\}$ is called a regular matrix pair. Discussions about
GSVD for regular and nonregular matrix pairs can be found in [35, 49] and [41,
50, 38], respectively. Write
$C_{q_{2}}=\mbox{diag}(c_{q_{1}+1},\dots,c_{q_{1}+q_{2}})$ with
$1>c_{q_{1}+1}\geq\cdots\geq c_{q_{1}+q_{2}}>0$ and
$S_{q_{2}}=\mbox{diag}(s_{q_{1}+1},\dots,s_{q_{1}+q_{2}})$ with
$0<s_{q_{1}+1}\leq\cdots\leq s_{q_{1}+q_{2}}<1$. Let
$c_{1}=\cdots=c_{q_{1}}=1,\ c_{q_{1}+q_{2}+1}=\cdots=c_{r}=0$ and
$s_{1}=\cdots=s_{q_{1}}=0,\ s_{q_{1}+q_{2}+1}=\cdots=s_{r}=1$. Write
$X=(x_{1},\ldots,x_{n})$, $P_{A}=(p_{A,1},\ldots,p_{A,m})$ and
$P_{L}=(p_{L,1},\ldots,p_{L,p})$. We call the tuple
$(c_{i},s_{i},x_{i},p_{A,i},p_{L,i})$ the $i$-th nontrivial GSVD components,
and the $i$-th largest generalized singular value is $\gamma_{i}:=c_{i}/s_{i}$
satisfying $c_{i}^{2}+s_{i}^{2}=1$, where $1\leq i\leq r$. In this paper, we
consider the nontrivial GSVD components and their computations.

Despite its remarkable capabilities, computing the GSVD poses significant
challenges. Early computational approaches for the GSVD were built upon
adaptations of algorithms designed for the SVD; for small-scale matrices,
there are several such numerical algorithms for full GSVD computation [53, 40,
5]. Recent development on stable computation of the CS decomposition (CSD)
[51] provides another alternative for small-scale GSVD computation. For large
and sparse problems, obtaining the full GSVD may not be feasible, yet it is
often necessary to compute only a subset of GSVD components relevant to
practical applications. Typically, this refers to certain extreme GSVD
components, which are those with the largest or smallest corresponding
generalized singular values, or interior GSVD components, which are those
whose corresponding desired generalized singular values are closest to a
specified target.

For large-scale GSVD computation, the first step is usually tranforming it as
an equivalent generalized eigendecomposition (GED) problem [4] or CSD problem.
The joint bidiagonalization (JBD) method proposed by Zha [54] can be used to
compute a few extreme GSVD components, which is essentially an indirect
procedure for CSD of the Q-factor of a regular $\\{A,L\\}$ (i.e. the Q matrix
in the QR factorization). This method relies on a JBD process that iteratively
reduces $\\{A,L\\}$ to bidiagonal forms simultaneously. Jia and Li [25, 26]
made a detailed analysis for the numerical behavior of the JBD method and the
convergence behavior for extreme GSVD components in finite precision
arithmetic. They proposed the semi-orthogonalization strategy and design a
partial reorthogonalization procedure to maintain regular convergence of the
computed approximate GSVD components. Subsequently, Li [33] analyzed the
influence of computational errors resulting from inaccurate inner iterations
in JBD on the convergence and final accuracy of computed GSVD components and
proposed a modified version of the JBD method. Recently, Alvarruiz et al. [2]
developed a thick restart technique for JBD to compute a partial GSVD,
enabling the storage and computation cost further saved. On the other hand,
the Jacobi–Davidson type algorithms [21] are capable of computing a few
interior GSVD components. A representative work is the Jacobi–Davidson GSVD
(JDGSVD) method proposed by Hochstenbach [20], which formulates the GSVD of
$\\{A,L\\}$ as a GED problem of an augmented symmetric matrix pair. This
method is further analyzed and developed in several subsequent work; see e.g.
[45, 23, 24]. As the development of contour integration technique for
eigenvalue problems of finding interior eigenvalues [47, 44], recently a
contour integral-based algorithm has been adopted for interior GSVD components
computation [36].

Apart from regarding the GSVD as an equivalent CSD or GED, there is very
little work on understanding and analyzing GSVD from other perspectives. In
[10], the authors studied the GSVD using a variational formulation analogous
to that of the SVD, providing a new understanding of the generalized singular
vectors. Recently, by treating $(A^{\top},L^{\top})^{\top}$ (more precisely,
the range space of it) as a point in the real Grassmann manifold
$\mathrm{Gr}(m+p,r)$ — the manifold composed of all $r$-dimensional subspaces
of $\mathbb{R}^{m+p}$ — the authors in [15] interpreted a modified form of
GSVD as a coordinate representation of $(A^{\top},L^{\top})^{\top}$. For
developing practical GSVD algorithms, however, these new perspectives on GSVD
are not enough. It would be beneficial to understand the GSVD from the
viewpoint of SVD so that many existing algorithms for large-scale SVD are
possible to be adapted for large-scale GSVD computation. One well-known result
is that $\\{\gamma_{i}\\}$ are the singular values of $AL^{{\dagger}}$ if $L$
has full column rank [54], where “${\dagger}$” is the Moor-Penrose
pseudoinverse. But for noninvertible and nonsquare $L$, generally the GSVD of
$\\{A,L\\}$ is not related to the SVD of $AL^{{\dagger}}$; this issue becomes
much more complicated for nonregular matrix pairs.

In this paper, we provide a new understanding of GSVD from the viewpoint of
SVD. This new perspective relies on the theory of singular value expansion
(SVE) for linear compact operators [16, §2.2], which is essentially the SVD if
the compact operator is a matrix between two Euclidean spaces. Denote by
$\mathcal{R}(\cdot)$ and $\mathcal{N}(\cdot)$ the range space and null space
of a matrix, respectively. By defining the positive semidefinite matrix
$M=A^{\top}A+L^{\top}L$, we first investigate the structure of trivial and
nontrivial GSVD components $x_{i}$, showing that those trivial $\\{x_{i}\\}$
form a basis for $\mathcal{N}(M)$ and any nontrivial $x_{i}$ belongs to the
coset $\bar{x}_{i}+\mathcal{N}(M)$, where $\bar{x}_{i}\in\mathcal{R}(M)$ is a
nontrivial GSVD component. Then we consider the nontrivial
$x_{i}\in\mathcal{R}(M)$ and other corresponding GSVD components. By
introducing a linear operator $\mathcal{A}$ induced by $\\{A,L\\}$ between two
finite-dimensional Hilbert spaces, where a non-Euclidean inner-product is
used, we show that the SVE of $\mathcal{A}$ is just the nontrivial GSVD
components of $A$, i.e. the first decomposition in (1.1a). Similarly, we
introduce a linear operator $\mathcal{L}$ induced by $\\{A,L\\}$ and show that
the SVE of $\mathcal{L}$ is just the nontrivial GSVD components of $L$. This
result reveals that the nontrivial part of the GVSD of $\\{A,L\\}$ is nothing
but the SVEs of the two linear operators induced by $\\{A,L\\}$. Combined with
the trivial GSVD components $\\{x_{i}\\}$, it completely characterizes the
structure of GSVD for any matrix pair with the same number of columns.

As a direct application of the above result, we propose a new iterative method
for computing several extreme nontrivial GSVD components of $\\{A,L\\}$. This
iterative method is a natural extension of the Golub-Kahan bidiagonalization
(GKB) method for large-scale SVD computation [18], which iteratively reduces a
matrix to a bidiagonal form by a Lanczos-type iterative process. There are
several variants and extensions for the standard GKB of a single matrice,
which are proposed to solve large-scale generalized least squares problems [6,
3], saddle point problems [14], or regularization of inverse problems [11, 12,
31, 32, 34]; most of them are named the generalized Golub-Kahan
bidiagonalization (gGKB) while some have other different names. As a natural
analogy of the standard GKB for SVD computation, we develop an operator-type
GKB for linear operators $\mathcal{A}$ and $\mathcal{L}$ to approximate their
SVE components, which is also named the gGKB process. We derive matrix-form
recursive relations for this operator-type GKB so that it can be used in
practical computations. Moreover, this approach offers a unified and general
treatment for extending the standard GKB, which can be used to derive nearly
all of the aforementioned gGKB processes.

Using the GSVD of $\\{A,L\\}$, we study several basic properties of the
proposed gGKB process. Due to the correspondence of GSVD and SVE, the gGKB
method can approximate well the extreme nontrivial GSVD components of
$\\{A,L\\}$, resulting in the gGKB_GSVD algorithm. We derive a relative
residual norm and its sharp upper bound for the computed nontrivial GSVD
components, which is a good measure of the approximating accuracy and can be
used in a stopping criterion for practical computations. Several preliminary
results about the convergence and accuracy of gGKB_GSVD in exact arithmetic
are provided, showing the effectiveness of this method.

The paper is organized as follows. In Section 2, we review several basic
properties of the GSVD. In Section 3, we introduce two linear operators
induced by $\\{A,L\\}$ to characterize the structure of GSVD by the SVEs of
them. In Section 4 we propose the new gGKB process and study its basic
properties. In Section 5, we propose the gGKB_GSVD algorithm for computing
nontrivial GSVD components. Numerical experiments are presented in Section 6,
and concluding remarks follow in Section 7.

##  2 GSVD, SVD and Golub-Kahan bidiagonalization

We review several basic properties of the GSVD of $\\{A,L\\}$ presented in
Theorem 1.1. The nontrivial generalized singular values of $\\{A,L\\}$ in
descending order are

$\underbrace{\infty,\dots,\infty}_{q_{1}},\ \
\underbrace{c_{q_{1}+1}/s_{q_{1}+1},\dots,c_{q_{1}+q_{2}}/s_{q_{1}+q_{2}}}_{q_{2}},\
\ \underbrace{0,\dots,0}_{q_{3}}.$ (2.1)

We remark that the three numbers $q_{1}$, $q_{2}$, $q_{3}$ are uniquely
determined by the property of $\\{A,L\\}$, and some of them may be zero in
certain instances. The nontrivial GSVD components are linked by the vector-
form relations

$\begin{cases}Ax_{i}=c_{i}p_{A,i}\\\ Lx_{i}=s_{i}p_{L,i}\\\
s_{i}A^{T}p_{A,i}=c_{i}L^{T}p_{L,i}\end{cases}$ (2.2)

for $i=1,\dots,r$. For those trivial GSVD components, it holds that

$Ax_{i}=\mathbf{0},\ \ Lx_{i}=\mathbf{0},\ \ A^{T}p_{A,i}=\mathbf{0},\ \
L^{T}p_{L,i}=\mathbf{0}$ (2.3)

for $i=r+1,\dots,n$. The following result describes the structure of trivial
and nontrivial GSVD components $\\{x_{i}\\}$.

###### Proposition 2.1

Let $M=A^{\top}A+L^{\top}L$ and partition $X$ as  $\left(\vbox{\vbox{\halign{
\hfil$#$\hfil\kern 2\p@\kern\@leftbwd& \thinspace\hfil$#$\hfil&&
\quad\hfil$#$\hfil\cr X_{1}$\hfil\kern 2.0pt\kern 6.66669pt&X_{2}&\crcr\kern
2.0pt\cr r$\hfil\kern 2.0pt\kern 6.66669pt&n-r\cr} }\kern-27.34329pt}\kern
0.0pt\right)\kern-6.66669pt\kern 0.0pt$ . Then
$\mathcal{R}(X_{2})=\mathcal{N}(M)$. Moreover, let

$\bar{X}=(\bar{X}_{1}\ \ X_{2}),\ \ \
\bar{X}_{1}=(\mathcal{P}_{\mathcal{R}(M)}x_{1},\dots,\mathcal{P}_{\mathcal{R}(M)}x_{r}).$
(2.4)

Then it holds

$A=P_{A}C_{A}\bar{X}^{-1},\ \ \ L=P_{L}S_{L}\bar{X}^{-1}.$ (2.5)

###### Proof.

It is clear that $\mathcal{N}(M)=\mathcal{N}(A)\cap\mathcal{N}(L)$, and using
th GSVD of $\\{A,L\\}$, it is easy to obtain that
$\mathcal{N}(A)\cap\mathcal{N}(L)=\mathcal{R}(X_{2})$. Thus, we have
$\mathcal{R}(X_{2})=\mathcal{N}(M)$. Using the relation that
$\mathcal{P}_{\mathcal{R}(M)}x_{i}=x_{i}-\mathcal{P}_{\mathcal{N}(M)}x_{i}$,
for $1\leq i\leq r$ we have

$A\mathcal{P}_{\mathcal{R}(M)}x_{i}=Ax_{i}-A\mathcal{P}_{\mathcal{N}(M)}x_{i}=Ax_{i},\
\
L\mathcal{P}_{\mathcal{R}(M)}x_{i}=Lx_{i}-L\mathcal{P}_{\mathcal{N}(M)}x_{i}=Lx_{i}.$

Using the above two relations, it is easy to verify (2.5). ∎

Since
$\mathrm{dim}(\mathcal{N}(M))=n-r=\mathrm{rank}(\\{x_{i}\\}_{i=r+1}^{n})$, it
follows that $\\{x_{i}\\}_{i=r+1}^{n}$ forms a basis for $\mathcal{N}(M)$.
Proposition 2.1 also indicates that for any $x_{i}$ with $1\leq i\leq r$,
$\mathcal{P}_{\mathcal{R}(M)}x_{i}$ is also an $i$-th generalized singular
vector of $\\{A,L\\}$. Therefore, any nontrivial $x_{i}$ can be decomposed
into two components, one being
$\mathcal{P}_{\mathcal{R}(M)}x_{i}\in\mathcal{R}(M)$ and the other being an
arbitrary vector in $\mathcal{N}(M)$, which means that any nontrivial $x_{i}$
belongs to the coset $\bar{x}_{i}+\mathcal{N}(M)$, where
$\bar{x}_{i}\in\mathcal{R}(M)$ is the $i$-th nontrivial GSVD component. In
particular, we can focus on those nontrivial $x_{i}$ in $\mathcal{R}(M)$,
which results in the new form of GSVD (2.5). In the subsequent part, we always
consider this form of GSVD by requiring

$x_{i}\in\mathcal{R}(M)\ \ \mathrm{for}\ \ 1\leq i\leq r.$ (2.6)

There exists a direct relationship between the SVD and GSVD for a matrix pair
with a special property. If $L$ has full column rank, it follows from (1.1)
that $r=n$ and $q_{1}=0$, leading to

$AL^{{\dagger}}=P_{A}C_{A}X^{-1}[(P_{L}S_{L})X^{-1}]^{{\dagger}}=P_{A}C_{A}X^{-1}X(P_{L}S_{L})^{{\dagger}}=P_{A}(C_{A}S_{L}^{{\dagger}})P_{L}^{\top},$
(2.7)

where we have used the property that
$(M_{1}M_{2})^{{\dagger}}=M_{2}^{{\dagger}}M_{1}^{{\dagger}}$ if $M_{1}$ has
full column rank and $M_{2}$ has full row rank. Therefore, the SVD of
$AL^{{\dagger}}$ is $P_{A}(C_{A}S_{L}^{{\dagger}})P_{L}^{\top}$ with the
singular values be $\\{\gamma_{i}\\}_{i=1}^{n}$.

For the above case, one can compute the SVD of $AL^{{\dagger}}$ to get the
GSVD components 111For numerical computations, it is not recommended to
explicitly form $AL^{{\dagger}}$ due to its numerical unstability, especially
when $L$ is close to column rank-deficient.. For a large-scale matrix, the GKB
process is a basic routine for computing a few extreme SVD components. At each
iteration, it reduces the matrix to a lower-order bidiagonal matrix and
generates two Krylov subspaces. The Rayleigh-Ritz procedure is then exploited
to approximate extreme SVD components using the bidiagonal matrix and Krylov
subspaces [4].

In the following part of the paper, we characterize the GSVD from the
viewpoint of SVD for any matrix pair with the same number of columns. Then we
generalize the GKB process so that it can be used to compute extreme GSVD
components.

##  3 Characterizing GSVD by singular value expansion of linear operators

We first discuss linear operators between two finite-dimensional Hilbert
spaces. Then we study the GSVD of $\\{A,L\\}$ using the singular value
expansion of linear operators. Note that Section 3.1 is quite general without
requiring $M=A^{\top}A+L^{\top}L$.

###  3.1 Linear operator induced by matrices

Suppose $G\in\mathbb{R}^{m\times m}$ is symmetric positive definite. It is
obviously that $\langle u,u^{\prime}\rangle_{G}:=u^{\top}Gu^{\prime}$ defines
an inner product on $\mathbb{R}^{m}$ such that
$(\mathbb{R}^{m},\langle\cdot,\cdot\rangle_{G})$ is an $m$-dimensional Hilbert
space. On the other hand, for any symmetric positive semidefinite matrix
$M\in\mathbb{R}^{n\times n}$ with $\mathrm{rank}(M)=r$, the bilinear form
$\langle v,v^{\prime}\rangle_{M}:=v^{\top}Mv^{\prime}$ is not a well-defined
inner product on $\mathbb{R}^{n}$ if $r<n$. In this case, we consider the
inner product on the subspace $\mathcal{R}(M)$.

###### Lemma 3.1

The bilinear form $\langle v,v^{\prime}\rangle_{M}:=v^{\top}Mv^{\prime}$ for
any $v,v^{\prime}\in\mathcal{R}(M)$ is an inner product, and
$(\mathcal{R}(M),\langle\cdot,\cdot\rangle_{M})$ is an $r$-dimensional Hilbert
space.

###### Proof.

The statement $\mathrm{dim}(\mathcal{R}(M))=\mathrm{rank}(M)$ is a basic
property. We need to show that $\langle\cdot,\cdot\rangle_{M}$ is indeed an
inner product, i.e., it is a symmetric and positive bilinear form on
$\mathcal{R}(M)\times\mathcal{R}(M)$. We only need to prove the positiveness.
To see it, let $v\in\mathcal{R}(M)$ such that $\langle
v,v\rangle_{M}=v^{\top}Mv=0$. It follows $v\in\mathcal{N}(M)$. Note that
$\mathcal{R}(M)\cap\mathcal{N}(M)=\\{\mathbf{0}\\}$ since $M$ is symmetric,
which leads to $v=\mathbf{0}$. ∎

Given a matrix $A\in\mathbb{R}^{m\times n}$. Define the linear map

$\mathcal{A}:(\mathcal{R}(M),\langle\cdot,\cdot\rangle_{M})\rightarrow(\mathbb{R}^{m},\langle\cdot,\cdot\rangle_{G}),\
\ \ v\mapsto Av,$ (3.1)

where $v$ and $Av$ are column vectors under the canonical bases of
$\mathbb{R}^{n}$ and $\mathbb{R}^{m}$. Let $W_{r}\in\mathbb{R}^{n\times r}$
and $Z\in\mathbb{R}^{m\times m}$ be two matrices whose columns constitute
orthonormal bases of $(\mathcal{R}(M),\langle\cdot,\cdot\rangle_{M})$ and
$(\mathbb{R}^{m},\langle\cdot,\cdot\rangle_{G})$, respectively, i.e.
$W_{r}^{\top}MW_{r}=I_{r}$ and $Z^{\top}GZ=I_{m}$. Then we have the
commutative diagram:

${(\mathcal{R}(M),\langle\cdot,\cdot\rangle_{M})}$${(\mathbb{R}^{m},\langle\cdot,\cdot\rangle_{G})}$${(\mathbb{R}^{r},\langle\cdot,\cdot\rangle_{2})}$${(\mathbb{R}^{m},\langle\cdot,\cdot\rangle_{2})\
,}$$\scriptstyle{\mathcal{A}}$$\scriptstyle{\pi_{1}}$$\scriptstyle{\pi_{2}}$$\scriptstyle{[\mathcal{A}]}$
(3.2)

where $[\mathcal{A}]$ denotes the matrix representation of $\mathcal{A}$ under
bases $W_{r}$ and $Z$, and $\pi_{1}$ and $\pi_{2}$ are two linear maps:

$\displaystyle\pi_{1}:(\mathbb{R}^{r},\langle\cdot,\cdot\rangle_{2})\rightarrow(\mathcal{R}(M),\langle\cdot,\cdot\rangle_{M}),\
\ \ x\mapsto W_{r}x,$ (3.3)
$\displaystyle\pi_{2}:(\mathbb{R}^{m},\langle\cdot,\cdot\rangle_{2})\rightarrow(\mathbb{R}^{m},\langle\cdot,\cdot\rangle_{G}),\
\ \ \ \ \ y\mapsto Zy.$ (3.4)

Note that $\pi_{1}$ and $\pi_{2}$ are two isomorphism maps such that
$(\mathcal{R}(M),\langle\cdot,\cdot\rangle_{M})\cong(\mathbb{R}^{r},\langle\cdot,\cdot\rangle_{2})$
and
$(\mathbb{R}^{m},\langle\cdot,\cdot\rangle_{G})\cong(\mathbb{R}^{m},\langle\cdot,\cdot\rangle_{2})$.
Since $\mathcal{A}$ is a bounded linear operator, we can define its adjoint

$\mathcal{A}^{*}:(\mathbb{R}^{m},\langle\cdot,\cdot\rangle_{G})\rightarrow(\mathcal{R}(M),\langle\cdot,\cdot\rangle_{M}),\
\ \ u\mapsto\mathcal{A}^{*}u$ (3.5)

by the relation

$\langle\mathcal{A}v,u\rangle_{G}=\langle\mathcal{A}^{*}u,v\rangle_{M}$ (3.6)

for any $v\in(\mathcal{R}(M),\langle\cdot,\cdot\rangle_{M})$ and
$u\in(\mathbb{R}^{m},\langle\cdot,\cdot\rangle_{G})$. We have the following
corresponding commutative diagram:

${(\mathcal{R}(M),\langle\cdot,\cdot\rangle_{M})}$${(\mathbb{R}^{m},\langle\cdot,\cdot\rangle_{G})}$${(\mathbb{R}^{r},\langle\cdot,\cdot\rangle_{2})}$${(\mathbb{R}^{m},\langle\cdot,\cdot\rangle_{2})\
,}$$\scriptstyle{\mathcal{A}^{*}}$$\scriptstyle{\pi_{1}}$$\scriptstyle{\pi_{2}}$$\scriptstyle{[\mathcal{A}^{*}]}$
(3.7)

where $[\mathcal{A}^{*}]$ is the matrix representation of $\mathcal{A}^{*}$
under bases $W_{r}$ and $Z$. The following result gives the matrix-forms of
$[\mathcal{A}]$ and $[\mathcal{A}^{*}]$.

###### Lemma 3.2

The matrix representations of $\mathcal{A}$ and $\mathcal{A}^{*}$ under bases
$W_{r}$ and $Z$ are

$[\mathcal{A}]=Z^{-1}AW_{r},\ \ \ [\mathcal{A}^{*}]=W_{r}^{\top}A^{\top}GZ.$
(3.8)

###### Proof.

By the commutative diagram (3.2), we have
$\mathcal{A}\circ\pi_{1}(x)=\pi_{2}\circ[\mathcal{A}](x)$ for any
$x\in\mathbb{R}^{r}$, which is equivalent to $AW_{r}x=Z[\mathcal{A}]x$ Thus,
we have $AW_{r}=Z[\mathcal{A}]$, leading to $[\mathcal{A}]=Z^{-1}AW_{r}$.
Similarly, by the commutative diagram (3.7), we have

$\mathcal{A}^{*}\circ\pi_{2}(y)=\pi_{1}\circ[\mathcal{A}^{*}](y)\ \
\Leftrightarrow\ \ \mathcal{A}^{*}Zy=W_{r}[\mathcal{A}^{*}]y$

for any $y\in\mathbb{R}^{m}$, which leads to

$\mathcal{A}^{*}Z=W_{r}[\mathcal{A}^{*}].$ (3.9)

From the definition of $\mathcal{A}^{*}$, we have

$\langle\mathcal{A}\circ\pi_{1}(x),\pi_{2}(y)\rangle_{G}=\langle\pi_{1}(x),\mathcal{A}^{*}\circ\pi_{2}(y)\rangle_{M}$

for any $x\in\mathbb{R}^{r}$ and $y\in\mathbb{R}^{m}$, which can also be
written as

$(AW_{r}x)^{\top}GZy=(W_{r}x)^{\top}M(\mathcal{A}^{*}Zy)\ \ \Leftrightarrow\ \
x^{\top}W_{r}^{\top}A^{\top}GZy=x^{\top}W_{r}^{\top}M\mathcal{A}^{*}Zy.$

Thus, we have

$W_{r}^{\top}A^{\top}GZ=W_{r}^{\top}M\mathcal{A}^{*}Z$ (3.10)

Combining (3.9) and (3.10) and using $W_{r}^{\top}MW_{r}=I_{r}$, we finally
obtain

$[\mathcal{A}^{*}]=W_{r}^{\top}MW_{r}[\mathcal{A}^{*}]=W_{r}^{\top}M\mathcal{A}^{*}Z=W_{r}^{\top}A^{\top}GZ.$

The proof is completed. ∎

The following result will be used throughout the paper.

###### Lemma 3.3

If $\mathcal{R}(W_{r})=\mathcal{R}(M)$ and $W_{r}^{\top}MW_{r}=I_{r}$, then
the Moor-Penrose pseudoinverse of $M$ can be expressed as

$M^{{\dagger}}=W_{r}W_{r}^{\top}.$

###### Proof.

Let $\bar{M}=W_{r}W_{r}^{\top}$. We only need to verify the following four
identities:

$\displaystyle M\bar{M}M=M,\ \ \ (M\bar{M})^{\top}=M\bar{M},$
$\displaystyle\bar{M}M\bar{M}=\bar{M},\ \ \ (\bar{M}M)^{\top}=\bar{M}M.$

The third identity is the easiest to verify:
$\bar{M}M\bar{M}=W_{r}W_{r}^{\top}MW_{r}W_{r}^{\top}=W_{r}W_{r}^{\top}=\bar{M}$.
Suppose the compact-form eigenvalue decomposition of $M$ is

$M=P_{r}\Lambda_{r}P_{r}^{\top},\ \ \
\Lambda_{r}=\operatorname{diag}(\lambda_{1},\dots,\lambda_{r}),$

where $\lambda_{1}\geq\cdots\geq\lambda_{r}>0$ and
$P_{r}\in\mathbb{R}^{n\times r}$ with $2$-orthonormal columns. Since
$\mathcal{R}(W_{r})=\mathcal{R}(M)=\mathcal{R}(P_{r})$, there exist
$D\in\mathbb{R}^{r\times r}$ such that $W_{r}=P_{r}D$. It follows that
$I_{r}=W_{r}^{\top}MW_{r}=D^{\top}\Lambda_{r}D.$ Therefore, it follows that

$M\bar{M}M=MP_{r}DD^{\top}P_{r}^{\top}M=P_{r}\Lambda_{r}DD^{\top}\Lambda_{r}P_{r}^{\top}=P_{r}\Lambda_{r}P_{r}^{\top}=M,$

since $\Lambda_{r}DD^{\top}=D^{\top}\Lambda_{r}D=I_{r}$. Similarly, we have

$\displaystyle
M\bar{M}=P_{r}\Lambda_{r}P_{r}^{\top}P_{r}DD^{\top}P_{r}^{\top}=P_{r}\Lambda_{r}DD^{\top}P_{r}^{\top}=P_{r}P_{r}^{\top}=(M\bar{M})^{\top},$
$\displaystyle\bar{M}M=P_{r}DD^{\top}P_{r}^{\top}M=P_{r}DD^{\top}\Lambda_{r}P_{r}^{\top}=P_{r}P_{r}^{\top}=(\bar{M}M)^{\top}.$

Now all the four identities have been verified. ∎

###  3.2 Characterizing GSVD by singular value expansion

In this subsection, we consider the simpler case that $G=I_{m}$, since it has
direct connections with the GSVD of a matrix pair. For notational simplicity,
let $\mathcal{X}=(\mathcal{R}(M),\langle\cdot,\cdot\rangle_{M})$ and
$\mathcal{Y}=(\mathbb{R}^{m},\langle\cdot,\cdot\rangle_{2})$. For the linear
compact operator

$\mathcal{A}:\mathcal{X}\rightarrow\mathcal{Y},\ \ \ v\mapsto Av$ (3.11)

between the two Hilbert spaces $\mathcal{X}$ and $\mathcal{Y}$, where $v$ and
$Av$ are column vectors under the canonical bases of $\mathbb{R}^{n}$ and
$\mathbb{R}^{m}$, it has the singular value expansion (SVE) with finite terms;
see e.g. [28, §15.4]. Here we use the terminology “SVE” instead of “SVD” to
distinguish it from the SVD of a matrix. The theory of SVE for $\mathcal{A}$
states that there exist positive scalars
$\sigma_{1}\geq\cdots\geq\sigma_{d}>0$, two orthonormal systems
$\\{f_{i}\\}_{i=1}^{d}\subseteq\mathcal{X}$ and
$\\{h_{i}\\}_{i=1}^{d}\subseteq\mathcal{Y}$ such that

$\mathcal{A}f_{i}=\sigma_{i}h_{i},\ \ \ \mathcal{A}^{*}h_{i}=\sigma_{i}f_{i},$
(3.12)

and any $v\in\mathcal{X}$ has the expansion

$v=v_{0}+\sum_{i=1}^{d}\langle v,f_{i}\rangle_{M}f_{i}$ (3.13)

with some $v_{0}\in\mathrm{ker}(\mathcal{A})$, and

$\mathcal{A}v=\sum_{i=1}^{d}\sigma_{i}\langle v,f_{i}\rangle_{M}h_{i},$ (3.14)

where $d=\mathrm{dim}(\mathrm{im}(\mathcal{A}))$. Here we use
$\mathrm{ker}(\cdot)$ and $\mathrm{im}(\cdot)$ to denote the kernel and image
of a linear operator, respectively, to distinguish them from the null space
$\mathcal{N}(\cdot)$ and range space $\mathcal{R}(\cdot)$ of a matrix.

The following result provides more details about the SVE of $\mathcal{A}$.

###### Theorem 3.1

For any $A\in\mathbb{R}^{m\times n}$ and symmetric positive semidefinite
matrix $M\in\mathbb{R}^{n\times n}$ with rank $r$, define the linear operator
$\mathcal{A}$ as (3.11). Then there exist an $M$-orthonormal matrix
$F\in\mathbb{R}^{n\times r}$, a $2$-orthonormal matrix
$H\in\mathbb{R}^{m\times m}$ and a diagonal matrix
$\Sigma\in\mathbb{R}^{m\times r}$, such that for any $v\in\mathcal{X}$ and
$u\in\mathcal{Y}$, it holds that

$\mathcal{A}v=H\Sigma F^{\top}Mv,\ \ \
\mathcal{A}^{*}u=F\Sigma^{\top}H^{\top}u$ (3.15)

under the canonical bases of $\mathbb{R}^{n}$ and $\mathbb{R}^{m}$.

###### Proof.

Let $\mathcal{X}_{1}=\mathrm{span}\\{f_{i}\\}_{i=1}^{d}$. We first prove
$\mathcal{X}=\mathrm{ker}(\mathcal{A})\oplus\mathcal{X}_{1}$. Noticing (3.13),
we only need to prove
$\mathrm{ker}(\mathcal{A})\cap\mathcal{X}_{1}=\\{\mathbf{0}\\}$. Let
$v=\sum_{j=1}^{d}\mu_{j}f_{j}\in\mathrm{ker}(\mathcal{A})\cap\mathcal{X}_{1}$.
By (3.14), it follows
$\mathbf{0}=\mathcal{A}v=\sum_{j=1}^{n}\sigma_{i}\mu_{i}h_{i}$, leading to
$\sigma_{i}\mu_{i}=0$ for $1\leq i\leq d$. Since $\sigma_{i}>0$, we have
$\mu_{i}=0$ for $1\leq i\leq d$, thereby $v=\mathbf{0}$. Then we prove
$\mathrm{ker}(\mathcal{A})\perp_{M}\mathcal{X}_{1}$, where $\perp_{M}$ is the
orthogonal relation in $\mathcal{X}$. For any $v\in\mathrm{ker}(\mathcal{A})$
and any $f_{i}$, by (3.12) we have
$f_{i}=\sigma_{i}^{-1}\mathcal{A}^{*}h_{i}$, thereby

$\langle v,f_{i}\rangle_{M}=\langle
v,\sigma_{i}^{-1}\mathcal{A}^{*}h_{i}\rangle_{M}=\langle\mathcal{A}v,\sigma_{i}^{-1}h_{i}\rangle_{2}=\langle\mathbf{0},\sigma_{i}^{-1}h_{i}\rangle_{2}=0.$

Note $\mathrm{dim}(\mathcal{X}_{1})=d$. Therefore, we can find $r-d$
$M$-orthonormal vectors in $\mathrm{ker}(\mathcal{A})$ that are $M$-orthogonal
to each $f_{i}$. Denote these vectors by $\\{f_{d+1},\dots,f_{r}\\}$. Then
$\\{f_{i}\\}_{i=1}^{r}$ constitute a complete orthonormal basis for
$\mathcal{X}$. From (3.14) we have
$\mathrm{im}(\mathcal{A})=\mathrm{span}\\{h_{i}\\}_{i=1}^{d}=:\mathcal{Y}_{1}$.
Using the relation
$\mathrm{ker}(\mathcal{A}^{*})=\mathrm{im}(\mathcal{A})^{\perp}=\mathcal{Y}_{1}^{\perp}$,
where the orthogonality is taken in
$(\mathbb{R}^{m},\langle\cdot,\cdot\rangle_{2})$, there exist $m-d$
2-orthonormal
$\\{h_{d+1},\dots,h_{m}\\}\subseteq\mathrm{ker}(\mathcal{A}^{*})$ such that
$\\{h_{i}\\}_{i=1}^{m}$ constitute a complete orthonormal basis for
$\mathcal{Y}$.

Therefore, for any $v\in\mathcal{X}$, it can be written as
$v=\sum_{i=1}^{r}\langle v,f_{i}\rangle_{M}f_{i}$, and thereby

$\mathcal{A}v=\sum_{i=1}^{r}\langle
v,f_{i}\rangle_{M}\mathcal{A}f_{i}=\sum_{i=1}^{r}\sigma_{i}\langle
v,f_{i}\rangle_{M}h_{i},$

where we define $\sigma_{d+1}=\cdots=\sigma_{r}=0$. Similarly, for any
$u\in\mathcal{Y}$ with the expansion $u=\sum_{i=1}^{m}\langle
u,h_{i}\rangle_{2}h_{i}$, it holds

$\mathcal{A}^{*}u=\sum_{i=1}^{m}\langle
u,h_{i}\rangle_{M}\mathcal{A}^{*}h_{i}=\sum_{i=1}^{r}\sigma_{i}\langle
u,h_{i}\rangle_{2}f_{i}.$

Let the matrices $F=(f_{1},\dots,f_{r})$, $H=(h_{1},\dots,h_{m})$ and
$\Sigma=\begin{pmatrix}\Sigma_{d}&\\\
&\mathbf{0}\end{pmatrix}\in\mathbb{R}^{m\times r}$ with
$\Sigma_{d}=\operatorname{diag}(\sigma_{1},\dots,\sigma_{d})$. Then (3.15) is
just the matrix-form of the above two identities. ∎

One can verify that (3.12), (3.13) and (3.14) can be derived from (3.15).
Therefore, the two relations in (3.15) describe completely the SVE of
$\mathcal{A}$. In the following part, we use the notation

$\mathcal{A}\sim H\Sigma F^{\top}$ (3.16)

to denote the SVE of $\mathcal{A}$. From the proof of Theorem 3.1, we have the
following basic relations for the four important subspaces:

$\displaystyle\mathrm{ker}(\mathcal{A})=\mathrm{span}\\{f_{i}\\}_{i=d+1}^{r},\
\ \ \ \ \mathrm{im}(\mathcal{A})=\mathrm{span}\\{h_{i}\\}_{i=1}^{d},$ (3.17a)
$\displaystyle\mathrm{ker}(\mathcal{A}^{*})=\mathrm{span}\\{h_{i}\\}_{i=d+1}^{m},\
\ \ \mathrm{im}(\mathcal{A}^{*})=\mathrm{span}\\{f_{i}\\}_{i=1}^{d}.$ (3.17b)

From the theory of SVE for linear compact operators, if the multiplicity of
$\sigma_{i}$ is $1$, then the corresponding $f_{i}$ and $h_{i}$ are uniquely
determined (at most differing by a sign). If the multiplicity of $\sigma_{i}$
is $m_{i}>1$, then there are $m_{i}$ corresponding linear independent
$\\{f_{i}\\}$ and $\\{h_{i}\\}$, respectively, which are $M$-orthonormal and
$2$-orthonormal.

Based on Theorem 3.1, now we can use the SVE of $\mathcal{A}$ to characterize
the GSVD of $\\{A,L\\}$. Remember that we consider those
$x_{i}\in\mathcal{R}(M)$ for $1\leq i\leq r$.

###### Theorem 3.2

Let the GSVD of $\\{A,L\\}$ be (1.1) and let $\mathcal{A}$ be defined as
(3.11) with $M=A^{\top}A+L^{\top}L$. Partition $P_{A}$ and $X$ as

$P_{A}=\vtop{\hbox{$\left(\vbox{\vbox{\halign{ \hfil$#$\hfil\kern
2\p@\kern\@leftbwd& \thinspace\hfil$#$\hfil&& \quad\hfil$#$\hfil\cr
P_{A1}$\hfil\kern 2.0pt\kern 6.66669pt&P_{A2}&P_{A3}&m\crcr\kern 2.0pt\cr
q_{1}$\hfil\kern 2.0pt\kern 6.66669pt&q_{2}&m-q_{1}-q_{2}\cr}
}\kern-29.39662pt}\kern 0.0pt\right)\kern-6.66669pt\kern 0.0pt$}\kern
29.39662pt}\ ,\ \ \ X=\vtop{\hbox{$\left(\vbox{\vbox{\halign{
\hfil$#$\hfil\kern 2\p@\kern\@leftbwd& \thinspace\hfil$#$\hfil&&
\quad\hfil$#$\hfil\cr X_{1}$\hfil\kern 2.0pt\kern
6.66669pt&X_{2}&X_{3}&X_{4}&n\crcr\kern 2.0pt\cr q_{1}$\hfil\kern 2.0pt\kern
6.66669pt&q_{2}&q_{3}&n-r\cr} }\kern-29.28773pt}\kern
0.0pt\right)\kern-6.66669pt\kern 0.0pt$}\kern 29.28773pt}$ (3.18)

and let $\widetilde{X}_{1}=(X_{1}\ X_{2}\ X_{3})$. Then the SVE of
$\mathcal{A}$ has the form

$\mathcal{A}\sim P_{A}\Sigma_{A}\widetilde{X}_{1}^{\top},$ (3.19)

and it holds that

$\displaystyle\mathrm{ker}(\mathcal{A})=\mathcal{R}(X_{3}),\ \ \ \ \ \ \ \ \ \
\mathrm{im}(\mathcal{A})=\mathcal{R}((P_{A1}\ P_{A2})),$ (3.20a)
$\displaystyle\mathrm{ker}(\mathcal{A}^{*})=\mathcal{R}(P_{A3}),\ \ \ \ \ \ \
\mathrm{im}(\mathcal{A}^{*})=\mathcal{R}((X_{1}\ X_{2})).$ (3.20b)

###### Proof.

Using the GSVD of $\\{A,L\\}$, we have

$\displaystyle
A^{\top}A+L^{\top}L=X^{-\top}\left(\begin{pmatrix}\Sigma_{A}^{\top}\Sigma_{A}&\\\
&\mathbf{0}\end{pmatrix}+\begin{pmatrix}\Sigma_{L}^{\top}\Sigma_{L}&\\\
&\mathbf{0}\end{pmatrix}\right)X^{-1}=X^{-\top}\begin{pmatrix}I_{r}&\\\
&\mathbf{0}\end{pmatrix}X^{-1},$

which leads to $\mathrm{rank}(M)=r$ and

$\displaystyle\begin{pmatrix}I_{r}&\\\
&\mathbf{0}\end{pmatrix}=\begin{pmatrix}\widetilde{X}_{1}^{\top}\\\
X_{4}^{\top}\end{pmatrix}M\begin{pmatrix}\widetilde{X}_{1}&X_{4}\end{pmatrix}=\begin{pmatrix}\widetilde{X}_{1}^{\top}M\widetilde{X}_{1}&\widetilde{X}_{1}^{\top}MX_{4}\\\
X_{4}^{\top}M\widetilde{X}_{1}&X_{4}^{\top}MX_{4}\end{pmatrix}.$

Therefore, we have $\widetilde{X}_{1}^{\top}M\widetilde{X}_{1}=I_{r}$. Note
that $\mathcal{R}(\widetilde{X}_{1})\subseteq\mathcal{R}(M)$. It follows that
$\widetilde{X}_{1}$ is an $M$-orthonormal basis of
$(\mathcal{R}(M),\langle\cdot,\cdot\rangle_{M})$, thereby we obtain from (3.3)
that $M^{{\dagger}}=\widetilde{X}_{1}\widetilde{X}_{1}^{\top}$. Notice from
(1.1) that

$A(\widetilde{X}_{1}\ X_{4})=P_{A}(\Sigma_{A}\ \mathbf{0})\Rightarrow
A\widetilde{X}_{1}=P_{A}\Sigma_{A}.$

Thus, we have
$AM^{{\dagger}}=A\widetilde{X}_{1}\widetilde{X}_{1}^{\top}=P_{A}\Sigma_{A}\widetilde{X}_{1}^{\top}$.
For any $v\in(\mathcal{R}(M),\langle\cdot,\cdot\rangle_{M})$, it holds

$\displaystyle\mathcal{A}v=\mathcal{A}\mathcal{P}_{\mathcal{R}(M)}v=AM^{{\dagger}}Mv=P_{A}\Sigma_{A}\widetilde{X}_{1}^{\top}Mv.$

Using the commutative diagram (3.7) and noticing $G=Z=I_{m}$ for the current
case, we have for any $u\in(\mathbb{R}^{m},\langle\cdot,\cdot\rangle_{2})$
that

$\mathcal{A}^{*}u=\pi_{1}\circ[\mathcal{A}^{*}](u)=\widetilde{X}_{1}(\widetilde{X}_{1}^{\top}A^{\top})u=M^{{\dagger}}A^{\top}u=(AM^{{\dagger}})^{\top}u=\widetilde{X}_{1}\Sigma_{A}^{\top}P_{A}^{\top}u,$
(3.21)

where we have used $[\mathcal{A}^{*}]=\widetilde{X}_{1}^{\top}A^{\top}$ by
Lemma 3.2. This proves that the SVE of $\mathcal{A}$ has the form
$P_{A}\Sigma_{A}\widetilde{X}_{1}^{\top}$.

From the SVE of $\mathcal{A}$ we have
$\mathrm{dim}(\mathrm{im}(\mathcal{A}))=q_{1}+q_{2}$. Using the relations
(3.17), it follows that $\mathrm{im}(\mathcal{A})=\mathcal{R}((P_{A1}\
P_{A2}))$. Since $P_{A}$ is a 2-orthogonal matrix, we then have
$\mathrm{ker}(\mathcal{A}^{*})=\mathcal{R}(P_{A3})$. The other two relations
can also be verified easily. ∎

Corresponding to Theorem 3.2, we have the following result.

###### Theorem 3.3

Define $\mathcal{L}$ as

$\mathcal{L}:(\mathcal{R}(M),\langle\cdot,\cdot\rangle_{M})\rightarrow(\mathbb{R}^{p},\langle\cdot,\cdot\rangle_{2}),\
\ \ v\mapsto Lv,$ (3.22)

where $v$ and $Lv$ are column vectors under the canonical bases of
$\mathbb{R}^{n}$ and $\mathbb{R}^{p}$. Partition $P_{L}$ as

$P_{L}=\vtop{\hbox{$\left(\vbox{\vbox{\halign{ \hfil$#$\hfil\kern
2\p@\kern\@leftbwd& \thinspace\hfil$#$\hfil&& \quad\hfil$#$\hfil\cr
P_{L1}$\hfil\kern 2.0pt\kern 6.66669pt&P_{L2}&P_{L3}&p\crcr\kern 2.0pt\cr
p-q_{2}-q_{3}$\hfil\kern 2.0pt\kern 6.66669pt&q_{2}&q_{3}\cr}
}\kern-29.42773pt}\kern 0.0pt\right)\kern-6.66669pt\kern 0.0pt$}\kern
29.42773pt}.$ (3.23)

Then the SVE of $\mathcal{L}$ has the form

$\mathcal{L}\sim P_{L}\Sigma_{L}\widetilde{X}_{1}^{\top},$ (3.24)

and it holds that

$\displaystyle\mathrm{ker}(\mathcal{L})=\mathcal{R}(X_{1}),\ \ \ \ \ \ \ \ \ \
\mathrm{im}(\mathcal{L})=\mathcal{R}((P_{L2}\ P_{L3})),$ (3.25a)
$\displaystyle\mathrm{ker}(\mathcal{L}^{*})=\mathcal{R}(P_{L1}),\ \ \ \ \ \ \
\ \mathrm{im}(\mathcal{L}^{*})=\mathcal{R}((X_{2}\ X_{3})).$ (3.25b)

Proposition 2.1 together with Theorem 3.2 and Theorem 3.3 characterizes
completely the GSVD of $\\{A,L\\}$ based on the SVEs of linear operators
$\mathcal{A}$ and $\mathcal{L}$. Particularly, they show that the nontrivial
part $\widetilde{X}_{1}$ is the common right SVE components of $\mathcal{A}$
and $\mathcal{L}$, while $P_{A}$ and $P_{L}$ are the left SVE components of
$\mathcal{A}$ and $\mathcal{L}$, respectively. Moreover, the relations (3.20)
and (3.25) use the image spaces and kernel spaces of $\mathcal{A},\
\mathcal{A}^{*}$ and $\mathcal{L},\ \mathcal{L}^{*}$ to describe the structure
of each GSVD blocks and give a new explanation of the three numbers $q_{1}$,
$q_{2}$ and $q_{3}$ in (1.1).

Based on the SVE characterization of GSVD, we can expect to modify those
algorithms for large-scale SVD computation for large-scale GSVD computation.
To compute nontrivial extreme GSVD components, we generalize the standard GKB
process from the viewpoint of linear operators.

##  4 Generalizing the Golub-Kahan bidiagonalization

In this section, the generalization of GKB is quite general without requiring
$M=A^{\top}A+L^{\top}L$, and we follow the notations and assumptions in
Section 3.1. For the linear operator in (3.1), the iterative process of GKB
can be described as follows; see [9] for discussions about GKB for linear
compact operators. Choosing a nonzero vector
$b\in(\mathbb{R}^{m},\langle\cdot,\cdot\rangle_{G})$, the basis recursive
relations are

$\begin{cases}\beta_{1}u_{1}=b,\\\
\alpha_{i}v_{i}=\mathcal{A}^{*}u_{i}-\beta_{i}v_{i-1},\\\
\beta_{i+1}u_{i+1}=\mathcal{A}v_{i}-\alpha_{i}u_{i},\end{cases}$ (4.1)

where $u_{i}\in(\mathbb{R}^{m},\langle\cdot,\cdot\rangle_{G})$ and
$v_{i}\in(\mathcal{R}(M),\langle\cdot,\cdot\rangle_{M})$, and $\alpha_{i}$ and
$\beta_{i}$ are positive scalars such that $\|v_{i}\|_{M}=\|u_{i}\|_{G}=1$.
Note that $v_{0}:=\mathbf{0}$ for the initial step.

For the purpose of practical computation, we need to derive a matrix-form
expression of the recursive relations. Using the isomorphisms $\pi_{1}$ and
$\pi_{2}$, denote $u_{i}$ and $v_{i}$ by $u_{i}=Zy_{i}$ and $v_{i}=W_{r}x_{i}$
with $y_{i}\in\mathbb{R}^{m}$ and $x_{i}\in\mathbb{R}^{r}$ for any $i\geq 1$.
Then we have

$\displaystyle\mathcal{A}v_{i}=\mathcal{A}\circ\pi_{1}(x_{i})=\pi_{2}\circ[\mathcal{A}]x_{i}=ZZ^{-1}AW_{r}x_{i}=Av_{i},$
$\displaystyle\mathcal{A}^{*}u_{i}=\mathcal{A}^{*}\circ\pi_{2}(y_{i})=\pi_{1}\circ[\mathcal{A}^{*}]y_{i}=W_{r}W_{r}^{\top}A^{\top}GZy_{i}=M^{{\dagger}}A^{\top}Gu_{i}.$

Therefore, the last two recursions in (4.1) can be written in the matrix-
vector forms

$\begin{cases}\alpha_{i}v_{i}=M^{{\dagger}}A^{\top}Gu_{i}-\beta_{i}v_{i-1}\\\
\beta_{i+1}u_{i+1}=Av_{i}-\alpha_{i}u_{i}.\end{cases}$ (4.2)

Using (4.2), the GKB of $\mathcal{A}$ can be proceeded step by step. We name
the above iterative process the generalized Golub-Kahan bidiagonalization
(gGKB). The pseudocode of gGKB is shown in Algorithm 1.

###### Remark 4.1

If $G$ is also positive semidefinite, define the linear operator
$\mathcal{A}:(\mathcal{R}(M),\langle\cdot,\cdot\rangle_{M})\rightarrow(\mathcal{R}(G),\langle\cdot,\cdot\rangle_{G})$
be $v\mapsto Av$. In this case, a similar gGKB process can be obtained. A
slight difference is that the initial vector should satisfy
$b\in\mathcal{R}(G)$.

Algorithm 1 Generalized Golub-Kahan bidiagonalization (gGKB)

Input: $A\in\mathbb{R}^{m\times n}$, $M\in\mathbb{R}^{n\times n}$,
$G\in\mathbb{R}^{m\times m}$, $b\in\mathbb{R}^{m}$

1:Initialize: let $\beta_{1}=\|b\|_{G}$, $u_{1}=b/\beta_{1}$,

2:Compute $\bar{s}=A^{\top}Gu_{1}$, $s=M^{{\dagger}}\bar{s}$,

3:$\alpha_{1}=\|s\|_{M}$, $v_{1}=s/\alpha_{1}$

4:for $i=1,2,\dots,k,$ do

5: $q=Av_{i}-\alpha_{i}u_{i}$,

6: $\beta_{i+1}=\|q\|_{G}$, $u_{i+1}=q/\beta_{i+1}$;

7: $\bar{s}=A^{\top}Gu_{i+1}$, $s=M^{{\dagger}}\bar{s}-\beta_{i+1}v_{i}$

8: $\alpha_{i+1}=\|s\|_{M}$, $v_{i+1}=s/\alpha_{i+1}$

9:end for

Output: $\\{\alpha_{i},\beta_{i}\\}_{i=1}^{k+1}$,
$\\{u_{i},v_{i}\\}_{i=1}^{k+1}$

For large-scale matrices, we can not directly compute $M^{{\dagger}}$. In this
case, using the relation

$M^{{\dagger}}\bar{s}=\mathop{\mathrm{argmin}}_{s\in\mathbb{R}^{n}}\|Ms-\bar{s}\|_{2},$
(4.3)

we compute $M^{{\dagger}}\bar{s}$ by iteratively solving the above least
squares problems. This can be done efficiently by using the LSQR algorithm
[42]. In this case, gGKB has the nested inner-outer iteration structure.

If $G=I_{m}$ and $M=I_{n}$, the gGKB becomes the standard GKB. If $G=I_{m}$
and $M$ is invertible, the gGKB is equivalent to the generalizations of GKB
prosed in [3, 11, 12, 32, 31] with different forms. The following result
describes the basic property of gGKB, very similar to that of the standard
GKB.

###### Proposition 4.1

For each $v_{i}$ generated by gGKB, it holds that $v_{i}\in\mathcal{R}(M)$.
The group of vectors $\\{v_{i}\\}_{i=1}^{k}$ is an $M$-orthonormal basis of
the Krylov subspace

$\mathcal{K}_{k}(M^{{\dagger}}A^{\top}GA,M^{{\dagger}}A^{\top}Gb)=\mathrm{span}\\{(M^{{\dagger}}A^{\top}GA)^{i}M^{{\dagger}}A^{\top}Gb\\}_{i=0}^{k-1},$
(4.4)

and $\\{u_{i}\\}_{i=1}^{k}$ is a $G$-orthonormal basis of the Krylov subspace

$\mathcal{K}_{k}(AM^{{\dagger}}A^{\top}G,b)=\mathrm{span}\\{(AM^{{\dagger}}A^{\top}G)^{i}b\\}_{i=0}^{k-1}.$
(4.5)

###### Proof.

To get a better understanding of gGKB, we give two proofs.

#### The first proof.

We prove $v_{i}\in\mathcal{R}(M)$ by mathematical induction. First note
$\mathcal{R}(M^{{\dagger}})=\mathcal{R}(M)$ for the symmetric $M$. For $i=1$,
we have $\alpha_{1}=M^{{\dagger}}A^{\top}Gu_{1}\in\mathcal{R}(M)$. Suppose
$v_{i}\in\mathcal{R}(M)$ for $i\geq 1$. From the recursion (4.2) we obtain

$\alpha_{i+1}v_{i+1}=M^{{\dagger}}A^{\top}Gu_{i}-\beta_{i+1}v_{i}\in\mathcal{R}(M),$

leading to $v_{i+1}\in\mathcal{R}(M)$. To prove the second property, we
exploit the theory about the GKB of
$\mathcal{A}:(\mathcal{R}(M),\langle\cdot,\cdot\rangle_{M})\rightarrow(\mathbb{R}^{m},\langle\cdot,\cdot\rangle_{G})$
with starting vector $b$, which states that $\\{v_{i}\\}_{i=1}^{k}$ and
$\\{u_{i}\\}_{i=1}^{k}$ are $M$-orthonormal basis and $G$-orthonormal basis of
the two Krylov subspaces

$\displaystyle\mathcal{K}_{k}(\mathcal{A}^{*}\mathcal{A},\mathcal{A}^{*}b):=\mathrm{span}\\{(\mathcal{A}^{*}\mathcal{A})^{i}\mathcal{A}^{*}b\\}_{i=0}^{k-1},$
$\displaystyle\mathcal{K}_{k}(\mathcal{A}\mathcal{A}^{*},b):=\mathrm{span}\\{(\mathcal{A}\mathcal{A}^{*})^{i}b\\}_{i=0}^{k-1},$

respectively. For any $u\in\mathbb{R}^{m}=\pi_{2}(y)=Zy$, from the commutative
diagram (3.7) and using Lemma 3.2 and Lemma 3.3, we obtain

$\displaystyle\mathcal{A}^{*}u=\mathcal{A}^{*}\circ\pi_{2}(y)=\pi_{1}\circ[\mathcal{A}^{*}](y)=W_{r}(W_{r}^{\top}A^{\top}GZ)y=M^{{\dagger}}A^{\top}Gu.$

Therefore, we have

$(\mathcal{A}^{*}\mathcal{A})^{i}\mathcal{A}^{*}b=(M^{{\dagger}}A^{\top}GA)^{i}M^{{\dagger}}A^{\top}Gb,\
\ \ (\mathcal{A}\mathcal{A}^{*})^{i}b=(AM^{{\dagger}}A^{\top}G)^{i}b.$

#### The second proof.

Using the coordinates of $u_{i}$ and $v_{i}$ under bases $W_{r}$ and $Z$, we
can write the last two relations in (4.1) as

$\displaystyle\alpha_{i}W_{r}x_{i}=\mathcal{A}^{*}Zy_{i}-\beta_{i}W_{r}x_{i-1},\
\ \ \beta_{i+1}Zy_{i+1}=\mathcal{A}W_{r}x_{i}-\alpha_{i}Zy_{i},$

where $v_{i}=W_{r}x_{i}$ and $u_{i}=Zy_{i}$. Note that $Z^{\top}GZ=I_{m}$
implies $Z^{-1}=Z^{\top}G$. Letting $\bar{b}=Z^{-1}b$, multiplying from left
the first and second relations by $W_{r}^{\top}M$ and $Z^{-1}$, and using
(3.10), we obtain

$\begin{cases}\beta_{1}y_{1}=\bar{b}_{1}\\\
\alpha_{i}x_{i}=W_{r}^{\top}A^{\top}GZy_{i}-\beta_{i}x_{i-1}\\\
\beta_{i+1}y_{i+1}=Z^{\top}GAW_{r}x_{i}-\alpha_{i}y_{i}.\end{cases}$ (4.6)

Since $G^{\top}=G$, if follows that (4.6) is the standard GKB of matrix
$Z^{\top}GAW_{r}$ with starting vector $\bar{b}$. Therefore,
$\\{x_{i}\\}_{i=1}^{k}$ and $\\{y_{i}\\}_{i=1}^{k}$ are 2-orthonormal bases of
the two Krylov subspaces

$\displaystyle\mathrm{span}\\{((Z^{\top}GAW_{r})^{\top}Z^{\top}GAW_{r})^{i}(Z^{\top}GAW_{r})^{\top}\bar{b}\\}_{i=0}^{k-1},$
$\displaystyle\mathrm{span}\\{(Z^{\top}GAW_{r}(Z^{\top}GAW_{r})^{\top})^{i}\bar{b}\\}_{i=0}^{k-1},$

respectively. Note that

$\displaystyle\ \ \ \
W_{r}((Z^{\top}GAW_{r})^{\top}Z^{\top}GAW_{r})^{i}(Z^{\top}GAW_{r})^{\top}\bar{b}$
$\displaystyle=W_{r}(W_{r}^{\top}A^{\top}GZZ^{\top}GAW_{r})^{i}W_{r}^{\top}A^{\top}GZ\bar{b}$
$\displaystyle=W_{r}(W_{r}^{\top}A^{\top}GAW_{r})^{i}W_{r}^{\top}A^{\top}Gb$
$\displaystyle=(W_{r}W_{r}^{\top}A^{T}GA)^{i}W_{r}W_{r}^{\top}A^{\top}Gb$
$\displaystyle=(M^{{\dagger}}A^{\top}GA)^{i}M^{{\dagger}}A^{\top}Gb.$

We have $v_{i}=W_{r}x_{i}\in\mathcal{R}(M)$ and obtain (4.4). Similarly, we
can obtain (4.5). ∎

It is easy to verify that (4.6) is equivalent to (4.1). Note from (3.2) that
$[\mathcal{A}^{*}]=[\mathcal{A}]^{\top}$ since $Z^{-1}=Z^{\top}G$. Therefore,
the matrix representations of $\mathcal{A}$ and $\mathcal{A}^{*}$ are
$Z^{\top}GAW_{r}\in\mathbb{R}^{m\times r}$ and $(Z^{\top}GAW_{r})^{\top}$,
respectively, which maps a coordinate vector from $\mathbb{R}^{r}$ to
$\mathbb{R}^{m}$. In this sense, we can say that the recursive relation (4.6)
is the coordinate representation for the gGKB of $\mathcal{A}$ under bases
$W_{r}$ and $Z$.

From the first proof of Proposition 4.1, we have $s\in\mathcal{R}(M)$.
Therefore, for each $i\geq 1$, if
$s=M^{{\dagger}}A^{\top}Gu_{i}-\beta_{i}v_{i-1}\neq\boldsymbol{0}$, then
$\alpha_{i}=\|s\|_{M}\neq 0$. This indicates that even if $M$ is not positive
definite, the gGKB does not terminate as long as $s$ or
$q=Av_{i}-\alpha_{i}u_{i}$ is nonzero. Here “terminate” means that
$\alpha_{i}$ or $\beta_{i}$ equals zero at the current step. Suppose gGKB does
not terminate before the $k$-th iteration, i.e. $\alpha_{i}\beta_{i}\neq 0$
for $1\leq i\leq k$. Then the $k$-step gGKB process generates an
$M$-orthonormal matrix
$V_{k+1}=(v_{1},\dots,v_{k+1})\in\mathbb{R}^{n\times(k+1)}$ and a
$G$-orthonormal matrix
$U_{k+1}=(u_{1},\dots,u_{k+1})\in\mathbb{R}^{m\times(k+1)}$, satisfying the
relations

$\displaystyle\beta_{1}U_{k+1}e_{1}=b,$ (4.7a) $\displaystyle
AV_{k}=U_{k+1}B_{k},$ (4.7b) $\displaystyle
M^{{\dagger}}A^{\top}GU_{k+1}=V_{k}B_{k}^{T}+\alpha_{k+1}v_{k+1}e_{k+1}^{\top},$
(4.7c)

where $e_{1}$ and $e_{k+1}$ are the first and $(k+1)$-th columns of $I_{k+1}$,
and

$B_{k}=\begin{pmatrix}\alpha_{1}&&&\\\ \beta_{2}&\alpha_{2}&&\\\
&\beta_{3}&\ddots&\\\ &&\ddots&\alpha_{k}\\\
&&&\beta_{k+1}\end{pmatrix}\in\mathbb{R}^{(k+1)\times k}$ (4.8)

has full column rank. Note that it may occurs that $\beta_{k+1}=0$, which
means gGKB terminates just at the $k$-th step, and in this case
$v_{k+1}=\mathbf{0}$.

We emphasize that gGKB will eventually terminate at most $\min\\{m,r\\}$
steps, since the column rank of $U_{k}$ or $V_{k}$ can not exceed
$\min\\{m,r\\}$. Using the GSVD of $\\{A,L\\}$, we can give a detailed
description of the “terminate step” of gGKB, defined as

$k_{t}=\min\\{k:\alpha_{k+1}\beta_{k+1}=0\\}.$ (4.9)

For any closed subspace $\mathcal{G}$ of a Hilbert space, denote by
$\mathcal{P}_{\mathcal{G}}$ the projection operator onto $\mathcal{G}$. We
have the following result.

###### Theorem 4.1

Define the linear operator $\mathcal{A}$ as (3.11) with
$M=A^{\top}A+L^{\top}L$, where the GSVD of $\\{A,L\\}$ is as (1.1). Suppose
there are $l$ distinct positive $c_{i}$ in $\Sigma_{A}$ with $l$ subspaces
$\mathcal{G}_{1},\dots,\mathcal{G}_{l}$ spanned by the corresponding
$p_{A,i}$. Then $k_{t}$ equals to the number of nonzero elements in
$\\{\mathcal{P}_{\mathcal{G}_{1}}b,\dots,\mathcal{P}_{\mathcal{G}_{l}}b\\}$.

The following lemma is needed to prove this theorem.

###### Lemma 4.1

For any square matrix $C$ and a vector $v$, define the degree of $v$ with
respect to $C$ as

$\mathrm{deg}_{C}(v)=\min\\{k:\exists\ p\in\mathcal{P}_{k}\ \mathrm{s.t.}\
p(C)v=\mathbf{0}\\},$

where $\mathcal{P}_{k}$ is the set of all polynomials with degrees not bigger
than $k$. Then we have

$\mathrm{deg}_{AM^{{\dagger}}A^{\top}}(b)=\mathrm{deg}_{M^{{\dagger}}A^{\top}A}(M^{{\dagger}}A^{\top}b)=l.$
(4.10)

###### Proof.

First notice that $\mathrm{deg}_{C}(v)$ is nothing but the maximum rank of
$\\{C^{i}v\\}_{i=0}^{k}$ with respect to $k\geq 0$. By Theorem 3.2, we have
$AM^{{\dagger}}A^{\top}=A\widetilde{X}_{1}\widetilde{X}_{1}^{\top}A^{\top}$.
Using the relation $A\widetilde{X}_{1}=P_{A}\Sigma_{A}$, we have
$AM^{{\dagger}}A^{\top}=P_{A}(\Sigma_{A}\Sigma_{A}^{\top})P_{A}^{\top}$, which
is the eigenvalue decomposition of $AM^{{\dagger}}A^{\top}$. Thus, the
positive eigenvalues of $AM^{{\dagger}}A^{\top}$ are
$1,c_{q_{1}+1}^{2},\dots,c_{q_{1}+q_{2}}^{2}$ with the corresponding
eigenvectors be the columns of $(P_{A1}\ P_{A2})$, and the corresponding
eigenspaces are subspaces $\mathcal{G}_{1},\dots,\mathcal{G}_{l}$. Denote the
$l$ distinct positive eigenvalues by $\lambda_{1},\dots,\lambda_{l}$ and let
$G_{i}$ be those matrices with $2$-orthonormal columns spanning
$\mathcal{G}_{i}$ for $1\leq i\leq l$. Then we can write the eigenvalue
decomposition of $AM^{{\dagger}}A^{\top}$ as
$AM^{{\dagger}}A^{\top}=\sum_{i=1}^{l}\lambda_{i}G_{i}G_{i}^{\top}$, and we
have $\mathcal{P}_{\mathcal{G}_{i}}=G_{i}G_{i}^{\top}$. For each $j\geq 0$, it
follows that

$w_{j}:=(AM^{{\dagger}}A^{\top})^{j}b=\sum_{i=1}^{l}(\lambda_{i}G_{i}G_{i}^{\top})^{j}b=\sum_{i=1}^{l}\lambda_{i}^{j}G_{i}G_{i}^{\top}b,$

since $G_{i}$ are mutually $2$-orthogonal. Without loss of generality, suppose
there are $s$ nonzero elements in
$\\{\mathcal{P}_{\mathcal{G}_{1}}b,\dots,\mathcal{P}_{\mathcal{G}_{l}}b\\}$
and
$g_{i}=\mathcal{P}_{\mathcal{G}_{i}}b/\|\mathcal{P}_{\mathcal{G}_{i}}b\|_{2}\neq\mathbf{0}$
for $1\leq i\leq s$. Then $\\{g_{i}\\}_{i=1}^{s}$ are mutually 2-orthogonal,
and

$w_{j}=\sum_{i=1}^{s}\lambda_{i}^{j}g_{i}\|\mathcal{P}_{\mathcal{G}_{i}}b\|_{2}=\sum_{i=1}^{s}\lambda_{i}^{j}g_{i}(g_{i}^{\top}b),$

since

$g_{i}^{\top}b=(G_{i}G_{i}^{\top}b)^{\top}b/\|\mathcal{P}_{\mathcal{G}_{i}}b\|_{2}=\|G_{i}^{\top}b\|_{2}^{2}/\|G_{i}^{\top}b\|_{2}=\|\mathcal{P}_{\mathcal{G}_{i}}b\|_{2}.$

Thus, the rank of $\\{w_{j}\\}_{j=0}^{k}$ is at most $s$, leading to
$\mathrm{deg}_{AM^{{\dagger}}A^{\top}}(b)\leq s$. On the other hand, for
$1\leq k\leq s$, by setting $\bar{w}_{j}:=g_{j}(g_{j}^{\top}b)$, we have
$(w_{1}\dots,w_{k})=(\bar{w}_{1},\dots,\bar{w}_{s})T_{k}$, where

$T_{k}=\begin{pmatrix}1&\lambda_{1}&\cdots&\lambda_{1}^{k-1}\\\
1&\lambda_{2}&\cdots&\lambda_{2}^{k-1}\\\ \vdots&\vdots&\cdots&\vdots\\\
1&\lambda_{s}&\cdots&\lambda_{s}^{k-1}\end{pmatrix}=:\vtop{\hbox{$\left(\vbox{\vbox{\halign{
\hfil$#$\hfil\kern 2\p@\kern\@leftbwd& \thinspace\hfil$#$\hfil&&
\quad\hfil$#$\hfil\cr T_{k1}$\hfil\kern 2.0pt\kern 6.66669pt&k\crcr\kern
2.0pt\cr T_{k2}$\hfil\kern 2.0pt\kern 6.66669pt&\ \ s-k\cr k$\hfil\kern
2.0pt\kern 6.66669pt\crcr} }\kern-39.12221pt}\kern
0.0pt\right)\kern-6.66669pt\kern 0.0pt$}\kern 39.12221pt}\ .$

Since $\lambda_{i}\neq\lambda_{j}$ for $1\leq i\neq j\leq k$, the Vandermonde
matrix $T_{k1}$ is invertible, thereby $T_{k}$ has full column rank.
Therefore, the rank of $\\{w_{i}\\}_{i=1}^{k}$ is $k$ for $1\leq k\leq s$,
leading to $\mathrm{deg}_{AM^{{\dagger}}A^{\top}}(b)\geq s$. This proves
$\mathrm{deg}_{AM^{{\dagger}}A^{\top}}(b)=s$.

To prove $\mathrm{deg}_{M^{{\dagger}}A^{\top}A}(M^{{\dagger}}A^{\top}b)=s$, it
is sufficient to show

$\mathrm{rank}\left(\\{(M^{{\dagger}}A^{\top}A)^{j}M^{{\dagger}}A^{\top}b\\}_{i=0}^{k}\right)=\mathrm{rank}\left(\\{w_{j}\\}_{j=0}^{k}\right)$

for any $k\geq 0$. Notice that

$(M^{{\dagger}}A^{\top}A)^{j}M^{{\dagger}}A^{\top}b=M^{{\dagger}}A^{\top}(AM^{{\dagger}}A^{\top})^{j}b=M^{{\dagger}}A^{\top}w_{i}$

and

$M^{{\dagger}}A^{\top}=\widetilde{X}_{1}(A\widetilde{X}_{1})^{\top}=\widetilde{X}_{1}\Sigma_{A}^{\top}P_{A}^{\top}=(X_{1}\
X_{2}C_{q_{2}})(P_{A1}\ P_{A2})^{\top}.$

Let $\tilde{w}_{j}=(M^{{\dagger}}A^{\top}A)^{j}M^{{\dagger}}A^{\top}b$. It
follows that
$\mathrm{rank}\left(\\{\tilde{w}_{j}\\}_{j=0}^{k}\right)\leq\mathrm{rank}\left(\\{w_{j}\\}_{j=0}^{k}\right)$.
To prove the inverse inequality, suppose $\\{w_{j}\\}_{j=0}^{k}$ are
independent. We only need to show $\\{\tilde{w}_{j}\\}_{j=0}^{k}$ are
independent. If there exist real numbers $\mu_{0},\dots,\mu_{k}$ such that
$\sum_{j=0}^{k}\mu_{j}\tilde{w}_{j}=\mathbf{0}$. Then
$M^{{\dagger}}A^{\top}Wz=\mathbf{0}$, where $W=(w_{0},\dots,w_{k})$ has full
column rank and $z=(\mu_{0},\dots,\mu_{k})^{\top}$. By the expression of
$M^{{\dagger}}A^{\top}$, it follows that
$Wz\in\mathcal{N}(M^{{\dagger}}A^{\top})=\mathcal{R}(P_{A3})$. On the other
hand, from
$AM^{{\dagger}}A^{\top}=P_{A}(\Sigma_{A}\Sigma_{A}^{\top})P_{A}^{\top}$ we get
$Wz\in\mathcal{R}(W)\subseteq\mathcal{R}(AM^{{\dagger}}A^{\top})=\mathcal{R}((P_{A1}\
P_{A2}))$. Since $\mathcal{R}((P_{A1}\
P_{A2}))\cap\mathcal{R}(P_{A3})=\\{\mathbf{0}\\}$, we obtain
$Wz=\mathbf{0}\Rightarrow z=\mathbf{0}$, meaning that
$\\{\tilde{w}_{j}\\}_{j=0}^{k}$ are independent. This completes the proof. ∎

  * Proof of Theorem 4.1.

Suppose gGKB terminates at the $k_{t}$-th step. By Proposition 4.1, the rank
of $\\{u_{i}\\}_{i=1}^{k_{t}}$ is $k_{t}$, implying
$k_{t}\leq\mathrm{deg}_{AM^{{\dagger}}A^{\top}}(b)=s$ by Lemma 4.1. Then we
show $k_{t}\geq s$. Notice from the relations (4.1) and (4.2) that

$\displaystyle\alpha_{1}\beta_{1}v_{1}=M^{{\dagger}}A^{\top}b,$
$\displaystyle\alpha_{i+1}\beta_{i+1}v_{i+1}=M^{{\dagger}}A^{\top}Av_{i}-(\alpha_{i}^{2}+\beta_{i+1}^{2})v_{i}-\alpha_{i}\beta_{i}v_{i-1}$

for $1\leq i\leq k_{t}$, where we have used

$\displaystyle M^{{\dagger}}A^{\top}Av_{i}$
$\displaystyle=\alpha_{i}M^{{\dagger}}A^{\top}u_{i}+\beta_{i+1}M^{{\dagger}}A^{\top}u_{i+1}$
$\displaystyle=\alpha_{i}(\alpha_{i}v_{i}+\beta_{i}v_{i+1})+\beta_{i+1}(\alpha_{i+1}v_{i+1}+\beta_{i+1}v_{i}).$

Therefore, it follows that

$v_{i+1}=\frac{1}{\alpha_{i+1}\beta_{i+1}}\left(M^{{\dagger}}A^{\top}Av_{i}-(\alpha_{i}^{2}+\beta_{i+1}^{2})v_{i}-\alpha_{i}\beta_{i}v_{i-1}\right)$

for $1\leq i<k_{t}$. Combining with
$v_{1}=\frac{1}{\alpha_{1}\beta_{1}}M^{{\dagger}}A^{\top}b$, the above
recursion leads to

$v_{k+1}=\sum_{i=0}^{k}\xi_{i}(M^{{\dagger}}A^{\top}A)^{i}M^{{\dagger}}A^{\top}b,\
\ \ \xi_{k}=1/\Pi_{i=1}^{k+1}\alpha_{i}\beta_{i}\neq 0$

for $1\leq k<k_{t}$. Since
$\mathbf{0}=\alpha_{k_{t}+1}\beta_{k_{t}+1}v_{k_{t}+1}$ is a linear
combination of $v_{k_{t}}$ and $v_{k_{t}-1}$ with nonzero coefficients, the
above identity implies that $\alpha_{k_{t}+1}\beta_{k_{t}+1}v_{k_{t}+1}$ must
be a linear combination of
$\\{(M^{{\dagger}}A^{\top}A)^{i}M^{{\dagger}}A^{\top}b\\}_{i=0}^{k_{t}}$ with
nonzero coefficients, thereby
$\\{(M^{{\dagger}}A^{\top}A)^{i}M^{{\dagger}}A^{\top}b\\}_{i=0}^{k_{t}}$ is
linear dependent. By Lemma 4.1, it follows that $k_{t}\geq s$. Finally, we
obtain $k_{t}=s$. $\hfill\Box$

Just as the standard GKB can be employed to approximate extreme SVD
components, we will utilize gGKB to approximate nontrivial extreme GSVD
components.

##  5 GSVD computation by generalized Golub-Kahan bidiagonalization

We first show that gGKB can be used to approximate the SVE components. Then we
use Theorem 3.2 and Theorem 3.3 to relate these approximations to the
nontrivial GSVD components.

###  5.1 Computing nontrivial GSVD components by gGKB

Suppose gGKB does not terminate before the $k$-th step. Then the compact-form
SVD of $B_{k}$ can be written as

$B_{k}=Y_{k}\Theta_{k}H_{k}^{\top},\ \ \
\Theta_{k}=\mathrm{diag}\left(\theta_{1}^{(k)},\dots,\theta_{k}^{(k)}\right),\
\ \ \theta_{i}^{(k)}>\dots>\theta_{k}^{(k)}>0,$ (5.1)

where
$Y_{k}=\left(y_{1}^{(k)},\dots,y_{k}^{(k)}\right)\in\mathbb{R}^{(k+1)\times
k}$ and
$H_{k}=\left(h_{1}^{(k)},\dots,h_{k}^{(k)}\right)\in\mathbb{R}^{k\times k}$
are two 2-orthonormal matrices. The approximation to the SVE triplet
$\left(c_{i},p_{A,i},x_{i}\right)$ of $\mathcal{A}$ is defined as

$\left(\bar{c}_{i}^{(k)},\bar{p}_{A,i}^{(k)},\bar{x}_{i}^{(k)}\right):=\left(\theta_{i}^{(k)},U_{k+1}y_{i}^{(k)},V_{k}h_{i}^{(k)}\right).$
(5.2)

To measure the quality of this approximation, we give the following result.

###### Theorem 5.1

The approximate SVE triplet for $\mathcal{A}$ satisfies

$\displaystyle\mathcal{A}\bar{x}_{i}^{(k)}-\bar{c}_{i}^{(k)}\bar{p}_{A,i}^{(k)}$
$\displaystyle=0,$ (5.3a)
$\displaystyle\mathcal{A}^{*}\bar{p}_{A,i}^{(k)}-\bar{c}_{i}^{(k)}\bar{x}_{i}^{(k)}$
$\displaystyle=\alpha_{k+1}v_{k+1}e_{k+1}^{\top}y_{i}^{(k)}.$ (5.3b)

###### Proof.

Note that $\mathcal{A}v=Av$. The first relation can be verified using (4.7b):

$\mathcal{A}\bar{v}_{i}^{(k)}-\bar{\sigma}_{i}^{(k)}\bar{u}_{i}^{(k)}=AV_{k}h_{i}^{(k)}-\theta_{i}^{(k)}U_{k+1}y_{i}^{(k)}=U_{k+1}\left(B_{k}h_{i}^{(k)}-\theta_{i}^{(k)}y_{i}^{(k)}\right)=0.$

For the second relation, using (3.21) that
$\mathcal{A}^{*}u=M^{{\dagger}}A^{\top}u$, we obtain from (4.7c) that

$\displaystyle\mathcal{A}^{*}\bar{p}_{A,i}^{(k)}-\bar{c}_{i}^{(k)}\bar{x}_{i}^{(k)}$
$\displaystyle=M^{{\dagger}}A^{\top}U_{k+1}y_{i}^{(k)}-\theta_{i}^{(k)}V_{k}h_{i}^{(k)}$
$\displaystyle=\left(V_{k}B_{k}^{\top}+\alpha_{k+1}v_{k+1}e_{k+1}^{\top}\right)y_{i}^{(k)}-\theta_{i}^{(k)}V_{k}h_{i}^{(k)}$
$\displaystyle=V_{k}(B_{k}^{\top}y_{i}^{(k)}-\theta_{i}^{(k)}h_{i}^{(k)})+\alpha_{k+1}v_{k+1}e_{k+1}^{\top}y_{i}^{(k)}$
$\displaystyle=\alpha_{k+1}v_{k+1}e_{k+1}^{\top}y_{i}^{(k)}.$

The proof is completed. ∎

Therefore, the triplet
$\left(\bar{\sigma}_{i}^{(k)},\bar{p}_{A,i}^{(k)},\bar{x}_{i}^{(k)}\right)$
can be accepted as a satisfied SVE triplet at the iteration that
$\left|\alpha_{k+1}v_{k+1}e_{k+1}^{\top}y_{i}^{(k)}\right|$ is sufficiently
small. Using the connection between the SVE of $\mathcal{A}$ and the GSVD of
$\\{A,L\\}$ revealed by Theorem 3.2, the tuple

$\left(\bar{c}_{i}^{(k)},\bar{s}_{i}^{(k)},\bar{p}_{A,i}^{(k)},\bar{x}_{i}^{(k)}\right):=\left(\theta_{i}^{(k)},(1-(\theta_{i}^{(k)})^{2})^{1/2},U_{k+1}y_{i}^{(k)},V_{k}h_{i}^{(k)}\right)$
(5.4)

can be used as a good approximation to a GSVD component. To further measure
the quality of this approximation, note from (2.2) that

$s_{i}^{2}A^{\top}Ax_{i}=c_{i}^{2}L^{\top}Lx_{i},\ \ \ 1\leq i\leq r.$ (5.5)

This is a well-known basic relation for GSVD, which indicates that the
nontrivial generalized eigenvalues of the generalized eigenvalue problem
$A^{T}Ax=\lambda L^{T}Lx$ are $\\{\gamma_{i}^{2}\\}_{i=1}^{r}$ and the
corresponding generalized eigenvectors are $\\{x_{i}\\}_{i=1}^{r}$ [19, §8.7].
Now we can give the following result.

###### Theorem 5.2

The above approximate GSVD tuple for $\\{A,L\\}$ satisfies

$(\bar{s}_{i}^{(k)})^{2}A^{\top}A\bar{x}_{i}^{(k)}-(\bar{c}_{i}^{(k)})^{2}L^{\top}L\bar{x}_{i}^{(k)}=\alpha_{k+1}\beta_{k+1}Mv_{k+1}e_{k}^{\top}h_{i}^{(k)}.$
(5.6)

###### Proof.

First notice from (5.5) that

$\displaystyle
A^{\top}Ax_{i}=(c_{i}^{2}+s_{i}^{2})A^{\top}Ax_{i}=c_{i}^{2}(A^{\top}A+L^{\top}L)x_{i}=c_{i}^{2}Mx_{i},$
$\displaystyle
L^{\top}Lx_{i}=(c_{i}^{2}+s_{i}^{2})L^{\top}Lx_{i}=s_{i}^{2}(A^{\top}A+L^{\top}L)x_{i}=s_{i}^{2}Mx_{i}$

for $1\leq i\leq r$. Since $\widetilde{X}_{1}=(x_{1},\dots,x_{r})$ is an
$M$-orthonormal basis of $(\mathcal{R}(M),\langle\cdot,\cdot\rangle_{M})$, it
follows that $A^{\top}A\mathcal{R}(M)\subseteq\mathcal{R}(M)$ and
$L^{\top}L\mathcal{R}(M)\subseteq\mathcal{R}(M)$. Therefore, we have
$A^{\top}A\bar{x}_{i}^{(k)},L^{\top}L\bar{x}_{i}^{(k)}\in\mathcal{R}(M)$ due
to $\bar{x}_{i}^{(k)}=V_{k}h_{i}^{(k)}\in\mathcal{R}(M)$. By Theorem 5.1, we
have

$\displaystyle\ \ \ \
M^{{\dagger}}[(\bar{s}_{i}^{(k)})^{2}A^{\top}A\bar{x}_{i}^{(k)}-(\bar{c}_{i}^{(k)})^{2}L^{\top}L\bar{x}_{i}^{(k)}]=M^{{\dagger}}[A^{\top}A\bar{x}_{i}^{(k)}-(\bar{c}_{i}^{(k)})^{2}M\bar{x}_{i}^{(k)}]$
$\displaystyle=\theta_{i}^{(k)}M^{{\dagger}}A^{\top}U_{k+1}y_{i}^{(k)}-(\theta_{i}^{(k)})^{2}V_{k}h_{i}^{(k)}$
$\displaystyle=\theta_{i}^{(k)}\left(V_{k}B_{k}^{\top}+\alpha_{k+1}v_{k+1}e_{k+1}^{\top}\right)y_{i}^{(k)}-(\theta_{i}^{(k)})^{2}V_{k}h_{i}^{(k)}$
$\displaystyle=\alpha_{k+1}\beta_{k+1}v_{k+1}e_{k}^{\top}h_{i}^{(k)},$

where we have used $B_{k}^{\top}y_{i}^{(k)}=\theta_{i}^{(k)}h_{i}^{(k)}$ and
$B_{k}h_{i}^{(k)}=\theta_{i}^{(k)}y_{i}^{(k)}$. Multiplying the above equality
by $M$ and using $\mathcal{P}_{\mathcal{R}(M)}=MM^{{\dagger}}$, we finally
obtain (5.6). ∎

Combining Theorems 5.1 and 5.2, it is more proper to use the residual norm

$\|r_{i}^{(k)}\|_{2}:=\left(\|A\bar{x}_{i}^{(k)}-\bar{c}_{i}^{(k)}\bar{p}_{A,i}^{(k)}\|_{2}^{2}+\|(\bar{s}_{i}^{(k)})^{2}A^{\top}A\bar{x}_{i}^{(k)}-(\bar{c}_{i}^{(k)})^{2}L^{\top}L\bar{x}_{i}^{(k)}\|_{2}^{2}\right)^{1/2}$
(5.7)

to measure the quality of the approximate GSVD components of $A$. Since
$\|v_{k+1}\|_{M}=1$, it follows from (5.6) that

$\|r_{i}^{(k)}\|_{2}/\|(A^{\top},L^{\top})^{\top}\|_{2}\leq\alpha_{k+1}\beta_{k+1}|e_{k}^{\top}h_{i}^{(k)}|,$
(5.8)

because $\|M\|_{2}^{1/2}=\|(A^{\top},L^{\top})^{\top}\|_{2}$. The easily
computed quantity $\alpha_{k+1}\beta_{k+1}|e_{k}^{\top}h_{i}^{(k)}|$ is an
upper bound of the scaling-invariant relative residual norm
$\|r_{i}^{(k)}\|_{2}/\|(A^{\top},L^{\top})^{\top}\|_{2}$, which can be used in
a stopping criterion.

We present the pseudocode of the gGKB-based GSVD computation (computing the
GSVD components of $A$) in Algorithm 2. We remark that in order to approximate
the GSVD components of $L$, the gGKB of $\mathcal{L}$ should be used; the
spirit is the same as that for $\mathcal{A}$ and we omit it.

Algorithm 2 The gGKB-based GSVD computation (gGKB_GSVD)

Input: $A\in\mathbb{R}^{m\times n}$, $L\in\mathbb{R}^{p\times n}$,
$\textrm{tol}>0$

1:Initialize: choose a nonzero $b\in\mathbb{R}^{m}$; form
$M=A^{\top}A+L^{\top}L$

2:Compute $\beta_{1},\ \alpha_{1},\ u_{1},\ v_{1}$ by gGKB

3:for $i=1,2,\dots,k,$ do

4: Compute $\beta_{k+1},\ \alpha_{k+1},\ u_{k+1},\ v_{k+1}$ by gGKB; form
$B_{k},\ U_{k+1}$ and $V_{k}$

5: Compute the SVD of $B_{k}$ as (5.1)

6: if $\alpha_{k+1}\beta_{k+1}|e_{k}^{\top}h_{i}^{(k)}|<\mathrm{tol}$ then

7: Compute
$\left(\bar{c}_{i}^{(k)},\bar{s}_{i}^{(k)},\bar{p}_{A,i}^{(k)},\bar{x}_{i}^{(k)}\right):=\left(\theta_{i}^{(k)},(1-(\theta_{i}^{(k)})^{2})^{\frac{1}{2}},U_{k+1}y_{i}^{(k)},V_{k}h_{i}^{(k)}\right)$

8: end if

9:end for

Output: Approximate GSVD components
$\left(\bar{c}_{i}^{(k)},\bar{s}_{i}^{(k)},\bar{p}_{A,i}^{(k)},\bar{x}_{i}^{(k)}\right)$

###  5.2 Convergence and accuracy

We provide preliminary results about the convergence and accuracy of gGKB_GSVD
for GSVD computation. The following result demonstrates the good property of
gGKB_GSVD at the terminate step.

###### Theorem 5.3

Following the notations and assumptions of Theorem 4.1, then at the $k_{t}$-th
step, the gGKB_GSVD algorithm computes exactly $k_{t}$ GSVD components
corresponding to the nonzero elements in
$\\{\mathcal{P}_{\mathcal{G}_{1}}b,\dots,\mathcal{P}_{\mathcal{G}_{l}}b\\}$.

###### Proof.

By Theorems 5.1 and 5.2, at the terminate step of gGKB, it computes exact the
SVE components of $\mathcal{A}$, which are also the exact GSVD components of
$\\{A,L\\}$ by Theorem 3.2. Following the notations in the proof of Lemma 4.1,
we need to prove that the $s$ vectors $\bar{p}_{A,i}^{(s)}$ belong separately
to the invariant subspaces $\mathcal{G}_{1},\dots,\mathcal{G}_{s}$. Since
$\theta_{i}^{(s)}>0$ have different values and $\mathcal{G}_{i}$ are mutually
2-orthogonal, these $\bar{p}_{A,i}^{(s)}$ must belong to different invariant
subspaces. Therefore, we only need to prove
$\mathcal{P}_{\mathcal{G}}\bar{p}_{A,i}^{(s)}=\bar{p}_{A,i}^{(s)}$ for each
$1\leq i\leq s$, where
$\mathcal{G}=\mathcal{G}_{1}\oplus\cdots\oplus\mathcal{G}_{s}$. From the proof
of Lemma 4.1, we have
$\bar{p}_{A,i}^{(s)}\in\mathcal{K}_{s}(AM^{{\dagger}}A^{\top},b)=\mathrm{span}\\{w_{i}\\}_{i=0}^{s-1}$,
and

$\widetilde{W}_{s}:=(w_{0},\dots,w_{s-1})=\widetilde{G}_{s}\begin{pmatrix}g_{1}^{\top}b&&\\\
&\ddots&\\\
&&g_{1}^{\top}b\end{pmatrix}\begin{pmatrix}1&\lambda_{1}&\cdots&\lambda_{1}^{s-1}\\\
1&\lambda_{2}&\cdots&\lambda_{2}^{s-1}\\\ \vdots&\vdots&\cdots&\vdots\\\
1&\lambda_{s}&\cdots&\lambda_{s}^{s-1}\end{pmatrix}=:\widetilde{G}_{s}\Lambda_{s}T_{s},$

where $\widetilde{G}_{s}=(g_{1},\dots,g_{s})$. Since $g_{i}^{\top}b\neq 0$ and
$T_{s}$ is nonsingular, it follows that
$\mathcal{R}(\widetilde{W}_{s})=\mathcal{R}(\widetilde{G}_{s})$. Thus, we can
write $\bar{p}_{A,i}^{(s)}$ as $\bar{p}_{A,i}^{(s)}=\widetilde{G}_{s}z$ with a
nonzero $z\in\mathbb{R}^{s}$. Now we immediately obtain

$\mathcal{P}_{\mathcal{G}}\bar{p}_{A,i}^{(s)}=\sum_{i=1}^{s}\mathcal{P}_{\mathcal{G}_{i}}\bar{p}_{A,i}^{(s)}=\sum_{i=1}^{s}G_{i}G_{i}^{\top}(\widetilde{G}_{s}z)=(\widetilde{G}_{s}\widetilde{G}_{s}^{\top})\widetilde{G}_{s}z=\bar{p}_{A,i}^{(s)},$

which is the desired result. ∎

To investigate the convergence behavior of the approximations, we give the
following result to describe the convergence speed of the Ritz values
$\theta_{i}^{(k)}$.

###### Theorem 5.4

For any $1\leq i\leq q_{1}+q_{2}$, let

$\theta_{i}=\arccos\frac{|c_{i}p_{A,i}^{\top}b|}{\|\Sigma_{A}^{\top}P_{A}^{\top}b\|_{2}}.$
(5.9)

Then at the $k$-th iteration of gGKB_GSVD, it holds

$0\leq
c_{i}^{2}-(\theta_{i}^{(k)})^{2}\leq(c_{1}^{2}-c_{r}^{2})\left(\frac{\kappa_{i}^{(k)}\tan\theta_{i}}{C_{k-i}(1+2\gamma_{i})}\right),$
(5.10)

where $C_{j}(t)$ is the $j$-th Chebyshev polynomial

$C_{j}(t)=\frac{1}{2}\left[(t+\sqrt{t^{2}-1})^{k}+(t+\sqrt{t^{2}-1})^{-k}\right],\
\ |t|\geq 1$

and

$\gamma_{i}=\frac{c_{i}^{2}-c_{i+1}^{2}}{c_{i+1}^{2}-c_{r}^{2}},\ \ \
\kappa_{i}^{(k)}=\prod_{j=1}^{i-1}\frac{(\theta_{j}^{(k)})^{2}-c_{r}^{2}}{(\theta_{j}^{(k)})^{2}-c_{i}^{2}}\
(i>1),\ \ \ \kappa_{1}^{(k)}=1\ (i=1).$

###### Proof.

Using the relations (4.6) with $G=Z=I_{m}$ and $W_{r}=\widetilde{X}_{1}$, it
follows that the coordinate representation of gGKB is the standard GKB of
$A\widetilde{X}_{1}$ with starting vector $b$. This GKB process is equivalent
to the symmetric Lanczos process of
$(A\widetilde{X}_{1})^{\top}A\widetilde{X}_{1}\in\mathbb{R}^{r\times r}$ with
starting vector $(A\widetilde{X}_{1})^{\top}b$, which generates 2-orthogonal
vectors $\\{u_{i}\\}_{i=1}^{k}$ and the symmetric tridiagonal matrix
$B_{k}^{\top}B_{k}$; see e.g. [30]. Since
$A\widetilde{X}_{1}=P_{A}\Sigma_{A}$, it follows that
$(A\widetilde{X}_{1})^{\top}A\widetilde{X}_{1}=I_{r}(\Sigma_{A}^{\top}\Sigma_{A})I_{r}^{\top}$
is the eigenvalue decomposition of
$(A\widetilde{X}_{1})^{\top}A\widetilde{X}_{1}$, and
$(A\widetilde{X}_{1})^{\top}b=\Sigma_{A}^{\top}P_{A}^{\top}b$. Since the
$i$-th eigenvector of $(A\widetilde{X}_{1})^{\top}A\widetilde{X}_{1}$ is
$e_{i}$ and

$(A\widetilde{X}_{1})^{\top}b=\Sigma_{A}^{\top}P_{A}^{\top}b=\begin{pmatrix}(P_{A1}^{\top}b)^{\top}&(P_{A2}^{\top}b)^{\top}&\mathbf{0}\end{pmatrix}^{\top},$

the angle between $(A\widetilde{X}_{1})^{\top}b$ and $e_{i}$ for
$q_{1}+q_{2}+1\leq i\leq r$ is $\pi/2$, and for $1\leq i\leq q_{1}+q_{2}$ the
angle is expressed as (5.9). Notice that the eigenvalues of
$B_{k}^{\top}B_{k}$ are $(\theta_{i}^{(k)})^{2}$. Using the convergence theory
of the symmetric Lanczos process (see e.g. [46, Theorem 6.4]), we immediately
obtain (5.10). ∎

Theorem 5.4 indicates that the convergence rate of $\theta_{i}^{(k)}$
primarily depends on two factors: the clossness between $b$ and the
corresponding vector $p_{A,i}$ and the degree of separation of
$\theta_{i}^{(k)}$ from others. Therefore, usually we can expect rapid
convergence to the extreme and well-separated positive $c_{i}$. Note again
that the approximations will not converge to the GSVD components corresponding
to those zero $c_{i}$, since the angle between $(A\widetilde{X}_{1})^{\top}b$
and $e_{i}$ is $\pi/2$ for $q_{1}+q_{2}+1\leq i\leq r$. The convergence
behavior of $\bar{p}_{A,i}^{(k)}$ and $\bar{x}_{i}^{(k)}$ can also be
described similarly based on the convergence theory of the symmetric Lanczos
process, but the mathematical expressions are more complex. Interested readers
can refer to [46, §6.6]

We remark that all the aforementioned results are derived for the gGKB with
exact computations., i.e. we do not take into account rounding errors and
computational errors arising from iteratively solving (4.3). In the presence
of rounding errors, the Lanczos-type iterative process behaves very
differently from that in exact arithmetic. One well-known result is that the
orthogonality of $u_{i}$ and $v_{i}$ will gradually lost, which leads to a
delay of convergence of approximations and the appearance of spurious
convergent quantities [30]. Also, the inaccurate computation of
$M^{{\dagger}}\bar{s}$ may affect the final accuracy of the approximations.
These issues for gGKB_GSVD will be addressed in future work. We will
demonstrate several of them in the following numerical experiments.

##  6 Experimental results

We report some experimental results to demonstrate the performance of
gGKB_GSVD for computing nontrivial extreme GSVD components. All the
experiments are performed in MATLAB R2023b using double precision. The codes
are available at https://github.com/Machealb/gsvd_iter. For the starting
vector of gGKB for $A$ and $L$, we use the random vector
$b=\texttt{randn(m,1)}$ and $b=\texttt{randn(p,1)}$ with random seed
rng(2024), respectively.

#### Example 1

The matrix pair $\\{A,L\\}$ is constructed as follows. Set $m=n=p=1000$. Let
$C_{A}=\mathrm{diag}(\\{c_{i}\\}_{i=1}^{n})$ with
$c(1)=1,c(2)=0.95,c(3)=0.90$, $c(4:n-6)=\texttt{linspace(0.88,0.12,n-6)}$ and
$c(n-2)=0.1,c(n-1)=0.05,c(n)=0.01$, and let
$S_{L}=\mathrm{diag}(\\{s_{i}\\}_{i=1}^{n})$ with $s_{i}=(1-c_{i}^{2})^{1/2}$.
Then let $W$ be an orthogonal matrix by letting
$W=\texttt{gallery(`orthog',n,2)}$, and $D=\texttt{diag(linspace(1,100,n))}$.
Finally let $A=C_{A}W^{\top}D$ and $L=S_{L}W^{\top}D$. By the construction,
$\\{A,L\\}$ is a regular matrix pair, and the $i$-th generalized singular
value of $\\{A,L\\}$ are $c_{i}/s_{i}$, where the corresponding generalized
singular vectors are the $i$-th columns of $I_{n}$, $I_{n}$ and $D^{-1}W$.

(a)

(b)

(c)

(d)

Figure 6.1: Convergence and accuracy of approximations to $c_{i}$. Top:
convergence of Ritz values $\theta_{i}^{(k)}$ to largest/smallest $c_{i}$ by
gGKB_GSVD without reorthogonalization (left) and with full reorthogonalization
(right). Bottom: error curves (full reorthogonalization).

In this experiment, we use gGKB_GSVD to compute several largest and smallest
generalized singular values, and show the convergence behavior of Ritz values
$\theta_{i}^{(k)}$ where gGKB is performed with and without
reorthogonalization. This is a small-scale problem, thereby we directly
compute $M^{-1}$ for computing $M^{-1}\bar{s}$ at each iteration of gGKB.
Figure 6.1 shows the convergence of the first three largest and smallest Ritz
values, where in the top two subfigures the right vertical lines indicate the
values of $c_{i}$. There are four findings. (1) If no reorthogonalization is
used for gGKB, then as the iteration proceeds, the converged Ritz values may
suddenly jump up (also may jump down for converging to those smallest $c_{i}$)
to become a ghost and then converge to the next larger $c_{i}$; this
phenomenon leads to the appearance of spurious copies of computed $c_{i}$. (2)
If gGKB is performed with full reorthogonalization, the convergence of the
Ritz values remains regular, and the first three largest and smallest Ritz
values converge to the first three largest and smallest $c_{i}$, respectively.
(3) The final accuracy of the approximated $c_{i}$ is around
$\mathcal{O}(\mathbf{u})$, where $\mathbf{u}=2^{-53}\approx 10^{-16}$ is the
roundoff unit of double precision. (4) The convergence to those largest
$c_{i}$ is faster than the convergence to those smallest $c_{i}$; also, the
convergence to $c_{1},c_{2}$ and $c_{n},c_{n-1}$ is faster than that to
$c_{3}$ and $c_{n-2}$. This is because $c_{1},c_{2}$ are more well-separated
from others than $c_{3}$; the same reason applies to $c_{n},c_{n-1}$ and
$c_{n-2}$.

(a)

(b)

(c)

(d)

Figure 6.2: Convergence and accuracy of approximations to $s_{i}$. Top:
convergence of Ritz values $\theta_{i}^{(k)}$ to largest/smallest $s_{i}$ by
gGKB_GSVD without reorthogonalization (left) and with full reorthogonalization
(right). Bottom: error curves (full reorthogonalization).

We also test using the gGKB of $\mathcal{L}$ to approximate several largest
and smallest $s_{i}$. The convergence behavior of the Ritz values and error
curves are plotted in Figure 6.2. In addition to the findings closely
resembling those depicted in Figure 6.1, there are two additional insights.
First, we find that the first three smallest Ritz values converge to
$s_{2},s_{3},s_{4}$ instead of $s_{1},s_{2},s_{3}$. This is because to
$s_{1}=0$, which can not be converged upon by Ritz values, as revealed by
Theorem 5.3 and Theorem 5.4. Therefore, we should use gGKB of $\mathcal{A}$ to
compute the generalized singular values with value $\infty$. Second, we find
from the bottom two subfigures that those smallest $s_{i}$ can be approximated
more quickly than the largest ones, due to their well-separated locations.
Notice that those smallest $s_{i}$ correspond to the largest generalized
singular values. Therefore, for computing GSVD components corresponding to
those smallest generalized singular values, using gGKB_GSVD based on the gGKB
of $\mathcal{L}$ may be a better choice.

#### Example 2

The matrix $A$ named well1850 is taken from the SuiteSparse matrix collection
[13], and the matrix $L$ is set as

$L=\begin{pmatrix}1.1&-1&&&\\\ &1.1&-1&&\\\ &&\ddots&\ddots&\\\ &&&1.1&-1\\\
\end{pmatrix}\in\mathbb{R}^{(n-1)\times n}.$

This is a regular matrix pair. We use the Matlab build-in function gsvd.m to
compute the full GSVD of $\\{A,L\\}$ as the baseline of comparison.

(a)

(b)

Figure 6.3: Error curves of the approximate GSVD components by gGKB_GSVD and
relative residual norm with its upper bound. Left: approximations to the 1-st
GSVD components. Right: approximations to the n-th GSVD components.

In this experiment, we test the performance of gGKB_GSVD for computing the
first and $n$-th GSVD components of $A$ for a matrix pair with nonsquare
matrices. We only show the results for the gGKB of $\mathcal{A}$ and omit the
similar results for the gGKB of $\mathcal{L}$. Full reorthogonalization is
used and $M^{-1}$ is computed directly. The errors for the approximated
generalized singular vectors are measured by $\sin\angle(x,y)$ between two
vectors. We also plot the variation of the relative residual norm and its
upper bound $\alpha_{k+1}\beta_{k+1}|e_{k}^{\top}h_{i}^{(k)}|$. Figure 6.3
shows that gGKB_GSVD can approximate very well the two group extreme GSVD
components, with final accuracy around $\mathcal{O}(\mathbf{u})$. The upper
bound $\alpha_{k+1}\beta_{k+1}|e_{k}^{\top}h_{i}^{(k)}|$ exhibits a nearly
identical decreasing trend as the relative residual norm. Therefore, it is a
highly suitable quantity to be employed in the stopping criterion. We also
observe that the convergence to the first GSVD components is faster than the
convergence to the and $n$-th.

#### Example 3

The matrix pair $\\{A,L\\}$ is constructed as follows. Set $m=n=p=1000$ and
set $r=900$. Let $C_{A}=(\Sigma_{A},\mathbf{0})$ with
$C_{A}=\mathrm{diag}(\\{c_{i}\\}_{i=1}^{r})$, where $c(1)=0.99,c(2)=0.98$,
$c(3:r-4)=\texttt{linspace(0.96,0.06,r-4)}$ and $c(r-1)=0.04,c(r-1)=0.02$, and
let $S_{L}=(\Sigma_{L},\mathbf{0})$ with
$\Sigma_{L}=\mathrm{diag}(\\{s_{i}\\}_{i=1}^{r})$ and
$s_{i}=(1-c_{i}^{2})^{1/2}$. Then let $W=\texttt{gallery(`orthog',n,2)}$ be an
orthogonal matrix and $D=\texttt{diag(linspace(1,10,n))}$. Finally let
$A=C_{A}W^{\top}D$ and $L=S_{L}W^{\top}D$. By the construction, we have
$\mathrm{rank}((A^{\top},L^{\top})^{\top})=r<n$, and the nontrivial GSVD
components are $c_{i},\ s_{i}$ and the $i$-th columns of $I_{n}$, $I_{n}$ and
$D^{-1}W$ for $1\leq i\leq r$. For each nontrivial $x_{i}$, we compute
$\mathcal{P}_{\mathcal{R}(M)}x_{i}=MM^{{\dagger}}x_{i}$ to get the
corresponding right generalized singular vector belonging to $\mathcal{R}(M)$.
We use gGKB_GSVD to compute $x_{i}$ that belongs to $\mathcal{R}(M)$.

(a)

(b)

Figure 6.4: Error curves of the approximate GSVD components by gGKB_GSVD and
relative residual norm with its upper bound, where
$\mathrm{rank}((A^{\top},L^{\top})^{\top})=r<n$. Left: approximations to the
1-st GSVD components. Right: approximations to the r-th GSVD components.

In this experiment, we test the performance of gGKB_GSVD for computing the
first and $r$-th GSVD components for a nonregular matrix pair. We directly
compute $M^{{\dagger}}$ and use full reorthogonalization for gGKB. Figure 6.4
shows very good performance of the algorithm: (1) the two extreme GSVD
components $(c_{i},p_{A,i},x_{i})$ for $i=1,r$ can be approximated with final
accuracy around $\mathcal{O}(\mathbf{u})$; (2) the relative residual norm and
its upper bound $\alpha_{k+1}\beta_{k+1}|e_{k}^{\top}h_{i}^{(k)}|$ follow
nearly identical decreasing curves, with both eventually stabilizing at a
level around $\mathcal{O}(\mathbf{u})$. Again, we observe that the convergence
to the first GSVD components is faster than the convergence to the $r$-th.

(a)

(b)

(c)

(d)

Figure 6.5: Accuracy of computed GSVD components by gGKB_GSVD, where
$s=M^{{\dagger}}\bar{s}$ at each gGKB iteration is computed by solving
$\min_{s}\|Ms-\bar{s}\|_{2}$ using lsqr.m with different stopping tolerance
tol. Top: tol=$10^{-10}$. Top: tol=$10^{-8}$.

#### Example 4

The matrix pair $\\{A,L\\}$ is constructed as follows. Set $m=n=p=5000$. Let
$C_{A}=\mathrm{diag}(\\{c_{i}\\}_{i=1}^{r})$ with $c(1)=0.99,c(2)=0.97$,
$c(3:n-4)=\texttt{linspace(0.95,0.15,n-4)}$ and $c(n-1)=0.1,c(n)=0.05$. Let
$S_{L}=\mathrm{diag}(\\{s_{i}\\}_{i=1}^{n})$ with $s_{i}=(1-c_{i}^{2})^{1/2}$.
Then let $W=\texttt{gallery(`orthog',n,2)}$ be an orthogonal matrix and
$D=\texttt{diag(linspace(1,10,n))}$. Finally let $A=C_{A}W^{\top}D$ and
$L=S_{L}W^{\top}D$. By the construction, $\\{A,L\\}$ is a regular matrix pair,
and the $i$-th GSVD components are $c_{i},\ s_{i}$ and $i$-th columns of
$I_{n}$, $I_{n}$ and $D^{-1}W$.

We use this experiment to demonstrate the impact of inaccuracy in the
computation of $M^{{\dagger}}\bar{s}$ on the final accurate of the approximate
GSVD components. We use the Matlab build-in function lsqr.m to solve (4.3)
iteratively with stopping tolerance $\texttt{tol}=10^{-10},10^{-8}$ at each
iteration of gGKB, respectively. Figure 6.5 shows the decrease of relative
errors of the first and $n$-th approximate GSVD components with the two
stopping tolerances. We observe that the computational accuracy of
$M^{{\dagger}}\bar{s}$ significantly affects the final accuracy of both the
generalized singular values and vectors. As the computational accuracy
deteriorates, so does the final accuracy of the computed GSVD components.
Further theoretical investigation into this issue should be conducted in
future research.

##  7 Conclusion and outlook

Based on the theory of singular value expansion (SVE) of linear compact
operators, we have provided a new understanding of the GSVD of $\\{A,L\\}$
with $A\in\mathbb{R}^{m\times n}$ and $L\in\mathbb{R}^{p\times n}$. By
defining the positive semidefinite matrix $M=A^{\top}A+L^{\top}L$, we have
shown that: (1) the trivial GSVD components $\\{x_{i}\\}$ form a basis for
$\mathcal{N}(M)$ and any nontrivial $x_{i}$ belongs to the coset
$\bar{x}_{i}+\mathcal{N}(M)$, where $\bar{x_{i}}\in\mathcal{R}(M)$ is a
nontrivial GSVD component; (2) the nontrivial GSVD components of $A$ and $L$
are just the SVEs of the linear operators
$\mathcal{A}:(\mathcal{R}(M),\langle\cdot,\cdot\rangle_{M})\rightarrow(\mathbb{R}^{m},\langle\cdot,\cdot\rangle_{2}),\
v\mapsto Av$ and
$\mathcal{L}:(\mathcal{R}(M),\langle\cdot,\cdot\rangle_{M})\rightarrow(\mathbb{R}^{p},\langle\cdot,\cdot\rangle_{2}),\
v\mapsto Lv$, respectively. As a direct application of this result, we have
developed an operator-type Golub-Kahan bidiagonalization (GKB) for
$\mathcal{A}$ and $\mathcal{L}$, leading to a novel generalized GKB (gGKB)
process. We have used the GSVD of $\\{A,L\\}$ to study basic properties of
gGKB and proposed the gGKB_GSVD algorithm to compute several nontrivial
extreme GSVD components of large-scale matrix pairs. Preliminary results about
convergence and accuracy of gGKB_GSVD for GSVD computation have been provided,
and numerical experiments are presented to demonstrate the effectiveness of
this method.

The idea of this paper offers potential directions for developing new
algorithms for large-scale GSVD computation. Note that the SVE of
$\mathcal{A}$ or $\mathcal{L}$ can be treated as a “weighted” SVD, where the
weight matrix $M$ induces a non-Euclidean inner product. Therefore, existing
SVD algorithms based on Krylov subspace projection may be modified to
approximate the SVE and consequently, the nontrivial GSVD components.

## References

  * [1] O. Alter, P. O. Brown, and D. Botstein.  Generalized singular value decomposition for comparative analysis of genome-scale expression data sets of two different organisms.  Proceedings of the National Academy of Sciences, 100(6):3351–3356, 2003.
  * [2] F. Alvarruiz, C. Campos, and J. E. Roman.  Thick-restarted joint Lanczos bidiagonalization for the GSVD.  Journal of Computational and Applied Mathematics, 440:115506, 2024\.
  * [3] M. Arioli.  Generalized Golub–Kahan bidiagonalization and stopping criteria.  SIAM Journal on Matrix Analysis and Applications, 34(2):571–592, 2013.
  * [4] Z. Bai, J. Demmel, J. Dongarra, A. Ruhe, and H. van der Vorst.  Templates for the solution of algebraic eigenvalue problems: a practical guide.  SIAM, 2000.
  * [5] Z. Bai and J. W. Demmel.  Computing the generalized singular value decomposition.  SIAM Journal on Scientific Computing, 14(6):1464–1486, 1993.
  * [6] S. J. Benbow.  Solving generalized least-squares problems with LSQR.  SIAM Journal on Matrix Analysis and Applications, 21(1):166–177, 1999.
  * [7] K. Bhuyan, S. Singh, and P. Bhuyan.  Application of generalized singular value decomposition to ionospheric tomography.  In Annales Geophysicae, volume 22, pages 3437–3444. Copernicus Publications Göttingen, Germany, 2004.
  * [8] Å. Björck.  Numerical Methods for Least Squares Problems.  SIAM, Philadelphia, 1996.
  * [9] N. A. Caruso and P. Novati.  Convergence analysis of LSQR for compact operator equations.  Linear Algebra Appl., 583:146–164, 2019.
  * [10] M. T. Chu, R. E. Funderlic, and G. H. Golub.  On a variational formulation of the generalized singular value decomposition.  SIAM Journal on Matrix Analysis and Applications, 18(4):1082–1092, 1997.
  * [11] J. Chung and A. K. Saibaba.  Generalized hybrid iterative methods for large-scale Bayesian inverse problems.  SIAM Journal on Scientific Computing, 39(5):S24–S46, 2017.
  * [12] J. Chung, A. K. Saibaba, M. Brown, and E. Westman.  Efficient generalized Golub–Kahan based methods for dynamic inverse problems.  Inverse Problems, 34(2):024005, 2018.
  * [13] T. A. Davis and Y. Hu.  The university of florida sparse matrix collection.  ACM Trans. Math. Software (TOMS), 38(1):1–25, 2011.
  * [14] A. Dumitrasc, C. Kruse, and U. Ruede.  Generalized Golub-Kahan bidiagonalization for nonsymmetric saddle point systems.  arXiv preprint arXiv:2310.06952, 2023.
  * [15] A. Edelman and Y. Wang.  The GSVD: Where are the ellipses?, matrix trigonometry, and more.  SIAM Journal on Matrix Analysis and Applications, 41(4):1826–1856, 2020.
  * [16] H. W. Engl, M. Hanke, and A. Neubauer.  Regularization of Inverse Problems.  Kluwer Academic Publishers, 2000.
  * [17] L. M. Ewerbring and F. T. Luk.  Canonical correlations and generalized SVD: applications and new algorithms.  Journal of computational and applied mathematics, 27(1-2):37–52, 1989.
  * [18] G. Golub and W. Kahan.  Calculating the singular values and pseudo-inverse of a matrix.  Journal of the Society for Industrial and Applied Mathematics, Series B: Numerical Analysis, 2(2):205–224, 1965.
  * [19] G. H. Golub and C. F. Van Loan.  Matrix Computations.  The Johns Hopkins University Press, Baltimore, 4th edition, 2013.
  * [20] M. Hochstenbach.  A Jacobi–Davidson type method for the generalized singular value problem.  Linear Algebra and its Applications, 431(3-4):471–487, 2009.
  * [21] M. E. Hochstenbach and Y. Notay.  The Jacobi–Davidson method.  GAMM-Mitteilungen, 29(2):368–382, 2006.
  * [22] P. Howland, M. Jeon, and H. Park.  Structure preserving dimension reduction for clustered text data based on the generalized singular value decomposition.  SIAM Journal on Matrix Analysis and Applications, 25(1):165–179, 2003.
  * [23] J. Huang and Z. Jia.  On choices of formulations of computing the generalized singular value decomposition of a large matrix pair.  Numerical Algorithms, 87:689–718, 2021.
  * [24] J. Huang and Z. Jia.  Two harmonic Jacobi–Davidson methods for computing a partial generalized singular value decomposition of a large matrix pair.  Journal of Scientific Computing, 93(2):41, 2022.
  * [25] Z. Jia and H. Li.  The joint bidiagonalization process with partial reorthogonalization.  Numer. Algor., 88(2):965–992, 2021.
  * [26] Z. Jia and H. Li.  The joint bidiagonalization method for large GSVD computations in finite precision.  SIAM Journal on Matrix Analysis and Applications, 44(1):382–407, 2023.
  * [27] B. Kågström.  The generalized singular value decomposition and the general ($A-\lambda B$)-problem.  BIT Numerical Mathematics, 24:568–583, 1984.
  * [28] R. Kress.  Linear Integral Equations.  Springer New York, NY, 2013.
  * [29] S. Kuo, W. Yeih, and Y. Wu.  Applications of the generalized singular-value decomposition method on the eigenproblem using the incomplete boundary element formulation.  Journal of Sound and Vibration, 235(5):813–845, 2000.
  * [30] R. M. Larsen.  Lanczos bidiagonalization with partial reorthogonalization.  DAIMI Report Series, (537), 1998.
  * [31] H. Li.  A preconditioned Krylov subspace method for linear inverse problems with general-form Tikhonov regularization.  arXiv preprint arXiv:2308.06577, 2023.
  * [32] H. Li.  Subspace projection regularization for large-scale Bayesian linear inverse problems.  arXiv preprint arXiv:2310.18618, 2023.
  * [33] H. Li.  The joint bidiagonalization of a matrix pair with inaccurate inner iterations.  SIAM Journal on Matrix Analysis and Applications, 45(1):232–259, 2024.
  * [34] H. Li, J. Feng, and F. Lu.  Scalable iterative data-adaptive RKHS regularization.  arXiv preprint arXiv:2401.00656, 2024.
  * [35] R.-C. Li.  Bounds on perturbations of generalized singular values and of associated subspaces.  SIAM J. Matrix Anal. Appl., 14(1):195–234, 1993.
  * [36] Y. Liu, X. Shan, and M. Shao.  A contour integral-based algorithm for computing generalized singular values.  arXiv preprint arXiv:2401.00121, 2023.
  * [37] K. Nakamura, K. Nakadai, and G. Ince.  Real-time super-resolution sound source localization for robots.  In 2012 IEEE/RSJ International Conference on Intelligent Robots and Systems, pages 694–699. IEEE, 2012.
  * [38] C. C. Paige.  A note on a result of Sun Ji-Guang: Sensitivity of the CS and GSV decompositions.  SIAM J. Numer. Anal., 21(1):186–191, 1984.
  * [39] C. C. Paige.  The general linear model and the generalized singular value decomposition.  Linear Algebra and its applications, 70:269–284, 1985.
  * [40] C. C. Paige.  Computing the generalized singular value decomposition.  SIAM Journal on Scientific and Statistical Computing, 7(4):1126–1146, 1986.
  * [41] C. C. Paige and M. A. Saunders.  Towards a generalized singular value decomposition.  SIAM J. Numer. Anal., 18(3):398–405, 1981.
  * [42] C. C. Paige and M. A. Saunders.  LSQR: An algorithm for sparse linear equations and sparse least squares.  ACM Trans. Math. Software, 8:43–71, 1982.
  * [43] C. H. Park and H. Park.  A relationship between linear discriminant analysis and the generalized minimum squared error solution.  SIAM Journal on Matrix Analysis and Applications, 27(2):474–492, 2005.
  * [44] P. T. Peter Tang and E. Polizzi.  FEAST as a subspace iteration eigensolver accelerated by approximate spectral projection.  SIAM Journal on Matrix Analysis and Applications, 35(2):354–390, 2014.
  * [45] A. Refahi Sheikhani and S. Kordrostami.  New iterative methods for generalized singular-value problems.  Mathematical Sciences, 11(4):257–265, 2017.
  * [46] Y. Saad.  Numerical methods for large eigenvalue problems: revised edition.  SIAM, 2011.
  * [47] T. Sakurai and H. Sugiura.  A projection method for generalized eigenvalue problems using numerical integration.  Journal of computational and applied mathematics, 159(1):119–128, 2003.
  * [48] J. M. Speiser and C. Van Loan.  Signal processing computations using the generalized singular value decomposition.  In Real-Time Signal Processing VII, volume 495, pages 47–57. SPIE, 1984.
  * [49] J.-G. Sun.  Perturbation analysis for the generalized singular value problem.  SIAM J. Numer. Anal., 20(3):611–625, 1983.
  * [50] J.-g. Sun.  Perturbation theorems for generalized singular values.  J. Comput. Math., 1(3):233–242, 1983.
  * [51] B. D. Sutton.  Stable computation of the CS decomposition: Simultaneous bidiagonalization.  SIAM Journal on Matrix Analysis and Applications, 33(1):1–21, 2012\.
  * [52] C. F. Van Loan.  Generalizing the singular value decomposition.  SIAM J. Numer. Anal., 13(1):76–83, 1976.
  * [53] C. F. Van Loan.  Computing the CS and the generalized singular value decompositions.  Numer. Math., 46(4):479–491, 1985.
  * [54] H. Zha.  Computing the generalized singular values/vectors of large sparse or structured matrix pairs.  Numer. Math., 72(3):391–417, 1996.