Block Orthogonalization Methods

• Block orthogonalization is a technique that enforces orthogonality on groups of vectors or subspaces rather than individual directions, enabling structured and efficient computations.
• It encompasses methods like block QR factorizations, Householder transforms, and randomized sketching, which improve numerical stability and reduce communication overhead.
• The approach finds applications in optimization, tensor computations, and distributed algorithms, addressing issues with finite-precision arithmetic and hardware constraints.

Block orthogonalization denotes a family of constructions in which orthogonality is imposed, maintained, or exploited on groups of vectors, subspaces, or structured updates rather than on single directions. In the literature considered here, the term covers block $QR$ factorizations of column partitions, orthogonalization of selected two-dimensional subspaces inside coordinate descent, row-block orthogonal updates on the Stiefel manifold, orthogonalization kernels for sets of Tensor Train tensors, and communication-aware orthogonalization of Krylov blocks on distributed architectures (Barlow et al., 2011, Dreier et al., 2022, Yuan, 2023). Across these settings, the recurring objective is to remove redundancy inside a block, preserve exact or approximate orthogonality, and expose short recurrences or low-communication updates that are not available in one-vector-at-a-time formulations.

1. Formal viewpoints and mathematical models

A standard linear-algebraic formulation is the project-and-normalize problem: given an orthogonal basis $Q\in\mathbb R^{n\times k}$ and a new block $X\in\mathbb R^{n\times s}$, compute a block $U$ such that $\vspan(Q,U)=\vspan(Q,X)$ and $Q^TU=0$. In the hardware-aware $PQR$ formulation, this is written as

$$\begin{bmatrix} Q & X \end{bmatrix} = \begin{bmatrix} Q & U \end{bmatrix} \begin{bmatrix} I & P\ 0 & N \end{bmatrix},$$

where $P$ is the projection matrix and $N$ is the normalizer (Dreier et al., 2022). This unifies what Gram–Schmidt traditionally treats as two separate steps.

For block $Q\in\mathbb R^{n\times k}$0 itself, a matrix $Q\in\mathbb R^{n\times k}$1 is partitioned into column blocks

$Q\in\mathbb R^{n\times k}$2

and factored as $Q\in\mathbb R^{n\times k}$3, where $Q\in\mathbb R^{n\times k}$4 has orthonormal columns and $Q\in\mathbb R^{n\times k}$5 is block upper triangular (Barlow et al., 2011). The same block structure appears in symmetric saddle point systems

$Q\in\mathbb R^{n\times k}$6

where block Gram–Schmidt is used to obtain a block $Q\in\mathbb R^{n\times k}$7 factorization of $Q\in\mathbb R^{n\times k}$8 and then solve $Q\in\mathbb R^{n\times k}$9 by triangular back substitution (Okulicka-Dłużewska et al., 2013).

A different but equally central model appears on the Stiefel manifold. In the row-block coordinate method OBCD, the constrained problem is

$X\in\mathbb R^{n\times s}$0

and each iteration updates only $X\in\mathbb R^{n\times s}$1 selected rows through

$X\in\mathbb R^{n\times s}$2

This preserves $X\in\mathbb R^{n\times s}$3 exactly while acting only on a local row block (Yuan, 2023).

At a more structural level, polynomial block Krylov subspaces are identified with spaces of matrix polynomials. With starting block $X\in\mathbb R^{n\times s}$4,

$X\in\mathbb R^{n\times s}$5

and the induced matrix-valued inner product is

$X\in\mathbb R^{n\times s}$6

Under a no-deflation assumption, this gives an isometric isomorphism between block Krylov spaces and spaces of matrix polynomials; the rational and extended cases similarly lead to matrix-valued rational functions and Laurent matrix polynomials (Rinelli et al., 16 May 2026). This identifies block orthogonalization with the orthogonalization of matrix-valued functions.

2. Canonical algorithmic patterns

The dominant algorithmic paradigms differ mainly in how they represent a block and how they enforce orthogonality.

Paradigm	Representative mechanism	Source
Reorthogonalized block Gram–Schmidt	Two block CGS passes with local $X\in\mathbb R^{n\times s}$7	(Barlow et al., 2011)
Block Householder	Generalized Householder map $X\in\mathbb R^{n\times s}$8	(He et al., 16 Feb 2026)
Small-block orthogonalization	Gram–Schmidt on a selected two-column block	(Jin et al., 2022)
Row-block orthogonal update	Left multiplication of selected rows by $X\in\mathbb R^{n\times s}$9	(Yuan, 2023)
Tensor block kernels	CGS, MGS, CGS2, MGS2, Gram, Householder in TT format	(Coulaud et al., 2022)

The reorthogonalized block classical Gram–Schmidt algorithm BCGS2 is the canonical two-pass scheme. For a current block $U$0 and previously computed basis $U$1, one block CGS pass forms

$U$2

The second pass repeats the same process on the tentative block and combines the coefficients; in exact arithmetic the result has the form $U$3 (Barlow et al., 2011). In saddle point problems, the same pattern is specialized to the two-block partition $U$4, with the first and second blocks orthogonalized by BCGS or BCGS2 and each local block factorized by thin Householder $U$5 (Okulicka-Dłużewska et al., 2013).

Householder-based block orthogonalization appears in a more specialized two-stage form when an orthonormal block $U$6 is already available and a new block $U$7 must be orthogonalized against it. The generalized Householder transformation

$U$8

maps a canonical orthonormal block to $U$9, after which only the trailing $\vspan(Q,U)=\vspan(Q,X)$0 block of $\vspan(Q,U)=\vspan(Q,X)$1 is $\vspan(Q,U)=\vspan(Q,X)$2-factorized (He et al., 16 Feb 2026). A notable feature is that the method uses only the square submatrix $\vspan(Q,U)=\vspan(Q,X)$3, not a full reorthogonalization of $\vspan(Q,U)=\vspan(Q,X)$4.

A minimal block construction appears in the greedy double subspaces coordinate descent method for overdetermined least squares. At iteration $\vspan(Q,U)=\vspan(Q,X)$5, two active columns $\vspan(Q,U)=\vspan(Q,X)$6 and $\vspan(Q,U)=\vspan(Q,X)$7 are selected, and the second is orthogonalized against the first by

$\vspan(Q,U)=\vspan(Q,X)$8

The associated coordinate-space direction is orthogonalized in the same way, and the update becomes a true projection onto the intersection of two selected hyperplanes rather than two independent coordinate corrections (Jin et al., 2022).

In Tensor Train format, block orthogonalization is realized through six kernels acting on a set $\vspan(Q,U)=\vspan(Q,X)$9: TT-CGS, TT-MGS, TT-CGS2, TT-MGS2, TT-Gram, and TT-Householder. Each method inserts TT-rounding

$Q^TU=0$0

to control rank growth during repeated additions and subtractions (Coulaud et al., 2022).

3. Stability, conditioning, and reorthogonalization

The central numerical issue in block orthogonalization is loss of orthogonality under finite precision. For BCGS2, the classical finite-precision target is

$Q^TU=0$1

with explicit dependence on local block conditions (Barlow et al., 2011). In the saddle point setting, BCGS2 with Householder inner $Q^TU=0$2 is proved backward stable under a mild assumption on the matrix $Q^TU=0$3, whereas plain BCGS can exhibit severe loss of orthogonality and poor residual behavior (Okulicka-Dłużewska et al., 2013).

Recent low-synchronization analyses sharpen this picture. One line of work shows that a reorthogonalized BCGS variant with one synchronization point per block column can attain $Q^TU=0$4 loss of orthogonality under the explicit condition

$Q^TU=0$5

while a two-synchronization version needs only

$Q^TU=0$6

An adaptive strategy combines them to use as few synchronization points as possible under the less restrictive condition $Q^TU=0$7 (Carson et al., 2024). A parallel analysis, however, traces stability degradation step by step as synchronizations are removed and concludes that the one-synchronization-per-block-column variant cannot in general be guaranteed stable in practice; in that framework, degradation begins already when passing from four synchronization points to three (Carson et al., 2024). Taken together, these results indicate that low-synchronization stability depends delicately on the exact algebraic reformulation, not only on the synchronization count.

For TT tensors, the same hierarchy familiar from matrix computations reappears with the TT-rounding tolerance $Q^TU=0$8 replacing the unit roundoff. The paper’s summary table reports TT-CGS and TT-Gram with $Q^TU=0$9-type behavior, while TT-CGS2, TT-MGS2, and TT-Householder achieve $PQR$0-level orthogonality loss (Coulaud et al., 2022). This suggests that reorthogonalization and Householder-type constructions remain the most robust kernels even after low-rank compression.

Two-stage Householder orthogonalization provides a different stability model. Its error bounds depend on the conditioning of the small matrix $PQR$1, not on $PQR$2, and suitable choices of the auxiliary unitary $PQR$3 make $PQR$4 uniformly well-conditioned. In the recommended second choice,

$PQR$5

so

$PQR$6

with

$PQR$7

This is the paper’s basis for describing the method as unconditionally stable (He et al., 16 Feb 2026).

4. Communication-avoiding and hardware-aware designs

Communication cost is a principal driver of modern block orthogonalization research. In $PQR$8-step GMRES, the original BCGS2 with CholQR2 orthogonalization requires five synchronizations every $PQR$9 steps, whereas BCGS-PIP2 reduces this to two synchronizations every $$\begin{bmatrix} Q & X \end{bmatrix} = \begin{bmatrix} Q & U \end{bmatrix} \begin{bmatrix} I & P\ 0 & N \end{bmatrix},$$0 steps by using block Pythagorean inner products to fuse inter-block and intra-block information (Yamazaki et al., 2024). The same paper introduces a two-stage scheme with one synchronization every $$\begin{bmatrix} Q & X \end{bmatrix} = \begin{bmatrix} Q & U \end{bmatrix} \begin{bmatrix} I & P\ 0 & N \end{bmatrix},$$1 steps in stage 1 and one synchronization every $$\begin{bmatrix} Q & X \end{bmatrix} = \begin{bmatrix} Q & U \end{bmatrix} \begin{bmatrix} I & P\ 0 & N \end{bmatrix},$$2 steps in stage 2, while keeping the step size $$\begin{bmatrix} Q & X \end{bmatrix} = \begin{bmatrix} Q & U \end{bmatrix} \begin{bmatrix} I & P\ 0 & N \end{bmatrix},$$3 small for stability. On up to 192 NVIDIA V100 GPUs on Summit, this two-stage scheme reduced orthogonalization time and total time-to-solution by factors of up to $$\begin{bmatrix} Q & X \end{bmatrix} = \begin{bmatrix} Q & U \end{bmatrix} \begin{bmatrix} I & P\ 0 & N \end{bmatrix},$$4 and $$\begin{bmatrix} Q & X \end{bmatrix} = \begin{bmatrix} Q & U \end{bmatrix} \begin{bmatrix} I & P\ 0 & N \end{bmatrix},$$5 over the original $$\begin{bmatrix} Q & X \end{bmatrix} = \begin{bmatrix} Q & U \end{bmatrix} \begin{bmatrix} I & P\ 0 & N \end{bmatrix},$$6-step GMRES; the original $$\begin{bmatrix} Q & X \end{bmatrix} = \begin{bmatrix} Q & U \end{bmatrix} \begin{bmatrix} I & P\ 0 & N \end{bmatrix},$$7-step solver had already achieved speedups of $$\begin{bmatrix} Q & X \end{bmatrix} = \begin{bmatrix} Q & U \end{bmatrix} \begin{bmatrix} I & P\ 0 & N \end{bmatrix},$$8 and $$\begin{bmatrix} Q & X \end{bmatrix} = \begin{bmatrix} Q & U \end{bmatrix} \begin{bmatrix} I & P\ 0 & N \end{bmatrix},$$9 over standard GMRES (Yamazaki et al., 2024).

The hardware-aware $P$0 framework pushes this idea further by combining distinct kernels at different hardware levels. CholeskyQR on the structured matrix $P$1 recovers BCGS-PIP, with

$P$2

so both $P$3 and $P$4 are obtained in one synchronization (Dreier et al., 2022). TSQR is then lifted to the same $P$5 setting, producing TreeTSPQR and FlatTSPQR. The paper’s recommendation is explicit: use TSPQR with Householder locally inside a node, and BCGS-PIP+ for the inter-node reduction stage, because Householder is locally stable and cache-friendly while MPI-level reductions are better handled by optimized all-reduce patterns (Dreier et al., 2022).

Random sketching has recently been used to stabilize the most fragile step of block orthogonalization in $P$6-step GMRES. In RandCholQR, a sketch

$P$7

is first orthogonalized by Householder $P$8 in the sketch space, yielding a triangular factor $P$9; the original block is then transformed by $N$0 before a final CholQR in the original space (Yamazaki et al., 20 Mar 2025). The resulting one-stage randomized method achieves overall orthogonality error of order machine precision whenever the corresponding block vectors are numerically full rank, and the Perlmutter experiments report enhanced numerical stability without a significant increase in execution time (Yamazaki et al., 20 Mar 2025).

5. Optimization, least squares, and learning

In optimization, block orthogonalization often appears not as basis construction but as a structured update rule. The least-squares method GDSCD is a small-block example tailored to highly coherent columns. With

$N$1

the method orthogonalizes a selected two-column block by Gram–Schmidt and obtains the descent identity

$N$2

The amplification factor $N$3 is precisely the gain due to block orthogonalization when the chosen columns are highly coherent (Jin et al., 2022). In the reported experiments, this effect is decisive: for a $N$4 consistent system with $N$5, GCD does not converge within the iteration cap, 2SGS uses $N$6 iterations and $N$7 seconds, while GDSCD uses only $N$8 iterations and $N$9 seconds; for a $Q\in\mathbb R^{n\times k}$00 consistent system with $Q\in\mathbb R^{n\times k}$01, 2SGS uses $Q\in\mathbb R^{n\times k}$02 iterations and $Q\in\mathbb R^{n\times k}$03 seconds, whereas GDSCD uses $Q\in\mathbb R^{n\times k}$04 iterations and $Q\in\mathbb R^{n\times k}$05 seconds (Jin et al., 2022).

OBCD uses exact row-block orthogonal transformations to stay feasible on the Stiefel manifold. With block size $Q\in\mathbb R^{n\times k}$06, each update globally solves a small nonsmooth optimization over $Q\in\mathbb R^{n\times k}$07, and the method converges to $Q\in\mathbb R^{n\times k}$08-block-$Q\in\mathbb R^{n\times k}$09 stationary points with ergodic rate $Q\in\mathbb R^{n\times k}$10; under the KL inequality it also admits non-ergodic convergence-rate statements (Yuan, 2023). For $Q\in\mathbb R^{n\times k}$11, the active block is parameterized explicitly by rotations and reflections,

$Q\in\mathbb R^{n\times k}$12

showing that the local orthogonalization step ranges over the full $Q\in\mathbb R^{n\times k}$13, not only $Q\in\mathbb R^{n\times k}$14 (Yuan, 2023).

In machine learning, the same logic reappears in gradient space. One line of work orthogonalizes layerwise component gradients by SVD,

$Q\in\mathbb R^{n\times k}$15

and then uses $Q\in\mathbb R^{n\times k}$16 in place of the raw gradient in SGDM, Adam, or LARS (Tuddenham et al., 2022). The blocks are filters, columns, or other layer components, and the method is best viewed as component-wise block orthogonalization in update space rather than parameter space. A distributed optimizer perspective appears in MuonBP, where model-parallel gradient matrices are usually orthogonalized per shard and only periodically reassembled for full orthogonalization; the theory requires two distinct stepsizes for blockwise and full steps, and on an 8B model with eight-way tensor parallelism and ZeRO optimizer state sharding the method yields an $Q\in\mathbb R^{n\times k}$17 throughput increase over Muon with no degradation in performance (Khaled et al., 19 Oct 2025).

6. Orthogonalization-free and generalized alternatives

A notable modern trend is to avoid explicit block orthogonalization altogether while retaining its effect at the level of the final subspace or distribution. In spectral clustering, four orthogonalization-free block methods evolve an unconstrained feature matrix $Q\in\mathbb R^{n\times k}$18 by minimizing either

$Q\in\mathbb R^{n\times k}$19

or

$Q\in\mathbb R^{n\times k}$20

or by triangularized fixed-point variants of these objectives (Pang et al., 2023). The $Q\in\mathbb R^{n\times k}$21-based methods converge to features of the form $Q\in\mathbb R^{n\times k}$22 or $Q\in\mathbb R^{n\times k}$23, while the $Q\in\mathbb R^{n\times k}$24-based methods converge to weighted features $Q\in\mathbb R^{n\times k}$25 or $Q\in\mathbb R^{n\times k}$26. In this sense, explicit orthogonalization is replaced by a nonconvex landscape whose minimizers are already the desired invariant subspace up to orthogonal or sign transformations (Pang et al., 2023).

A more radical orthogonalization bypass appears in quantum projection DPP sampling. The standard preprocessing cost is a global $Q\in\mathbb R^{n\times k}$27 decomposition $Q\in\mathbb R^{n\times k}$28, but the proposed alternative uses only column normalization

$Q\in\mathbb R^{n\times k}$29

then corrects the resulting non-orthogonal state preparation by rejection sampling or amplitude amplification (Fanuel et al., 7 Mar 2025). The acceptance probability is exactly

$Q\in\mathbb R^{n\times k}$30

so orthogonalization is traded for a determinant-weighted correction rather than approximated locally or blockwise. The paper is explicit that this is not a block orthogonalization method in the numerical-linear-algebra sense; the correction is global, not blockwise (Fanuel et al., 7 Mar 2025).

At the highest level of abstraction, the structural correspondence between block Krylov subspaces and orthogonal matrix polynomials shows that block orthogonalization can be understood as orthogonalization of matrix-valued functions under induced inner products. In the unitary case, the associated measure is supported on the unit circle, and Szegő and CMV recurrences transfer to block Krylov bases (Rinelli et al., 16 May 2026). This suggests that many apparently different block orthogonalization algorithms are instances of a common multiplication-by-$Q\in\mathbb R^{n\times k}$31 recurrence in a matrix-valued orthogonal function space.

A recurring limitation across the literature is that the strongest guarantees are usually conditional. Reorthogonalized block Gram–Schmidt depends on local block conditioning; small-block methods such as GDSCD only capture pairwise interactions; OBCD’s stronger optimality notions rely on exact small-block solves; TT methods inherit sensitivity to TT-rounding accuracy; and orthogonalization-free alternatives recover invariant subspaces or equivalent features, not necessarily a canonical orthonormal basis (Jin et al., 2022, Coulaud et al., 2022, Yuan, 2023, Pang et al., 2023). The persistent theme is therefore not the elimination of orthogonality, but the search for formulations in which orthogonality is enforced where it is most valuable and relaxed where it is most expensive.