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). **N1 ✅\*** largely subsumed by the B1 block-FD
> Hessian; the remaining FD-step/tol coupling note folds into **N2**. **N2** (tolerances.md)
> and **N7** (conditioning diagnostic) **open** — scheduled S2/S4 in
> [`finding-orchestration.md`](finding-orchestration.md). 298/298 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`:
```cpp
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`:
```cpp
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.