Skip to content

fix(#130): make TP evidence norms underflow/overflow-safe - #136

Merged
marcohost33-maker merged 25 commits into
mainfrom
chatgpt/issue-130-scaled-norms
Sep 10, 2026
Merged

marcohost33-maker merged 25 commits into
mainfrom
chatgpt/issue-130-scaled-norms

Conversation

@marcohost33-maker

Copy link
Copy Markdown
Owner

Summary

Closes #130 on top of PR #121 (base=pr107-fix) as a focused numerical-evidence correction.

trace_preservation_defect() used ordinary Euclidean/Frobenius norms. Their internal sum of squares can underflow to 0.0 or overflow to inf for finite matrices even when the mathematical norm is representable. Because those values feed the trace-preservation applicability decision, arithmetic failure could become scientific evidence.

Fix

  • add numerics.norms.scaled_euclidean_norm();
  • use one shared norm contract for both ||vec(I)^H L||_2 and ||L||_F;
  • scale by an exact power of two before squaring, then restore with ldexp;
  • handle real and imaginary components separately so selecting the scale does not itself overflow on a finite complex value;
  • return inf only when an input component or the true norm is non-representable, and propagate NaN explicitly.

The design follows LAPACK xLASSQ's numerical contract: scaled sum-of-squares accumulation to avoid both overflow and underflow. Reference: https://www.netlib.org/lapack/explore-html/d8/d76/group__lassq_gae8f40b0a34771b4f2d9c863de3af7be5.html

Measurement vs decision policy

This PR deliberately removes the old coupling where a naive Frobenius overflow was used as an implicit refusal signal. A finite, representable norm must be measured as finite. Whether the resulting zero-mode evidence is usable is a separate structural decision.

The historical +-1e308 counterexample remains fail-closed: after the Frobenius norm is correctly measured at about sqrt(2)*1e308, the operator passes the relative TP applicability comparison because its trace defect is tiny relative to its scale, but the zero-mode band swallows the entire measured spectrum. band_discriminates() therefore refuses certification. Safety is retained at the layer that actually has the evidence to make that decision.

Tests

New tests/test_scaled_norms.py covers:

  • ordinary tiny values around 1e-200;
  • subnormal values around 1e-320;
  • large values whose naive squares overflow but whose true norm is finite;
  • a genuinely non-representable true norm -> inf;
  • complex values near float64 range limits;
  • scale equivariance across 500 orders of magnitude;
  • explicit zero/NaN/inf semantics.

tests/test_pr121_review_round23.py is updated so the former overflow-as-policy test now requires the representable norm to be finite and requires the structural certificate to remain unresolved. Existing tiny non-TP and healthy-GKSL controls remain.

Scope boundary

Change-control note

As with #134, CHANGELOG/CITATION edits are left for the integration point rather than adding another conflict to the large #121 stack. Final integration must record this results-correctness change with exact-head CI evidence.

@chatgpt-codex-connector

Copy link
Copy Markdown

You have reached your Codex usage limits for code reviews. You can see your limits in the Codex usage dashboard.
To continue using code reviews, you can upgrade your account or add credits to your account and enable them for code reviews in your settings.

Copy link
Copy Markdown
Owner Author

Current exact-head verification on d4d635f34b0a4fa94eae926693dd7abd626d1403 found a deterministic formatting failure, not a scientific-test failure: reusable CI stops at Ruff W292 because three rewritten files lack a terminal newline (src/liouscope/numerics/linalg.py, tests/test_pr121_review_round20.py, tests/test_pr121_review_round23.py). Quality Contract, QuTiP and CodeQL are green. This is being treated as a real post-write integrity defect, not a CI flake. Do not call #136 verified until the EOF defects are corrected and exact-head CI is rerun.

Copy link
Copy Markdown
Owner Author

Exact-head verification is complete on 59f861a571c5ee06d3c5aa2917079fe10eb22641. The EOF repair was byte-preserving: exactly one terminal newline was appended to each of the three W292 files, the one-shot repair workflow self-deleted, and the external marker deletion produced the clean candidate head. On that exact head, CI (Python 3.10–3.14), Quality Contract, QuTiP cross-checks, CodeQL, and ci-reusable-pilot all completed successfully. This verifies the #130 slice itself; it does not by itself promote the underlying pr107-fix stack.

marcohost33-maker and others added 3 commits September 3, 2026 00:40
…underflowed

Three findings of the round-1 external review on PR #147.

[gls.py:227] When the whitened residuals carry a finite, non-zero RMS below
the smallest positive float64, ``log_rss`` and ``log_sigma`` are both finite
and the profile likelihood is computable -- but ``exp(log_sigma)`` underflows
to 0.0, and this branch returned success=False, log_likelihood=nan,
likelihood_degenerate=True. ``_fit_with_model`` then scored the model inf and
dropped it from selection. That is the ABSOLUTE SCALE BOUNDARY issue #135 set
out to remove, put back one layer down: the same curve in different rate units
either is or is not an AICc candidate.

Reproduced at 6ddd053 before any change, on the reviewer's construction: one
minimum-subnormal residual (5e-324) among 64 otherwise-zero points gives
log_rss = -1488.8801438427624, log_sigma = -746.519513463061 and a profile
log-likelihood of +47686.4367955108, while exp(-746.52) underflows because the
smallest subnormal is exp(-744.44).

The fit now keeps its finite log-space likelihood and stays selectable; only
the scale-dependent evidence is withheld. GLSFitOutput gains
``scale_unavailable`` (additive, default False). ``sigma`` is NaN there,
deliberately NOT 0.0: ``_ar1_resample`` consumes it as the innovation standard
deviation, and 0.0 would draw identical replicates -- a zero-width confidence
interval, which is the failure mode of an uncertainty pipeline rather than a
conservative one. ``parametric_bootstrap`` refuses on the new flag rather than
on ``success``, so "no interval" and "no estimate" stay distinguishable, and
``compute_relaxation_layer`` reports the CI as NaN through its existing
RuntimeError handler.

[aicc.py:40, twice] The likelihood replacement is a methodology change that
can reorder AICc, change the selected model and change the reported rate, and
the commit recorded it in neither place AGENTS.md requires. CITATION.cff now
carries the pending correction (AGENTS.md:121-123) and CHANGELOG.md an
[Unreleased] "Changed" entry for the likelihood itself plus a "Fixed" entry
for the underflow repair above (AGENTS.md:195-197). Both state what the change
is NOT: only the likelihood path was corrected, model selection is not
scale-free in general, and per-model likelihoods are not manifest fields, so
the behavioural change is invisible in input_hash -- archived analyses whose
selected model matters must be re-run.

Discrimination proven per GUARD LINE, 3/3, run DK-20260902T224741-cc2ae3b913bd
(ledger Vero/Data/diskriminierung_runs.jsonl): the underflow detection, the
NaN sigma and the bootstrap refusal each turn their test red at an assertion
when removed, green again after restore. The first attempt scored ROT DURCH
ABSTURZ because the test asserted the warning via ``pytest.warns``, whose
failure is not attributable to a guard line; it now records warnings and
asserts on them.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01FPKDUXBV4qXpWX8ZtpqgAZ
…refused by arithmetic

Round-2 external review of PR #147.

- _types.py / relaxation.py: FitResult carries `scale_unavailable` and
  `_fit_with_model` copies it. Without it the persisted report showed only
  `bca_ci_beta = (nan, nan)` -- the value ANY bootstrap or jackknife failure
  produces -- and the reason lived in a warning the artefact drops.
- aicc.py: the standardised RSS enters as `-0.5 * RSS`, so its representable
  range reaches `2 * float64.max`. The guard compared against `log(max)` and
  returned -inf for the octave above it; the bound now carries `+ log(2)` and
  the exponential is taken after subtracting `log(2)` where that is needed.
  Below that octave the arithmetic is untouched, so no existing likelihood
  moves by even one ulp (pinned by a bit-identity test).

Beleglauf folgt (volle Suite, rund 5 min).

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01FPKDUXBV4qXpWX8ZtpqgAZ
marcohost33-maker and others added 3 commits September 3, 2026 02:34
…OverflowError

The mutation run reported ROT DURCH ABSTURZ:OverflowError for this control:
without the guard, math.exp raises before any assertion runs, so the death was
not attributable and the case did not count as a proof. The outcome is captured
and its TYPE checked with isinstance instead.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01FPKDUXBV4qXpWX8ZtpqgAZ
…n two logs

Thread PRRT_kwDOSXoNQ86evTEG, src/liouscope/fitting/aicc.py:90.

Round 2 widened the bound to log_max + log(2) and the boundary was still
decided in LOG SPACE, where at this magnitude it cannot be decided.

Measured, sigma=1 with a single residual sqrt(float_max)*sqrt(2):
  log_rss          = 710.475860073944
  log_max + log(2) = 710.4758600739439   -- exactly ONE ulp lower
so the guard returned -inf, although the quantity the formula actually forms
is 1.7976931348623155e+308 = 0.9999999999999999 * float_max, finite.

The obvious smaller repair does not work either, and the code and a test now
say so: exp(log_rss - log(2)) raises OverflowError on this same input,
because one ulp of a logarithm near 710 is a factor of ~1e-16 in the value --
exactly the width of the decision being made.

So the halved sum is MATERIALISED from the residuals: halving each term
before accumulating keeps every partial sum inside float64 whenever the true
half-RSS is, and float64 addition itself decides the boundary. There is no
rounded threshold left to be one ulp wrong about, and an overflow is now a
measurement rather than a prediction.

Verified after the change:
  reviewer's fixture      -> -1.7976931348623155e+308 (was -inf)
  4x float_max            -> -inf
  1e3x float_max          -> -inf
  64 residuals overflowing only in the SUM -> -inf   (a per-term bound would
                                                      have missed this one)
  ordinary/sigma=None/empty/nan/inf/all-zero inputs unchanged
The lower octave keeps 0.5 * exp(...) untouched, so no previously computable
likelihood moves by one ulp -- the promise round 2 made.

Baseline pr147-work @ cee633f: 1094 passed, exit 0.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01FPKDUXBV4qXpWX8ZtpqgAZ
…teness

Follow-up inside the same finding, found by the mutation run
DK-20260903T074723-246b5d512bf2, which reported the finiteness check BLIND.

Measuring why exposed a latent state collapse rather than dead code:
  half = inf -> without the check the subtraction already yields -inf
                (measured), so the check is provably equivalent there
  half = nan -> without the check the function returns nan; WITH a blanket
                isfinite() check it would return -inf
and those are different answers. This function returns NaN by contract for
"model-selection likelihood unavailable"; -inf means "dropped from
selection". A guard that silently turns the first into the second is the
same substitution this whole review wave is about.

NaN is unreachable on this path -- a NaN residual is caught by the
isnan(log_rss) gate above, and every term 0.5*s*s is non-negative -- so the
condition now states the invariant it actually means.

M7 stays BLIND after the change and now provably so, rather than by
accident. 19/19 in tests/test_issue135_likelihood_scale.py, ruff clean,
mypy clean.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01FPKDUXBV4qXpWX8ZtpqgAZ
@marcohost33-maker

Copy link
Copy Markdown
Owner Author

Status 2026-09-07 — conflicting, and it is the largest unique body of work in the stack

mergeStateStatus  DIRTY   (mergeable = CONFLICTING against pr107-fix)
head 59f861a5 -> check-runs: 13, test-matrix jobs: 5, non-success: 0
vs pr107-fix@bac5fe0d : 13 commits missing here, 16 unique to this PR

Two things follow from those numbers.

1. The green matrix is real but earned on a superseded tree. Those 13 runs come from the now-closed exact-head verification PR #137, which carried this SHA against main. So this is not one of the CI-blind PRs — #131 and #147 are. But the tree it certifies is 13 commits behind pr107-fix, and this PR conflicts with its own base, so the rollup describes something that cannot currently be merged anywhere.

2. Sixteen unique commits is the most in the stack, and three other PRs sit on top of this branch: #139 and #147 both declare chatgpt/issue-130-scaled-norms as base, and #149 carries #139's head. A conflicting base propagates upward — resolving it here is the unblocking action for four PRs, not one.

No issue band added: the base is pr107-fix, not main, so Closes #130 would be structurally inert. It belongs in the body only once this PR targets main.

What Marco has to decide

Order matters more than usual here: #121 -> then rebase this branch -> then #139/#147/#149 on top. Rebasing this one first, before #121 lands, would mean doing the conflict work twice.

Read-only measurement. No push, no merge, no rebase, no re-run, nothing resolved.

…e-safe-likelihood

fix(#135): make Gaussian profile likelihood scale-safe
Base automatically changed from pr107-fix to main September 8, 2026 14:32
Two conflicts, both resolved in favour of the numerically correct side.

src/liouscope/numerics/linalg.py -- trace_preservation_defect. Both sides
solve the SAME problem (norms that survive under-/overflow) differently:
main keeps np.linalg.norm plus two one-directional underflow rescues and
deliberately leaves overflow saturated, because the round-21 refusal of a
non-finite reference scale used `fro == inf` as its trigger. #136 replaces
both norms with the xLASSQ-style scaled_euclidean_norm and refuses the
round-20 counterexample structurally, via
trace_preservation_componentwise_error == 1.

Kept #136. The deciding measurement is arithmetic, not preference: for the
round-20 operator diag(1,2,3,4) with L[0,1]=+1e308, L[3,1]=-1e308 the true
Frobenius norm is sqrt(2)*1e308 = 1.4142135623730950643e308, verified
against exact rational arithmetic at 60 decimal digits. It is REPRESENTABLE
in float64; np.linalg.norm returning inf is a defect, not a policy.

tests/test_pr121_review_round24.py -- consequence, not cosmetics.
test_the_one_directional_overflow_policy_is_untouched asserted
`not isfinite(fro)` on that exact matrix, while #136's
tests/test_pr121_review_round23.py asserts `isfinite(fro)` on the byte-
identical matrix. No implementation can satisfy both. The round-24 test
pinned the MECHANISM (an overflow) rather than the property it existed to
protect (the operator stays refused). Rewrote it to pin the property:
cert.applicable is False still holds, now carried by the componentwise
backward error. Measured, not assumed.

CHANGELOG.md -- additive bullet lists under one `### Fixed` heading; union.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01VG6CeFJVgem6n4Tsam1onb
marcohost33-maker added a commit that referenced this pull request Sep 10, 2026
One conflict, identical to the one in #136: linalg.py trace_preservation_defect.
#139 branched off #136 at 59f861a, so it carries the same scaled-norm side.
Resolved the same way and verified byte-identical: the resolved file hashes
b545eca04914e2c4d11502c6d5da52224812cb44e21a16dee83957161dbf11a9, the same
blob as in aef4ce9. Justification is recorded there (exact-rational oracle:
the round-20 operator's true Frobenius norm is representable, np.linalg.norm
returning inf was the defect, and the protected property -- the operator
stays inapplicable -- is carried by the componentwise backward error).

tests/test_pr121_review_round24.py taken from aef4ce9 for the same reason.

Verified the #113 claim rather than trusting it. The body says the primitive
removes ONE left null vector and that further traceless zero modes of a
degenerate stationary manifold correctly remain in L0. Measured:

  amplitude damping (unique)   dim ker(L)=1 -> 0 zero modes left in L0
  pure dephasing (2-dim ker)   dim ker(L)=2 -> 1 zero mode  left in L0
  two decoupled blocks (d=4)   dim ker(L)=4 -> 3 zero modes left in L0

k_reduced = k_full - 1 throughout, and the nonzero spectrum is preserved
exactly in every case (12 of 12 for d=4). The code does not pretend
uniqueness. The shipped test only covers the 2-dim case; the d=4 case above
is a stronger control that the claim survives.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01VG6CeFJVgem6n4Tsam1onb
@marcohost33-maker
marcohost33-maker merged commit 38f8652 into main Sep 10, 2026
13 checks passed
@marcohost33-maker
marcohost33-maker deleted the chatgpt/issue-130-scaled-norms branch September 10, 2026 18:27
marcohost33-maker added a commit that referenced this pull request Sep 10, 2026
…iring it

Baseline 33e2602: 1201 passed, exit 0,
coverage 94.35%.

Reproduced the P2 finding before touching the code. For `d = 16`, first column
carrying `+2.5e307` in the first eight trace rows and `-2.5e307` in the
remaining eight:

  vec(I)^H @ L, entry 0                  = (inf+nanj)   [overflow + invalid]
  trace_preservation_defect              = nan, scale = 1e308
  trace_preservation_componentwise_error = 0.0
  restrict_to_traceless                  -> ValueError, "not representable as
                                            finite float64"

This is the mirror image of the finding repaired in 5c151f9. That one was
fail-OPEN on an invalid generator; this is fail-CLOSED on a valid one, and the
scale-safe componentwise gate that reads the operator correctly never runs,
because the normwise gate raises first.

Checked against exact rational arithmetic rather than against another float
computation, following the reasoning PR #136 used for the Frobenius norm: the
column sum is EXACTLY zero as a Fraction, and the Frobenius norm is exactly
`4 * 2.5e307 = 1e308`, inside float64 range. The generator is legal and fully
representable; only an intermediate product overflows.

Worth recording because it is the same class one level up: the first oracle
written for this check was an ordinary sorted Python `sum`, and it overflowed
to `-inf` itself. An oracle computed the same way as the thing it judges is
not an oracle.

Root cause is the SHARED helper `linalg.py:331`, `vec_i.conj() @ L_super`, on
which both certified eigensolver paths also depend -- so the repair belongs
there and not as a second copy in traceless.py.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01VG6CeFJVgem6n4Tsam1onb
marcohost33-maker added a commit that referenced this pull request Sep 11, 2026
…139)

* feat(#113): add traceless-space restriction primitive

* chore(#113): record the reproduction of review finding B1 before repairing it

Baseline 82a30b0: 1195 passed, exit 0
(Python 3.12.10, numpy 2.5.3, scipy 1.18.1).

Reproduced the P1 review finding on `restrict_to_traceless` before touching
the code, so the repair is evidenced rather than believed. For the 4x4 input
`L[1, 0] = 1e300`, `L[0, 2] = 1`:

  global defect     = 1.000000e+00
  operator_scale    = 1.000000e+300
  global ratio      = 1.000000e-300   (default tp_rtol = 1e-10)  -> gate passes
  componentwise err = 1.000000e+00                               -> gate should fail
  RESULT: a TracelessRestriction is returned, no fail-closed.

Sharper than the report: the returned restriction carries its own refutation.
`invariance_defect = 7.071068e-01` and `reconstruction_defect = 7.071068e-01`
are order one, i.e. the module MEASURES that the traceless subspace is not
invariant under this generator and draws no consequence from it. The normwise
ratio 1e-300 is not a wrong reading of the trace equation for column 2; it is
no reading of it at all, because the unrelated 1e300 entry sets the scale that
divides it. The gate carries zero bits here.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01VG6CeFJVgem6n4Tsam1onb

* fix(#139): gate the traceless restriction on the componentwise trace equation

PR #139 review, finding B1. `restrict_to_traceless` admitted a generator on
the normwise ratio `||q^H L|| / ||L||` alone. That ratio divides ONE number by
the norm of the whole operator, so an entry taking part in no violated trace
equation dilutes a violation elsewhere until the quotient is round-off. The
gate then carries no information about the equation it is supposed to be
about.

Measured before the change, on the 4x4 input `L[1, 0] = 1e300`, `L[0, 2] = 1`:
global ratio 1e-300 (accepted at the default `tp_rtol = 1e-10`), while the
trace equation for column 2 has relative error 1. The restriction handed back
carried its own refutation: `invariance_defect = 7.07e-01` and
`reconstruction_defect = 7.07e-01`. The module measured that the traceless
subspace was NOT invariant under the generator it had just accepted -- and the
invariance of that subspace is the single assumption the construction rests on.

Issue #130 already built the correct reading, and both certified eigensolver
paths already gate on it (`linalg.py:1116-1117` and `linalg.py:1392-1393`).
This uses THE SAME comparison rather than a second, slightly different one:
`trace_preservation_componentwise_error` is already dimensionless, so it is
compared directly against `tp_rtol` with no scale factor, and a non-finite
reading is refused instead of compared. Two gates that measure one quantity
differently are a defect class of their own; a test asserts the agreement.

`TracelessRestriction` now carries `trace_componentwise_error` so the reading
that decided admission is auditable rather than recomputed by callers.

Not a behaviour change for legal input, measured: the componentwise error is
0.0 for the amplitude-damping and dephasing fixtures at every rate scale
tested (1e-100 .. 1e100), and 7.8e-18 / 2.3e-17 / 1.8e-17 for the stiff jump
network at fast rates 1e8 / 1e10 / 1e12 -- all far below the default
tolerance.

Positive controls, because a gate that refuses everything also passes every
negative test:
- the refusal fixture with only the offending entry removed, 1e300 and all,
  is still accepted (`trace_componentwise_error == 0.0`, invariance defect 0);
- the physical stiff generator is accepted at the default tolerance AND
  refused once `tp_rtol` drops below its measured nonzero reading, which is
  what shows the gate is on the code path rather than vacuously passing.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01VG6CeFJVgem6n4Tsam1onb

* docs(#113): record the traceless-restriction methodology in the citation record

PR #139 review, finding B2. AGENTS.md:121-123 requires that changes touching
results or methodology update `CITATION.cff`, and the relevant
`MANIFEST_SCHEMA` version if the run manifest contract changes. The commit
that introduced `numerics/traceless.py` left both untouched.

Written to the convention CITATION.cff states for itself (added 2026-08-09):
`abstract` describes the surface of the CITED release only, and capability
sentences for merged-but-unreleased work belong in CHANGELOG.md under
[Unreleased] plus a "pending for the next cut" note, to be moved into the
abstract when the release is cut. Putting it in `abstract` now would credit
the published v0.5.0 artifact with a method it does not contain.

`MANIFEST_SCHEMA` is deliberately NOT bumped: the run manifest contract is
unchanged. No diagnostic is routed through `restrict_to_traceless` -- the only
callers in the tree are its tests -- so no reported number moves and no
manifest field is added or altered. The entry says so explicitly, because the
neighbouring pending notes are corrections that DO change numerical results
and a reader must not have to guess which kind this is. When a diagnostic is
later routed through the restriction, that change is the one to describe as a
correction.

Verified locally: .github/scripts/check_claim_safety.py passes (22 public-
facing markdown files), .github/scripts/check_workflow_hardening.py passes,
and the added CITATION.cff lines are comments only, so the CFF mapping itself
is byte-unchanged.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01VG6CeFJVgem6n4Tsam1onb

* chore(#139): record the reproduction of review finding P2 before repairing it

Baseline 33e2602: 1201 passed, exit 0,
coverage 94.35%.

Reproduced the P2 finding before touching the code. For `d = 16`, first column
carrying `+2.5e307` in the first eight trace rows and `-2.5e307` in the
remaining eight:

  vec(I)^H @ L, entry 0                  = (inf+nanj)   [overflow + invalid]
  trace_preservation_defect              = nan, scale = 1e308
  trace_preservation_componentwise_error = 0.0
  restrict_to_traceless                  -> ValueError, "not representable as
                                            finite float64"

This is the mirror image of the finding repaired in 5c151f9. That one was
fail-OPEN on an invalid generator; this is fail-CLOSED on a valid one, and the
scale-safe componentwise gate that reads the operator correctly never runs,
because the normwise gate raises first.

Checked against exact rational arithmetic rather than against another float
computation, following the reasoning PR #136 used for the Frobenius norm: the
column sum is EXACTLY zero as a Fraction, and the Frobenius norm is exactly
`4 * 2.5e307 = 1e308`, inside float64 range. The generator is legal and fully
representable; only an intermediate product overflows.

Worth recording because it is the same class one level up: the first oracle
written for this check was an ordinary sorted Python `sum`, and it overflowed
to `-inf` itself. An oracle computed the same way as the thing it judges is
not an oracle.

Root cause is the SHARED helper `linalg.py:331`, `vec_i.conj() @ L_super`, on
which both certified eigensolver paths also depend -- so the repair belongs
there and not as a second copy in traceless.py.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01VG6CeFJVgem6n4Tsam1onb

* fix(#139): sum the trace equations without overflowing on the way to zero

PR #139 review, P2. `trace_preservation_defect` assembled `vec(I)^H L` with an
ordinary matrix product. That product accumulates in ordinary units, so a
column whose trace equation cancels EXACTLY can still leave the float64 range
on the way there. Measured for d=16 with `+2.5e307` in eight trace rows of
column 0 and `-2.5e307` in the other eight: the product returned `inf+nanj`,
the defect came back `nan` beside a perfectly representable scale of 1e308,
and `restrict_to_traceless` refused a legal generator. Checked against exact
rational arithmetic, not against another float computation: the column sum is
exactly zero and the Frobenius norm is exactly `4 * 2.5e307 = 1e308`.

This is the mirror image of the finding repaired in 5c151f9 -- that one was
fail-open on an invalid generator, this is fail-closed on a valid one -- and
the scale-safe componentwise gate that reads this operator correctly (0.0)
never got to run, because the normwise gate raised first.

Repaired at the SHARED helper rather than in `traceless.py`, because both
certified eigensolver paths depend on the same assembly. New primitive
`norms.scaled_column_sums` follows the contract the rest of that module
already uses: power-of-two scaling, accumulate, shift back.

Three decisions in it were measured rather than assumed:

- The scale is PER COLUMN. One global scale also stops the overflow, but a
  matrix holding a column near 1e300 beside one near 1e-300 then flushes the
  small column to zero -- turning a real trace defect into an apparent exact
  zero. Verified: the global variant returns 0 for a column whose sum is
  1e-300. That is the same failure in the opposite, silent direction.
- Accumulation is `math.fsum`, matching `scaled_cancellation_ratio`, which
  reads the same trace equations. A vectorised variant benchmarked 2.4x to
  11.9x faster and returned 0 for the column `[1, 1e16, -1e16]`, whose exact
  sum is 1: ordinary summation does not round such a defect away, it deletes
  it, and a deleted defect reads as an exactly trace-preserving generator.
  Measured cost of keeping fsum: 0.07 ms at d=4 against a 0.26 ms eigensolve,
  falling to 0.30% of the eigensolve at d=16.
- Scaling is a precondition, not a refinement: unscaled `math.fsum` on those
  columns raises `OverflowError` rather than returning a sentinel.

CHANGES A REPORTED NUMBER, deliberately. `trace_defect` is now the correctly
rounded exact sum of the represented coefficients, where the matrix product
returned whatever its summation order happened to produce. The certificate
verdicts do not move at any sane tolerance -- the affected readings are at
rounding level -- but the number itself appears in persisted certificate audit
metadata. `test_stiff_generator_is_exactly_trace_preserving` asserted
`defect == 0.0` for the stiff fixture and is renamed to
`..._is_trace_preserving_to_rounding`: adjudicated with `fractions.Fraction`,
the exact sum of that generator's represented coefficients is 1.158640e-11 in
column 10, i.e. 2.29e-17 relative. The old equality was a property of the
summation route, not of the generator, and the test's actual claim -- that the
lost zero mode is the solver's doing and not a malformed input -- is unchanged.

A 1-D argument used to reach the product and return a number that described no
superoperator; it is now reported through the same unevaluable channel the
non-square case already uses.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01VG6CeFJVgem6n4Tsam1onb

* test(#139): make the P2 anchor fail as a verdict, not as an escaping error

The discrimination run classified this node ROT DURCH ABSTURZ rather than as
evidence: under the mutation, `restrict_to_traceless` raised and the test
died of ValueError instead of an assertion. A test that dies of an unexpected
error type proves that something broke, not what.

The reading is now asserted at `trace_preservation_defect` first -- under the
mutation that is `nan`, so the failure is a statement about a number -- and
the integration call is wrapped so a refusal reports as `pytest.fail`.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01VG6CeFJVgem6n4Tsam1onb

* fix(#139): give the real and imaginary parts their own scale

Found by turning the P2 finding against its own repair, before pushing it.

`scaled_column_sums` scaled per column, which closes the underflow between a
huge column and a tiny one. Real and imaginary parts of ONE column still
shared a scale, inherited from `_finite_component_scale` -- and that is
correct for a norm, where nothing cancels, but wrong for a sum. Measured on
the column [1e300+1e-300j, -1e300, 0], whose sum is 1e-300j: the real parts
annihilate, the answer is the imaginary part 600 orders of magnitude below the
scale that was chosen, and the function returned exactly 0j.

That is the same failure class as the one being repaired, one level down, and
in the fail-OPEN direction: a real trace defect reported as an exact
cancellation. Each component now carries its own power-of-two scale.

Cost after the change, measured: 0.095 ms at d=4 against a 0.339 ms
eigensolve, 0.20% of the eigensolve at d=16.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01VG6CeFJVgem6n4Tsam1onb

* docs(#139): record the overflow repair and the number it changes

AGENTS.md:121-123 requires the citation record to follow methodology and
results changes. This one changes both a method and a reported number, so it
is written as a correction rather than an addition, in the same shape the
neighbouring pending notes use.

CITATION.cff gets a 'pending for the next cut' comment -- the abstract still
describes released v0.5.0 only, per the convention that file states for itself
-- and CHANGELOG.md [Unreleased] carries the detail: what overflowed, what the
exact rational adjudication says, that `trace_defect` now reads 1.158640e-11
instead of 0.0 on the stiff fixture, and that certificate verdicts do not move
because the affected readings are at rounding level.

MANIFEST_SCHEMA is deliberately not bumped, and the reason is stated rather
than assumed: `trace_defect` reaches persisted reports through
`zero_mode_certificate` (_types.py:70) but is not a run-manifest field, so the
manifest contract is untouched and `input_hash` does not move.

Verified: check_claim_safety.py passes (22 public-facing markdown files),
check_workflow_hardening.py passes, and the added CITATION.cff lines are
comments only, so the CFF mapping is byte-unchanged.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01VG6CeFJVgem6n4Tsam1onb

* fix(types): out explizit annotieren -- mypy-Gate fiel nur unter Python 3.10

Die CI-Matrix meldete auf `test (ubuntu-latest, 3.10)`:

  norms.py:193: error: Incompatible types in assignment (expression has type
  "ndarray[tuple[int, ...], dtype[Any]]", variable has type
  "ndarray[tuple[int], dtype[Any]]")

Ursache: `out = np.zeros(cols, dtype=complex)` wird als EINDIMENSIONAL
inferiert (`tuple[int]`); die spaetere Neuzuweisung
`out = accumulated[0] + 1j * accumulated[1]` liefert den allgemeineren
Shape-Typ `tuple[int, ...]`. Die explizite Annotation `np.ndarray` deckt
beide Formen und entspricht der Rueckgabe-Signatur der Funktion.

Warum das lokal nicht auffiel: die mypy-Konfiguration zielt fest auf
`python_version = "3.12"` (pyproject.toml:178, damit mypy die NumPy->=2.5-Stubs
mit PEP-695-`type`-Statement parsen kann). Unter 3.10 installiert die Matrix
eine aeltere NumPy, deren Stubs den Shape-Typ enger aufloesen. Der Unterschied
liegt also in der STUB-Version, nicht in der Sprachversion -- und ist ohne ein
3.10-venv lokal nicht sichtbar.

Belegt lokal (Python 3.12): mypy exit 0 ueber 55 Dateien; gezielte Tests
(-k "norm or traceless") 94 passed, exit 0. Der 3.10-Beleg steht aus und
kommt aus der CI dieses Pushes -- Stand hier ehrlich "lokal gruen, 3.10 offen".

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01VG6CeFJVgem6n4Tsam1onb

* chore(#139): reproduce the third finding and judge the approach before patching

Baseline e48c879.

REPRODUCED. `scaled_column_sums([1e300, -1e300, 1e-300])` returns 0, where the
exact sum and unscaled `math.fsum` both return 1e-300. It propagates:
`trace_preservation_defect` reports 0.0 for a column whose exact defect is
1e-300, and `restrict_to_traceless(..., tp_rtol=0)` therefore ACCEPTS a
generator that is not exactly trace preserving. Fail-open.

One of my own suspicions was wrong and is recorded as such: I expected
`scaled_cancellation_ratio` (#130) to share the defect, since it shares the
scaling. It does not. The exact ratio there is 1e-300/2e300 = 5e-601, which is
itself zero in float64, so its 0.0 is the correct answer. Only the normwise
path is affected, exactly as the review says.

JUDGEMENT, because this is the third level of one class in one function and
the pattern is the finding: ONE power-of-two scale per column and component
cannot carry this, structurally. A scale SHIFTS the window 2^-1074..2^1024, it
does not widen it. Every addend below `max * 2^-1074` becomes zero before fsum
sees it. Measured: a gap of 2^-1000 survives, 2^-1074 is lost. Chosen too
large the remainder dies, too small it overflows -- once a column spans the
range AND the dominant terms cancel, no choice is right. A fourth patch on
this construction would earn a fifth.

REPLACEMENT, and it is simpler than what it replaces: `math.fsum` is already
exact over the full range; its only failure is intermediate overflow, which is
order-dependent (measured: eight +2.5e307 then eight -2.5e307 raises, the same
values alternating do not). So try fsum first -- the normal path, exact, with
no scaling and hence no underflow -- and fall back to two exponent bands split
at 2^512 only when it raises. Above the split, scaling by 2^-1024 lands every
term in [2^-512, 1): no underflow, and k terms cannot overflow. Below it, k
terms sum below k*2^512 < 2^1024 with no scaling at all. Two bands suffice
provably, and no addend is ever sacrificed.

Validated before building, against exact rational arithmetic: ten edge cases
exact, including eight +2.5e307, eight -2.5e307 and a 1e-320 remainder, which
neither the current code nor bare fsum can do; and 4000 random cases with
exponents from -1060 to 1020 and forced cancellation, with zero deviations.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01VG6CeFJVgem6n4Tsam1onb

* fix(#139): sum exactly instead of choosing a scale that cannot exist

PR #139, third review round on this function, and the third level of one class
in it: overflow between columns, then a shared scale for real and imaginary
parts, now underflow within a component. The pattern is the finding, so the
construction is replaced rather than patched again.

REPRODUCED FIRST. `scaled_column_sums([1e300, -1e300, 1e-300])` returned 0
where the exact sum and unscaled `math.fsum` both return 1e-300;
`trace_preservation_defect` therefore reported 0.0, and
`restrict_to_traceless(..., tp_rtol=0)` -- the strictest request a caller can
make -- accepted a generator that is not trace preserving.

WHY THE APPROACH COULD NOT BE PATCHED. A power-of-two scale SHIFTS the
representable window 2**-1074 .. 2**1024; it does not widen it. Every addend
below `max * 2**-1074` is flushed to zero before the accumulation starts, and
choosing the scale smaller reintroduces the overflow the scaling exists to
prevent. Measured: a gap of 2**-1000 survives, 2**-1074 is lost. Once a column
spans the range and its dominant terms cancel, no scale is correct.

WHAT REPLACES IT. `math.fsum` is already exact over the whole range -- its one
failure is an intermediate that leaves it, and that is order-dependent (eight
+2.5e307 then eight -2.5e307 raises; the same values alternating do not). So
`_overflow_safe_fsum` calls it plainly first, with no scaling and hence no
underflow, and falls back to two exponent bands split at 2**512 only when it
raises. Above the split a shift by 2**-1024 lands every term in [2**-512, 1),
far above the subnormal floor, and k such terms sum below k. Below the split k
terms sum below k * 2**512 < 2**1024 and need no scaling at all. Neither band
can overflow or underflow, so two bands suffice provably and no addend is ever
sacrificed. Overflow in the recombination is the honest answer.

Validated against exact rational arithmetic, which cannot share the failure
mode the way a float oracle can: ten edge cases exact -- including eight
+2.5e307, eight -2.5e307 and a 1e-320 remainder in ONE column, which neither
the old construction nor bare fsum can compute -- and 500 randomised columns in
the suite plus 4000 in the pre-build probe, exponents -1060 to 1020 with forced
cancellation, zero deviations.

Cheaper as well as more correct, measured with warm-up: 0.0245 ms at d=4
against a 0.179 ms eigensolve (13.6%, down from 28%), 1.43% at d=8, 0.25% at
d=16. The band fallback costs 2.7x the normal path and is entered only on
overflow.

Also removes the seam that made the 3.10 type gate fail in e48c879: `out` is
no longer reassigned to a differently shaped array, only indexed, so its shape
type cannot drift between NumPy stub versions.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01VG6CeFJVgem6n4Tsam1onb

* test(#139): pin the known limit of the band fallback, found by self-probe

Turned the finding against its own repair, as in the previous round, and this
time the first probe was itself the problem: it reported zero deviations while
entering the band fallback once in 6000 cases. A probe that does not reach the
path it judges measures nothing. With a guaranteed-overflow kernel the picture
is different and honest:

  plain path : 20000 randomised columns, 0 deviations (exactly rounded)
  band path  : 6713 of 20000 columns, 3 deviations, each one ulp in the last
               place of a result of order 1e182 or larger

So there IS a fourth level -- but it is not the same class. The class this
function fights deletes information and reads as an exact zero; double rounding
in the recombination moves the last bit. The distinction was measured, not
argued: 19965 further columns constructed to cancel INSIDE the band path
produced no spurious zero and no sign flip, so no defect is deleted and none
changes direction.

Not repaired here, deliberately: removing the last ulp needs the exact
remainder `math.fsum` discards, i.e. a Shewchuk accumulator returning its
partial sums, which is a different construction and not what this review asked
for. The limit is documented at the function and pinned by a test that asserts
one-ulp accuracy, no fabricated zero and no sign flip, and that fails if the
fixture ever stops reaching the fallback.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01VG6CeFJVgem6n4Tsam1onb

* docs(#139): correct the record to the construction that now exists

The CHANGELOG and CITATION entries from the previous round described a
per-column, per-component scaling that the third review round replaced. A
record that describes a construction the code no longer has is a defect of its
own, so both are corrected to what is there: exact fsum per column and
component, with an exponent-banded fallback entered only on overflow.

The known one-ulp limit of that fallback is recorded with it, measured rather
than asserted, so the entry does not claim more than the code delivers.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01VG6CeFJVgem6n4Tsam1onb

* fix(#139): abstain instead of guessing when the bands cannot be recombined

PR #139, fifth review round. Two exponent bands rounded independently cannot
carry a residual that survives only across their boundary. Reproduced exactly
as reported, for `[7e307]*3 + [-7e307]*3 + [2**513] + [-2**511]*4 + [1e-300]`
whose exact sum is 1e-300: the bands are `+2**513` and `-2**513 + 1e-300`, the
low band rounds correctly to exactly `-2**513` because the remainder is far
below its last place, and the two then annihilate to 0.0. Used as the twelve
trace-row coefficients of a finite 144x144 operator, both trace-preservation
readings read zero and `restrict_to_traceless(..., tp_rtol=0)` accepted an
operator that is not trace preserving.

This is NOT a fifth approximation. The exact answer needs an accumulator that
keeps Shewchuk partial sums across the band boundary until one final rounding,
which is a different construction and does not belong in this pull request.
What belongs here is refusing to guess: whether a band sum was rounded is
cheap to establish -- re-sum the band against its own negated total and see
what is left -- and the recombination is checked by its fast-two-sum residual.
When all three are exact the result is the exact sum; otherwise the function
returns `nan` and the existing fail-closed channel refuses the operator.

A result whose accuracy the function cannot vouch for is reported as UNKNOWN,
not as ZERO. That distinction is the whole subject of this pull request.

The contract this pins is stronger than the one-ulp bound it replaces: every
number returned is now the exact sum, and everything else is an abstention.
The test is rewritten accordingly and asserts both that the fallback is
reached and that it does not abstain on everything -- a check that abstains
always would pass every negative test and measure nothing.

Measured, because an abstention that costs legitimate input is not acceptable:
across the trace-related suite the fallback is entered 10 times in 10280 calls
and abstains once. Physical stiff generators never reach it at all -- their
trace equations do not overflow `math.fsum` for fast rates up to 1e12 -- which
a positive control now asserts directly rather than assuming.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01VG6CeFJVgem6n4Tsam1onb

* fix(#139): drop a guard that abstained without cause

The discrimination run classified the recombination guard BLIND, and chasing
that down showed the guard was the problem rather than the test. In all 6 of
109 fallback columns where it alone decided the outcome, the answer it refused
would have been CORRECT.

The reason is straightforward once stated: at that point both band sums are
already established as exact, so a single IEEE-754 addition of them returns
the correctly rounded exact total by definition. There is nothing left to
verify, and refusing anyway costs real answers.

Abstaining without cause is its own defect. A function that refuses everything
passes every negative test and is still broken, which is exactly the trap the
positive control in this suite exists to catch.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01VG6CeFJVgem6n4Tsam1onb

* docs(#139): record the abstention instead of the one-ulp bound it replaced

The previous entries promised a fallback accurate to one ulp. After the fifth
review round that is no longer what the code does, and the new contract is the
stronger one: every value the fallback returns is the exact sum, and anything
it cannot vouch for is reported as unmeasurable so the gates refuse it.

Recorded with the numbers that bound the cost: the fallback is entered 10
times in 10280 calls across the trace-related suite and abstains once, and
physical stiff generators never reach it.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01VG6CeFJVgem6n4Tsam1onb

---------

Co-authored-by: Claude Opus 5 (1M context) <noreply@anthropic.com>
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

fix(numerics): make trace-preservation norms underflow/overflow-safe

1 participant