Skip to content

GPU Speedup Scaling Model

Jump to bottom
s-sasaki-earthsea-wizard edited this page May 16, 2026 · 1 revision

GPU Speedup Scaling Model

Reference page for the empirical observation across correct_topography and invert_network bench scenes that speedup tracks K (= number of valid pixels to invert) much more strongly than D (= num_date) or P (= per-pixel design matrix width). This page works out the structural reason from the per-pixel solve cost up, fits coefficients against measured wall times, and gives a quick "expected speedup for a new scene" lookup.

Empirical data behind this page:

Algorithm-side discussion of why the per-pixel solver choice matters at all is in gels QR vs. Normal Equation + Cholesky.

TL;DR

K is the batch dimension — GPU parallelism amortises the per-pixel Python + scipy dispatch overhead that CPU pays K times in a for loop. D and P are per-solve internal dimensions — they appear in both walls with similar prefactors and don't move the speedup ratio much.

The speedup-vs-scene model is

$$ \sigma(K, D, P) ;\approx; \frac{K,(\alpha + \gamma,D P^2)}{F + K,(\beta + \delta,D P) + \varepsilon,K D P^2} $$

with three regimes:

Regime Condition σ behaviour Example
(A) Overhead-bound K small, $K\alpha \gg K\gamma DP^2$, $F \gg K\delta DP$ $\sigma \approx K\alpha / (F + K\beta)$ — K-only, D drops out Tier 1 small-K cluster: σ = 2.4–2.6× regardless of D ∈ [24, 114]
(B) Transition K small, D large enough that $K\delta DP \sim F$ K cancels; σ depends on D, P, prefactor ratios; can decrease in D SanFranBay (K=326k, D=333) drops to 1.76×
(C) Compute-bound K large, $K\gamma DP^2 \gg K\alpha$ and $K\delta DP \gg F$ $\sigma \to \gamma P / \delta$ asymptote (≈ 10× for correct_topography) Galapagos (K = 3.4M): σ = 6.15×, approaching asymptote

Notation

Symbol Meaning Tier 1 + Galapagos range
$K$ Number of valid pixels to invert (image_shape × valid_mask count) 226k – 3.4M
$D$ num_date — length of the per-pixel time series 24 – 333
$P$ Per-pixel design matrix width. For correct_topography: $P = 1 + (\text{polyOrder} + 1) + #\text{stepDates}$ (i.e. 1 geometric column + polynomial + step columns). For invert_network: $P \approx D - 1$ 4 – 6 (correct_topography)

Per-pixel solve work — correct_topography (normal-equation Cholesky)

For each valid pixel both CPU and GPU paths solve

$$ G_0 \in \mathbb{R}^{D \times P}, \quad N = G_0^\top G_0 \in \mathbb{R}^{P \times P}, \quad y = G_0^\top d \in \mathbb{R}^P, \quad L L^\top = N $$

then back-solve $L u = y$, $L^\top x = u$, and reconstruct $\hat d = G_0 x$ for the ts_cor / ts_res outputs.

Per-pixel flop count:

step flops
$N = G_0^\top G_0$ $D P^2$ ← dominant for $D \gg P$
$y = G_0^\top d$ $D P$
Cholesky on $N$ $P^3 / 3$
Triangular solves $2 P^2$
ts_cor / ts_res reconstruction $D P$ each

For the Tier 1 ranges ($D \in [24, 333]$, $P \in [4, 6]$): per-pixel flops are in the $\sim 400$ – $\sim 5{,}300$ range. Crucially this is linear in D, quadratic in P.

Wall-time decomposition

CPU per-pixel scipy loop

CPU iterates over the K pixels in Python; each pixel calls scipy.linalg.lstsq:

$$ T_{\text{cpu}}(K, D, P) ;\approx; K \cdot \bigl(\alpha + \gamma \cdot D P^2 \bigr) $$

  • $\alpha \approx 10\text{–}40,\mu\text{s}$/pixel — D-independent Python + scipy dispatch overhead. The Python for loop pays this per iteration regardless of D.
  • $\gamma \approx 0.1\text{–}1,\text{ns}$/flop — single-core scalar compute on the dominant $D P^2$ flop count.

The relative size of $\alpha$ vs $\gamma D P^2$ determines whether CPU is "interpreter-bound" or "compute-bound". At Tier 1 values $\gamma D P^2 \approx 40,\text{ns}$ – $1.6,\mu\text{s}$/pixel, so most scenes are interpreter-bound (α dominates).

GPU batched path

GPU does one batched solve across all K pixels simultaneously:

$$ T_{\text{gpu}}(K, D, P) ;\approx; F ;+; K,(\beta + \delta \cdot D P) ;+; \varepsilon \cdot K D P^2 $$

  • $F \approx 3,\text{s}$ — Framework fixed cost: CUDA init, timeseries.h5 + geometryGeo.h5 reads, mask building, kernel launches. K-independent.
  • $\beta \sim \mu\text{s}$/pixel — CPU-side bookkeeping (mask application, index construction) before the GPU dispatch.
  • $\delta \sim \text{ns}$/(D,P)-element — memory traffic: the (K, D, P) G_batch tensor is built on CPU via numpy, transferred via PCIe (Gen4 ~32 GB/s host→device), then read by GPU kernels; outputs (K, D) ts_cor / ts_res flow back.
  • $\varepsilon \approx 10^{-12},\text{s/flop}$ — GPU peak compute. The $\varepsilon K D P^2$ term is sub-second for every scene measured; never the dominant cost on the GPU side.

Crucially: D enters the GPU wall dominantly through $\delta$ (memory bandwidth) rather than through $\varepsilon$ (compute), because the (K, D, P) tensor at K = 326k, D = 333 is already 1.65 GiB at float32 — moving that data dominates seconds-scale wall time, not the 60–100 ms of actual matmul work.

Why K is the batch dimension and D, P are not

The Python for loop on CPU pays the interpreter cost $\alpha$ for every pixel — K times in total. The GPU batched dispatch fuses all K pixels into one kernel call; the K-loop cost is amortised across the batch and effectively replaced by the framework constant $F$.

D and P, by contrast, are inside the per-pixel solve:

  • A larger D means each per-pixel lstsq call processes more data, on both sides. CPU pays $\gamma D P^2$ per pixel; GPU pays $\delta D P$ per pixel (memory) + $\varepsilon D P^2$ per pixel (compute, hidden). D appears in both with similar coefficients → speedup ratio is roughly D-invariant.
  • A larger P means more columns of the per-pixel design matrix. CPU flops scale $P^2$; GPU memory scales $P$ (and compute $P^2$). P is the only way to get the GPU compute term to dominate over its memory term, but for the practical range $P \in [4, 10]$ this never happens at small K because $\varepsilon$ stays sub-second.

The clean way to see this: in the K-cancelled limit (regime B / C),

$$ \sigma_\infty(D, P) ;=; \frac{\alpha + \gamma D P^2}{\beta + \delta D P} $$

K disappears. The GPU advantage that scales with K — the $K\alpha / (F + K\beta)$ term — is gone. What remains is the per-pixel compute-vs-memory ratio, which is bounded above by $\gamma P / \delta$ (taking $D \to \infty$) — an asymptote dictated by the per-flop CPU cost vs the per-byte GPU memory cost, both intrinsic to the hardware

  • implementation.

Speedup formula in the three regimes

Regime (A) — small K, small D: overhead-bound

Conditions: $K\alpha \gg K\gamma DP^2$ (interpreter > compute on CPU) and $F \gg K\delta DP$ (framework > per-element memory on GPU).

$$ \sigma ;\approx; \frac{K\alpha}{F + K\beta} $$

D disappears. Speedup scales only with K. This is the cluster Kuju / SanFranSF / Fernandina(SSD) sit in: σ ≈ 2.4–2.6× regardless of D ∈ [24, 114] or P ∈ [4, 6].

Regime (B) — small K, large D: transition

Conditions: $K$ still too small for $K\gamma DP^2 \gg K\alpha$ on CPU, but D is large enough that $K\delta DP \gtrsim F$ on GPU.

K cancels and the formula collapses to $\sigma_\infty(D, P)$. Because the GPU memory-traffic term grows linearly in D while interpreter overhead $\alpha$ on CPU is D-flat, σ can actually decrease as D grows in this regime: the GPU wall is rising faster than the CPU wall.

This is SanFranBay (K = 326k, D = 333, P = 4): σ drops to 1.76×.

Regime (C) — large K: compute-bound

Conditions: $K\gamma DP^2 \gg K\alpha$ on CPU and $K\delta DP \gg F$ on GPU.

Same K-cancelled formula as regime (B):

$$ \sigma ;\to; \sigma_\infty(D, P) ;=; \frac{\alpha + \gamma D P^2}{\beta + \delta D P} $$

But now the $\gamma D P^2$ term dominates $\alpha$ in the numerator and $\delta D P$ dominates $\beta$ in the denominator, so

$$ \sigma ;\to; \frac{\gamma D P^2}{\delta D P} ;=; \frac{\gamma P}{\delta} $$

— independent of D, linear in P, and bounded by the per-flop CPU vs per-byte GPU coefficient ratio. For correct_topography at P = 4 on RTX 5080 this evaluates to ≈ 10× (report_galapagos.md §2).

Galapagos (K = 3.4M, P = 4) at σ = 6.15× is within striking distance of this asymptote.

Empirical fit — five-scene Tier 1 + Galapagos data

correct_topography, fork commit 60c9e3ac, local SSD, RTX 5080 / PyTorch 2.11+cu128:

Scene K D P $T_{\text{cpu}}$ (s) $T_{\text{gpu}}$ (s) σ Regime
Kuju (ALOS / ROI_PAC) 226k 24 4 9.14 3.61 2.53× (A)
SanFranSF (ARIA) 260k 114 4 9.32 3.84 2.43× (A)
Fernandina (ISCE2, SSD) 270k 98 6 11.15 4.36 2.56× (A)
SanFranBay (GMTSAR) 326k 333 4 14.53 8.23 1.76× (B)
Galapagos (ISCE2) 3,400k 98 4 163.83 26.66 6.15× (C)

Approximate coefficient fit (single-significant-figure, ±30%):

  • $\alpha \approx 40,\mu\text{s}$/pixel — from Kuju: $T_{\text{cpu}}/K \approx 9.14,\text{s} / 226\text{k} \approx 40,\mu\text{s}$ at small D where the $\gamma D P^2$ term is negligible.
  • $F \approx 3,\text{s}$ — from Kuju: $T_{\text{gpu}}$ at smallest scene is 3.6 s, mostly framework overhead.
  • $\delta \approx 12,\text{ns}$/element — from SanFranBay: ($T_{\text{gpu}} - F$) / (K D P) ≈ 5.2 s / 4.3e8 ≈ 12 ns.
  • $\gamma \approx 0.3,\text{ns}$/flop — from Galapagos: $T_{\text{cpu}} / (K D P^2) \approx 164$ s / 5.3e9 ≈ 0.3 ns.

Cross-check by plugging back into the σ formula:

Scene Predicted σ Measured σ
Kuju $\approx 9.0 / (3.0 + 0.26) = 2.8$× 2.53×
SanFranSF $\approx 10.5 / (3.0 + 1.4) = 2.4$× 2.43×
SanFranBay $\approx 13.5 / (3.0 + 5.2) = 1.65$× 1.76×
Galapagos $\approx 152 / (3.0 + 16.0) = 8.0$× 6.15×

Qualitative story holds; absolute prefactor estimates are ±30% slop-level (real $T_{\text{gpu}}$ has additional terms like h5 read cost that this model lumps into $F$ but which actually scale weakly with K + D).

When does GPU dispatch pay off? — practical lookup

correct_topography on RTX 5080 / SSD, ROM-style estimates:

K Expected σ Note
< 100k ≲ 2× K too small to amortise $F$; CPU is interpreter-bound and the K-batch advantage doesn't make it to wall time
100k – 300k 2 – 2.5× Regime (A); D and P don't matter much in this range
300k – 1M 2.5 – ≈ 4× K-axis under-sampled by current data; first crossing into regime (B)/(C); a Tier 2 scene would pin this
1M – 10M 4 – ≈ 10× GPU clearly worthwhile, approaching the $\gamma P / \delta$ asymptote
> 10M ≈ 10× Diminishing returns past Galapagos scale

Storage modifier

The GPU path is more I/O sensitive than the CPU path because its compute portion is sub-second; the timeseries.h5 read becomes the dominant share of $F$ when storage is slow. Measured on Fernandina (K = 270k, report_fernandina.md §9):

Storage $T_{\text{cpu}}$ (s) $T_{\text{gpu}}$ (s) σ
local SSD 11.15 4.36 2.56×
CIFS / NAS 7.45 7.64 0.97×

Same scene, same hardware. NAS storage erases the GPU win at K = 270k. Prefer SSD for the GPU path whenever practical.

Why invert_network reaches 36× while correct_topography tops at ~10×

Both steps share the same K-as-batch-axis structure. The asymptote difference comes from the per-pixel flop count:

correct_topography (Galapagos) invert_network (Galapagos)
Design matrix $G_0 \in \mathbb{R}^{D \times P} = 98 \times 4$ $A \in \mathbb{R}^{\text{num_ifg} \times (D-1)} = 475 \times 97$
Per-pixel flops $D P^2 \approx 1{,}600$ $\text{num_ifg} \times (D-1)^2 \approx 4{,}500{,}000$
Ratio 1× ~2,800× more per-pixel work

In the regime-(C) asymptote $\sigma_\infty = \gamma P / \delta$, replacing $P$ with the larger effective design-matrix width $(D-1)$ for invert_network lifts the asymptote ceiling by roughly $(D-1) / P \approx 24\times$ at Galapagos. Measured invert_network Galapagos σ = 36.4× wall (44.4× internal) (report_large_scene.md) is consistent with this lift.

This is also why the algorithm-side optimisation (gels → batched Cholesky in gels QR vs. Normal Equation + Cholesky) mattered so much for invert_network at K = 3.4M and matters comparatively less for correct_topography: per-pixel work for invert_network is large enough to expose the QR launch-overhead pathology, while correct_topography's smaller per-pixel work hides launch overhead inside framework overhead anyway.

How to project speedup for a new step

Apply the same decomposition:

  1. Count per-pixel flops of the dominant kernel ($D P^2$ for normal-equation Cholesky; different for other algorithms).
  2. Identify the dominant memory-traffic term on GPU: typically the shape of the largest per-pixel intermediate that has to flow through PCIe + GPU memory hierarchy.
  3. Estimate the per-pixel CPU interpreter overhead $\alpha$: usually in the tens-of-µs range when calling scipy or numpy from a Python for-loop.
  4. Find the asymptote $\sigma_\infty \approx \gamma \times (\text{per-pixel-flops}) / \delta \times (\text{per-pixel-bytes})$ using the hardware coefficients ($\gamma \approx 0.3,\text{ns}$/flop, $\delta \approx 12,\text{ns}$/element on the system measured here).
  5. Find the K crossover where $K \alpha / F \approx 1$: this is the K at which speedup starts climbing out of regime (A). For our system $F \approx 3,\text{s}$, $\alpha \approx 40,\mu\text{s}$, crossover at $K \approx 75\text{k}$. Below this, GPU dispatch is overhead-bound regardless of D or P.
  6. Verify the asymptote is much higher than the crossover-region speedup ($\approx 2$× in our case). If not, GPU dispatch won't pay off at any practical K, and the step is not worth GPU-porting.

Related pages

Clone this wiki locally