ConformalLabpp/doc/reviewer/numerical-stability-audit-2026-05-31.md

Numerical-Stability Audit — ConformalLabpp

Date: 2026-05-31 Auditor: External reviewer (Claude Opus 4.8) Scope: Tolerance hierarchy, conditioning, catastrophic cancellation, and the hard-coded "magic" constants across code/include/. Focus: Floating-point robustness of a numerical library — distinct from port-faithfulness (java-port-audit.md) and from internal error propagation (test-coverage-error-handling-audit-2026-05-31.md).

Status legend: 🔴 Critical · 🟡 Important · 🔵 Polish

✅ Resolution status (2026-05-31, Session 1): N3 ✅ (selectable clamp mode HardJava|SmoothBarrier), N4/N6 ✅ (named constexpr constants), N5 ✅ (Kahan stable triangle area), N7 ✅ (S2 — min_ldlt_pivot / sparse_qr_fallback_used in NewtonResult). N1 ✅* largely subsumed by the B1 block-FD Hessian; the remaining FD-step/tol coupling note folds into N2. N2 (tolerances.md) open — scheduled S4 in finding-orchestration.md. 301/301 tests green.

Relationship to existing audits: java-port-audit.md already verified several per-formula numerical edge cases (degenerate triangles #1, overflow centring #10). This audit takes the cross-cutting view those miss: how the tolerances relate to each other, where constants are hard-coded vs. centralized, and where the formulas can lose precision. No finding here duplicates java-port-audit.

---

Summary table

ID	Sev	Title	Location
N1	🔴	FD-Hessian step (1e-5) sets an accuracy floor (~1e-10) too close to the Newton tol (1e-8) — can stall quadratic convergence	newton_solver.hpp:408,569
N2	🟡	Tolerance hierarchy is undocumented and uncoordinated — 7+ independent literals with no stated relationship	many headers
N3	🟡	Discontinuous clamp b<0 → 0.01 breaks gradient/Hessian continuity at the domain boundary	hyper_ideal_functional.hpp:179,274
N4	🟡	Magic constants (-30.0, 0.01, 1e-15, 1.0-1e-15) are unnamed, undocumented, and scattered	spherical_functional.hpp, spherical_geometry.hpp, hyper_ideal_functional.hpp
N5	🟡	Cotangent/area formula is cancellation-prone for needle/cap triangles	euclidean_hessian.hpp:81-93
N6	🔵	constants.hpp holds only π/2π — all tolerances live as bare literals	constants.hpp
N7	🔵	No conditioning diagnostic on the Hessian before the linear solve	newton_solver.hpp

---

N1 — 🔴 FD-Hessian step size fights the convergence tolerance

Evidence

• newton_hyper_ideal and newton_inversive_distance build the Hessian by finite differences with hess_eps = 1e-5 (newton_solver.hpp:408, :569).
• Both solvers converge against tol = 1e-8 on ‖G‖∞ (:406, :567).

Problem

A central-difference Hessian with step h = 1e-5 has error O(h²) + O(ε_machine/h). With h = 1e-5 the truncation term is ~1e-10 and the round-off term is 1e-16/1e-5 = 1e-11; the achievable Hessian accuracy floor is therefore **1e-10–1e-11**. That is only ~2–3 decades below the convergence target tol = 1e-8. Near the solution, where Newton should converge quadratically, an inaccurate Hessian degrades the step to linear convergence and can stall the iteration just above tol — exactly the regime test_newton_phase9a.cpp exercises.

This is also why the line-search "stall" exit (FINDING-H1 / error-handling audit) can trigger on otherwise well-posed problems: the FD Hessian, not the geometry, is the limiting factor.

Fix

• Primary: replace the FD Hessian with the analytic / block-FD paths (cross-ref api-performance-audit B1 — the HyperIdeal block-FD already exists and is more accurate as well as faster).
• Until then: document the coupling, and consider a richer FD stencil or Richardson extrapolation if FD must stay. A relative step h = 1e-6·(1+|x_i|) is more robust than an absolute 1e-5.

Acceptance criteria

• The FD-step/tol relationship is documented, or the analytic/block-FD Hessian replaces FD so the floor drops well below tol.
• A test demonstrates quadratic (not stalled) convergence near the solution.

---

N2 — 🟡 Tolerance hierarchy is undocumented and uncoordinated

Evidence (the full set of independent literals)

Constant	Value	Location	Meaning
Newton tol	1e-8	newton_solver.hpp (all 5 solvers)	convergence on ‖G‖∞
FD hess_eps	1e-5	newton_solver.hpp:408,569, hyper_ideal_hessian.hpp:67	Hessian FD step
gradient_check eps/tol	1e-5 / 1e-4	hyper_ideal_functional.hpp:474-475	test-only FD check
Armijo c1	1e-4	newton_solver.hpp:149	sufficient-decrease
sparsity drop	1e-15	newton_solver.hpp:599	triplet prune threshold
Gauss-Bonnet tol	1e-8	gauss_bonnet.hpp:134,154	χ residual
spherical clamp	1e-15	spherical_geometry.hpp:34	asin domain guard

Problem

These seven thresholds were each chosen locally (several "to match the Java FunctionalTest"), with no documented relationship. Yet they interact: N1 shows hess_eps vs tol; the sparsity drop 1e-15 interacts with hess_eps (FD values below 1e-15 are pruned, but FD noise is ~1e-11, so the prune threshold is effectively never the limiting factor — possibly dead). A reviewer cannot tell which are load-bearing and which are arbitrary.

Fix

Add a short doc/math/tolerances.md (or a section in constants.hpp) that states each tolerance, its role, and the constraints between them (e.g. "FD floor must be ≪ Newton tol", "Armijo c1 ∈ [1e-4, 1e-1]"). Centralize them as named constexpr.

Acceptance criteria

• Every numerical threshold is named, documented, and its relationship to the others stated.

---

N3 — 🟡 Discontinuous clamp breaks gradient/Hessian continuity

Evidence

hyper_ideal_functional.hpp:179-181 and :274-276:

if (v1b && b1 < 0.0) b1 = 0.01;   // negative scale → snap to 0.01

Problem

This snaps any negative b to a constant 0.01. The map b ↦ max(b, 0.01) is continuous but not differentiable at b = 0, and the snap-to-constant makes the returned angles locally independent of b for b < 0 — so the analytic gradient and the FD Hessian disagree across the boundary, and Newton steps that cross b = 0 get an inconsistent local model. The comment says it "mirrors the Java original", but the Java code is a reference, not a correctness proof: a hard clamp inside the function being differentiated is a known source of Newton stalls.

Fix

• Prefer keeping iterates in the feasible region via the line search (reject steps that drive b < b_min) rather than clamping inside the energy evaluation.
• If clamping must stay, use a smooth barrier / softplus so the derivative is continuous, and make b_min a named constant (N4).

Acceptance criteria

• The energy/gradient used by Newton is C¹ across the feasibility boundary, or the line search keeps iterates strictly feasible so the clamp never fires.

---

N4 — 🟡 Unnamed magic constants

Evidence

• spherical_functional.hpp: m.lambda0[e] = -30.0; (very-short-edge floor)
• spherical_geometry.hpp:34: if (half >= 1.0) half = 1.0 - 1e-15;
• hyper_ideal_functional.hpp: 0.01 clamp (N3), b < 0.0 thresholds
• newton_solver.hpp:599: if (std::abs(val) > 1e-15)

Problem

Each is a bare literal with at most an inline comment. -30.0 (≈ exp(-15), i.e. edge length ~3e-7) and 1.0 - 1e-15 (the asin/acos domain guard) encode real numerical reasoning that is invisible to the next maintainer and impossible to tune consistently.

Fix

Promote to named constexpr in constants.hpp with a one-line rationale each, e.g. constexpr double LOG_LENGTH_FLOOR = -30.0; // exp(-15): edge below ~3e-7 treated as 0.

Acceptance criteria

• No bare numeric tolerance/floor literals remain in the math headers; each is a documented named constant.

---

N5 — 🟡 Cotangent/area formula is cancellation-prone

Evidence

euclidean_hessian.hpp:81-93:

const double t12 = -l12 + l23 + l31;
const double t23 = +l12 - l23 + l31;
const double t31 = +l12 + l23 - l31;
...
const double denom2 = 2.0 * std::sqrt(t12 * t23 * t31 * l123);   // = 8·Area

Problem

For a needle triangle (one edge ≈ sum of the other two), one of t12/t23/t31 is the difference of nearly-equal lengths → catastrophic cancellation, and the area via sqrt of the product loses precision. The cotangent weights then carry large relative error, which feeds directly into the Hessian and the linear solve. The <= 0 guard catches the fully-degenerate case but not the ill-conditioned-but-valid one.

Fix

Use a numerically stable area formula (Kahan's formula for triangle area from sorted side lengths) instead of the raw product-under-sqrt. This is a well-known, drop-in improvement for cotangent-Laplacian assembly.

Acceptance criteria

• A sliver-triangle test shows the cotangent weights match a high-precision (e.g. long double / Kahan) reference to acceptable relative error.

---

N6 — 🔵 constants.hpp holds only π

Evidence

constants.hpp defines exactly PI and TWO_PI. Every tolerance from N2 lives elsewhere as a literal.

Fix

Extend constants.hpp (or a sibling tolerances.hpp) to be the single source of truth for the N2/N4 constants. Tie to N2's documentation.

---

N7 — 🔵 No conditioning diagnostic before the linear solve

Evidence

solve_with_fallback (newton_solver.hpp:65) tries LDLT then SparseQR but never reports the conditioning of H.

Problem

On an ill-conditioned-but-non-singular Hessian, LDLT "succeeds" and returns a low-accuracy step with no signal. The sparse_qr_fallback_used flag only fires on outright LDLT failure, not on near-singularity.

Fix

Optionally expose an estimate (e.g. smallest pivot magnitude from LDLT, or a cheap condition estimate) in the diagnostics, so callers/tests can flag near-singular solves. Low priority; pairs naturally with the FINDING-I1 status enum (error-handling audit).

---

What is already good

• constants.hpp centralizes π at 30 digits — the right instinct, just not applied to tolerances (N6).
• The spherical asin domain guard (spherical_geometry.hpp:34) exists at all — many implementations forget it (it just needs a named constant, N4).
• Central differences (not forward) are used for FD gradients/Hessians — the correct O(h²) choice.
• java-port-audit.md #10 already flagged the inversive-distance edge-length overflow centring — consistent with this audit's spirit; coordinate the two.

Suggested order

1. N1 (FD step vs tol) — highest correctness impact; largely subsumed by the api-performance-audit B1 analytic/block-FD work.
2. N3 (discontinuous clamp) — fixes a real Newton-stall source.
3. N5 (stable area) — drop-in robustness for the Euclidean Hessian.
4. N2 + N4 + N6 (document + centralize tolerances/constants) — do together.
5. N7 (conditioning diagnostic) — polish, pairs with I1.