reg_field/reg_flow: diag_bending omits fold-back term under folding boundary conditions

[balbasty]

Summary

diag_bending's (and diag_membrane's) computed diagonal does not equal the true M[i][i] of the materialized operator matrix under a folding boundary condition (e.g. DCT2): the diag_* family computes b00 - Σ m_k·f_k and omits the fold-back term +Σ_{k: n(k)=i} m_k that should apply when a folded neighbour lands back on the centre voxel.

Measured discrepancy: ~13.6 (2D, DCT2) / ~49.4 (3D, DCT2) for diag_bending; ~11.0 (2D, DCT2) / ~1.1 (3D?) for diag_membrane. Exact under DFT (diff ~3.6e-15), where no folding onto the centre occurs for the grid sizes tested.

Scope

This is pre-existing and general — it affects the plain (unweighted) diag_bending/diag_membrane, not something introduced by the RLS/JRLS work or by the fastfields/fastfields-kernels#38 self-adjointness fix. diag_bending_rls/_jrls inherit the identical mismatch (no better, no worse) from the plain operator.

Is this actually a bug?

Possibly intentional — diag_* is documented/used as a preconditioner diagonal for relaxation solvers, where an approximate diagonal (ignoring the fold-back correction) is a common and acceptable simplification, not a correctness requirement in the same sense as self-adjointness. Needs someone with more context on the original design intent to decide whether:

1. This is a deliberate approximation (document it as such and close), or
2. diag_* is meant to expose the exact matrix diagonal (e.g. for direct solvers or diagnostics) and needs the fold-back term added.

Reproduction

Materialize the full matvec_bending/matvec_membrane operator (unit impulses, no weight map) under DCT2 boundary on a small grid, compare M[i][i] against diag_bending(i)/diag_membrane(i) at boundary voxels where a neighbour folds back onto the centre.

Surfaced during the investigation for fastfields/fastfields-kernels#34 (bending RLS/JRLS self-adjointness); not folded into that fix (#38) since it's unrelated to self-adjointness and pre-dates the RLS/JRLS work entirely.

Workstream: claude-jitfields-to-fastfields

[balbasty]

Triage for the 0.0.dev0 release: confirmed non-blocking, staying open/blocked pending the design decision.

Traced where diag_bending/diag_membrane are actually consumed (fastfields-cpu-impl/reg_field.h):

• relax_bending_ (the implementation behind the shipping field_relax op) computes the residual via matvec_bending<isub> — the exact, self-adjoint operator (fixed under regularisers: fix bending RLS/JRLS cross-term self-adjointness bug fastfields-kernels#38) — and only feeds diag_bending<set> in as the local block fed to PosDef::relax_ for that sweep's update.
• That's a stationary (Richardson/block-Jacobi-style) iteration: x ← x + M⁻¹(g − Ax) with A = the exact operator (via matvec) and M = the local block built from diag_bending. The fixed point of that iteration is Ax=g regardless of what M is, as long as M is nonsingular and the iteration converges — M only affects convergence rate/stability, not the limit. So an approximate (fold-back-omitting) M doesn't corrupt the converged result of field_relax; it can only make convergence slower/less stable in the folding-boundary case.
• field_diag is also exposed standalone in the binding, but its own docstring already frames it as "Diagonal (preconditioner) of the field regulariser operator" — consistent with the approximation being the documented contract, not exact-diagonal semantics.
• test_field_relax_reduces_residual (fastfields-dlpack tests) validates convergence via the matvec-based residual, not via diag_bending, which matches this reading and is presumably why it passes despite this issue being open.

So this is confined to case (b) in your own framing — it's the preconditioner path, not the exact-diagonal-for-direct-solvers path, and it does not corrupt any path that ships and is used at runtime today (no direct solver consumes diag_* as an exact diagonal in this codebase yet). I did not attempt to resolve the design-intent question (deliberate approximation vs. exact-diagonal support) — that's still yours to decide. Leaving the bug/blocked labels as-is.

---

Generated by Claude Code