591a4f5Reduce the Laplace-Beltrami matvec to 6 transforms per species per iterationDan Fortunato 1# Reducing spherical-harmonic transforms in $\Delta_\Gamma$
3**Summary.** Algorithm 4 costs 12 transforms per matvec (8 syntheses, 4 analyses). A flux-form
4reformulation, with weights chosen so that every analyzed field is smooth on $S^2$, evaluates the
5same operator in **6 transforms** (4 syntheses, 2 analyses). Notation follows `algos.pdf`.
7---
9## 1. Where the current cost comes from
11| Algorithm 4 line | Work | Transforms |
12|---|---|---|
13| 1 | $\partial_\theta u,\ \partial_\varphi u$ | 2 $\mathcal{S}$ |
14| 5 | analysis of 3 Cartesian components of $\nabla_\Gamma u$ | 3 $\mathcal{A}$ |
15| 6 | $\partial_\theta$ and $\partial_\varphi$ of each of those 3 components | 6 $\mathcal{S}$ |
16| 8 | final analysis | 1 $\mathcal{A}$ |
17| | | **12** |
19Two sources of waste:
211. The gradient is carried in **ambient $\mathbb{R}^3$ components** — 3 fields for an intrinsically
22 2-dimensional object.
232. Line 6 takes **both** derivatives of **each** component, where the divergence needs only one
24 derivative of each of two fluxes.
26There is also a possible free win independent of everything below: Algorithm 1 as specified returns
27all five derivatives. The Laplacian path needs only $\partial_\theta u$ and $\partial_\varphi u$;
28the second-derivative and mixed-derivative machinery is used exclusively by Algorithm 3 (curvature).
29If `surface_screened_laplacian` calls Algorithm 1 wholesale, it is doing 5 syntheses where 2 suffice
30at lines 1 and 6.
32---
34## 2. The constraint that shapes the solution
36The Cartesian design in Algorithm 4 exists to avoid pole singularities, and it is correct to do so.
37The relevant property is smoothness **as a scalar function on $S^2$**, since that is what controls
38SH coefficient decay and hence whether $\mathcal{A}$ is meaningful.
40| Quantity | Smooth on $S^2$? |
41|---|---|
42| $\partial_\varphi u$ | yes — exactly band-limited, eq. (2.2) |
43| $\sin\theta\,\partial_\theta u$ | yes — exactly band-limited, eq. (2.4) |
44| $\partial_\theta u$ | **no** — bounded, but $\varphi$-dependent limit at the poles |
45| $g_{\theta\theta},\ g^{\theta\theta},\ V_\theta,\ V_\varphi$ | **no** in general |
46| $(\nabla_\Gamma u)_x,\ (\nabla_\Gamma u)_y,\ (\nabla_\Gamma u)_z$ | yes |
48Concretely, for the ellipsoid $X = (a\sin\theta\cos\varphi,\ b\sin\theta\sin\varphi,\ c\cos\theta)$,
49$|X_\theta|^2 \to a^2\cos^2\varphi + b^2\sin^2\varphi$ as $\theta\to 0$: no limit exists.
51Algorithm 4 never analyzes anything in the "no" rows — the non-smooth quantities appear only as
52pointwise grid factors. **Any replacement must preserve this property.** The naive flux form
53$P = \sqrt{g}\,(g^{\theta\theta}u_\theta + g^{\theta\varphi}u_\varphi)$,
54$Q = \sqrt{g}\,(g^{\varphi\theta}u_\theta + g^{\varphi\varphi}u_\varphi)$ does not:
55on the round sphere with $u = x$, $Q = -\sin\varphi$, which is not a function on $S^2$.
57### The correct weighting
59$P$ and $Q$ are $\sqrt{g}$ times the contravariant components of $G := \nabla_\Gamma u$. Using
60$\det[X_\theta, X_\varphi, n] = -\sqrt{g}$:
62$$P = -\,G\cdot(X_\varphi \times n), \qquad \sin\theta\,Q = -\,G\cdot\big(n \times \sin\theta\,X_\theta\big).$$
64Every factor on the right is smooth on $S^2$: $G$ smooth, $n$ smooth, $X_\varphi$ smooth, and
65$\sin\theta\,X_\theta$ smooth because it is exactly band-limited by the recurrence already
66implemented. Hence $P$ and $\tilde{Q} := \sin\theta\,Q$ are analyzable, on the same footing and for
67the same structural reason as the Cartesian gradient components.
69The cross products are the smoothness certificate only — they are not needed in the code.
71---
73## 3. Precompute (once per surface, grid space)
75Replaces `_precompute_metric_quantities()`. From the embedding coefficients $\hat{X}^m_\ell$, obtain
76$X_\varphi$ and $\sin\theta\,X_\theta$ componentwise via Algorithm 1, where the latter is the
77**undivided** output of Algorithm 1 line 4, i.e. $\mathcal{S}(v^m_\ell)$ with no $/\sin\theta$.
78Then, pointwise:
80$$\tilde{g}_{\theta\theta} := |\sin\theta\,X_\theta|^2, \qquad
81 \tilde{g}_{\theta\varphi} := (\sin\theta\,X_\theta)\cdot X_\varphi, \qquad
82 g_{\varphi\varphi} := |X_\varphi|^2$$
84$$J := \frac{\sqrt{\tilde{g}_{\theta\theta}\,g_{\varphi\varphi} - \tilde{g}_{\theta\varphi}^{\,2}}}{\sin^2\theta}
85 \qquad\text{so that } \sqrt{\det g} = J\sin\theta$$
87(The radicand is $\sin^2\theta\det g = J^2\sin^4\theta$, so the square root is $J\sin^2\theta$ — hence
88$\sin^2\theta$, not $\sin\theta$, in the denominator.)
90Store four scalar grid arrays:
92$$p_1 = \frac{g_{\varphi\varphi}}{J\sin^2\theta}, \qquad
93 p_2 = -\frac{\tilde{g}_{\theta\varphi}}{J\sin^2\theta}, \qquad
94 q_2 = \frac{\tilde{g}_{\theta\theta}}{J\sin^2\theta}, \qquad
95 r = \frac{1}{J\sin^2\theta}$$
97All four are bounded: the $\sin^2\theta$ denominators cancel against vanishing numerators
98($\tilde{g}_{\theta\varphi} = O(\sin^2\theta)$, $g_{\varphi\varphi} = O(\sin^2\theta)$), the same
99finite limits the current $V_\theta, V_\varphi$ have. Note $p_1, p_2, q_2$ are bounded where
100$V_\varphi = O(1/\sin\theta)$ is not.
102**Three scalar arrays replace the six components of $V_\theta, V_\varphi$.** The surface
103representation is unchanged: $X_\theta, X_\varphi$ still come componentwise from $\hat{X}^m_\ell$
104via Algorithm 1.
106---
108## 4. The per-matvec algorithm
110Input $\{u^m_\ell\}$; output $\{(\Delta_\Gamma u)^m_\ell\}$.
112| # | Step | Transforms |
113|---|---|---|
114| 1 | $v^m_\ell \leftarrow \alpha^+(\ell-1,m)u^m_{\ell-1} + \alpha^-(\ell+1,m)u^m_{\ell+1}$; $A \leftarrow \mathcal{S}(v^m_\ell)$ | $\mathcal{S}$ |
115| 2 | $B \leftarrow \mathcal{S}(im\,u^m_\ell)$ | $\mathcal{S}$ |
116| 3 | $P \leftarrow p_1 A + p_2 B$, $\tilde{Q} \leftarrow p_2 A + q_2 B$ — pointwise | — |
117| 4 | $\hat{P} \leftarrow \mathcal{A}(P)$, $\hat{\tilde{Q}} \leftarrow \mathcal{A}(\tilde{Q})$ | 2 $\mathcal{A}$ |
118| 5 | $s^m_\ell \leftarrow \alpha^+(\ell-1,m)\hat{P}^m_{\ell-1} + \alpha^-(\ell+1,m)\hat{P}^m_{\ell+1} + im\,\hat{\tilde{Q}}^m_\ell$ | — |
119| 6 | $\Delta_\Gamma u \leftarrow r \cdot \mathcal{S}(s^m_\ell)$ | $\mathcal{S}$ |
120| 7 | $\{(\Delta_\Gamma u)^m_\ell\} \leftarrow \mathcal{A}(\Delta_\Gamma u)$; zero $\ell \ge L-2$ | $\mathcal{A}$ |
122**4 syntheses + 2 analyses = 6**, versus 12.
124- $A$ and $B$ are exactly $\sin\theta\,\partial_\theta u$ and $\partial_\varphi u$.
125- Steps 1 and 5 use the **same** precomputed $\alpha^\pm$ table; step 5 is the adjoint-style reuse
126 of the shift already implemented for step 1. Adding the two flux contributions in coefficient
127 space before synthesizing is what saves the final pair of transforms.
128- The **only** division by $\sin\theta$ anywhere is folded into $p_1, p_2, q_2, r$ at precompute
129 time. The per-matvec path contains none.
0d99c91Differentiate the phi flux in grid spaceDan Fortunato 131### Implemented variation (2026-08-05): the $\varphi$-flux never needs the Legendre basis
133Step 4's analysis of $\tilde{Q}$ exists only so step 5 can apply $\partial_\varphi$ — but
134$\partial_\varphi$ is diagonal in the Fourier index, so the implementation differentiates
135$\tilde{Q}$ on the grid instead: FFT each latitude row, multiply mode $m$ by $im$ (zeroing
136$m \ge L-2$ to mirror the top-degree filter; the Fourier analysis stage truncates $m > m_{\max}$
137for free), inverse FFT. **5 Legendre transforms + one Legendre-free FFT derivative**, versus 6.
138The caveat is that the grid route skips $\tilde Q$'s band projection in $\ell$; measured, this
139does not bite — the band-edge spectra are identical to the 6-transform route's (the $m$ mask and
140the final analysis's projection contain it), the Algorithm-4 A/B agreement is unchanged
141($3.6\times10^{-4}$ after 20 steps at $L=63$), and the step gets ~8% faster at $L=255$
142(~2% at $L=127$, where transform batching had already amortized most of what this removes).
591a4f5Reduce the Laplace-Beltrami matvec to 6 transforms per species per iterationDan Fortunato 144---
146## 5. Numerical trade-off
148Both schemes contain two powers of $\sin\theta$ division in total. What differs is **placement**.
150- **Algorithm 4** spends them in separate stages, one before line 5 and one after. The intervening
151 analysis suppresses the polar spike: Gauss–Legendre weights give $w_1 = O(L^{-2})$ at the polar
152 ring, so a grid error of $\varepsilon L$ there contributes
153 $\sim L^{-2}\cdot L^{1/2}\cdot \varepsilon L = \varepsilon L^{-1/2}$ to any coefficient. Stage 2
154 then starts from clean coefficients and incurs a *fresh* $\varepsilon L$. The two amplifications
155 never multiply. Net grid-space relative error: $\varepsilon L$.
156- **The new scheme** has no division at all through step 5, then pays for both powers at once in
157 $r = O(L^2)$ at step 6 — one event, with no intervening analysis to break it in half. Net
158 grid-space relative error: $\varepsilon L^2$.
160The mechanism: with $N \approx L+1$ Gauss–Legendre nodes, $1 - x_1 = O(N^{-2})$ so
161$\sin\theta_1 = O(N^{-1})$. Since $s = J\sin^2\theta\,\Delta_\Gamma u$ is $O(L^{-2})$ at the polar
162ring but $O(1)$ over the bulk, and synthesis commits roundoff scaled by the field's *global* size at
163every node alike, multiplying by $r \sim L^2$ recovers the signal and inflates the noise.
165| | divisions | placement | grid-space relative error |
166|---|---|---|---|
167| Algorithm 4 | $\sin\theta$, $\sin\theta$ | separated by $\mathcal{A}$ | $\varepsilon L$ |
168| Six-transform | $\sin^2\theta$ | all at the end | $\varepsilon L^2$ |
170**This likely does not reach the returned coefficients.** Step 7's analysis suppresses the spike
171exactly as line 5 does today: $L^{-2}\cdot L^{1/2}\cdot\varepsilon L^2 = \varepsilon L^{1/2}$,
172comparable to the ordinary $\varepsilon\sqrt{L}$ accumulation of a transform pair — and the new
173scheme runs half as many transforms, lowering that baseline. Inside the implicit solve, GMRES sees
174only coefficients, so the extra power should be invisible.
176It matters only if grid values of $\Delta_\Gamma u$ are consumed directly: a nonlinear reaction
177term, max-norm diagnostics, or an adaptive error estimator.
179Algorithm 1 line 7 already divides by $\sin^2\theta$, so the code is exposed to $\varepsilon L^2$
180today — just on the second-derivative path, which the Laplacian never touches.
182---
184## 5a. float32 / WebGPU
186Target is WebGPU, which is float32-only: $\varepsilon = 2^{-24} \approx 6\times10^{-8}$ (spacing
187$2^{-23} \approx 1.2\times10^{-7}$). No float64 fallback exists on device. All estimates in §5 are
188linear in $\varepsilon$, so they scale directly:
190| | $\varepsilon\sqrt{L}$ (coeffs) | $\varepsilon L$ (Alg. 4 grid) | $\varepsilon L^2$ (new, grid) |
191|---|---|---|---|
192| $L=64$ | $5\times10^{-7}$ | $4\times10^{-6}$ | $2\times10^{-4}$ |
193| $L=128$ | $7\times10^{-7}$ | $8\times10^{-6}$ | $1\times10^{-3}$ |
194| $L=256$ | $1\times10^{-6}$ | $1.5\times10^{-5}$ | $4\times10^{-3}$ |
196**Coefficient space is fine** (~$10^{-6}$), which is the floor a float32 iterative solve sits at
197anyway. **Grid space is not**: 0.1–0.4% relative on the polar rings at $L\ge128$. For a
198reaction–diffusion solver this matters only if grid-space $\Delta_\Gamma u$ is consumed outside the
199matvec. If the IMEX splitting evaluates $f(u)$ from $u$ on the grid (typical), it never is.
201**Two float32-specific arguments in favour of the new scheme:**
203- Baseline SHT roundoff accumulates per transform ($\sim\varepsilon\sqrt{L}$ to $\varepsilon L$
204 each). Running 6 transforms instead of 12 halves that accumulation. On the coefficient-space error
205 GMRES actually sees, this plausibly outweighs the polar term — the new scheme may be *more*
206 accurate end-to-end in float32. Not asserted without measurement.
207- 3 weight arrays instead of 6 halves per-matvec texture/buffer traffic. On GPU that is often the
208 real bottleneck, independent of arithmetic.
210**Mitigations available without float64:**
212- **Pairwise or blocked summation in the Legendre sum over $\ell$.** The single highest-value
213 float32 change, and it benefits the existing code too. See §5b.
214- **Double-float (`f32x2`) arithmetic** for the pointwise steps 3 and 6 if needed — cheap, no
215 transforms involved. Does not help with transform roundoff, which is the dominant term, so try
216 summation order first.
217- **CPU precompute in float64.** JS `Number` is float64, so §3 can run on the CPU regardless of
218 WebGPU's limits, with float32 weights uploaded. Cost is CPU-side SHTs plus upload, paid once per
219 surface update; viable if the surface evolves slowly or is prescribed analytically, likely too
220 slow if the metric is rebuilt every timestep. Per the correction below, this is probably
221 unnecessary.
222- **Cap $L$.** All the error terms grow with $L$; float32 sets a practical ceiling that float64
223 would not.
225---
227## 5b. Summation order in the Legendre transform
229This is orthogonal to the 12→6 change, applies equally to the current code, and in float32 is
230probably worth more than the transform-count reduction. Do it first and independently, so its effect
231can be measured on its own.
233**Why.** Every $\varepsilon L$ and $\varepsilon L^2$ in §5 rides on the per-transform roundoff
234floor, and in float32 that floor is set by *how the sums are accumulated*, not by the mathematics.
235For each $(m, \theta_i)$ the synthesis evaluates
237$$u^m(\theta_i) = \sum_{\ell=|m|}^{L} u^m_\ell\,\bar P^m_\ell(\cos\theta_i),$$
239an $O(L)$-term sum. Error growth by accumulation strategy, for an $N$-term sum:
241| Strategy | Worst case | Typical (random signs) |
242|---|---|---|
243| Sequential | $\varepsilon N$ | $\varepsilon\sqrt{N}$ |
244| Pairwise / tree | $\varepsilon\log_2 N$ | $\varepsilon\sqrt{\log_2 N}$ |
245| Kahan compensated | $\varepsilon$ (+ $O(\varepsilon^2 N)$) | $\varepsilon$ |
247At $L=256$ in float32 that is the difference between $\sim1.5\times10^{-5}$ and $\sim5\times10^{-7}$
248per transform — more than an order of magnitude, for no change in operation count.
250**On GPU this may already be partly free.** A workgroup tree reduction over $\ell$ *is* pairwise
251summation. The failure mode is a serial `for` loop over $\ell$ inside a single thread, which is the
252natural way to write the shader if each thread owns one $(m,\theta_i)$ pair and is exactly the
253$\varepsilon N$ row above. Check which shape the kernel has before assuming anything.
255**Where it applies.**
257- Synthesis $\mathcal{S}$: the sum over $\ell$, as above. The $\varphi$-direction FFT is already
258 tree-structured and needs no attention.
259- Analysis $\mathcal{A}$: the quadrature sum over latitude nodes $\theta_i$ carries the identical
260 problem and the identical fix. It also matters more here, because this is the step relied on in
261 §5 to suppress the polar spike — a noisy quadrature sum weakens exactly the mechanism the
262 six-transform scheme depends on.
264**Practical notes.**
266- Blocked summation (accumulate in blocks of 8–32, then combine) captures most of the pairwise
267 benefit with a simpler kernel and better register behaviour than a full tree.
268- Kahan costs ~4 flops per term and is usually bandwidth-hidden on GPU; worth benchmarking rather
269 than assuming it is too expensive.
270- For $m>0$ near the poles, $\bar P^m_\ell(\cos\theta)$ spans many orders of magnitude across $\ell$.
271 Summing smallest-magnitude-first helps, and is nearly free here because the terms are already
272 roughly ordered by $\ell$.
273- Standard stable recurrences for $\bar P^m_\ell$ (and guarding their under/overflow in float32's
274 narrower exponent range) are a separate prerequisite — no summation strategy rescues inaccurate
275 Legendre values.
277**Measurement.** Transform a band-limited field forward then back and compare to the input, in
278float32, sweeping $L\in\{64,128,256\}$. Sequential accumulation shows error growing roughly linearly
279in $L$; pairwise shows near-flat growth. This isolates the transform floor from everything else in
280§7 and should be run before the validation gate there, since it sets the baseline that gate is
281measured against.
283### Measured (2026-08-04, Dawn/Metal, `scripts/sht-accuracy.ts`)
285The sweep was run and the summation-order changes tried. Outcome: **withdrawn — the floor here is
286not summation-limited.**
288| $L$ | grid | rel-$L_2$ roundtrip | worst degree |
289|---|---|---|---|
290| 63 | 64×128 | $3.4\times10^{-6}$ | $\ell=62$: $4.5\times10^{-6}$ |
291| 127 | 128×256 | $4.7\times10^{-6}$ | $\ell=110$: $6.4\times10^{-6}$ |
292| 255 | 256×512 | $1.1\times10^{-5}$ | $\ell=246$: $1.4\times10^{-5}$ |
294- **The analysis side already sums pairwise.** The quadrature over latitudes is a workgroup
295 tree/subgroup reduction (`leg_analys`); only the synthesis has the serial per-thread $\ell$-loop.
296- **Kahan is unavailable on WebGPU in practice.** Dawn/Metal compiles WGSL with fast-math: a probe
297 kernel evaluates $((10^8 + 1) - 10^8) - 1$ to $0$, so the compensation folds away and Kahan
298 compiles to plain summation (bit-identical results, verified).
299- **Blocked summation (B=16) in the synthesis $\ell$-loop moved nothing**: $1.128\times10^{-5}
300 \to 1.128\times10^{-5}$ at $L=255$ (low digits shift, confirming the reordering was live), while
301 costing ~5% per round trip at $L=255$. Reverted.
302- **Diagnosis:** the worst error concentrates at the top degrees — the signature of the Legendre
303 *recurrence* error (chains of length $\sim\ell$), not of $\ell$-uniform accumulation noise. This
304 is the "standard stable recurrences are a separate prerequisite" caveat above: the floor is set
305 by the accuracy of the $\bar P^m_\ell$ values themselves, and no summation strategy touches it.
306- The measured floor ($\sim\varepsilon L^{0.85}$, $1.1\times10^{-5}$ at $L=255$) is what the §7
307 validation gate should be read against.
309---
311## 6. Code changes
313| Location | Change |
314|---|---|
315| `src/surface_gradient/partial_derivatives` | Expose $\mathcal{S}(v^m_\ell)$ **pre-division** (flag or separate entry point). Needed by both the precompute and step 1. |
316| `SurfaceDiffOperator._precompute_metric_quantities()` | Return `p1, p2, q2, r` instead of `V_theta, V_phi`, per §3. |
317| `src/surface_screened_laplacian::surface_screened_laplacian()` | Replace body with §4. Both `for i in {x,y,z}` loops disappear. |
318| `SurfaceDiffOperator._precompute_curvature()`, Algorithm 3 | **Unchanged.** Still needs $X_{\theta\theta}, X_{\theta\varphi}, X_{\varphi\varphi}$ and the full Algorithm 1. |
319| `src/timestepping::make_implicit_op()`, Algorithm 5 | **Unchanged.** Only what line 8 calls changes. |
320| `src/real_embedding.py` | **Unchanged.** |
322The deprecated `SurfaceDiffOperator` methods for $\Delta_\Gamma$ and $(I + c\Delta_\Gamma)$ are the
323natural place to keep the old path as a reference implementation for the validation below.
325---
327### Correction: precompute conditioning
329An earlier draft claimed the polar relative error in $\tilde g_{\theta\theta}$ is $\varepsilon L^2$,
330making float64 precompute essential. That was wrong by a factor of $L$, in the safe direction.
331$\tilde g_{\theta\theta}$ is not synthesized directly; it is the square of $\sin\theta\,X_\theta$,
332which *is* synthesized, is $O(\sin\theta)$ at the poles, and carries absolute error $\varepsilon$ —
333so relative error $\varepsilon L$, preserved (up to a factor 2) by squaring. Same for
334$g_{\varphi\varphi}$ and $\tilde g_{\theta\varphi}$. The determinant combination is $O(\sin^4\theta)$
335and so are both of its terms, so there is no extra cancellation generically; $J$ inherits
336$\sim\varepsilon L$, i.e. $\sim10^{-5}$ at $L=128$ in float32. Acceptable.
338Caveat: this assumes the difference is not small compared to its terms, which fails if $X_\theta$
339becomes nearly parallel to $X_\varphi$ (near-degenerate parametrization). Worth a runtime check on
340$\det g$ if the surface can deform that far.
342Precompute error is also a *fixed* perturbation, identical every matvec, so it perturbs which
343operator is being solved but injects no noise into the Krylov space — GMRES converges normally.
345---
347## 7. Validation, in order
349At $\varepsilon=6\times10^{-8}$ there is no margin for the $\varepsilon\sqrt{L}$ suppression estimate
350in §5 to be off by an order of magnitude. Step 2 is a **gate**, not a confirmation.
3521. **Smoothness check (do this first).** On a deformed, non-axisymmetric surface, form $P$ and
353 $\tilde{Q}$ on the grid and compare their SH coefficient decay against $(\nabla_\Gamma u)_x$ from
354 the current code. Matching tails confirm both are genuinely smooth on $S^2$. If this fails,
355 nothing else is worth doing. Run this in float64 on CPU — it is a mathematical check, not a
356 precision one.
3572. **Coefficient-space diff in float32 at production $L$**, against a float64 CPU reference
358 implementation of Algorithm 4. Landing near $10^{-6}$ means the suppression argument holds.
359 Landing near $10^{-4}$ means the polar spike is surviving the analysis and the polar rings need
360 separate handling.
3613. **Sweep $L \in \{64,128,256\}$** and fit the growth exponent of (2). Flat-ish confirms
362 suppression; growth like $L^2$ means it is not working.
3634. **Grid-space max-norm diff near the poles.** Expect the extra power of $L$ here. If only this
364 grows and (2) stays flat, the scheme is fine for use inside the implicit solve.
3655. **GMRES iteration count and final achieved residual**, float32, versus the current code. The
366 operator is the same, so iterations should be unchanged; a stall above tolerance that does not
367 occur in float64 indicates the matvec noise floor is binding.
369---
371## 8. Suggestions considered and withdrawn
373- **Splitting $r$ across steps 3 and 6** to keep $\sin^1$ scaling. Does not work: $P/\sin\theta$ is
374 not smooth on $S^2$ (round sphere, $u = x$: $P = \sin\theta\cos\theta\cos\varphi$, so
375 $P/\sin\theta \to \cos\varphi$ at the pole). That $\sin\theta$ must stay in $r$. Algorithm 4 *can*
376 split because its intermediate — the Cartesian gradient — is smooth; that smoothness is precisely
377 what the six extra transforms buy.
378- **Weighting by $J$ to get an SPD operator and use PCG.** Avoiding the $1/\sqrt{g}$ makes
379 $\mathcal{L}$ self-adjoint, but the mass term becomes multiplication by $J$, costing its own
380 synthesis/analysis pair. A wash on transform count; worth it only if the CG properties themselves
381 are wanted. Discrete symmetry would also hold only to quadrature accuracy unless products are
382 dealiased (3/2 rule).
383- **Mixed precision (float64 for step 6's synthesis and the $r$ multiply only).** Unavailable:
384 WebGPU is float32-only. Superseded by the mitigations in §5a.
385- If the $\varepsilon L^2$ ever does bind, the remaining remedies are a shifted or uniform-in-$\theta$
386 latitude grid (removing the $O(L^{-2})$ node clustering) or a separate local formula for the polar
387 rings — both more work than the transform savings justify without a specific reason. In float32
388 the ceiling on $L$ may bind first and be the cheaper accommodation.
390## 9. Related: external vector transforms
392Steps 1–2 and 4–5 together are a vector/spin-weighted spherical harmonic transform. If SHTns
393(`spat_to_SHsphtor`, `SHsphtor_to_spat`) or SPHEREPACK (`gradgs`, `divgs`) can be linked, the
394gradient and divergence each become a single library call, the pole divisions are handled internally,
395and the hand-rolled $\alpha^\pm$ recurrences on this path are no longer needed.
397## 10. Beyond transform count
399The other lever is iteration count rather than cost per iteration. Since $M^{-1}A$ approaches
400multiplication by $1/J$ at high $\ell$, folding a mean or smoothed $J$ into the preconditioner could
401reduce GMRES iterations by more than any of the above reduces transforms.
e4d6a3bPrecondition with the operator's symbol; project the correction onto the bandDan Fortunato 402
403### Measured (2026-08-05, Richardson iteration, fp32)
405Implemented, with two corrections the measurements forced.
407**The right constant is the minimax over the symbol, not over $J$.** The high-$\ell$ per-mode
408factor is governed by the full principal symbol: in the orthonormal frame the symbol matrix is
409$S = (1/J)\begin{pmatrix} p_1 & p_2 \\ p_2 & q_2\end{pmatrix}$, whose eigenvalues $\mu(x)$ are the
410inverse squared principal stretches — direction matters. Preconditioning with $\lambda/\hat J$
411contracts every mode and direction iff $\hat J\mu \in (0,2)$, so
412$$\hat J = 2/(\mu_{\min} + \mu_{\max}), \qquad \text{rate} = (\mu_{\max}-\mu_{\min})/(\mu_{\max}+\mu_{\min}) < 1.$$
413The det-based mean of $J$ (this section's original suggestion; $\mu$'s geometric mean, exact only
414for conformal surfaces) is insufficient: on the shipped ellipsoid it leaves directional
415high-degree bands with amplification $> 1$ — patterns went qualitatively high-frequency at
416moderate settings and diverged as niter or $L$ grew. With the symbol-based constant every
417niter/geometry combination in the test sweep converges (peanut: $\mu \in [0.44, 6.2]$, plain rate
4185.2, preconditioned rate 0.87).
420**The correction must be band-projected.** Algorithm 5's "zero $\ell \ge L-2$" is load-bearing:
421without applying the same mask to the correction $d\Delta u$, the top two degrees iterate toward
422the *undiffused* right-hand side — each Richardson iteration strips more of their implicit
423diffusion, at species-dependent rates, manufacturing a spurious Turing band at the band edge
424(observed on the round sphere: top-degree energy growing $\sim 3\%$/step at $L=127$, 8 iterations).
426**Payoff shape:** on mildly deformed surfaces one iteration already reaches the $\sim10^{-4}$ fp32
427accumulation floor, so iteration counts do not drop — the speculation above does not hold at fp32.
428The gain is reach and correctness: stiff geometries and high niter/$L$ combinations that
429previously diverged (or silently shifted the pattern's wavelength) now converge with a
430resolution-independent spectrum.