Lamé flow regulariser is not self-adjoint at reflecting boundaries (critical)

[balbasty]

Part of #25. Critical / tricky — needs review before merge.

The linear-elastic (shears/div) flow operator was not self-adjoint under reflecting boundaries (DCT2/DST*): the cross-channel block was computed as D_c D_d with the field's natural boundary reads instead of D_c^T D_d, so ⟨Lx,y⟩ ≠ ⟨x,Ly⟩ except under DFT — breaking the SPD assumption of CG/relaxation solvers (and torch's self-adjoint backward).

Fix (fastfields-kernels/regularisers/flow/{2d,3d}.h, matvec_all + matvec_lame): read the transposed axis through the adjoint boundary via a new bound::transpose() (DCT1↔DST1, DCT2↔DST2; DFT/Zero/Replicate unchanged), and correct a corner sign. Interior stencil unchanged.

Verified by new operator-symmetry tests in fastfields-cpu-lib (2D + 3D, all penalty combos, DCT2/DFT/mixed) and torch gradcheck.

Upstream (jitfields shares the buggy kernel): balbasty/jitfields#5. fastfields intentionally diverges to the correct operator.

[balbasty]

Review outcome (via extended symmetry tests).

Reviewing the fix across all reflecting boundaries surfaced a real scope limit:

• Self-adjoint (verified, 2D+3D): DCT2, DST2, DFT, Zero — ⟨Lx,y⟩ == ⟨x,Ly⟩ to round-off.
• NOT self-adjoint: DCT1 / DST1 (whole-sample reflect). Their forward and adjoint use different reflection centres (reflect_{N-1} vs reflect_{N+1}), so a single companion-boundary read cannot reproduce D^T; a same-index transpose map (as used for DCT2↔DST2) can't make it SPD.

Decision: keep the fix (exact for the half-sample family + DFT/Zero, which covers the boundaries flow regularisation actually uses — DCT2/Neumann and DFT), and document DCT1/DST1 as a known non-SPD limitation rather than overclaim. Recorded at the bound::transpose() definition (fastfields-kernels/bounds.h) and in the cpu-lib symmetry test. A proper DCT1/DST1 adjoint would need a distinct (non-companion) discretisation — tracked as future work under #25.

---

Generated by Claude Code