Generalized Orthogonalization Algorithm

• Generalized orthogonalization algorithms are methods that replace the Euclidean inner product with application-specific constraints to define orthogonality.
• They extend classical QR, Gram–Schmidt, and Householder techniques by incorporating projections, retractions, and lifting to preserve model structure.
• These methods improve numerical stability and efficiency in varied applications such as optimization, tensor networks, statistical inference, and machine learning.

Generalized orthogonalization algorithm denotes a class of procedures that extend classical orthogonalization beyond Euclidean column orthogonality to structured constraints such as $A^\top S A = S$, $X^\top B X = I_p$, Tensor Train gauge conditions, orthogonal Procrustes systems with multiple unknown orthogonal factors, and projection onto application-specific consistent subspaces (Saraeb, 2024, Shustin et al., 2021, Takeda et al., 30 Jun 2026, Benitez et al., 2024). Across these formulations, the central operation is the construction of an orthogonal or orthogonality-constrained representative by projection, factorization, lifting, retraction, or gauge fixing, with the objective of preserving an ambient model structure while improving conditioning, identifiability, stability, or invariance.

1. Core mathematical patterns

A generalized orthogonalization problem replaces the standard Euclidean constraint $Q^\top Q = I$ by a geometry adapted to the application. In bilinear-form preservation, one seeks matrices $A$ satisfying $A^\top S A = S$ for a fixed invertible symmetric or skew-symmetric $S$, with the associated $S$-bilinear form $\langle x,y\rangle_S = x^\top S y$ (Saraeb, 2024). In optimization on the generalized Stiefel manifold, the feasible set is

$$\mathrm{St}_B(n,p)=\{X\in\mathbb{R}^{n\times p}\mid X^\top B X = I_p\},$$

with tangent-space projections and $B$-orthonormalizing retractions replacing standard QR (Shustin et al., 2021). In Tensor Train (TT) methods, orthogonality is imposed on left- or right-matricizations of TT cores,

$X^\top B X = I_p$0

thereby fixing the TT gauge and producing an orthonormal parameterization (Takeda et al., 30 Jun 2026).

Other variants operate on structured linear spaces rather than on matrices with orthonormal columns. For pairwise comparison matrices, the orthogonalization is performed after logarithmic transformation $X^\top B X = I_p$1, where reciprocity becomes skew-symmetry and consistency becomes additive consistency. The target is the subspace

$X^\top B X = I_p$2

and the consistent approximation is the orthogonal projection of $X^\top B X = I_p$3 onto $X^\top B X = I_p$4 under a generalized Frobenius inner product $X^\top B X = I_p$5 (Benitez et al., 2024). In orthogonalization-based hypothesis testing, the target is a conditional-mean orthogonality relation,

$X^\top B X = I_p$6

obtained after splitting the data by additive symmetric noise (Dharamshi et al., 29 Jun 2026).

These formulations share a common structural feature: orthogonality is defined relative to a chosen inner product, manifold, gauge, or moment condition. This suggests that generalized orthogonalization is less a single algorithm than a reusable design principle specialized to the geometry of the model.

2. Extensions of QR, Gram–Schmidt, and Householder

A central line of development generalizes classical orthogonalization kernels while preserving the numerical guarantees associated with Householder and QR methods. In two-stage orthogonalization, the task is to orthogonalize a matrix $X^\top B X = I_p$7 against an existing orthonormal block $X^\top B X = I_p$8, first by removing components in $X^\top B X = I_p$9, then by orthonormalizing the residual. The generalized Householder transformation

$Q^\top Q = I$0

maps $Q^\top Q = I$1 to $Q^\top Q = I$2 whenever $Q^\top Q = I$3, and is used to build a two-stage Householder orthogonalization that is unconditionally stable and requires orthogonalizing only the square submatrix $Q^\top Q = I$4 rather than the full $Q^\top Q = I$5 (He et al., 16 Feb 2026). With the QR-based choice $Q^\top Q = I$6 from $Q^\top Q = I$7, the construction yields $Q^\top Q = I$8 and $Q^\top Q = I$9; the polar-based choice gives $A$0 (He et al., 16 Feb 2026).

Compact WY representations supply another generalization of Householder orthogonalization. A product of reflectors $A$1 is represented as

$A$2

which reduces synchronization from $A$3 to $A$4 and enables Level-3 BLAS application $A$5 (Ishigami et al., 2012). In inverse iteration for symmetric tridiagonal eigenproblems, replacing modified Gram–Schmidt by compact WY orthogonalization produced reported speedups ranging from $A$6 to $A$7, while preserving Householder-level orthogonality (Ishigami et al., 2012).

The hardware-aware PQR framework recasts orthogonalization as

$A$8

and composes algorithms across the hardware hierarchy rather than fixing a single kernel (Dreier et al., 2022). The framework combines CholeskyQR or BCGS-PIP reductions, TSQR variants, and Householder kernels so that local subproblems fit in cache while inter-node reductions exploit optimized collectives. Stability is tracked through

$A$9

with $A^\top S A = S$0 behavior for Householder, $A^\top S A = S$1 for BMGS, and $A^\top S A = S$2 for BCGS and BCGS-PIP; the reiterated BCGS-PIP+ is stable under $A^\top S A = S$3 (Dreier et al., 2022).

For bilinear-form-preserving orthogonalization, generalized Gram–Schmidt and polar decomposition are carried out in the $A^\top S A = S$4-inner product. The $A^\top S A = S$5-polar construction forms

$A^\top S A = S$6

so that $A^\top S A = S$7 (Saraeb, 2024). For symmetric $A^\top S A = S$8, a commuting-block Schur construction generates Haar-random elements on $A^\top S A = S$9 by diagonalizing $S$0, drawing independent Haar blocks on the eigenspace multiplicities, and setting $S$1 (Saraeb, 2024). For skew-symmetric $S$2, the analogous compact subgroup is generated through unitary blocks and the real embedding

$S$3

which yields Haar sampling on $S$4 (Saraeb, 2024).

3. Lifted, manifold, and polar formulations

A different strand generalizes orthogonalization by lifting nonconvex orthogonality constraints into convex or Riemannian formulations. In the generalized orthogonal Procrustes setting

$S$5

the problem is homogenized to $S$6, lifted through the block Gram matrix $S$7, and relaxed to the semidefinite program

$S$8

For generic $S$9 with $S$0, the relaxation returns a rank-$S$1 solution and recovers the orthogonal factors exactly up to the homogenization symmetry (Zhang et al., 2015). In the $S$2-matrix generalization, exact recovery holds for generic $S$3 if $S$4, and under perturbations $S$5, the rounded estimates satisfy $S$6 (Zhang et al., 2015).

On the generalized Stiefel manifold, orthogonalization is implemented by retractions such as the $S$7-polar map

$S$8

or a $S$9-QR based on Cholesky factorization of $\langle x,y\rangle_S = x^\top S y$0 (Shustin et al., 2021). The corresponding tangent projection under the standard $\langle x,y\rangle_S = x^\top S y$1-metric is

$\langle x,y\rangle_S = x^\top S y$2

and randomized preconditioning is introduced by replacing the natural metric $\langle x,y\rangle_S = x^\top S y$3 with $\langle x,y\rangle_S = x^\top S y$4, often built from a sketch $\langle x,y\rangle_S = x^\top S y$5 through

$\langle x,y\rangle_S = x^\top S y$6

For CountSketch with $\langle x,y\rangle_S = x^\top S y$7, the sketch satisfies

$\langle x,y\rangle_S = x^\top S y$8

hence $\langle x,y\rangle_S = x^\top S y$9 (Shustin et al., 2021). This metric approximation controls the Riemannian Hessian conditioning in CCA and FDA, where orthogonalization enters through $$\mathrm{St}_B(n,p)=\{X\in\mathbb{R}^{n\times p}\mid X^\top B X = I_p\},$$0-retractions rather than Euclidean QR (Shustin et al., 2021).

Polar orthogonalization also reappears in iterative approximation schemes. For $$\mathrm{St}_B(n,p)=\{X\in\mathbb{R}^{n\times p}\mid X^\top B X = I_p\},$$1, the target orthogonal factor is

$$\mathrm{St}_B(n,p)=\{X\in\mathbb{R}^{n\times p}\mid X^\top B X = I_p\},$$2

and the classical third-order Newton–Schulz iteration

$$\mathrm{St}_B(n,p)=\{X\in\mathbb{R}^{n\times p}\mid X^\top B X = I_p\},$$3

is replaced by Chebyshev-optimized odd polynomials in CANS (Grishina et al., 12 Jun 2025). In the degree-3 setting, the interval error obeys $$\mathrm{St}_B(n,p)=\{X\in\mathbb{R}^{n\times p}\mid X^\top B X = I_p\},$$4, with $$\mathrm{St}_B(n,p)=\{X\in\mathbb{R}^{n\times p}\mid X^\top B X = I_p\},$$5 (Grishina et al., 12 Jun 2025). The resulting approximate orthogonalization is used both in Muon-style optimizers and as a practical retraction alternative on the Stiefel manifold (Grishina et al., 12 Jun 2025).

Orthogonalization is also fused with matrix preconditioning in Pro-KLShampoo. Restricting one Kronecker factor to a spike-and-flat family,

$$\mathrm{St}_B(n,p)=\{X\in\mathbb{R}^{n\times p}\mid X^\top B X = I_p\},$$6

the method orthogonalizes the complement via

$$\mathrm{St}_B(n,p)=\{X\in\mathbb{R}^{n\times p}\mid X^\top B X = I_p\},$$7

and an exact identity shows that this recovers the algebraic form of full KL-Shampoo’s complement whitening (Sun et al., 7 May 2026). On GPT-2 124M/350M and LLaMA 134M/450M, the reported wallclock savings to reach matched loss levels were $$\mathrm{St}_B(n,p)=\{X\in\mathbb{R}^{n\times p}\mid X^\top B X = I_p\},$$8, $$\mathrm{St}_B(n,p)=\{X\in\mathbb{R}^{n\times p}\mid X^\top B X = I_p\},$$9, $B$0, and $B$1, respectively (Sun et al., 7 May 2026).

4. Tensor-network and streaming generalizations

In tensor methods, generalized orthogonalization is often expressed as gauge fixing on a tensor network. For an order-$B$2 tensor $B$3 in TT form,

$B$4

Online TT-ALS enforces left- and right-orthogonality incrementally on TT core matricizations (Takeda et al., 30 Jun 2026). Under orthonormal environments,

$B$5

the local ALS subproblem has the exact closed-form minimizer

$B$6

A deterministic single sweep first right-orthogonalizes through LQ factorizations and then left-orthogonalizes through QR, absorbing triangular factors into adjacent cores via gauge transformations that preserve the TT representation (Takeda et al., 30 Jun 2026).

This orthogonalization-centered construction yields explicit theoretical guarantees. The local objective is monotonically nonincreasing under each exact update, and the temporal error satisfies

$B$7

whenever $B$8 (Takeda et al., 30 Jun 2026). In the uniform setting, the dominant per-step complexity is

$B$9

which reduces the rank dependence from quadratic to linear relative to prior online TT methods such as TT-FOA with $X^\top B X = I_p$00 (Takeda et al., 30 Jun 2026).

The empirical results in that work make the role of orthogonalization unusually explicit. At $X^\top B X = I_p$01, batch TT-ALS and mini-batch variants ran out of memory, TT-FOA slowed to approximately $X^\top B X = I_p$02 s/frame with degraded error, whereas Online TT-ALS maintained $X^\top B X = I_p$03 at approximately $X^\top B X = I_p$04 s/frame (Takeda et al., 30 Jun 2026). On grayscale video streams with $X^\top B X = I_p$05, the method reported $X^\top B X = I_p$06, $X^\top B X = I_p$07, $X^\top B X = I_p$08, and $X^\top B X = I_p$09 at $X^\top B X = I_p$10 ms/frame (Takeda et al., 30 Jun 2026). In the ablation without QR/LQ orthogonalization, the reconstruction error drifted from $X^\top B X = I_p$11 at $X^\top B X = I_p$12 to $X^\top B X = I_p$13 at $X^\top B X = I_p$14, whereas with orthogonality it stayed near $X^\top B X = I_p$15 (Takeda et al., 30 Jun 2026).

A related tensor-oriented study compares CGS, MGS, CGS2, MGS2, Gram/Cholesky, and Householder kernels in TT format with TT-rounding at tolerance $X^\top B X = I_p$16 (Coulaud et al., 2022). The reported conclusion is that classical round-off bounds appear to persist with the unit round-off $X^\top B X = I_p$17 effectively replaced by the TT-rounding accuracy $X^\top B X = I_p$18: Householder, CGS2, and MGS2 stay near the rounding floor, MGS degrades roughly like $X^\top B X = I_p$19, and CGS or Gram behave like $X^\top B X = I_p$20 (Coulaud et al., 2022). In this setting, generalized orthogonalization is inseparable from rank control, since TT-rounding is both the dominant cost and the analog of finite-precision arithmetic.

5. Randomized, greedy, and first-order inference-oriented variants

Generalized orthogonalization also appears as a computational strategy for sparse approximation and Krylov subspace construction. In generalized orthogonal least-squares (GOLS), the standard OLS rule

$X^\top B X = I_p$21

is extended to a block rule that selects $X^\top B X = I_p$22 columns at each step, followed by recursive projector downdates

$X^\top B X = I_p$23

The implementation in the paper has dominant complexity $X^\top B X = I_p$24, and the associated OLS theory shows exact recovery in $X^\top B X = I_p$25 iterations with probability at least $X^\top B X = I_p$26 provided $X^\top B X = I_p$27 for Gaussian or Bernoulli designs, whereas the paper states no formal recovery theorem for GOLS itself (Hashemi et al., 2016).

Randomized Gram–Schmidt replaces high-dimensional inner products by sketch-space calculations. Given a sketch $X^\top B X = I_p$28, the algorithm computes $X^\top B X = I_p$29, solves a small least-squares problem in the sketch space to obtain the projection coefficients, and normalizes with the sketched norm $X^\top B X = I_p$30 (Balabanov et al., 2020). The resulting $X^\top B X = I_p$31 is exactly orthonormal in the sketched inner product and approximately orthonormal in $X^\top B X = I_p$32, with singular values controlled by the embedding distortion. In the reported complexity comparison, RGS uses one pass over $X^\top B X = I_p$33 per iteration, requires about half the high-dimensional work of CGS, and integrates into Arnoldi and GMRES, where the randomized GMRES residual is quasi-optimal over a slightly perturbed Krylov subspace (Balabanov et al., 2020).

In high-dimensional statistical estimation, orthogonalization is built into the dynamics of first-order methods through Onsager terms. The generalized first-order template introduces coefficients $X^\top B X = I_p$34 and $X^\top B X = I_p$35 so that the messages

$X^\top B X = I_p$36

become asymptotically orthogonal across time (Montanari et al., 2022). The paper proves a reduction from generalized first-order methods to AMP and then to an orthogonal AMP normal form, yielding scalar state-evolution lower bounds $X^\top B X = I_p$37 and $X^\top B X = I_p$38, with Bayes AMP or GAMP attaining those bounds under proportional asymptotics and under data distributions that extend beyond purely Gaussian designs to independent entries with bounded fourth moments (Montanari et al., 2022). Here orthogonalization is not a postprocessing step but the mechanism that decouples the iteration into effective scalar channels.

6. Debiasing, testing, and structured-decision applications

In machine learning models with non-linearities, generalized orthogonalization is used to remove linear explainability of predictions by sensitive variables after activation. For a monotone element-wise activation $X^\top B X = I_p$39, the corrected prediction-level transform in the GLM setting is

$X^\top B X = I_p$40

and the parameter-level correction can be written as the constrained problem

$X^\top B X = I_p$41

with MDMM updates on $X^\top B X = I_p$42 (Rügamer et al., 2024). For ReLU models, the paper states that residualizing features by $X^\top B X = I_p$43 suffices for the evaluation model with ReLU and $X^\top B X = I_p$44 loss, and the same projection principle extends to tensor-valued outputs through mode-1 multiplication by $X^\top B X = I_p$45 (Rügamer et al., 2024).

In orthogonalization-based hypothesis testing, the data are transformed by external symmetric noise into

$X^\top B X = I_p$46

and one computes under the null

$X^\top B X = I_p$47

The test then checks whether the residual moment

$X^\top B X = I_p$48

holds, with Wald-type calibration through $X^\top B X = I_p$49 (Dharamshi et al., 29 Jun 2026). The same construction extends to post-selection inference by selecting the null on $X^\top B X = I_p$50, computing a debiased statistic $X^\top B X = I_p$51, and using a conditional $X^\top B X = I_p$52 limit under the selected null (Dharamshi et al., 29 Jun 2026).

For pairwise comparisons, the objective is to find the closest consistent matrix to an inconsistent reciprocal matrix. In log-space, the consistent component is the projection of $X^\top B X = I_p$53 onto $X^\top B X = I_p$54 under

$X^\top B X = I_p$55

and the consistent multiplicative matrix is recovered by $X^\top B X = I_p$56 (Benitez et al., 2024). In the Frobenius case $X^\top B X = I_p$57, the projection has the closed form

$X^\top B X = I_p$58

while the orthogonal complement is characterized by $X^\top B X = I_p$59; for general $X^\top B X = I_p$60, the complement becomes $X^\top B X = I_p$61 (Benitez et al., 2024).

Across these variants, guarantees are tied to the geometry that defines orthogonality: exact recovery in semidefinite lifting (Zhang et al., 2015), unconditional stability in generalized Householder orthogonalization (He et al., 16 Feb 2026), provable $X^\top B X = I_p$62 validity in orthogonalization-based testing (Dharamshi et al., 29 Jun 2026), and monotonicity or state-evolution characterizations in TT and AMP settings (Takeda et al., 30 Jun 2026, Montanari et al., 2022). A plausible implication is that generalized orthogonalization is best understood as a transferable algorithmic motif: enforce the orthogonality relation native to the model, then exploit the resulting structure for numerical stability, statistical validity, or computational efficiency.