Generalized Golub-Kahan Bidiagonalization

• Generalized Golub-Kahan Bidiagonalization is a Krylov subspace method that extends classical bidiagonalization using weighted, non-Euclidean inner products for various inverse problems.
• It constructs paired basis sequences and reduces large-scale operators to low-dimensional bidiagonal or Hessenberg systems for efficient computation.
• Applications include Bayesian inversion, GSVD computation, saddle point systems, and tensor equations, offering robust regularization and noise mitigation.

Generalized Golub–Kahan bidiagonalization (GKB) denotes a family of Krylov subspace processes that extend classical Golub–Kahan bidiagonalization beyond the Euclidean, single-matrix least-squares setting. Across the recent literature, these extensions are used to handle weighted inner products, prior- and noise-weighted inverse problems, generalized singular value decomposition (GSVD) computation, saddle point systems, and tensor equations. A unifying feature is the construction of paired basis sequences and a small projected bidiagonal or Hessenberg problem, so that large-scale structure is accessed through matrix-vector products, weighted normalizations, and low-dimensional regularized subproblems rather than explicit factorization of the full operator (Chung et al., 2016).

1. Operator viewpoint and defining structure

In its most abstract formulation, generalized GKB replaces the Euclidean geometry of standard bidiagonalization by weighted or Hilbert-space inner products. A representative setting introduces a positive semidefinite matrix $M$ on the domain and a matrix $G$ on the data space, so that $A$ is viewed as an operator

$$A : (R(M), \langle \cdot,\cdot\rangle_M) \to (\mathbb{R}^m,\langle \cdot,\cdot\rangle_G),$$

with adjoint

$A^*u = M^\dagger A^T G u.$

The resulting generalized GKB recursion is

$$\beta_1 u_1 = b,\qquad \alpha_i v_i = A^*u_i - \beta_i v_{i-1},\qquad \beta_{i+1}u_{i+1} = Av_i - \alpha_i u_i,$$

with $\|u_i\|_G=\|v_i\|_M=1$ (Li, 2024).

This operator formulation is broad enough to cover several distinct algorithmic traditions. In Bayesian inverse problems, $M$ is induced by the prior covariance and $G$ by the noise covariance. In GSVD computation, $M=A^TA+L^TL$ induces the geometry required for matrix pairs $G$0. In saddle point systems, $G$1 and $G$2 encode blockwise energy norms. A plausible implication is that “generalized GKB” is best understood not as a single fixed algorithm, but as a structural template: alternating primal and dual recurrences under non-Euclidean orthogonality constraints.

A concise taxonomy used in the recent literature is as follows.

Variant	Core modification	Representative use
Weighted/Hilbert-space gGKB	Replace Euclidean adjoint and orthogonality by $G$3- and $G$4-inner products	GSVD components, weighted least squares
gen-GK for Bayesian inversion	Use $G$5- and $G$6-orthogonality	Prior-conditioned inverse problems
Generalized GKB for saddle points	Couple block solves and weighted recurrences	CRAIG/nsCRAIG-type solvers
Tensor GKB variants	Replace matrix multiplication by Einstein product or t-product	Image, video, multilinear inverse problems

2. Core bidiagonalization recurrences and projected problems

A central concrete instance appears in Bayesian inverse problems of the form

$G$7

To avoid explicit computation of $G$8 or $G$9, the variable transformation

$A$0

converts the MAP problem into

$A$1

with $A$2. The generalized GKB process then generates

$A$3

$A$4

together with the orthogonality relations

$A$5

The associated right Krylov space is

$A$6

and each iteration requires only matrix-vector products with $A$7, $A$8, $A$9, and solves with $$A : (R(M), \langle \cdot,\cdot\rangle_M) \to (\mathbb{R}^m,\langle \cdot,\cdot\rangle_G),$$0 (Chung et al., 2016).

The algorithmic purpose of these recurrences is projection. At step $$A : (R(M), \langle \cdot,\cdot\rangle_M) \to (\mathbb{R}^m,\langle \cdot,\cdot\rangle_G),$$1, one seeks $$A : (R(M), \langle \cdot,\cdot\rangle_M) \to (\mathbb{R}^m,\langle \cdot,\cdot\rangle_G),$$2, where $$A : (R(M), \langle \cdot,\cdot\rangle_M) \to (\mathbb{R}^m,\langle \cdot,\cdot\rangle_G),$$3 solves the reduced problem

$$A : (R(M), \langle \cdot,\cdot\rangle_M) \to (\mathbb{R}^m,\langle \cdot,\cdot\rangle_G),$$4

The full-space approximation is then

$$A : (R(M), \langle \cdot,\cdot\rangle_M) \to (\mathbb{R}^m,\langle \cdot,\cdot\rangle_G),$$5

This projected formulation explains why generalized GKB is frequently paired with hybrid regularization: the large-scale operator is compressed into a small bidiagonal model $$A : (R(M), \langle \cdot,\cdot\rangle_M) \to (\mathbb{R}^m,\langle \cdot,\cdot\rangle_G),$$6, while regularization is imposed only after projection.

The same projected-least-squares pattern persists in tensor formulations. In the Einstein-product setting, the global tensor GKB process produces a bidiagonal matrix $$A : (R(M), \langle \cdot,\cdot\rangle_M) \to (\mathbb{R}^m,\langle \cdot,\cdot\rangle_G),$$7 and a reduced Tikhonov problem in $$A : (R(M), \langle \cdot,\cdot\rangle_M) \to (\mathbb{R}^m,\langle \cdot,\cdot\rangle_G),$$8, while in the t-product setting the weighted tensor GKB process yields a reduced tensor system with small bidiagonal tensor $$A : (R(M), \langle \cdot,\cdot\rangle_M) \to (\mathbb{R}^m,\langle \cdot,\cdot\rangle_G),$$9 (Guide et al., 2020, Reichel et al., 2021).

3. Regularization, hybridization, and inverse-problem interpretation

The strongest concentration of generalized GKB applications is in large-scale inverse problems, where the method is used not only as a Krylov solver but as a regularization mechanism. In the generalized hybrid framework, the projected problem supports adaptive parameter choice by criteria such as GCV and the discrepancy principle, and the approach is reported to inherit the usual advantages of hybrid methods, including avoidance of semi-convergence and automatic estimation of the regularization parameter (Chung et al., 2016).

This role becomes clearer when the method is compared with standard LSQR. For the transformed Bayesian problem, gen-LSQR iterates are equivalent to LSQR iterates applied to the prior-conditioned Tikhonov problem

$A^*u = M^\dagger A^T G u.$0

but the generalized GKB formulation achieves this without forming $A^*u = M^\dagger A^T G u.$1 (Chung et al., 2016). This equivalence is operationally significant in settings where $A^*u = M^\dagger A^T G u.$2 is dense, defined on irregular grids, or only accessible through matrix-vector multiplication, including Matérn-class covariances.

The same projection principle underlies later decomposition methods for mixed-prior inverse problems. In the sdHybr/FGGK framework, generalized GKB is combined with flexible Krylov updates so that a solution $A^*u = M^\dagger A^T G u.$3 with smooth and sparse components can be approximated by

$A^*u = M^\dagger A^T G u.$4

and the projected problem incorporates both $A^*u = M^\dagger A^T G u.$5 and $A^*u = M^\dagger A^T G u.$6: $A^*u = M^\dagger A^T G u.$7 Within that reduced model, UPRE, DP, and WGCV are used for adaptive parameter selection (Chung et al., 2022).

The regularizing effect of GKB-based methods has also been analyzed from the residual side. For LSQR, LSMR, and CRAIG, residuals can be written explicitly in terms of bidiagonalization vectors, and the coefficients in those linear combinations quantify propagated noise. For CRAIG, the residual is a multiple of a single bidiagonalization vector, which makes its norm directly interpretable as a regularization indicator (Hnětynková et al., 2016). Complementary LSQR theory for discrete ill-posed problems with multiple singular values shows that, for severely and moderately ill-posed problems, LSQR attains 2-norm filtering best possible regularized solutions at semi-convergence, and that this behavior is governed by the low-rank approximations and Ritz values generated by Golub–Kahan bidiagonalization (Jia, 2020).

4. GSVD, filtered expansions, and generalized spectral approximation

A second major line of development treats generalized GKB as a computational vehicle for GSVD. For a matrix pair $A^*u = M^\dagger A^T G u.$8, the 2024 GSVD study characterizes the GSVD through singular value expansions of two induced operators on finite-dimensional Hilbert spaces and then derives a generalized GKB process tailored to compute nontrivial extreme GSVD components (Li, 2024).

In this setting, with

$A^*u = M^\dagger A^T G u.$9

gGKB generates $$\beta_1 u_1 = b,\qquad \alpha_i v_i = A^*u_i - \beta_i v_{i-1},\qquad \beta_{i+1}u_{i+1} = Av_i - \alpha_i u_i,$$0-orthonormal right basis vectors and $$\beta_1 u_1 = b,\qquad \alpha_i v_i = A^*u_i - \beta_i v_{i-1},\qquad \beta_{i+1}u_{i+1} = Av_i - \alpha_i u_i,$$1-orthonormal left basis vectors, produces a bidiagonal $$\beta_1 u_1 = b,\qquad \alpha_i v_i = A^*u_i - \beta_i v_{i-1},\qquad \beta_{i+1}u_{i+1} = Av_i - \alpha_i u_i,$$2, and then extracts approximate GSVD triplets from the SVD

$$\beta_1 u_1 = b,\qquad \alpha_i v_i = A^*u_i - \beta_i v_{i-1},\qquad \beta_{i+1}u_{i+1} = Av_i - \alpha_i u_i,$$3

The approximations are

$$\beta_1 u_1 = b,\qquad \alpha_i v_i = A^*u_i - \beta_i v_{i-1},\qquad \beta_{i+1}u_{i+1} = Av_i - \alpha_i u_i,$$4

The paper gives residual-based quality certificates and a computable upper bound involving $$\beta_1 u_1 = b,\qquad \alpha_i v_i = A^*u_i - \beta_i v_{i-1},\qquad \beta_{i+1}u_{i+1} = Av_i - \alpha_i u_i,$$5 (Li, 2024).

The inverse-problem literature had already connected generalized GKB to GSVD filtering. For the Bayesian MAP estimator, one can write a filtered GSVD expansion

$$\beta_1 u_1 = b,\qquad \alpha_i v_i = A^*u_i - \beta_i v_{i-1},\qquad \beta_{i+1}u_{i+1} = Av_i - \alpha_i u_i,$$6

and the corresponding gen-LSQR iterates act as truncated or polynomial-filtered GSVD approximants (Chung et al., 2016).

This spectral perspective clarifies an important point: generalized GKB is not merely a weighted least-squares solver. It is also a structure-preserving reduction mechanism for generalized spectral objects. That includes GSVD components of $$\beta_1 u_1 = b,\qquad \alpha_i v_i = A^*u_i - \beta_i v_{i-1},\qquad \beta_{i+1}u_{i+1} = Av_i - \alpha_i u_i,$$7, generalized matrix functions of rectangular matrices, and rational-Krylov generalizations in which the projected matrices become quasiseparable rather than bidiagonal (Arrigo et al., 2015, Casulli et al., 2021). In the latter case, the literature describes short recurrences that recover classical GKB when all poles are at infinity, which suggests a continuous transition from polynomial to rational generalized bidiagonalization (Casulli et al., 2021).

5. Saddle point systems and generalized block formulations

Generalized GKB has also become a notable solver framework for saddle point systems. For systems of the form

$$\beta_1 u_1 = b,\qquad \alpha_i v_i = A^*u_i - \beta_i v_{i-1},\qquad \beta_{i+1}u_{i+1} = Av_i - \alpha_i u_i,$$8

the Craig variant of generalized GKB operates in inner products induced by $$\beta_1 u_1 = b,\qquad \alpha_i v_i = A^*u_i - \beta_i v_{i-1},\qquad \beta_{i+1}u_{i+1} = Av_i - \alpha_i u_i,$$9 and $\|u_i\|_G=\|v_i\|_M=1$0, and, after an augmented Lagrangian modification

$\|u_i\|_G=\|v_i\|_M=1$1

can exhibit convergence that depends only weakly on model size in structural mechanics applications (Arioli et al., 2018). Parallel studies implemented in PETSc show that the resulting inner–outer GKB strategy scales effectively when the inner solves with $\|u_i\|_G=\|v_i\|_M=1$2 are handled by direct or multigrid-accelerated iterative methods (Kruse et al., 2020).

Inexactness in the inner solves does not invalidate the method, but it changes the accuracy model. The inner–outer GKB analysis shows that later inner solves may be relaxed because their perturbations are filtered by the decay of the GKB update coefficients $\|u_i\|_G=\|v_i\|_M=1$3, and dynamic relaxation strategies can reduce the total number of inner iterations substantially (Darrigrand et al., 2022).

For nonsymmetric positive definite leading blocks, the classical symmetric equivalence with CG breaks down. The nonsymmetric extension constructs a generalized GKB process in which the right basis requires full orthogonalization and the reduced operator becomes upper Hessenberg. In this regime, the method is theoretically equivalent to applying FOM to the Schur complement equation rather than CG, and practical stopping criteria can be based on residual recurrences and energy-norm estimates for the primal variable (Dumitrasc et al., 2023).

The 2025 extension to generalized saddle point systems with nonzero $\|u_i\|_G=\|v_i\|_M=1$4 block

$\|u_i\|_G=\|v_i\|_M=1$5

pushes this line further. It develops CRAIG and nsCRAIG analogues for $\|u_i\|_G=\|v_i\|_M=1$6, proves equivalence with Schur-complement reduction combined with inner CG or FOM, and derives residual- and energy-based stopping criteria (Wang et al., 3 Sep 2025). One recurring theme across these block formulations is that generalized GKB offers low memory requirements relative to unrestarted GMRES or full-system Krylov methods, because the bidiagonalization exploits the block structure directly.

6. Tensor, multilinear, and application-driven generalizations

A distinct branch of the literature generalizes GKB from matrices to tensors. In the Einstein-product formulation, tensor equations

$\|u_i\|_G=\|v_i\|_M=1$7

are treated by a global Golub–Kahan bidiagonalization process with updates

$\|u_i\|_G=\|v_i\|_M=1$8

$\|u_i\|_G=\|v_i\|_M=1$9

followed by a reduced tensor Tikhonov problem. This construction is aimed at multilinear systems and discrete ill-posed problems, with applications to color image and video restoration (Guide et al., 2020).

In the t-product setting, weighted tensor GKB methods introduce SPD tensors $M$0 and $M$1 to represent colored noise and regularization. The large-scale tensor problem

$M$2

is reduced by a weighted bidiagonalization process to a small projected problem, and the discrepancy principle is used to determine both $M$3 and the number of bidiagonalization steps (Reichel et al., 2021). A related tensor GKB line based on general tensor operators and Tikhonov regularization formulates reduced problems in the small bidiagonal tensor basis and applies discrepancy-principle parameter choice with quadrature support (Beik et al., 2019).

Other application-driven generalizations extend the GKB idea rather than only its ambient algebra. In inverse Cauchy problems, a CCBM-GKB method combines a coupled complex boundary reformulation with an infinite-dimensional GKB process and Givens rotations, yielding a discrepancy-principle regularization method faster than classic Landweber in the reported experiments (Gong et al., 2022). In X-ray CT with unmatched forward and back projectors, AB- and BA-GKB are proposed as preconditioned forms of GKB; these preserve a bidiagonal projected structure and are reported to be more resilient to semiconvergence than GMRES-type alternatives in the tested problems (Alsubhi, 30 Jun 2026). In image deblurring with Kronecker-product approximations, an enlarged-GKB strategy builds a projected space beyond the desired rank and then truncates, providing low-rank approximations for reordered blur matrices that compete favorably with randomized SVD in the reported tests (Alsubhi et al., 2024).

Taken together, these developments show that generalized GKB has become a methodological bridge between classical bidiagonalization and a wider class of weighted, block-structured, generalized-spectral, and tensor-structured problems. The literature does not present a single terminal form of the method; rather, it presents a common bidiagonalization paradigm adapted to the geometry, algebra, and regularization requirements of each problem class.