Skip to content

gels QR vs Normal Equation Cholesky

Jump to bottom
Syota_Sasaki edited this page May 12, 2026 · 7 revisions

gels QR vs. Normal Equation + Cholesky

Reference page for the Phase 2 design decision (issue #4) to switch the per-pixel WLS solver in estimate_timeseries_batch from cusolver's gels driver to batched Cholesky on the normal-equation form.

The two approaches are mathematically distinct and differ in how well their underlying GPU library implementations parallelise across pixels. Both axes matter; this page separates them so design trade-offs do not get confused.

Problem setting

For each pixel k, solve the over-determined weighted least-squares system

$$ x_k ;=; \arg\min_{x}; \big| W_k^{1/2},(G,x - y_k) \big|_2^{,2} $$

with shapes (FernandinaSenDT128 example):

Symbol Shape Meaning
G (num_pair, num_unknown) = (288, 97) Design matrix (shared across pixels for SBAS velocity formulation)
W_k (num_pair, num_pair) diagonal Per-pixel weight
y_k (num_pair,) Per-pixel observed phase differences
x_k (num_unknown,) Per-pixel unknowns (per-interval velocity or per-date phase)

The batch is over n pixels per chunk (auto-sized from free VRAM, typically n ≈ 1.9 × 10⁴).

Approach A: gels (QR-based)

Math

Factor A = QR via Householder reflections, where Q is orthogonal and R is upper triangular. Then solve R x = Q^T y by back-substitution. R has size (num_unknown, num_unknown) once truncated to the active rows.

cusolver implementation (cusolverDnSgels)

  • Each pixel is treated as an independent problem.
  • For each pixel, num_unknown = 97 Householder reflections are applied sequentially.
  • Each reflection requires updating all remaining columns: this shows up in the trace as ormtr_gemv_c (apply Q) and ormtr_gerc (rank-1 update) kernels.
  • "Batching" is implemented by iterating problems on the host; inside the GPU each problem is sequential.

Launch count

For n pixels:

$$ \text{launches} = n \times \text{num-unknown} \times 2 + n $$

For our shape: 19,403 × 97 × 2 + 19,403 ≈ 1.88 million launches per chunk. Empirical confirmation in report_profile.md.

Approach B: Normal Equation + Cholesky

Math

  1. Form the normal equation H_k x_k = b'_k, where H_k = G^T W_k G (symmetric positive definite if W_k^{1/2} G is full rank) and b'_k = G^T W_k y_k.
  2. Cholesky-factor H_k = L_k L_k^T (lower triangular L_k).
  3. Forward solve L_k z_k = b'_k, then back-solve L_k^T x_k = z_k.

torch / cusolver implementation

  • torch.bmm on shape (n, num_unknown, num_pair) × (n, num_pair, num_unknown) builds all H_k in one batched GEMM kernel.
  • torch.linalg.cholesky on a (n, num_unknown, num_unknown) tensor invokes cusolverDnSpotrfBatched, which factors all n matrices in one kernel launch.
  • torch.cholesky_solve issues 2-3 batched triangular-solve kernels for the forward + back substitution.

Launch count

For n pixels: roughly 5 kernel launches, regardless of n.

Why the launch-count gap is structural, not a library quirk

Both NVIDIA cusolver and AMD's rocSOLVER, and academic libraries like MAGMA, expose batched Cholesky but limit batched QR to small or specialised cases. The reason is in the algorithmic structure, not in implementation effort.

Cholesky is batch-parallel friendly

For SPD H, the Cholesky factor is computed by scanning columns left-to-right, with each column producing a triangular update of the trailing block. Within one matrix, this is sequential in column index but otherwise data-local. Across batch elements, each problem is fully independent and can be handed to its own thread block. Single-kernel batched Cholesky scales to thousands of problems per launch with negligible cross-batch synchronisation.

QR Householder is not

Each Householder reflection applies an outer-product update to all remaining columns of the current problem. The amount of work per column varies with the reflection index, and the reflection vector itself depends on the column being eliminated. Trying to batch many problems into one kernel produces:

  • Warp divergence: different batch elements may have already passed through pivoting decisions or numerical truncation thresholds at different reflection indices.
  • Irregular memory access: the active column varies per batch element, defeating coalesced loads.
  • Wasted SM occupancy: a fixed thread-block layout has to assume the worst case work per reflection.

cusolver's choice to expose gels as per-problem (with batch implemented on the host) is therefore a structural consequence of the QR algorithm, not an arbitrary library design.

A "tile-based batched QR" approach (used by MAGMA's batched_dgeqrf) recovers some parallelism, but with substantially more bookkeeping and is not what cusolver exposes via torch.linalg.lstsq.

Speed comparison (FernandinaSenDT128, RTX 5080)

gels (current) Normal-eq + Cholesky (predicted)
Math flops/px ~5.4 Mflops ~3 Mflops
Per-chunk kernel launches ~1.88 M ~5
Per-chunk GPU compute time 10.4 s ~6 s
Per-chunk launch overhead ~10 s ≪ 1 s
Per-chunk wall 19.9 s ≈ 6 s (~3× speedup)

Numbers in column 1 are profile-measured (one chunk of 19,403 px on RTX 5080); see report_profile.md for the methodology and full kernel breakdown. Numbers in column 2 are predicted from the launch-count and flop-count differences; empirical validation is part of issue #4 step 2.

The dominant contribution to predicted speedup is launch-count reduction (~370,000×), not the modest 0.6× math-flop reduction.

Numerical risk: condition-number squaring

The flip side of the normal-equation form is mathematically unavoidable:

$$ \text{cond}!\big(G^T W G\big) = \text{cond}!\big(W^{1/2} G\big)^{2} $$

This is independent of float precision or library implementation.

  • QR is built from orthogonal transforms, which preserve the condition number.
  • The normal equation explicitly multiplies G^T W G, which squares the condition number even before any factorisation begins.

For float32 (≈ 7-8 significant digits) the rule of thumb is:

cond(W^{1/2} G) Normal-eq cond float32 fitness
≤ 10² ≤ 10⁴ Safe (≤ 4 digits lost)
~10³ ~10⁶ Marginal (~6 digits lost)
≥ 10⁴ ≥ 10⁸ Too tight; expect numerical breakdown

SBAS network design matrices typically have cond(G) ≤ 10³ for well-connected networks, but specific datasets can violate this:

  • Sparsely-connected networks (few ifgrams per date)
  • Reference-date pathologies (date with very few connections used as the reference)
  • Ill-conditioned baseline geometry (long temporal or perpendicular baselines clustered together)

The risk is therefore data-dependent, not architectural.

Empirical validation gate (issue #4 step 2)

Per-pixel RMS of the new path against the gels reference path on FernandinaSenDT128 must satisfy

$$ \text{RMS}\big(x_{\text{Cholesky}} - x_{\text{gels}}\big) ;\le; 10^{-4} \times \text{signal scale} $$

with zero NaN / Inf and inv_quality ranking preserved on a sample set of pixels.

  • Gate passes → replace the body of estimate_timeseries_batch with the normal-equation Cholesky path; run the SSD bench harness; produce report_normal_eq.md.
  • Gate fails → pivot options:
    • Mixed precision: build H and factor L in float64, accept the extra memory but recover the lost digits.
    • mask_all_net GPU-path branch: precompute G⁺ once for fully-connected pixels (reduces ~70% of pixels in FernandinaSenDT128 to a single batched matmul), keep gels for the rest.
    • Re-open the previously out-of-scope mixed-precision discussion from Phase 1.

Summary table

Layer gels (current) Cholesky on normal-eq (proposed) Origin of difference
Math flop count per pixel 5.4 Mflops 3 Mflops Algorithm structure
Library batch parallelism Per-pixel iterative True batched Cholesky permits per-block parallelism, QR does not (algorithmic, not arbitrary)
Per-chunk kernel launches 1.88 M ~5 Direct consequence of the row above
Numerical conditioning cond(G) preserved (orthogonal QR) cond(G)² (squared by H = G^T W G) Pure mathematics; independent of precision and implementation
Float32 fitness for SBAS Always safe Conditional on cond(G) ≤ ~10³ Data-dependent; validated empirically

The two columns marked "algorithmic / mathematical" would not change even with a perfect rewrite of cusolver. The columns marked "library implementation" are choices that follow from the math but are still implementation-layer concerns.

References

Clone this wiki locally