From 021c28d6c0e517653a2324981bf6b45da97b4ff9 Mon Sep 17 00:00:00 2001 From: Marco Hostettler Date: Tue, 1 Sep 2026 20:21:31 +0200 Subject: [PATCH 01/22] feat(numerics): add scaled sum-of-squares norm primitive --- src/liouscope/numerics/norms.py | 64 +++++++++++++++++++++++++++++++++ 1 file changed, 64 insertions(+) create mode 100644 src/liouscope/numerics/norms.py diff --git a/src/liouscope/numerics/norms.py b/src/liouscope/numerics/norms.py new file mode 100644 index 00000000..d16b0dde --- /dev/null +++ b/src/liouscope/numerics/norms.py @@ -0,0 +1,64 @@ +"""Overflow/underflow-safe Euclidean norm primitives. + +The implementation follows the same numerical contract as LAPACK's xLASSQ +family: scale values before squaring, accumulate in the scaled domain, then +restore the scale. LiouScope uses a power-of-two scale so the rescaling itself +does not introduce a reciprocal overflow for subnormal complex values and does +not round the mantissas merely to choose numerical units. +""" + +from __future__ import annotations + +import numpy as np + + +def scaled_euclidean_norm(values: np.ndarray) -> float: + """Return ``sqrt(sum(abs(values)**2))`` without spurious under/overflow. + + This is the vector 2-norm and, when ``values`` is a matrix, its Frobenius + norm. For finite input it returns a finite number whenever the true norm is + representable as ``float64``. It returns ``inf`` only when an infinite + input component or the mathematical norm itself is not representable, and + propagates ``NaN`` rather than fabricating evidence from corrupted input. + + Complex values are treated as their real and imaginary components, so the + calculation is mathematically ``sqrt(sum(re**2 + im**2))``. We avoid + computing ``abs(z)`` before selecting the scale because a finite complex + value near the top of the floating-point range can overflow in the modulus. + + The power-of-two scaling is an xLASSQ-style specialization: if ``m`` is the + largest finite component magnitude and ``m = f * 2**e`` with + ``0.5 <= f < 1``, shifting every component by ``2**(-e)`` places the + largest component below one. Squaring cannot overflow, at least one square + is order one, and restoring with ``ldexp`` overflows iff the true norm does. + """ + arr = np.asarray(values) + if arr.size == 0: + return 0.0 + + real = np.asarray(np.real(arr), dtype=float) + imag = np.asarray(np.imag(arr), dtype=float) + if np.any(np.isnan(real)) or np.any(np.isnan(imag)): + return float("nan") + if np.any(np.isinf(real)) or np.any(np.isinf(imag)): + return float("inf") + + component_max = float( + max( + float(np.max(np.abs(real))), + float(np.max(np.abs(imag))), + ) + ) + if component_max == 0.0: + return 0.0 + + exponent = int(np.frexp(component_max)[1]) + scaled_real = np.ldexp(real, -exponent) + scaled_imag = np.ldexp(imag, -exponent) + sumsq = float( + np.sum(scaled_real * scaled_real, dtype=float) + + np.sum(scaled_imag * scaled_imag, dtype=float) + ) + scaled_norm = float(np.sqrt(sumsq)) + with np.errstate(over="ignore", under="ignore"): + return float(np.ldexp(scaled_norm, exponent)) From 87027f3dd4287d0d43bf747f7aa0f6f3bec387ed Mon Sep 17 00:00:00 2001 From: Marco Hostettler Date: Tue, 1 Sep 2026 20:22:10 +0200 Subject: [PATCH 02/22] test(#130): pin xLASSQ-style norm extremes and equivariance --- tests/test_scaled_norms.py | 59 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 59 insertions(+) create mode 100644 tests/test_scaled_norms.py diff --git a/tests/test_scaled_norms.py b/tests/test_scaled_norms.py new file mode 100644 index 00000000..d6937f57 --- /dev/null +++ b/tests/test_scaled_norms.py @@ -0,0 +1,59 @@ +"""Numerical-contract tests for the scaled sum-of-squares primitive (#130).""" + +from __future__ import annotations + +import numpy as np +import pytest + +from liouscope.numerics.norms import scaled_euclidean_norm + + +def test_tiny_finite_values_do_not_underflow_to_zero() -> None: + values = np.array([1.0e-200, -2.0e-200, 3.0e-200]) + got = scaled_euclidean_norm(values) + + assert got > 0.0 + assert got == pytest.approx(np.sqrt(14.0) * 1.0e-200, rel=2.0e-15) + + +def test_subnormal_values_keep_a_nonzero_norm() -> None: + values = np.array([1.0e-320, -2.0e-320, 3.0e-320]) + got = scaled_euclidean_norm(values) + + assert got > 0.0 + assert got == pytest.approx(3.7416e-320, rel=2.0e-4) + + +def test_large_finite_values_stay_finite_when_true_norm_is_representable() -> None: + values = np.array([1.0e308, -1.0e308, 1.0e308]) + got = scaled_euclidean_norm(values) + + assert np.isfinite(got) + assert got == pytest.approx(np.sqrt(3.0) * 1.0e308, rel=2.0e-15) + + +def test_nonrepresentable_true_norm_fails_as_infinity_not_a_false_finite_value() -> None: + values = np.array([1.0e308, -1.0e308, 1.0e308, -1.0e308]) + assert scaled_euclidean_norm(values) == float("inf") + + +def test_complex_components_follow_the_same_scaled_contract() -> None: + got = scaled_euclidean_norm(np.array([1.0e308 + 1.0e308j])) + + assert np.isfinite(got) + assert got == pytest.approx(np.sqrt(2.0) * 1.0e308, rel=2.0e-15) + + +@pytest.mark.parametrize("scale", [1.0e-250, 1.0e-100, 1.0, 1.0e100, 1.0e250]) +def test_norm_is_scale_equivariant_over_representable_ranges(scale: float) -> None: + base = np.array([1.0, -2.0, 3.0, -4.0]) * 1.0e-10 + reference = scaled_euclidean_norm(base) + got = scaled_euclidean_norm(scale * base) + + assert got == pytest.approx(abs(scale) * reference, rel=3.0e-15) + + +def test_zero_nan_and_inf_semantics_are_explicit() -> None: + assert scaled_euclidean_norm(np.zeros(4)) == 0.0 + assert np.isnan(scaled_euclidean_norm(np.array([1.0, np.nan]))) + assert scaled_euclidean_norm(np.array([1.0, np.inf])) == float("inf") From 1f94a4151c072f881d17953f607e3f5bbbd78b15 Mon Sep 17 00:00:00 2001 From: Marco Hostettler Date: Tue, 1 Sep 2026 20:24:07 +0200 Subject: [PATCH 03/22] fix(#130): use scaled norms for trace-preservation evidence --- src/liouscope/numerics/linalg.py | 79 ++++++++------------------------ 1 file changed, 18 insertions(+), 61 deletions(-) diff --git a/src/liouscope/numerics/linalg.py b/src/liouscope/numerics/linalg.py index 36bac452..72aa1892 100644 --- a/src/liouscope/numerics/linalg.py +++ b/src/liouscope/numerics/linalg.py @@ -27,6 +27,7 @@ ZERO_MODE_APOSTERIORI_MARGIN, ZERO_MODE_EPS_FACTOR, ) +from .norms import scaled_euclidean_norm @dataclass(frozen=True, slots=True) @@ -252,23 +253,27 @@ def operator_zero_tolerance( def trace_preservation_defect(L_super: np.ndarray) -> tuple[float, float]: - """Return ``(||vec(I)^H L||_2, ||L||_F)`` for a vectorised generator. + """Return robust ``(||vec(I)^H L||_2, ||L||_F)`` for a generator. Every GKSL generator is trace preserving, which in column-stacking convention is the exact statement ``vec(I)^H L = 0``. The quotient of the two returned numbers is therefore a dimensionless, unit-invariant measure of how far ``L`` is from being a legal generator. - Both norms square their entries before summing, so a generator expressed - in small enough -- but still perfectly normal -- rate units loses its - reference scale to underflow: for ``diag([0, -1e-200, -2e-200, -3e-200])`` - every square falls below ``5e-324`` and the pair comes back ``(0.0, 0.0)``. - The applicability gate then reads ``0 > tp_rtol * tiny``, which is False, - so a demonstrably non-trace-preserving operator becomes certificate- - applicable -- while the SAME operator at ``1e-150`` is correctly refused. - A pure change of rate unit must not decide whether an operator is a legal - generator, so the norms are recomputed on a scaled copy when the reference - scale has been lost (round-23 review, PR #121). + Issue #130: both quantities are evidence, so their arithmetic must be + correct independently of the policy that consumes them. Ordinary + ``np.linalg.norm`` squares before summing and can therefore return the + false sentinels ``0.0`` (underflow) or ``inf`` (overflow) even when the + mathematical norm of finite input is representable. Both norms now use + :func:`scaled_euclidean_norm`, an xLASSQ-style power-of-two scaled + sum-of-squares primitive. It returns a finite norm whenever the true + float64 result is representable and returns infinity only when the input or + the true norm is genuinely non-representable. + + This deliberately separates measurement from certificate policy. A very + large but finite operator must not be rejected merely because a naive norm + overflowed; structural guards such as ``band_discriminates`` decide whether + the resulting zero-mode evidence actually discriminates. """ L_super = np.asarray(L_super) n = L_super.shape[0] @@ -276,56 +281,8 @@ def trace_preservation_defect(L_super: np.ndarray) -> tuple[float, float]: if d * d != n: return float("nan"), float("nan") vec_i = np.eye(d, dtype=complex).reshape(-1, order="F") - defect = float(np.linalg.norm(vec_i.conj() @ L_super)) - fro = float(np.linalg.norm(L_super, ord="fro")) - # The rescue is deliberately ONE-DIRECTIONAL: it repairs underflow and - # leaves overflow alone. ``fro = inf`` is not an accident to be worked - # around but the input to the round-21 refusal of a non-finite reference - # scale -- the guard that closed the round-20 counterexample (cancelling - # +-1e308 entries, spectrum {1,2,3,4}, certified as entirely stationary). - # Computing that Frobenius norm by scaling would make it finite again at - # ~1.4e308, readmit the operator, and silently reopen a hole whose repair - # CI had to prove across five interpreter versions. Underflow and overflow - # look symmetric here and are not: only one of them has a documented - # decision behind it. - if fro == 0.0 and L_super.size: - # ``max(|Re|, |Im|)`` rather than ``max|.|``: the modulus of a finite - # complex entry can itself overflow, and this quantity cannot. It is - # within sqrt(2) of the true maximum magnitude, which is all that is - # needed -- it only selects the exponent to shift by. - m = float(np.max(np.maximum(np.abs(np.real(L_super)), np.abs(np.imag(L_super))))) - if m > 0.0: - # Scaled by a POWER OF TWO, not by ``m`` itself. Two reasons, the - # second of which cost a regression before it was measured: - # - # 1. ``ldexp`` shifts the exponent and leaves every mantissa bit - # alone, so the round trip is EXACT -- dividing by ``m`` and - # multiplying back introduces two roundings into a number whose - # whole purpose is to be compared against a tolerance. - # 2. ``L / m`` is COMPLEX division, and for a subnormal divisor - # NumPy's algorithm forms an intermediate reciprocal that - # overflows: measured for ``diag([0, -1e-310, -2e-310, - # -3e-310])``, ``L / m`` came back ``[nan, -inf, -inf, -inf]`` - # and both norms became NaN. The round-21 guard then refused the - # operator -- fail-closed, but on a value nothing had measured, - # and a healthy GKSL generator at that scale was refused with it - # (``applicable`` went True -> False for valid input). Repairing - # a fail-open must not install a fail-closed defect in its place. - # - # ``frexp`` returns ``m = mantissa * 2**exp`` with the mantissa in - # [0.5, 1), so shifting by ``-exp`` puts every entry in a range - # where neither the squares below nor the reconstruction can leave - # the normal range. A genuinely zero operator has ``m == 0`` and - # keeps ``(0.0, 0.0)``, which is the correct answer, not a lost one. - exp = int(np.frexp(m)[1]) - scaled = np.ldexp(np.real(L_super), -exp) + 1j * np.ldexp( - np.imag(L_super), -exp - ) - defect = float( - np.ldexp(float(np.linalg.norm(vec_i.conj() @ scaled)), exp) - ) - fro = float(np.ldexp(float(np.linalg.norm(scaled, ord="fro")), exp)) - return defect, fro + defect_vector = vec_i.conj() @ L_super + return scaled_euclidean_norm(defect_vector), scaled_euclidean_norm(L_super) def band_discriminates( From c2c31a4f2286a1c7e2ed79230a53643d4dbcc75f Mon Sep 17 00:00:00 2001 From: Marco Hostettler Date: Tue, 1 Sep 2026 20:24:55 +0200 Subject: [PATCH 04/22] test(#130): separate large-norm measurement from certificate policy --- tests/test_pr121_review_round23.py | 46 ++++++++++++++---------------- 1 file changed, 22 insertions(+), 24 deletions(-) diff --git a/tests/test_pr121_review_round23.py b/tests/test_pr121_review_round23.py index 4f1c67e5..0e8817ee 100644 --- a/tests/test_pr121_review_round23.py +++ b/tests/test_pr121_review_round23.py @@ -20,10 +20,10 @@ sit below the eigensolver's backward error. Every guard here is paired with a positive control: a guard that refuses -everything discriminates nothing and must not be able to pass this file. The -overflow twin of finding 12 has its own test for the opposite reason -- it -pins a refusal that must NOT be repaired, because the round-21 guard depends -on it. +everything discriminates nothing and must not be able to pass this file. Issue +#130 extends finding 12 symmetrically: a mathematically representable large norm +must be measured as finite, while the structural zero-mode certificate remains +responsible for refusing evidence that does not discriminate. """ from __future__ import annotations @@ -57,7 +57,7 @@ def _non_tp_operator() -> np.ndarray: # --------------------------------------------------------------------------- -# Finding 12: use underflow-safe norms for the TP gate. +# Finding 12: use underflow/overflow-safe norms for the TP gate. # --------------------------------------------------------------------------- @@ -110,35 +110,33 @@ def test_a_healthy_generator_in_the_same_rate_units_still_certifies() -> None: assert cert.certified, f"a GKSL generator at scale {scale:g} was not certified" -def test_the_overflow_direction_is_deliberately_not_rescued() -> None: - """The round-21 refusal must SURVIVE the round-23 repair. +def test_representable_large_norm_is_measured_then_rejected_structurally() -> None: + """Measurement arithmetic must not carry the certificate's policy. - The two directions look symmetric and are not. ``||L||_F = inf`` for the - round-20 counterexample (cancelling +-1e308 entries, spectrum {1,2,3,4}) is - the input to the round-21 guard, which refuses a non-finite reference - scale. Its true Frobenius norm is about 1.4e308 and therefore perfectly - representable, so a scaled computation would hand back a finite number, - readmit the operator and reopen a hole that took five interpreter versions - of CI to close. This test goes red the moment someone makes the rescue - symmetric. + The round-20 counterexample has two cancelling ``+-1e308`` entries. Its + true Frobenius norm is about ``sqrt(2)*1e308`` and is therefore finite in + float64, even though a naive sum of squares overflows. Issue #130 requires + that finite measurement to stay finite. Safety is preserved one layer + later: the zero-mode band contains the whole measured spectrum, so + ``band_discriminates`` leaves the certificate unresolved instead of + treating arithmetic overflow as a surrogate refusal signal. """ - # The round-20 operator, entry for entry (see - # tests/test_pr121_review_round20.py::_overflowing_norm_operator). op = np.diag([1.0, 2.0, 3.0, 4.0]).astype(complex) op[0, 1] = 1.0e308 op[3, 1] = -1.0e308 assert np.all(np.isfinite(op)), "the input must be FINITE for this to be the case" - with np.errstate(over="ignore"): - defect, fro = trace_preservation_defect(op) + + defect, fro = trace_preservation_defect(op) assert defect == pytest.approx(np.sqrt(17.0)), "a real, order-one trace defect" - assert not np.isfinite(fro), ( - "the overflowing reference scale was made finite; the round-21 guard " - "that refuses a non-finite scale now has nothing to refuse" - ) + assert np.isfinite(fro), "a representable Frobenius norm must not overflow" + assert fro == pytest.approx(np.sqrt(2.0) * 1.0e308, rel=2.0e-15) + with np.errstate(over="ignore"), warnings.catch_warnings(): warnings.simplefilter("ignore") cert = certified_eigvals(op)[-1] - assert not cert.applicable and not cert.certified + assert cert.applicable + assert not cert.certified + assert not cert.resolved def test_the_repair_reaches_into_the_subnormal_range() -> None: From c5492eac307b1f6b911c3a3bf445d91e9855f401 Mon Sep 17 00:00:00 2001 From: Marco Hostettler Date: Tue, 1 Sep 2026 20:25:20 +0200 Subject: [PATCH 05/22] test(#130): make tiny-norm assertions discriminating --- tests/test_scaled_norms.py | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/tests/test_scaled_norms.py b/tests/test_scaled_norms.py index d6937f57..fa499464 100644 --- a/tests/test_scaled_norms.py +++ b/tests/test_scaled_norms.py @@ -13,7 +13,10 @@ def test_tiny_finite_values_do_not_underflow_to_zero() -> None: got = scaled_euclidean_norm(values) assert got > 0.0 - assert got == pytest.approx(np.sqrt(14.0) * 1.0e-200, rel=2.0e-15) + # Compare after dividing by the known input scale: pytest.approx's default + # absolute tolerance would otherwise dwarf a 1e-200 expected value and let + # a wildly inaccurate non-zero result pass. + assert got / 1.0e-200 == pytest.approx(np.sqrt(14.0), rel=2.0e-15) def test_subnormal_values_keep_a_nonzero_norm() -> None: @@ -21,7 +24,9 @@ def test_subnormal_values_keep_a_nonzero_norm() -> None: got = scaled_euclidean_norm(values) assert got > 0.0 - assert got == pytest.approx(3.7416e-320, rel=2.0e-4) + # Subnormal quantisation limits the attainable relative accuracy here, so + # compare the dimensionless ratio rather than using an absolute tolerance. + assert got / 1.0e-320 == pytest.approx(np.sqrt(14.0), rel=2.0e-4) def test_large_finite_values_stay_finite_when_true_norm_is_representable() -> None: From f8dca7ad9d358885ffc659c1a155ecf7393d0c72 Mon Sep 17 00:00:00 2001 From: Marco Hostettler Date: Tue, 1 Sep 2026 20:27:46 +0200 Subject: [PATCH 06/22] test(#130): make scale-equivariance assertion dimensionless --- tests/test_scaled_norms.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/tests/test_scaled_norms.py b/tests/test_scaled_norms.py index fa499464..90eba316 100644 --- a/tests/test_scaled_norms.py +++ b/tests/test_scaled_norms.py @@ -54,8 +54,12 @@ def test_norm_is_scale_equivariant_over_representable_ranges(scale: float) -> No base = np.array([1.0, -2.0, 3.0, -4.0]) * 1.0e-10 reference = scaled_euclidean_norm(base) got = scaled_euclidean_norm(scale * base) + expected = abs(scale) * reference - assert got == pytest.approx(abs(scale) * reference, rel=3.0e-15) + assert got > 0.0 + # Compare the ratio to one so tiny expected values cannot pass through + # pytest.approx's default absolute tolerance. + assert got / expected == pytest.approx(1.0, rel=3.0e-15) def test_zero_nan_and_inf_semantics_are_explicit() -> None: From 45ba26048e485fe6be963c44c9b0f9f89b8eeee0 Mon Sep 17 00:00:00 2001 From: Marco Hostettler Date: Tue, 1 Sep 2026 20:38:19 +0200 Subject: [PATCH 07/22] feat(#130): add scale-safe componentwise cancellation metric --- src/liouscope/numerics/norms.py | 100 +++++++++++++++++++++++--------- 1 file changed, 74 insertions(+), 26 deletions(-) diff --git a/src/liouscope/numerics/norms.py b/src/liouscope/numerics/norms.py index d16b0dde..067d8012 100644 --- a/src/liouscope/numerics/norms.py +++ b/src/liouscope/numerics/norms.py @@ -1,36 +1,53 @@ -"""Overflow/underflow-safe Euclidean norm primitives. +"""Overflow/underflow-safe norm and cancellation primitives. -The implementation follows the same numerical contract as LAPACK's xLASSQ -family: scale values before squaring, accumulate in the scaled domain, then -restore the scale. LiouScope uses a power-of-two scale so the rescaling itself -does not introduce a reciprocal overflow for subnormal complex values and does -not round the mantissas merely to choose numerical units. +The Euclidean norm implementation follows the same numerical contract as +LAPACK's xLASSQ family: scale values before squaring, accumulate in the scaled +domain, then restore the scale. LiouScope uses a power-of-two scale so the +rescaling itself does not introduce a reciprocal overflow for subnormal complex +values and does not round mantissas merely to choose numerical units. """ from __future__ import annotations +import math + import numpy as np +def _finite_component_scale(values: np.ndarray) -> tuple[np.ndarray, np.ndarray, int] | None: + """Return real/imag parts shifted so their largest component is O(1). + + ``None`` denotes an exact all-zero array. NaN/inf are deliberately left to + public callers because their desired sentinel semantics differ by operation. + """ + arr = np.asarray(values) + real = np.asarray(np.real(arr), dtype=float) + imag = np.asarray(np.imag(arr), dtype=float) + component_max = float( + max( + float(np.max(np.abs(real))) if real.size else 0.0, + float(np.max(np.abs(imag))) if imag.size else 0.0, + ) + ) + if component_max == 0.0: + return None + exponent = int(np.frexp(component_max)[1]) + return np.ldexp(real, -exponent), np.ldexp(imag, -exponent), exponent + + def scaled_euclidean_norm(values: np.ndarray) -> float: """Return ``sqrt(sum(abs(values)**2))`` without spurious under/overflow. This is the vector 2-norm and, when ``values`` is a matrix, its Frobenius - norm. For finite input it returns a finite number whenever the true norm is - representable as ``float64``. It returns ``inf`` only when an infinite + norm. For finite input it returns a finite number whenever the true norm is + representable as ``float64``. It returns ``inf`` only when an infinite input component or the mathematical norm itself is not representable, and propagates ``NaN`` rather than fabricating evidence from corrupted input. Complex values are treated as their real and imaginary components, so the - calculation is mathematically ``sqrt(sum(re**2 + im**2))``. We avoid + calculation is mathematically ``sqrt(sum(re**2 + im**2))``. We avoid computing ``abs(z)`` before selecting the scale because a finite complex value near the top of the floating-point range can overflow in the modulus. - - The power-of-two scaling is an xLASSQ-style specialization: if ``m`` is the - largest finite component magnitude and ``m = f * 2**e`` with - ``0.5 <= f < 1``, shifting every component by ``2**(-e)`` places the - largest component below one. Squaring cannot overflow, at least one square - is order one, and restoring with ``ldexp`` overflows iff the true norm does. """ arr = np.asarray(values) if arr.size == 0: @@ -43,18 +60,10 @@ def scaled_euclidean_norm(values: np.ndarray) -> float: if np.any(np.isinf(real)) or np.any(np.isinf(imag)): return float("inf") - component_max = float( - max( - float(np.max(np.abs(real))), - float(np.max(np.abs(imag))), - ) - ) - if component_max == 0.0: + scaled = _finite_component_scale(arr) + if scaled is None: return 0.0 - - exponent = int(np.frexp(component_max)[1]) - scaled_real = np.ldexp(real, -exponent) - scaled_imag = np.ldexp(imag, -exponent) + scaled_real, scaled_imag, exponent = scaled sumsq = float( np.sum(scaled_real * scaled_real, dtype=float) + np.sum(scaled_imag * scaled_imag, dtype=float) @@ -62,3 +71,42 @@ def scaled_euclidean_norm(values: np.ndarray) -> float: scaled_norm = float(np.sqrt(sumsq)) with np.errstate(over="ignore", under="ignore"): return float(np.ldexp(scaled_norm, exponent)) + + +def scaled_cancellation_ratio(values: np.ndarray) -> float: + """Return ``abs(sum(values)) / sum(abs(values))`` scale-safely. + + This is a componentwise relative backward-error primitive for a scalar + linear equation assembled by cancellation. It is in ``[0, 1]`` for finite + input, is exactly scale invariant under representable non-zero rescaling, + and returns zero for an all-zero equation. + + Power-of-two scaling prevents both overflow and subnormal underflow before + summation. ``math.fsum`` is used separately on real and imaginary parts so + the numerator measures the represented equation rather than ordinary + left-to-right summation noise. The denominator is accumulated from stable + ``hypot`` magnitudes in the same scaled domain, so unrelated large entries + outside the equation cannot dilute its defect. + """ + arr = np.asarray(values) + if arr.size == 0: + return 0.0 + real = np.asarray(np.real(arr), dtype=float) + imag = np.asarray(np.imag(arr), dtype=float) + if np.any(~np.isfinite(real)) or np.any(~np.isfinite(imag)): + return float("nan") + + scaled = _finite_component_scale(arr) + if scaled is None: + return 0.0 + scaled_real, scaled_imag, _ = scaled + sum_real = math.fsum(float(x) for x in scaled_real.ravel()) + sum_imag = math.fsum(float(x) for x in scaled_imag.ravel()) + numerator = math.hypot(sum_real, sum_imag) + denominator = math.fsum( + math.hypot(float(re), float(im)) + for re, im in zip(scaled_real.ravel(), scaled_imag.ravel(), strict=True) + ) + if denominator == 0.0: + return 0.0 + return float(min(1.0, numerator / denominator)) From 423ea8efe199646f0b576e3a2c83690817d9ea26 Mon Sep 17 00:00:00 2001 From: Marco Hostettler Date: Tue, 1 Sep 2026 20:39:56 +0200 Subject: [PATCH 08/22] fix(#130): gate TP certificates on componentwise backward error --- src/liouscope/numerics/linalg.py | 794 +++++-------------------------- 1 file changed, 118 insertions(+), 676 deletions(-) diff --git a/src/liouscope/numerics/linalg.py b/src/liouscope/numerics/linalg.py index 72aa1892..b7580fa1 100644 --- a/src/liouscope/numerics/linalg.py +++ b/src/liouscope/numerics/linalg.py @@ -1,12 +1,10 @@ """General linear-algebra helpers used across diagnostics. -Includes: - -* :func:`eig_nonhermitian` --- wraps ``scipy.linalg.eig`` (LAPACK ``zgeev``), - the non-Hermitian eigensolver required by Liouvillians. Anchor D. -* :func:`is_hermitian`, :func:`is_density_matrix` --- input validation. -* :func:`support_check` --- enforces ``supp(rho_0) subset supp(rho_ss)`` - with ``eps = 1e-12`` regularisation (anchor J). +The module deliberately separates numerical measurements from scientific +certificate policy. In particular, trace-preservation norms are measured with +scaled arithmetic, while applicability of the structural zero-mode theorem is +also checked equation-by-equation so unrelated large matrix entries cannot hide +a local trace defect. """ from __future__ import annotations @@ -27,7 +25,7 @@ ZERO_MODE_APOSTERIORI_MARGIN, ZERO_MODE_EPS_FACTOR, ) -from .norms import scaled_euclidean_norm +from .norms import scaled_cancellation_ratio, scaled_euclidean_norm @dataclass(frozen=True, slots=True) @@ -50,21 +48,8 @@ def eig_nonhermitian( ) -> EigenDecomposition: """Return the eigendecomposition of a generic non-Hermitian matrix. - Always uses LAPACK ``zgeev`` (via ``scipy.linalg.eig``). Anchor D: never - fall back to ``eigh`` which assumes Hermiticity that Liouvillians do not - have. - - The promotion to ``complex128`` below is what makes the "always zgeev" - contract TRUE: unlike ``numpy.linalg`` (which computes in double - regardless of input dtype and only casts the result back), - ``scipy.linalg.eig`` dispatches by dtype and would run single-precision - ``cgeev`` on a ``complex64`` input — measured at ~30x the eigenvalue error - of the double solve on the same stored matrix. Every zero-mode/backward- - error tolerance downstream (issue #108) is calibrated against the double - solve, so the solver precision must not silently follow the caller's - storage dtype. Representation error already present in single-precision - INPUT data is the caller's data quality and is not (cannot be) undone - here. + Input is promoted to complex128 so the solver precision does not silently + follow a caller's complex64 storage dtype. """ A = np.asarray(A, dtype=complex) if A.ndim != 2 or A.shape[0] != A.shape[1]: @@ -85,30 +70,16 @@ def eig_nonhermitian( @dataclass(frozen=True, slots=True) class ZeroModeCertificate: - """Outcome of the structural zero-mode check on a computed spectrum. - - ``applicable`` is False when the operator is not trace preserving, in which - case no zero eigenvalue is guaranteed and the check says nothing. - """ + """Outcome of the structural zero-mode check on a computed spectrum.""" applicable: bool certified: bool solver: str residual: float - #: The RAW operator-derived band ``rtol * eps * ||L||_2`` that certification - #: was carried out against. REPORT AND DIAGNOSTICS ONLY -- it is NOT a - #: filter cutoff. After an a posteriori refinement (issue #113 second axis) - #: it is strictly LARGER than the tolerance actually applied, so filtering - #: by it discards the very slow mode the certificate rescued. Every filter - #: takes its cutoff from :meth:`zero_set_tolerance`. bound: float trace_defect: float zero_mode_count: int = 1 ambiguous_count: int = 0 - #: Tolerance the gap filters must apply. Equal to ``bound`` unless the a - #: posteriori certificate (issue #113 second axis) pulled a genuine slow - #: mode out of the band -- then it separates the refined zero set from that - #: mode, so D1/D3/D4 stop discarding physics the certificate has kept. zero_tolerance: float = float("nan") @property @@ -123,35 +94,7 @@ def zero_set_tolerance( atol: float | None = None, name: str = "eigenvalues", ) -> float: - """The ONE cutoff a zero-mode filter may apply to this spectrum. - - Round-17 review (PR #121). Five call sites had independently written - the same two-branch expression -- the certificate's band when it is - applicable, the radius-based proxy otherwise -- and four of them were - never migrated when the a posteriori refinement (issue #113 second - axis) made ``bound`` the wrong half of that branch. Duplicating a - DECISION at five sites is what made a four-fold miss possible, so the - decision is made here, once: - - * applicable certificate -> :attr:`applied_tolerance`: the raw band - ``bound``, unless the refinement pulled a genuine slow mode out of - it, in which case the refined split (never larger) is used, so the - rescued mode survives the filter it was rescued for; - * inapplicable certificate -> the radius-based proxy of - :func:`liouscope.numerics.scale.spectral_zero_tolerance`. Without - established trace preservation no zero mode is guaranteed and the - operator-norm band can exceed the entire spectrum of a strongly - non-normal input (measured: 2.2e3 against eigenvalues of order 1), - which would discard every mode and report a gapless D1; - * an explicit ``atol`` (the legacy pre-#108 absolute opt-in) still - wins over both. - - The spectrum is validated by ``spectral_zero_tolerance`` in every - branch, so an override cannot reintroduce silent acceptance of - corrupted solver output. - """ - # Deferred import: ``numerics.scale`` imports ``require_finite_square_2d`` - # from this module, so a module-level import would be circular. + """Return the one cutoff a zero-mode filter may apply to this spectrum.""" from .scale import spectral_zero_tolerance if atol is None and self.applicable: @@ -160,22 +103,7 @@ def zero_set_tolerance( @property def resolved(self) -> bool: - """True iff every mode is decidably zero or decidably non-zero. - - A count of zero modes greater than one is NOT by itself a problem: a - conserved quantity or symmetry sector gives a genuinely degenerate - stationary manifold whose extra zero modes sit at *machine zero*, and - the gap read off the complement is correct. - - The failure mode (issue #113) is different and is what - ``ambiguous_count`` measures: eigenvalues that fall inside the #108 - zero-mode tolerance while being far ABOVE the bare backward error - ``eps * ||L||``. Such a mode is neither machine-zero nor resolved -- - it lives inside the safety factor the tolerance carries, so the - arithmetic cannot say whether it is a genuine slow mode or noise. - When any exist, D1 would silently discard real physics and report the - next surviving (fast) mode instead. - """ + """True iff every mode is decidably zero or decidably non-zero.""" return self.applicable and self.certified and self.ambiguous_count == 0 def as_dict(self) -> dict[str, object]: @@ -204,47 +132,7 @@ def operator_zero_tolerance( rtol: float = ZERO_MODE_EPS_FACTOR, name: str = "L_super", ) -> float: - """Operator-derived zero-mode tolerance ``rtol * eps * ||L||_2``. - - Round-13 review (issue #112 follow-up): certification and downstream - filtering must use ONE scale. The certificate accepts a stationary - residual up to the true eigensolver backward error ``eps * ||L||_2``, - but :func:`liouscope.numerics.scale.spectral_zero_tolerance` -- which - only sees the spectrum -- can only use the spectral radius - ``max|lambda|`` as a proxy. For a strongly non-normal generator - ``||L||_2`` exceeds the radius by orders of magnitude (measured: 3.9e3 - on a 4x4 trace-preserving example), so a certified-resolved zero mode - with residual between the two thresholds survived the radius filter and - D1 reported ``~1e-12`` -- occasionally NEGATIVE -- instead of the true - gap ``1``. - - Every consumer that holds the OPERATOR (not just its spectrum) must - therefore filter on this scale -- but a consumer holding a - :class:`ZeroModeCertificate` must take its cutoff from - :meth:`ZeroModeCertificate.zero_set_tolerance` and NOT from the - certificate's ``bound``. Round-17 review (PR #121): the two WERE the - same expression when this paragraph was written, and stopped being so - with the a posteriori refinement (issue #113 second axis), which lowers - the applied tolerance below ``bound`` precisely when a genuine slow mode - was rescued from the band. Filtering by ``bound`` there discards that - mode again, one layer after the certificate saved it -- four call sites - did exactly that, following this paragraph. The radius-based - :func:`~liouscope.numerics.scale.spectral_zero_tolerance` remains the - correct fallback for spectrum-only call sites. Genuine slow modes - falling inside this coarser band are not silently swallowed: they are - flagged by the certificate's ambiguity split (``resolved=False``, - issue #113), which floors the verdict instead. - - Parameters - ---------- - L_super - The operator whose eigensolve backward error sets the scale. - rtol - Multiplier on ``eps * ||L||_2``, shared with the certificate and - with ``spectral_zero_tolerance``. - name - Argument name for the fail-closed validation message. - """ + """Operator-derived zero-mode tolerance ``rtol * eps * ||L||_2``.""" if not np.isfinite(rtol) or rtol < 0.0: raise ValueError(f"rtol must be finite and non-negative, got {rtol}") L_c = np.asarray(require_finite_square_2d(L_super, name=name), dtype=complex) @@ -252,65 +140,69 @@ def operator_zero_tolerance( return float(rtol * float(np.finfo(float).eps) * norm2) +def _hilbert_dimension_from_superoperator(L_super: np.ndarray) -> int | None: + """Return d for a d^2 x d^2 superoperator, else ``None``.""" + arr = np.asarray(L_super) + if arr.ndim != 2 or arr.shape[0] != arr.shape[1]: + return None + n = int(arr.shape[0]) + d = int(round(np.sqrt(n))) + return d if d * d == n else None + + def trace_preservation_defect(L_super: np.ndarray) -> tuple[float, float]: - """Return robust ``(||vec(I)^H L||_2, ||L||_F)`` for a generator. - - Every GKSL generator is trace preserving, which in column-stacking - convention is the exact statement ``vec(I)^H L = 0``. The quotient of the - two returned numbers is therefore a dimensionless, unit-invariant measure - of how far ``L`` is from being a legal generator. - - Issue #130: both quantities are evidence, so their arithmetic must be - correct independently of the policy that consumes them. Ordinary - ``np.linalg.norm`` squares before summing and can therefore return the - false sentinels ``0.0`` (underflow) or ``inf`` (overflow) even when the - mathematical norm of finite input is representable. Both norms now use - :func:`scaled_euclidean_norm`, an xLASSQ-style power-of-two scaled - sum-of-squares primitive. It returns a finite norm whenever the true - float64 result is representable and returns infinity only when the input or - the true norm is genuinely non-representable. - - This deliberately separates measurement from certificate policy. A very - large but finite operator must not be rejected merely because a naive norm - overflowed; structural guards such as ``band_discriminates`` decide whether - the resulting zero-mode evidence actually discriminates. + """Return robust ``(||vec(I)^H L||_2, ||L||_F)``. + + Both quantities are measurements and therefore use the same xLASSQ-style + scaled Euclidean norm. They remain finite whenever the true float64 norm is + representable; policy is handled separately by + :func:`trace_preservation_equation_defect` and the certificate gates. """ L_super = np.asarray(L_super) - n = L_super.shape[0] - d = int(round(np.sqrt(n))) - if d * d != n: + d = _hilbert_dimension_from_superoperator(L_super) + if d is None: return float("nan"), float("nan") vec_i = np.eye(d, dtype=complex).reshape(-1, order="F") defect_vector = vec_i.conj() @ L_super return scaled_euclidean_norm(defect_vector), scaled_euclidean_norm(L_super) +def trace_preservation_equation_defect(L_super: np.ndarray) -> float: + """Return the worst componentwise relative defect of ``vec(I)^H L = 0``. + + A global quotient ``||vec(I)^H L|| / ||L||`` is useful as a measurement, + but it is too weak as the precondition for an exact zero-mode theorem: an + unrelated entry of size ``1e308`` can make an order-one trace defect look + relatively tiny. Each column of ``vec(I)^H L = 0`` is therefore judged + against only the terms that actually participate in that scalar equation: + + ``eta_j = |sum_i L_(ii),j| / sum_i |L_(ii),j|``. + + The returned value is ``max_j eta_j``. This is a componentwise relative + backward-error criterion, analogous to LAPACK's BERR philosophy. The ratio + is evaluated with power-of-two scaling, so it is invariant to representable + rate-unit rescaling and remains defined in the subnormal range. + """ + L_super = np.asarray(L_super) + d = _hilbert_dimension_from_superoperator(L_super) + if d is None: + return float("nan") + if not np.all(np.isfinite(L_super)): + return float("nan") + identity_rows = np.arange(d, dtype=int) * (d + 1) + worst = 0.0 + for j in range(L_super.shape[1]): + ratio = scaled_cancellation_ratio(L_super[identity_rows, j]) + if not np.isfinite(ratio): + return float("nan") + worst = max(worst, ratio) + return float(worst) + + def band_discriminates( magnitudes: np.ndarray, zero_count: int, bound: float ) -> bool: - """Did the zero-mode band decide anything, or accept everything? - - Round-20 review (PR #121). ``residual <= bound`` certifies that a zero - mode was found. It says nothing about whether the band SEPARATED that mode - from the rest, and on the reviewer's counterexample it did not: a finite - 4x4 with cancelling ``+-1e308`` entries has ``||L||_2 = 1.4e308``, hence a - band of ``3.1e295``, and all four eigenvalues ``{1, 2, 3, 4}`` fall inside - it. The certificate came back ``certified=True, resolved=True`` with - ``zero_mode_count = 4`` -- a claim that the ENTIRE space is stationary, - for a generator whose spectrum is nothing of the kind. - - An acceptance region that contains every measured value has discriminated - nothing; passing it is not evidence. Note this is NOT a tightened - threshold -- tightening ``rtol`` would move the same band, not repair its - logic -- but a second, structural axis: the band must have had the - OPPORTUNITY to reject. - - The one legitimate way for a band to contain the whole spectrum is - ``bound == 0``: the exactly-zero generator, whose eigenvalues are exactly - zero and for which "everything is stationary" is the correct physics - measured exactly. A band of width zero accepts only exact zeros, so it - discriminates by construction and is excluded here. - """ + """Return whether the zero-mode acceptance band rejected anything.""" if bound <= 0.0: return True return zero_count < int(magnitudes.size) @@ -325,53 +217,15 @@ def certified_nonzero_modes( left_vectors: np.ndarray | None = None, margin: float = ZERO_MODE_APOSTERIORI_MARGIN, ) -> np.ndarray: - """Modes inside the zero band that are provably NOT the stationary mode. - - The magnitude band ``|lambda| <= rtol * eps * ||L||`` is a worst-case - backward-error argument: it assumes the eigensolver was as inaccurate as it - is allowed to be. On a well-conditioned generator it usually was not, and - the assumption is what discards genuine slow modes -- the - ``ZERO_MODE_AMBIGUITY_FACTOR`` comment records that no threshold on - ``|lambda|`` alone can separate the two populations, which are measured - only a factor of 2 apart. - - This is the second axis. For each in-band mode we compute the *a - posteriori* bound from the residuals actually attained, - - |lambda - lambda_hat| <~ max(||L x - lambda_hat x||, - ||L^H y - conj(lambda_hat) y||) / |y^H x| - - with ``||x|| = ||y|| = 1``. Were ``lambda`` exactly zero, that bound would - force ``|lambda_hat| <= bound``; so ``|lambda_hat| > margin * bound`` - certifies a genuine mode instead of merely scoring one. ``margin`` covers - the first-order caveat only -- see ``ZERO_MODE_APOSTERIORI_MARGIN`` for the - corpus measurement behind its value. - - Deliberately one-directional: the result can only move a mode OUT of the - zero set, never into it, so it cannot manufacture the false "unresolved" - verdict that the ambiguity split exists to avoid. - - Two invariants are enforced fail-closed: - - * fewer than two in-band modes -> nothing is certified. A trace-preserving - generator has at least one stationary mode, so a lone in-band mode is the - expected one and there is no decision to make. - * the smallest-magnitude in-band mode is never certified. Trace - preservation guarantees a zero eigenvalue exists; emptying the band would - contradict the very precondition under which this certificate applies. - """ + """Identify in-band modes that a posteriori residuals prove non-stationary.""" idx = np.flatnonzero(in_band) out = np.zeros(eigenvalues.shape, dtype=bool) if idx.size < 2: return out - # The stationary mode itself is never a candidate (see docstring). idx = idx[idx != idx[int(np.argmin(np.abs(eigenvalues[idx])))]] vr, vl, ref = right_vectors, left_vectors, eigenvalues if vr is None or vl is None: - # Own decomposition: the certificate is a property of the OPERATOR, so - # it is legitimate to certify a candidate ladder's spectrum from the - # operator's own eigenpairs -- but only where the two agree. try: ref, vl, vr = sla.eig(L_c, left=True, right=True) except (ValueError, sla.LinAlgError): @@ -382,15 +236,13 @@ def certified_nonzero_modes( lam = eigenvalues[i] j = int(np.argmin(np.abs(ref - lam))) mag = float(abs(ref[j])) - # Fail closed when the borrowed decomposition disagrees about this mode: - # a bound computed for a different eigenvalue certifies nothing. if abs(ref[j] - lam) > 0.1 * max(abs(lam), mag) or mag == 0.0: continue x = vr[:, j] / max(float(np.linalg.norm(vr[:, j])), tiny) y = vl[:, j] / max(float(np.linalg.norm(vl[:, j])), tiny) sep = float(abs(np.vdot(y, x))) if sep <= 0.0: - continue # defective pair: the first-order bound does not apply + continue res = max( float(np.linalg.norm(L_c @ x - ref[j] * x)), float(np.linalg.norm(L_c.conj().T @ y - np.conj(ref[j]) * y)), @@ -409,22 +261,7 @@ def refine_zero_band( right_vectors: np.ndarray | None = None, left_vectors: np.ndarray | None = None, ) -> tuple[np.ndarray, float]: - """Zero-mode membership plus the tolerance the gap filters must then apply. - - Returns ``(in_band, zero_tolerance)``. Without a certified non-stationary - mode both are exactly the pre-refinement values, so the healthy path is - unchanged bit for bit. - - The tolerance has to travel with the membership: D1/D3/D4 filter by - magnitude, so a mode this certificate rescues would be discarded again one - layer later if they kept reading the raw ``bound``. It is placed at the - geometric mean of the two populations -- the scale-free midpoint, which is - the right notion of "between" for quantities spanning decades. - - Fail-closed on inversion: if some mode kept as zero is LARGER than the - smallest certified one, no single threshold can express the split, and the - refinement is abandoned wholesale rather than applied half-way. - """ + """Return refined zero membership and its corresponding filter tolerance.""" nonzero = certified_nonzero_modes( L_c, eigenvalues, @@ -441,117 +278,41 @@ def refine_zero_band( return magnitudes <= bound, bound if lo <= 0.0: return in_band, hi * 0.5 - # ROUND-22 REVIEW (PR #121). ``sqrt(lo * hi)`` forms the PRODUCT first, - # and that product underflows to exactly 0 whenever ``lo * hi`` drops - # below ~5e-324 -- while ``lo``, ``hi`` and their geometric mean are all - # perfectly representable. The tolerance then becomes 0.0, the strict - # ``|lambda| > tol`` filters keep the numerical stationary residual as a - # physical mode, and D1 reports it as the gap. Since ``lo * hi`` scales as - # c^2 under a uniform rate rescale ``L -> cL`` while the refinement itself - # is scale free, the failure is purely a units artefact: the same physics - # in smaller rate units silently changes the answer. - # - # ``sqrt(lo) * sqrt(hi)`` is the same number mathematically and forms each - # factor separately, so it can neither underflow (sqrt of the smallest - # subnormal is ~2.2e-162) nor overflow (sqrt of the largest float is - # ~1.3e154, and the product of two such factors is at most the largest - # float again). return in_band, float(np.sqrt(lo) * np.sqrt(hi)) +def _tp_certificate_applicable( + L_c: np.ndarray, + tp_defect: float, + fro: float, + tp_rtol: float, +) -> bool: + """Return whether trace preservation is established strongly enough. + + The historical normwise criterion is retained for compatibility and broad + scale checking. The componentwise equation criterion is the load-bearing + addition from #130: it prevents unrelated large entries from diluting a + local violation of the exact trace equations. + """ + equation_defect = trace_preservation_equation_defect(L_c) + if not np.isfinite(tp_defect) or not np.isfinite(fro): + return False + if not np.isfinite(equation_defect): + return False + if equation_defect > tp_rtol: + return False + if tp_defect > tp_rtol * fro: + return False + return True + + def certified_eigvals( L_super: np.ndarray, *, rtol: float = ZERO_MODE_EPS_FACTOR, tp_rtol: float = 1.0e-10, ) -> tuple[np.ndarray, ZeroModeCertificate]: - """Eigenvalues of a Liouvillian, checked against an exact structural fact. - - Issue #112. For a trace-preserving generator ``vec(I)^H L = 0`` holds - *exactly*, so ``0`` is an exact eigenvalue and a correct eigensolve must - return some ``|lambda| <= rtol * eps * ||L||``. This is a **theorem, not a - threshold**: when the computed spectrum contains no such eigenvalue the - solve has demonstrably failed, and no choice of zero-mode tolerance - (relative or absolute, issue #108) can repair it -- the missing information - is not in the eigenvalues. - - Why this is needed - ------------------ - LAPACK's complex non-symmetric driver (``zgeev``) deflates a subdiagonal - when the Ahues-Tisseur test ``|h[i,i-1] h[i-1,i]| <= eps |h[i,i]| - |h[i-1,i-1] - h[i,i]|`` passes. On a STIFF generator the large diagonal - entries inflate that right-hand side: with a fast rate ``2.7e5`` alongside - slow rates ``~1e-5``, the bound reaches ``~1.6e-5`` and the entire slow - block deflates to its own DIAGONAL. Measured on a four-level classical jump - network, the exact zero mode vanished from the returned spectrum and a - spurious mode appeared in its place, so ``D1`` reported ``7.28e-6`` for a - generator whose true gap is ``1.07e-5``. - - That failure is invisible to conditioning: the per-eigenvalue condition - numbers of the wrong spectrum are ``1``-``25``, i.e. the solver reports the - wrong answer as well conditioned. Balancing does not fix it either - (measured). What does fix it is solving the same matrix by a route whose - deflation is not driven by the stiff diagonal. - - Repair ladder - ------------- - Candidates are tried in order and the FIRST one satisfying the certificate - WITHOUT ambiguous in-band modes is returned; the primary solve is - therefore byte-identical whenever it is already correct (measured: - 400/400 random four-level networks, stiff and well-conditioned alike, - keep the ``zgeev`` result). An ambiguous candidate does not end the - ladder (round-15 review): a later route can resolve the very modes it - cannot -- measured, ``zgeev`` certified a stiff network with one - ambiguous mode at ``6.3e-7`` while ``dgeev-real`` returned a clean - ``2.1e-15`` stationary mode. The best ambiguous candidate (fewest - ambiguous modes) is returned, fail-closed, only when every route stays - ambiguous. - - 1. ``zgeev`` -- the incumbent, via :func:`eig_nonhermitian`; - 2. ``dgeev`` on the real part, when ``L`` is real-valued (a different - LAPACK driver with a different deflation path); - 3. the complex Schur form ``zgees``, which deflates on the full Hessenberg - structure rather than eigenvalue-by-eigenvalue; - 4. an explicitly balanced ``zgeev``. - - If no candidate is certified, the eigenvalues of the best (smallest - residual) candidate are returned together with ``certified=False`` -- the - caller decides how loudly to fail. The function never silently substitutes - an uncertified spectrum for a certified one. - - Parameters - ---------- - L_super - ``d^2 x d^2`` vectorised generator. - rtol - Multiplier on the backward error ``eps * ||L||_2``, shared with - :func:`liouscope.numerics.scale.spectral_zero_tolerance`. - tp_rtol - Relative tolerance deciding whether ``L`` is trace preserving at all. - A non-trace-preserving operator has no guaranteed zero mode, so the - certificate is reported as ``applicable=False`` and the incumbent - spectrum is returned untouched. - - Returns - ------- - tuple - ``(eigenvalues, certificate)``. - """ - # ROUND-22 REVIEW (PR #121). A tolerance that is not a number cannot - # decide anything. With ``tp_rtol = NaN`` the trace-preservation test - # ``defect > tp_rtol * fro`` is False for EVERY operator (all NaN - # comparisons are), and with ``tp_rtol = inf`` its right-hand side is - # infinite -- so a demonstrably non-trace-preserving generator walks - # straight past the applicability gate and can come back - # ``applicable=True, certified=True``. Measured before this guard: - # ``certified_eigvals(diag([0,-1,-2,-3]), tp_rtol=nan)`` reported - # applicable and certified with a trace defect of 3.0. - # - # Round 21 closed exactly this shape one variable further along (a - # non-finite REFERENCE scale ``fro``) and left the PARAMETER beside it - # unchecked. The generalisation is the fix: every tolerance this layer - # compares against is validated the way ``operator_zero_tolerance`` and - # ``spectral_zero_tolerance`` already validate theirs. + """Eigenvalues checked against the exact trace-preserving zero-mode fact.""" if not np.isfinite(rtol) or rtol < 0.0: raise ValueError(f"rtol must be finite and non-negative, got {rtol}") if not np.isfinite(tp_rtol) or tp_rtol < 0.0: @@ -562,13 +323,6 @@ def certified_eigvals( norm2 = float(np.linalg.norm(L_c, 2)) if L_c.size else 0.0 bound = rtol * float(np.finfo(float).eps) * norm2 - # Round-17 review (PR #121). The primary solve must not be able to end the - # ladder. ``zgeev`` raising ``LinAlgError`` IS a nonconvergence -- exactly - # the failure the real-driver / Schur / balanced routes exist to repair -- - # and letting the exception propagate defeated the ladder precisely in the - # case it was built for. The error is carried instead and re-raised only - # if NO route produced a spectrum at all, so the caller still sees the - # original diagnosis when nothing worked. primary: np.ndarray | None primary_error: sla.LinAlgError | None try: @@ -576,44 +330,11 @@ def certified_eigvals( primary_error = None except sla.LinAlgError as exc: primary, primary_error = None, exc - if ( - not np.isfinite(tp_defect) - # ROUND-21: a non-finite REFERENCE scale is as unusable as a - # non-finite defect. ``defect <= tp_rtol * fro`` against an - # infinite ``fro`` admits every operator, so the band that - # follows is drawn from a scale that was never measured; which - # modes it then swallows is decided by the LAPACK build (2 of 4 - # on the CI runners, 4 of 4 locally). Trace preservation - # RELATIVE to the operator's own scale is not a statement that - # can be made about infinity, so it is refused here. - or not np.isfinite(fro) - # ROUND-23 REVIEW (PR #121), same finding one line further down. - # ``max(fro, tiny)`` floors the REFERENCE SCALE at a constant that - # is larger than the operator whenever the operator is subnormal, - # so the relative test stops being relative: measured for - # ``diag([0, -1e-320, -2e-320, -3e-320])``, defect 3e-320 against a - # floor-derived bound of 2.2e-318, hence 'trace preserving' for an - # operator whose relative defect is 0.80. Repairing the norms above - # without this leaves the rate-unit dependence intact below ~1e-317. - # The floor protected nothing: ``defect <= sqrt(d) * fro`` holds, so - # ``fro == 0`` forces ``defect == 0`` and the comparison ``0 > 0`` is - # False either way -- the exactly-zero generator stays applicable. - or tp_defect > tp_rtol * fro - ): - # No certificate applies without trace preservation, so there is - # nothing to repair TOWARDS and the ladder is not run: the primary - # failure is the honest answer. + + if not _tp_certificate_applicable(L_c, tp_defect, fro, tp_rtol): if primary is None: assert primary_error is not None raise primary_error - # ROUND-22 REVIEW (PR #121). ``zero_mode_count`` kept its dataclass - # default of 1 here, which asserts a stationary mode that nothing - # counted and nothing guarantees -- the whole point of - # ``applicable=False`` is that no zero eigenvalue is implied. Both - # certificate APIs reported ``zero_mode_count: 1`` for - # ``diag([1,2,3,4])``, a spectrum containing no zero at all, and the - # dict is persisted as audit metadata. 0 is the honest count: nothing - # was certified as stationary. return primary, ZeroModeCertificate( applicable=False, certified=False, @@ -625,31 +346,9 @@ def certified_eigvals( ) def _candidates() -> Iterator[tuple[str, np.ndarray]]: - # LAZY by design (twelfth-round review): the ladder stops at the first - # certified spectrum, and the healthy path is the primary solve -- an - # eager list would charge every caller three additional cubic - # decompositions that the loop then never looks at. if primary is not None: yield ("zgeev", primary) - # ``dgeev`` genuinely solves a DIFFERENT matrix unless L is EXACTLY - # real: ``np.allclose`` carries an absolute default atol, so a stiff - # complex generator in small rate units read as "real" and the repair - # silently deleted its Hamiltonian part (measured: D3 moved from 1e-7 - # to ~9.7e-6 under a pure L -> 1e-10 L rescale). Exact realness is the - # only scale-invariant criterion under which the route is valid. if not np.any(L_c.imag): - # GUARDED (round-16 review): an exactly real generator whose - # ``dgeev`` solve raises must not end the ladder. Without this - # suppression the exception propagates out of the generator and - # the two remaining repair routes below -- ``zgees-schur`` and - # ``balanced-zgeev`` -- never run, even when either could recover - # a certifiable spectrum. The route is a REPAIR attempt like the - # others; a repair that fails is a route that did not work, not a - # reason to abandon the remaining ones. - # - # The sibling ladder in ``certified_eig`` already guarded exactly - # this step. The omission here was the inconsistency, not the - # guard there. with contextlib.suppress(ValueError, sla.LinAlgError): yield ("dgeev-real", np.linalg.eigvals(L_c.real).astype(complex)) with contextlib.suppress(ValueError, sla.LinAlgError): @@ -658,28 +357,6 @@ def _candidates() -> Iterator[tuple[str, np.ndarray]]: balanced = sla.matrix_balance(L_c, permute=True)[0] yield ("balanced-zgeev", np.linalg.eigvals(balanced)) - # Split the zero-tolerance band (issue #113). Below the split a mode is - # machine-zero; above it -- but still inside the #108 tolerance, which - # carries a factor ``rtol`` on top -- it is neither machine-zero nor - # resolved. - # - # The factor is set from a MEASURED distribution, not from taste. Across - # 83 healthy generators (single/two-qubit families plus random GKSL, - # unique and degenerate steady states alike) the largest in-band - # |lambda| reaches 2.38 * eps*||L||, median 0.39. Unresolved slow modes - # were measured at 4.87 and 4.87e2. The two populations are therefore - # NOT cleanly separated -- the nearest pair is 2.38 against 4.87, about - # 2x -- so the split is placed at 30x, roughly an order of magnitude - # clear of the healthy maximum rather than midway between the two. - # - # That choice is deliberately conservative: a false NaN on a healthy - # generator destroys a correct analysis, while a missed marginal case - # leaves behaviour exactly as it was before this check existed. The cost - # is bounded reach -- above a spectral spread of ~1e14 the slow modes - # sink below the split and the defect is undetectable by any magnitude - # test, because double precision has genuinely lost the information. - # That residual is documented in issue #113 rather than papered over - # here. split = ZERO_MODE_AMBIGUITY_FACTOR * float(np.finfo(float).eps) * norm2 best: tuple[str, np.ndarray, float] | None = None ambiguous_best: tuple[str, np.ndarray, float, int, int, float] | None = None @@ -687,21 +364,10 @@ def _candidates() -> Iterator[tuple[str, np.ndarray]]: magnitudes = np.abs(ev) residual = float(np.min(magnitudes)) if ev.size else float("inf") if residual <= bound: - # Issue #113 second axis: a mode the a posteriori bound certifies - # as non-stationary leaves the zero set before it is counted, so a - # genuine slow mode no longer masquerades as a second zero mode and - # hands D1 the NEXT eigenvalue as the gap. - in_band, zero_tolerance = refine_zero_band( - L_c, ev, magnitudes, bound - ) + in_band, zero_tolerance = refine_zero_band(L_c, ev, magnitudes, bound) zero_count = int(np.count_nonzero(in_band)) ambiguous = int(np.count_nonzero(in_band & (magnitudes > split))) if not band_discriminates(magnitudes, zero_count, bound): - # ROUND-20: the band swallowed the whole spectrum, so this - # route proved nothing. Carried on as an UNcertified - # candidate, exactly like a route whose residual missed the - # band -- a later route on a different deflation path may - # still produce a spectrum the band can separate. if best is None or residual < best[2]: best = (name, ev, residual) continue @@ -716,16 +382,14 @@ def _candidates() -> Iterator[tuple[str, np.ndarray]]: zero_mode_count=zero_count, zero_tolerance=zero_tolerance, ) - # Round-15 review: an AMBIGUOUS candidate must not end the - # ladder -- a later route can resolve the very modes this one - # cannot (measured: zgeev certified a stiff network with one - # ambiguous mode at 6.3e-7 while dgeev-real returns a clean - # 2.1e-15 stationary mode). It is retained only as the - # fail-closed fallback when EVERY route stays ambiguous, keyed - # by the fewest ambiguous modes (ties: ladder order). if ambiguous_best is None or ambiguous < ambiguous_best[4]: ambiguous_best = ( - name, ev, residual, zero_count, ambiguous, zero_tolerance + name, + ev, + residual, + zero_count, + ambiguous, + zero_tolerance, ) continue if best is None or residual < best[2]: @@ -745,7 +409,6 @@ def _candidates() -> Iterator[tuple[str, np.ndarray]]: zero_tolerance=zero_tolerance, ) if best is None: - # Every route raised (only reachable when the primary solve did too). assert primary_error is not None raise primary_error return best[1], ZeroModeCertificate( @@ -765,64 +428,7 @@ def certified_eig( rtol: float = ZERO_MODE_EPS_FACTOR, tp_rtol: float = 1.0e-10, ) -> tuple[EigenDecomposition, ZeroModeCertificate]: - """Like :func:`certified_eigvals`, but also returns left/right eigenvectors. - - Issue #112 follow-up. The eigenvalue-only certificate protects D1/D3/D4, but - the layers that consume eigen*vectors* -- D19's slowest-mode overlap above - all -- were still reading the uncertified decomposition. That matters more - than it sounds: on a stiff generator ``zgeev`` can return a *spurious* slow - mode, and D19 then measures the initial state's overlap with an eigenvector - that corresponds to no physical mode. A wrong overlap on that rung fires a - false A11/F4 Mpemba candidate -- the highest-priority branch of the - classifier, and exactly the failure mode #108 was fixed to prevent. - - Repair ladder, narrower than the eigenvalue-only one - ---------------------------------------------------- - Only routes that actually produce eigenvectors in the ORIGINAL basis are - usable here: - - 1. ``zgeev`` -- the incumbent; - 2. ``dgeev`` on the real part, when ``L`` is real-valued. - - The Schur and balanced routes are deliberately absent: ``zgees`` yields no - eigenvectors directly, and the balanced solve returns them in the balanced - basis, where back-transforming re-introduces exactly the scaling that made - the decomposition unreliable. A narrower ladder is the honest choice -- - it means some spectra that ``certified_eigvals`` can repair are NOT - repairable here, and the certificate then reports ``certified=False`` - rather than pretending otherwise. - - Beyond the eigenvalue certificate, a candidate is accepted only when - every consumed mode's left AND right eigen*vector* residuals pass the - per-mode relative gate ``r_j <= max(VECTOR_RESIDUAL_REL_MAX * |lambda_j|, - bound)`` (round-15 review) -- the vectors are what D19 and the Petermann - factors consume, and a small ``|lambda|`` does not vouch for them. When - the vectors fail on every route, the returned certificate carries - ``certified=False`` with ``residual`` set to the offending (vector) - residual, so the downstream warning shows a number that actually - exceeds the printed bound. Ambiguous candidates likewise no longer end - the ladder (same rule as :func:`certified_eigvals`). - - Returns - ------- - tuple - ``(EigenDecomposition with left_vectors set, certificate)``. - """ - # ROUND-22 REVIEW (PR #121). A tolerance that is not a number cannot - # decide anything. With ``tp_rtol = NaN`` the trace-preservation test - # ``defect > tp_rtol * fro`` is False for EVERY operator (all NaN - # comparisons are), and with ``tp_rtol = inf`` its right-hand side is - # infinite -- so a demonstrably non-trace-preserving generator walks - # straight past the applicability gate and can come back - # ``applicable=True, certified=True``. Measured before this guard: - # ``certified_eigvals(diag([0,-1,-2,-3]), tp_rtol=nan)`` reported - # applicable and certified with a trace defect of 3.0. - # - # Round 21 closed exactly this shape one variable further along (a - # non-finite REFERENCE scale ``fro``) and left the PARAMETER beside it - # unchecked. The generalisation is the fix: every tolerance this layer - # compares against is validated the way ``operator_zero_tolerance`` and - # ``spectral_zero_tolerance`` already validate theirs. + """Like :func:`certified_eigvals`, also returning left/right eigenvectors.""" if not np.isfinite(rtol) or rtol < 0.0: raise ValueError(f"rtol must be finite and non-negative, got {rtol}") if not np.isfinite(tp_rtol) or tp_rtol < 0.0: @@ -833,13 +439,6 @@ def certified_eig( norm2 = float(np.linalg.norm(L_c, 2)) if L_c.size else 0.0 bound = rtol * float(np.finfo(float).eps) * norm2 - # Round-17 review (PR #121). The primary solve must not be able to end the - # ladder. ``zgeev`` raising ``LinAlgError`` IS a nonconvergence -- exactly - # the failure the real-driver / Schur / balanced routes exist to repair -- - # and letting the exception propagate defeated the ladder precisely in the - # case it was built for. The error is carried instead and re-raised only - # if NO route produced a spectrum at all, so the caller still sees the - # original diagnosis when nothing worked. primary: EigenDecomposition | None primary_error: sla.LinAlgError | None try: @@ -847,41 +446,11 @@ def certified_eig( primary_error = None except sla.LinAlgError as exc: primary, primary_error = None, exc - if ( - not np.isfinite(tp_defect) - # ROUND-21: a non-finite REFERENCE scale is as unusable as a - # non-finite defect. ``defect <= tp_rtol * fro`` against an - # infinite ``fro`` admits every operator, so the band that - # follows is drawn from a scale that was never measured; which - # modes it then swallows is decided by the LAPACK build (2 of 4 - # on the CI runners, 4 of 4 locally). Trace preservation - # RELATIVE to the operator's own scale is not a statement that - # can be made about infinity, so it is refused here. - or not np.isfinite(fro) - # ROUND-23 REVIEW (PR #121), same finding one line further down. - # ``max(fro, tiny)`` floors the REFERENCE SCALE at a constant that - # is larger than the operator whenever the operator is subnormal, - # so the relative test stops being relative: measured for - # ``diag([0, -1e-320, -2e-320, -3e-320])``, defect 3e-320 against a - # floor-derived bound of 2.2e-318, hence 'trace preserving' for an - # operator whose relative defect is 0.80. Repairing the norms above - # without this leaves the rate-unit dependence intact below ~1e-317. - # The floor protected nothing: ``defect <= sqrt(d) * fro`` holds, so - # ``fro == 0`` forces ``defect == 0`` and the comparison ``0 > 0`` is - # False either way -- the exactly-zero generator stays applicable. - or tp_defect > tp_rtol * fro - ): + + if not _tp_certificate_applicable(L_c, tp_defect, fro, tp_rtol): if primary is None: assert primary_error is not None raise primary_error - # ROUND-22 REVIEW (PR #121). ``zero_mode_count`` kept its dataclass - # default of 1 here, which asserts a stationary mode that nothing - # counted and nothing guarantees -- the whole point of - # ``applicable=False`` is that no zero eigenvalue is implied. Both - # certificate APIs reported ``zero_mode_count: 1`` for - # ``diag([1,2,3,4])``, a spectrum containing no zero at all, and the - # dict is persisted as audit metadata. 0 is the honest count: nothing - # was certified as stationary. return primary, ZeroModeCertificate( applicable=False, certified=False, @@ -897,7 +466,6 @@ def certified_eig( ) def _candidates() -> Iterator[tuple[str, EigenDecomposition]]: - # Lazy + exact realness, for the same reasons as the eigenvalue ladder. if primary is not None: yield ("zgeev", primary) if not np.any(L_c.imag): @@ -913,36 +481,12 @@ def _candidates() -> Iterator[tuple[str, EigenDecomposition]]: ) def _vector_residual(decomp: EigenDecomposition) -> float: - """Largest OFFENDING per-mode eigenvector residual; 0.0 when all pass. - - Round-15 review: the eigenvalue certificate alone does not guarantee - the eigen*vectors* -- the quantities D19 and the Petermann factors - actually consume. Measured on a valid stiff classical network: - certified and resolved with ``bound = 3.2e-7``, yet the LEFT vector - of the selected slow mode had residual ``3.1e-4`` -- three orders of - magnitude beyond the certificate, i.e. not an eigenvector of ``L`` - in any usable sense (the right vectors were fine). - - The acceptance test is PER MODE and RELATIVE: - ``r_j <= max(VECTOR_RESIDUAL_REL_MAX * |lambda_j|, bound)`` for every - non-zero-band mode, with ``r_j`` the larger of the unit-normalised - left/right residuals. ``r_j / |lambda_j|`` is the first-order - relative error scale of anything computed from the pair, and unlike - a single operator-scale cutoff it separates the measured - populations: a legitimate ``dgeev-real`` repair carries a few - percent on its slow modes, while a corrupt ``zgeev`` decomposition - carries 22% to 2900% (see - :data:`liouscope._consts.VECTOR_RESIDUAL_REL_MAX` for the measured - calibration). Modes inside the zero band are not gated here: their - vectors span the stationary manifold and are not consumed through - this decomposition. - """ ev = decomp.eigenvalues if not ev.size: return 0.0 vr = decomp.right_vectors vl = decomp.left_vectors - assert vl is not None # every ladder route computes left vectors + assert vl is not None tiny = np.finfo(float).tiny res_r = np.linalg.norm(L_c @ vr - vr * ev[None, :], axis=0) / np.maximum( np.linalg.norm(vr, axis=0), tiny @@ -965,14 +509,11 @@ def _vector_residual(decomp: EigenDecomposition) -> float: tuple[str, EigenDecomposition, float, int, int, float] | None ) = None vector_failed: tuple[str, EigenDecomposition, float, float] | None = None + for name, decomp in _candidates(): magnitudes = np.abs(decomp.eigenvalues) residual = float(np.min(magnitudes)) if magnitudes.size else float("inf") if residual <= bound: - # Issue #113 second axis: a mode the a posteriori bound certifies - # as non-stationary leaves the zero set before it is counted, so a - # genuine slow mode no longer masquerades as a second zero mode and - # hands D1 the NEXT eigenvalue as the gap. in_band, zero_tolerance = refine_zero_band( L_c, decomp.eigenvalues, @@ -984,9 +525,6 @@ def _vector_residual(decomp: EigenDecomposition) -> float: zero_count = int(np.count_nonzero(in_band)) ambiguous = int(np.count_nonzero(in_band & (magnitudes > split))) if not band_discriminates(magnitudes, zero_count, bound): - # ROUND-20, same guard as the sibling ladder. The two loops - # having drifted apart is what produced the round-16 finding; - # this one is added to both in the same commit. if best is None or residual < best[2]: best = (name, decomp, residual) continue @@ -1003,25 +541,24 @@ def _vector_residual(decomp: EigenDecomposition) -> float: zero_mode_count=zero_count, zero_tolerance=zero_tolerance, ) - # Eigenvalues fine, eigenvectors demonstrably not: keep for - # the fail-closed record and try the next route. if vector_failed is None or vec_residual < vector_failed[3]: vector_failed = (name, decomp, residual, vec_residual) continue - # Round-15 review: an ambiguous candidate must not end the - # ladder (same rule as certified_eigvals). if ambiguous_best is None or ambiguous < ambiguous_best[4]: ambiguous_best = ( - name, decomp, residual, zero_count, ambiguous, zero_tolerance + name, + decomp, + residual, + zero_count, + ambiguous, + zero_tolerance, ) continue if best is None or residual < best[2]: best = (name, decomp, residual) if ambiguous_best is not None: - name, decomp, residual, zero_count, ambiguous, zero_tolerance = ( - ambiguous_best - ) + name, decomp, residual, zero_count, ambiguous, zero_tolerance = ambiguous_best return decomp, ZeroModeCertificate( applicable=True, certified=True, @@ -1034,10 +571,6 @@ def _vector_residual(decomp: EigenDecomposition) -> float: zero_tolerance=zero_tolerance, ) if vector_failed is not None: - # No route produced consumable eigenvectors. certified=False keeps - # ``resolved`` False, so every consumer withholds -- the eigenvalue - # residual alone must not be allowed to read as success when the - # vectors it vouches for are not eigenvectors. name, decomp, residual, vec_residual = vector_failed return decomp, ZeroModeCertificate( applicable=True, @@ -1049,7 +582,6 @@ def _vector_residual(decomp: EigenDecomposition) -> float: zero_mode_count=0, ) if best is None: - # Every route raised (only reachable when the primary solve did too). assert primary_error is not None raise primary_error return best[1], ZeroModeCertificate( @@ -1064,23 +596,7 @@ def _vector_residual(decomp: EigenDecomposition) -> float: def require_finite_square_2d(A: np.ndarray, *, name: str = "matrix") -> np.ndarray: - """Validate ``A`` as a finite, square, 2-D array at a public boundary. - - Returns ``np.asarray(A)`` unchanged when valid, else raises a - :class:`ValueError` that names the offending argument *and* the concrete - defect. This is a fail-closed guard for public entry points: without it a - ``NaN``/``inf``-laden Liouvillian flows into ``scipy.linalg.expm`` / ``svd`` - and surfaces as an opaque, location-blind error deep in LAPACK - (``"array must not contain infs or NaNs"`` / ``"SVD did not converge"``), - leaving the caller no clue which input was malformed. - - Parameters - ---------- - A - Candidate array. - name - Human-readable argument name used in error messages. - """ + """Validate ``A`` as a finite, square, 2-D array at a public boundary.""" arr = np.asarray(A) if arr.ndim != 2 or arr.shape[0] != arr.shape[1]: raise ValueError(f"{name} must be a square 2-D array, got shape {arr.shape}") @@ -1097,14 +613,7 @@ def require_finite_square_2d(A: np.ndarray, *, name: str = "matrix") -> np.ndarr def hermiticity_defect(A: np.ndarray) -> tuple[float, float]: - """Return ``(max|A - A^H|, max|A|)`` -- the absolute defect and its scale. - - Both are elementwise sup-norms, so their quotient is the *relative* - Hermiticity defect: a dimensionless number that is unchanged by a pure - change of units ``A -> cA``. Returning the pair (rather than the quotient) - keeps the zero-operator case in the caller's hands, where the right - semantics are known. - """ + """Return ``(max|A - A^H|, max|A|)``.""" A = np.asarray(A) defect = float(np.max(np.abs(A - A.conj().T))) if A.size else 0.0 scale = float(np.max(np.abs(A))) if A.size else 0.0 @@ -1117,55 +626,7 @@ def is_hermitian( *, rtol: float = EPS_HERMITICITY, ) -> bool: - """Return True iff ``A`` is Hermitian. - - The gate is SCALE-RELATIVE by default (issue #109): ``A`` is accepted iff - ``max|A - A^H| <= rtol * max|A|``. An operator that is not trace-normalised - -- a Hamiltonian, a jump operator -- carries a physical dimension, so an - ABSOLUTE tolerance makes the validation verdict depend on the caller's - choice of units, and it fails in both directions: - - * **fail-open** (the serious one) -- a genuinely non-Hermitian ``H`` with a - relative defect of ``1e-6`` passes a fixed ``1e-9`` gate as soon as - ``||H|| <~ 1e-3``. The resulting generator is not GKSL at all, and every - downstream diagnostic then describes dynamics that are not a quantum - channel; - * **false rejection** -- an exactly Hermitian ``H`` obtained from a - similarity transform carries round-off ``~eps*||H||``, which exceeds a - fixed ``1e-9`` once ``||H|| >~ 1e8``, so valid input is refused. - - ``rtol`` is deliberately ``EPS_HERMITICITY = 1e-9``, i.e. the historical - absolute value re-read as a relative one. At unit scale (``max|A| ~ 1``, - the regime of every existing fixture) the gate is therefore unchanged; only - operators far from unit scale move, which is exactly the defect. Note that - ``1e-9`` is many orders above round-off (``~1e-16``) *by design*: this is a - validation gate against malformed input, not a round-off discriminator. - - ``rtol=0`` semantics inside :func:`numpy.allclose` are not used here. - ``np.allclose``'s own ``rtol`` compares against ``|b|`` elementwise, which - for ``b = A^H`` makes the tolerance largest exactly where ``A`` is largest - and vanish where ``A`` is zero -- not a statement about the operator's - scale. The scale is taken once, from ``max|A|``. - - Parameters - ---------- - A - Candidate matrix. Non-square (or non-2-D) input returns ``False``. - atol - Legacy ABSOLUTE tolerance. When given, ``A`` is accepted iff - ``max|A - A^H| <= atol``, reproducing the pre-#109 behaviour verbatim. - This is the opt-in for callers whose operator is normalised to a known - scale -- density matrices, where the absolute reading is already - relative to ``tr(rho) = 1`` (see :func:`is_density_matrix`). - rtol - Relative tolerance on ``max|A|``; ignored when ``atol`` is given. - - Notes - ----- - Zero-operator semantics: ``max|A| == 0`` gives ``tol = 0``, so the all-zero - matrix (exactly Hermitian) passes and nothing else does. This mirrors the - documented zero-operator behaviour of :func:`liouscope.numerics.scale.rate_scale`. - """ + """Return True iff ``A`` is Hermitian, scale-relative by default.""" A = np.asarray(A) if A.ndim != 2 or A.shape[0] != A.shape[1]: return False @@ -1180,23 +641,12 @@ def is_density_matrix( atol_trace: float = EPS_TRACE, atol_psd: float = 1.0e-10, ) -> bool: - """Return True iff ``rho`` is a valid density matrix. - - Checks Hermiticity, unit trace and positive semi-definiteness. - - The Hermiticity gate stays ABSOLUTE here (issue #109). A density matrix is - trace-normalised, so ``atol`` is already implicitly relative to a fixed - scale ``tr(rho) = 1``; the rate-dimension argument that makes an absolute - tolerance wrong for a Hamiltonian does not apply. Reading it as relative to - ``max|rho|`` would instead *tighten* the gate by up to a factor ``d`` for - nearly maximally mixed states -- a behaviour change with no defect behind it. - """ + """Return True iff ``rho`` is a valid density matrix.""" rho = np.asarray(rho) if not is_hermitian(rho, atol=1.0e-9): return False if abs(np.trace(rho) - 1.0) > atol_trace: return False - # Hermitian => use eigvalsh for stability. evals = np.linalg.eigvalsh((rho + rho.conj().T) / 2) return bool(evals.min() >= -atol_psd) @@ -1207,14 +657,7 @@ def support_check( *, eps: float = EPS_SUPP, ) -> tuple[bool, np.ndarray]: - """Verify ``supp(rho_initial) subset supp(rho_steady)``. - - Returns a tuple ``(ok, regularised_steady)`` where ``regularised_steady`` - is ``rho_steady + eps * I / d`` shifted to avoid singularities when - computing the relative entropy ``D(rho || pi)``. - - Anchor J: this guards against NaNs when ``rho_steady`` is rank-deficient. - """ + """Verify ``supp(rho_initial) subset supp(rho_steady)``.""" rho_initial = np.asarray(rho_initial) rho_steady = np.asarray(rho_steady) d = rho_steady.shape[0] @@ -1225,7 +668,6 @@ def support_check( return True, rho_steady null_space = evecs_steady[:, null_mask] - # Projector onto null space P_null = null_space @ null_space.conj().T leakage = float(np.real(np.trace(P_null @ rho_initial @ P_null))) ok = leakage < eps From f749da6a777de1d1d0821d81acd86edfd6c90387 Mon Sep 17 00:00:00 2001 From: Marco Hostettler Date: Tue, 1 Sep 2026 20:44:51 +0200 Subject: [PATCH 09/22] fix(#130): restore full linalg review context before minimal TP gate --- src/liouscope/numerics/linalg.py | 794 ++++++++++++++++++++++++++----- 1 file changed, 676 insertions(+), 118 deletions(-) diff --git a/src/liouscope/numerics/linalg.py b/src/liouscope/numerics/linalg.py index b7580fa1..72aa1892 100644 --- a/src/liouscope/numerics/linalg.py +++ b/src/liouscope/numerics/linalg.py @@ -1,10 +1,12 @@ """General linear-algebra helpers used across diagnostics. -The module deliberately separates numerical measurements from scientific -certificate policy. In particular, trace-preservation norms are measured with -scaled arithmetic, while applicability of the structural zero-mode theorem is -also checked equation-by-equation so unrelated large matrix entries cannot hide -a local trace defect. +Includes: + +* :func:`eig_nonhermitian` --- wraps ``scipy.linalg.eig`` (LAPACK ``zgeev``), + the non-Hermitian eigensolver required by Liouvillians. Anchor D. +* :func:`is_hermitian`, :func:`is_density_matrix` --- input validation. +* :func:`support_check` --- enforces ``supp(rho_0) subset supp(rho_ss)`` + with ``eps = 1e-12`` regularisation (anchor J). """ from __future__ import annotations @@ -25,7 +27,7 @@ ZERO_MODE_APOSTERIORI_MARGIN, ZERO_MODE_EPS_FACTOR, ) -from .norms import scaled_cancellation_ratio, scaled_euclidean_norm +from .norms import scaled_euclidean_norm @dataclass(frozen=True, slots=True) @@ -48,8 +50,21 @@ def eig_nonhermitian( ) -> EigenDecomposition: """Return the eigendecomposition of a generic non-Hermitian matrix. - Input is promoted to complex128 so the solver precision does not silently - follow a caller's complex64 storage dtype. + Always uses LAPACK ``zgeev`` (via ``scipy.linalg.eig``). Anchor D: never + fall back to ``eigh`` which assumes Hermiticity that Liouvillians do not + have. + + The promotion to ``complex128`` below is what makes the "always zgeev" + contract TRUE: unlike ``numpy.linalg`` (which computes in double + regardless of input dtype and only casts the result back), + ``scipy.linalg.eig`` dispatches by dtype and would run single-precision + ``cgeev`` on a ``complex64`` input — measured at ~30x the eigenvalue error + of the double solve on the same stored matrix. Every zero-mode/backward- + error tolerance downstream (issue #108) is calibrated against the double + solve, so the solver precision must not silently follow the caller's + storage dtype. Representation error already present in single-precision + INPUT data is the caller's data quality and is not (cannot be) undone + here. """ A = np.asarray(A, dtype=complex) if A.ndim != 2 or A.shape[0] != A.shape[1]: @@ -70,16 +85,30 @@ def eig_nonhermitian( @dataclass(frozen=True, slots=True) class ZeroModeCertificate: - """Outcome of the structural zero-mode check on a computed spectrum.""" + """Outcome of the structural zero-mode check on a computed spectrum. + + ``applicable`` is False when the operator is not trace preserving, in which + case no zero eigenvalue is guaranteed and the check says nothing. + """ applicable: bool certified: bool solver: str residual: float + #: The RAW operator-derived band ``rtol * eps * ||L||_2`` that certification + #: was carried out against. REPORT AND DIAGNOSTICS ONLY -- it is NOT a + #: filter cutoff. After an a posteriori refinement (issue #113 second axis) + #: it is strictly LARGER than the tolerance actually applied, so filtering + #: by it discards the very slow mode the certificate rescued. Every filter + #: takes its cutoff from :meth:`zero_set_tolerance`. bound: float trace_defect: float zero_mode_count: int = 1 ambiguous_count: int = 0 + #: Tolerance the gap filters must apply. Equal to ``bound`` unless the a + #: posteriori certificate (issue #113 second axis) pulled a genuine slow + #: mode out of the band -- then it separates the refined zero set from that + #: mode, so D1/D3/D4 stop discarding physics the certificate has kept. zero_tolerance: float = float("nan") @property @@ -94,7 +123,35 @@ def zero_set_tolerance( atol: float | None = None, name: str = "eigenvalues", ) -> float: - """Return the one cutoff a zero-mode filter may apply to this spectrum.""" + """The ONE cutoff a zero-mode filter may apply to this spectrum. + + Round-17 review (PR #121). Five call sites had independently written + the same two-branch expression -- the certificate's band when it is + applicable, the radius-based proxy otherwise -- and four of them were + never migrated when the a posteriori refinement (issue #113 second + axis) made ``bound`` the wrong half of that branch. Duplicating a + DECISION at five sites is what made a four-fold miss possible, so the + decision is made here, once: + + * applicable certificate -> :attr:`applied_tolerance`: the raw band + ``bound``, unless the refinement pulled a genuine slow mode out of + it, in which case the refined split (never larger) is used, so the + rescued mode survives the filter it was rescued for; + * inapplicable certificate -> the radius-based proxy of + :func:`liouscope.numerics.scale.spectral_zero_tolerance`. Without + established trace preservation no zero mode is guaranteed and the + operator-norm band can exceed the entire spectrum of a strongly + non-normal input (measured: 2.2e3 against eigenvalues of order 1), + which would discard every mode and report a gapless D1; + * an explicit ``atol`` (the legacy pre-#108 absolute opt-in) still + wins over both. + + The spectrum is validated by ``spectral_zero_tolerance`` in every + branch, so an override cannot reintroduce silent acceptance of + corrupted solver output. + """ + # Deferred import: ``numerics.scale`` imports ``require_finite_square_2d`` + # from this module, so a module-level import would be circular. from .scale import spectral_zero_tolerance if atol is None and self.applicable: @@ -103,7 +160,22 @@ def zero_set_tolerance( @property def resolved(self) -> bool: - """True iff every mode is decidably zero or decidably non-zero.""" + """True iff every mode is decidably zero or decidably non-zero. + + A count of zero modes greater than one is NOT by itself a problem: a + conserved quantity or symmetry sector gives a genuinely degenerate + stationary manifold whose extra zero modes sit at *machine zero*, and + the gap read off the complement is correct. + + The failure mode (issue #113) is different and is what + ``ambiguous_count`` measures: eigenvalues that fall inside the #108 + zero-mode tolerance while being far ABOVE the bare backward error + ``eps * ||L||``. Such a mode is neither machine-zero nor resolved -- + it lives inside the safety factor the tolerance carries, so the + arithmetic cannot say whether it is a genuine slow mode or noise. + When any exist, D1 would silently discard real physics and report the + next surviving (fast) mode instead. + """ return self.applicable and self.certified and self.ambiguous_count == 0 def as_dict(self) -> dict[str, object]: @@ -132,7 +204,47 @@ def operator_zero_tolerance( rtol: float = ZERO_MODE_EPS_FACTOR, name: str = "L_super", ) -> float: - """Operator-derived zero-mode tolerance ``rtol * eps * ||L||_2``.""" + """Operator-derived zero-mode tolerance ``rtol * eps * ||L||_2``. + + Round-13 review (issue #112 follow-up): certification and downstream + filtering must use ONE scale. The certificate accepts a stationary + residual up to the true eigensolver backward error ``eps * ||L||_2``, + but :func:`liouscope.numerics.scale.spectral_zero_tolerance` -- which + only sees the spectrum -- can only use the spectral radius + ``max|lambda|`` as a proxy. For a strongly non-normal generator + ``||L||_2`` exceeds the radius by orders of magnitude (measured: 3.9e3 + on a 4x4 trace-preserving example), so a certified-resolved zero mode + with residual between the two thresholds survived the radius filter and + D1 reported ``~1e-12`` -- occasionally NEGATIVE -- instead of the true + gap ``1``. + + Every consumer that holds the OPERATOR (not just its spectrum) must + therefore filter on this scale -- but a consumer holding a + :class:`ZeroModeCertificate` must take its cutoff from + :meth:`ZeroModeCertificate.zero_set_tolerance` and NOT from the + certificate's ``bound``. Round-17 review (PR #121): the two WERE the + same expression when this paragraph was written, and stopped being so + with the a posteriori refinement (issue #113 second axis), which lowers + the applied tolerance below ``bound`` precisely when a genuine slow mode + was rescued from the band. Filtering by ``bound`` there discards that + mode again, one layer after the certificate saved it -- four call sites + did exactly that, following this paragraph. The radius-based + :func:`~liouscope.numerics.scale.spectral_zero_tolerance` remains the + correct fallback for spectrum-only call sites. Genuine slow modes + falling inside this coarser band are not silently swallowed: they are + flagged by the certificate's ambiguity split (``resolved=False``, + issue #113), which floors the verdict instead. + + Parameters + ---------- + L_super + The operator whose eigensolve backward error sets the scale. + rtol + Multiplier on ``eps * ||L||_2``, shared with the certificate and + with ``spectral_zero_tolerance``. + name + Argument name for the fail-closed validation message. + """ if not np.isfinite(rtol) or rtol < 0.0: raise ValueError(f"rtol must be finite and non-negative, got {rtol}") L_c = np.asarray(require_finite_square_2d(L_super, name=name), dtype=complex) @@ -140,69 +252,65 @@ def operator_zero_tolerance( return float(rtol * float(np.finfo(float).eps) * norm2) -def _hilbert_dimension_from_superoperator(L_super: np.ndarray) -> int | None: - """Return d for a d^2 x d^2 superoperator, else ``None``.""" - arr = np.asarray(L_super) - if arr.ndim != 2 or arr.shape[0] != arr.shape[1]: - return None - n = int(arr.shape[0]) - d = int(round(np.sqrt(n))) - return d if d * d == n else None - - def trace_preservation_defect(L_super: np.ndarray) -> tuple[float, float]: - """Return robust ``(||vec(I)^H L||_2, ||L||_F)``. - - Both quantities are measurements and therefore use the same xLASSQ-style - scaled Euclidean norm. They remain finite whenever the true float64 norm is - representable; policy is handled separately by - :func:`trace_preservation_equation_defect` and the certificate gates. + """Return robust ``(||vec(I)^H L||_2, ||L||_F)`` for a generator. + + Every GKSL generator is trace preserving, which in column-stacking + convention is the exact statement ``vec(I)^H L = 0``. The quotient of the + two returned numbers is therefore a dimensionless, unit-invariant measure + of how far ``L`` is from being a legal generator. + + Issue #130: both quantities are evidence, so their arithmetic must be + correct independently of the policy that consumes them. Ordinary + ``np.linalg.norm`` squares before summing and can therefore return the + false sentinels ``0.0`` (underflow) or ``inf`` (overflow) even when the + mathematical norm of finite input is representable. Both norms now use + :func:`scaled_euclidean_norm`, an xLASSQ-style power-of-two scaled + sum-of-squares primitive. It returns a finite norm whenever the true + float64 result is representable and returns infinity only when the input or + the true norm is genuinely non-representable. + + This deliberately separates measurement from certificate policy. A very + large but finite operator must not be rejected merely because a naive norm + overflowed; structural guards such as ``band_discriminates`` decide whether + the resulting zero-mode evidence actually discriminates. """ L_super = np.asarray(L_super) - d = _hilbert_dimension_from_superoperator(L_super) - if d is None: + n = L_super.shape[0] + d = int(round(np.sqrt(n))) + if d * d != n: return float("nan"), float("nan") vec_i = np.eye(d, dtype=complex).reshape(-1, order="F") defect_vector = vec_i.conj() @ L_super return scaled_euclidean_norm(defect_vector), scaled_euclidean_norm(L_super) -def trace_preservation_equation_defect(L_super: np.ndarray) -> float: - """Return the worst componentwise relative defect of ``vec(I)^H L = 0``. - - A global quotient ``||vec(I)^H L|| / ||L||`` is useful as a measurement, - but it is too weak as the precondition for an exact zero-mode theorem: an - unrelated entry of size ``1e308`` can make an order-one trace defect look - relatively tiny. Each column of ``vec(I)^H L = 0`` is therefore judged - against only the terms that actually participate in that scalar equation: - - ``eta_j = |sum_i L_(ii),j| / sum_i |L_(ii),j|``. - - The returned value is ``max_j eta_j``. This is a componentwise relative - backward-error criterion, analogous to LAPACK's BERR philosophy. The ratio - is evaluated with power-of-two scaling, so it is invariant to representable - rate-unit rescaling and remains defined in the subnormal range. - """ - L_super = np.asarray(L_super) - d = _hilbert_dimension_from_superoperator(L_super) - if d is None: - return float("nan") - if not np.all(np.isfinite(L_super)): - return float("nan") - identity_rows = np.arange(d, dtype=int) * (d + 1) - worst = 0.0 - for j in range(L_super.shape[1]): - ratio = scaled_cancellation_ratio(L_super[identity_rows, j]) - if not np.isfinite(ratio): - return float("nan") - worst = max(worst, ratio) - return float(worst) - - def band_discriminates( magnitudes: np.ndarray, zero_count: int, bound: float ) -> bool: - """Return whether the zero-mode acceptance band rejected anything.""" + """Did the zero-mode band decide anything, or accept everything? + + Round-20 review (PR #121). ``residual <= bound`` certifies that a zero + mode was found. It says nothing about whether the band SEPARATED that mode + from the rest, and on the reviewer's counterexample it did not: a finite + 4x4 with cancelling ``+-1e308`` entries has ``||L||_2 = 1.4e308``, hence a + band of ``3.1e295``, and all four eigenvalues ``{1, 2, 3, 4}`` fall inside + it. The certificate came back ``certified=True, resolved=True`` with + ``zero_mode_count = 4`` -- a claim that the ENTIRE space is stationary, + for a generator whose spectrum is nothing of the kind. + + An acceptance region that contains every measured value has discriminated + nothing; passing it is not evidence. Note this is NOT a tightened + threshold -- tightening ``rtol`` would move the same band, not repair its + logic -- but a second, structural axis: the band must have had the + OPPORTUNITY to reject. + + The one legitimate way for a band to contain the whole spectrum is + ``bound == 0``: the exactly-zero generator, whose eigenvalues are exactly + zero and for which "everything is stationary" is the correct physics + measured exactly. A band of width zero accepts only exact zeros, so it + discriminates by construction and is excluded here. + """ if bound <= 0.0: return True return zero_count < int(magnitudes.size) @@ -217,15 +325,53 @@ def certified_nonzero_modes( left_vectors: np.ndarray | None = None, margin: float = ZERO_MODE_APOSTERIORI_MARGIN, ) -> np.ndarray: - """Identify in-band modes that a posteriori residuals prove non-stationary.""" + """Modes inside the zero band that are provably NOT the stationary mode. + + The magnitude band ``|lambda| <= rtol * eps * ||L||`` is a worst-case + backward-error argument: it assumes the eigensolver was as inaccurate as it + is allowed to be. On a well-conditioned generator it usually was not, and + the assumption is what discards genuine slow modes -- the + ``ZERO_MODE_AMBIGUITY_FACTOR`` comment records that no threshold on + ``|lambda|`` alone can separate the two populations, which are measured + only a factor of 2 apart. + + This is the second axis. For each in-band mode we compute the *a + posteriori* bound from the residuals actually attained, + + |lambda - lambda_hat| <~ max(||L x - lambda_hat x||, + ||L^H y - conj(lambda_hat) y||) / |y^H x| + + with ``||x|| = ||y|| = 1``. Were ``lambda`` exactly zero, that bound would + force ``|lambda_hat| <= bound``; so ``|lambda_hat| > margin * bound`` + certifies a genuine mode instead of merely scoring one. ``margin`` covers + the first-order caveat only -- see ``ZERO_MODE_APOSTERIORI_MARGIN`` for the + corpus measurement behind its value. + + Deliberately one-directional: the result can only move a mode OUT of the + zero set, never into it, so it cannot manufacture the false "unresolved" + verdict that the ambiguity split exists to avoid. + + Two invariants are enforced fail-closed: + + * fewer than two in-band modes -> nothing is certified. A trace-preserving + generator has at least one stationary mode, so a lone in-band mode is the + expected one and there is no decision to make. + * the smallest-magnitude in-band mode is never certified. Trace + preservation guarantees a zero eigenvalue exists; emptying the band would + contradict the very precondition under which this certificate applies. + """ idx = np.flatnonzero(in_band) out = np.zeros(eigenvalues.shape, dtype=bool) if idx.size < 2: return out + # The stationary mode itself is never a candidate (see docstring). idx = idx[idx != idx[int(np.argmin(np.abs(eigenvalues[idx])))]] vr, vl, ref = right_vectors, left_vectors, eigenvalues if vr is None or vl is None: + # Own decomposition: the certificate is a property of the OPERATOR, so + # it is legitimate to certify a candidate ladder's spectrum from the + # operator's own eigenpairs -- but only where the two agree. try: ref, vl, vr = sla.eig(L_c, left=True, right=True) except (ValueError, sla.LinAlgError): @@ -236,13 +382,15 @@ def certified_nonzero_modes( lam = eigenvalues[i] j = int(np.argmin(np.abs(ref - lam))) mag = float(abs(ref[j])) + # Fail closed when the borrowed decomposition disagrees about this mode: + # a bound computed for a different eigenvalue certifies nothing. if abs(ref[j] - lam) > 0.1 * max(abs(lam), mag) or mag == 0.0: continue x = vr[:, j] / max(float(np.linalg.norm(vr[:, j])), tiny) y = vl[:, j] / max(float(np.linalg.norm(vl[:, j])), tiny) sep = float(abs(np.vdot(y, x))) if sep <= 0.0: - continue + continue # defective pair: the first-order bound does not apply res = max( float(np.linalg.norm(L_c @ x - ref[j] * x)), float(np.linalg.norm(L_c.conj().T @ y - np.conj(ref[j]) * y)), @@ -261,7 +409,22 @@ def refine_zero_band( right_vectors: np.ndarray | None = None, left_vectors: np.ndarray | None = None, ) -> tuple[np.ndarray, float]: - """Return refined zero membership and its corresponding filter tolerance.""" + """Zero-mode membership plus the tolerance the gap filters must then apply. + + Returns ``(in_band, zero_tolerance)``. Without a certified non-stationary + mode both are exactly the pre-refinement values, so the healthy path is + unchanged bit for bit. + + The tolerance has to travel with the membership: D1/D3/D4 filter by + magnitude, so a mode this certificate rescues would be discarded again one + layer later if they kept reading the raw ``bound``. It is placed at the + geometric mean of the two populations -- the scale-free midpoint, which is + the right notion of "between" for quantities spanning decades. + + Fail-closed on inversion: if some mode kept as zero is LARGER than the + smallest certified one, no single threshold can express the split, and the + refinement is abandoned wholesale rather than applied half-way. + """ nonzero = certified_nonzero_modes( L_c, eigenvalues, @@ -278,41 +441,117 @@ def refine_zero_band( return magnitudes <= bound, bound if lo <= 0.0: return in_band, hi * 0.5 + # ROUND-22 REVIEW (PR #121). ``sqrt(lo * hi)`` forms the PRODUCT first, + # and that product underflows to exactly 0 whenever ``lo * hi`` drops + # below ~5e-324 -- while ``lo``, ``hi`` and their geometric mean are all + # perfectly representable. The tolerance then becomes 0.0, the strict + # ``|lambda| > tol`` filters keep the numerical stationary residual as a + # physical mode, and D1 reports it as the gap. Since ``lo * hi`` scales as + # c^2 under a uniform rate rescale ``L -> cL`` while the refinement itself + # is scale free, the failure is purely a units artefact: the same physics + # in smaller rate units silently changes the answer. + # + # ``sqrt(lo) * sqrt(hi)`` is the same number mathematically and forms each + # factor separately, so it can neither underflow (sqrt of the smallest + # subnormal is ~2.2e-162) nor overflow (sqrt of the largest float is + # ~1.3e154, and the product of two such factors is at most the largest + # float again). return in_band, float(np.sqrt(lo) * np.sqrt(hi)) -def _tp_certificate_applicable( - L_c: np.ndarray, - tp_defect: float, - fro: float, - tp_rtol: float, -) -> bool: - """Return whether trace preservation is established strongly enough. - - The historical normwise criterion is retained for compatibility and broad - scale checking. The componentwise equation criterion is the load-bearing - addition from #130: it prevents unrelated large entries from diluting a - local violation of the exact trace equations. - """ - equation_defect = trace_preservation_equation_defect(L_c) - if not np.isfinite(tp_defect) or not np.isfinite(fro): - return False - if not np.isfinite(equation_defect): - return False - if equation_defect > tp_rtol: - return False - if tp_defect > tp_rtol * fro: - return False - return True - - def certified_eigvals( L_super: np.ndarray, *, rtol: float = ZERO_MODE_EPS_FACTOR, tp_rtol: float = 1.0e-10, ) -> tuple[np.ndarray, ZeroModeCertificate]: - """Eigenvalues checked against the exact trace-preserving zero-mode fact.""" + """Eigenvalues of a Liouvillian, checked against an exact structural fact. + + Issue #112. For a trace-preserving generator ``vec(I)^H L = 0`` holds + *exactly*, so ``0`` is an exact eigenvalue and a correct eigensolve must + return some ``|lambda| <= rtol * eps * ||L||``. This is a **theorem, not a + threshold**: when the computed spectrum contains no such eigenvalue the + solve has demonstrably failed, and no choice of zero-mode tolerance + (relative or absolute, issue #108) can repair it -- the missing information + is not in the eigenvalues. + + Why this is needed + ------------------ + LAPACK's complex non-symmetric driver (``zgeev``) deflates a subdiagonal + when the Ahues-Tisseur test ``|h[i,i-1] h[i-1,i]| <= eps |h[i,i]| + |h[i-1,i-1] - h[i,i]|`` passes. On a STIFF generator the large diagonal + entries inflate that right-hand side: with a fast rate ``2.7e5`` alongside + slow rates ``~1e-5``, the bound reaches ``~1.6e-5`` and the entire slow + block deflates to its own DIAGONAL. Measured on a four-level classical jump + network, the exact zero mode vanished from the returned spectrum and a + spurious mode appeared in its place, so ``D1`` reported ``7.28e-6`` for a + generator whose true gap is ``1.07e-5``. + + That failure is invisible to conditioning: the per-eigenvalue condition + numbers of the wrong spectrum are ``1``-``25``, i.e. the solver reports the + wrong answer as well conditioned. Balancing does not fix it either + (measured). What does fix it is solving the same matrix by a route whose + deflation is not driven by the stiff diagonal. + + Repair ladder + ------------- + Candidates are tried in order and the FIRST one satisfying the certificate + WITHOUT ambiguous in-band modes is returned; the primary solve is + therefore byte-identical whenever it is already correct (measured: + 400/400 random four-level networks, stiff and well-conditioned alike, + keep the ``zgeev`` result). An ambiguous candidate does not end the + ladder (round-15 review): a later route can resolve the very modes it + cannot -- measured, ``zgeev`` certified a stiff network with one + ambiguous mode at ``6.3e-7`` while ``dgeev-real`` returned a clean + ``2.1e-15`` stationary mode. The best ambiguous candidate (fewest + ambiguous modes) is returned, fail-closed, only when every route stays + ambiguous. + + 1. ``zgeev`` -- the incumbent, via :func:`eig_nonhermitian`; + 2. ``dgeev`` on the real part, when ``L`` is real-valued (a different + LAPACK driver with a different deflation path); + 3. the complex Schur form ``zgees``, which deflates on the full Hessenberg + structure rather than eigenvalue-by-eigenvalue; + 4. an explicitly balanced ``zgeev``. + + If no candidate is certified, the eigenvalues of the best (smallest + residual) candidate are returned together with ``certified=False`` -- the + caller decides how loudly to fail. The function never silently substitutes + an uncertified spectrum for a certified one. + + Parameters + ---------- + L_super + ``d^2 x d^2`` vectorised generator. + rtol + Multiplier on the backward error ``eps * ||L||_2``, shared with + :func:`liouscope.numerics.scale.spectral_zero_tolerance`. + tp_rtol + Relative tolerance deciding whether ``L`` is trace preserving at all. + A non-trace-preserving operator has no guaranteed zero mode, so the + certificate is reported as ``applicable=False`` and the incumbent + spectrum is returned untouched. + + Returns + ------- + tuple + ``(eigenvalues, certificate)``. + """ + # ROUND-22 REVIEW (PR #121). A tolerance that is not a number cannot + # decide anything. With ``tp_rtol = NaN`` the trace-preservation test + # ``defect > tp_rtol * fro`` is False for EVERY operator (all NaN + # comparisons are), and with ``tp_rtol = inf`` its right-hand side is + # infinite -- so a demonstrably non-trace-preserving generator walks + # straight past the applicability gate and can come back + # ``applicable=True, certified=True``. Measured before this guard: + # ``certified_eigvals(diag([0,-1,-2,-3]), tp_rtol=nan)`` reported + # applicable and certified with a trace defect of 3.0. + # + # Round 21 closed exactly this shape one variable further along (a + # non-finite REFERENCE scale ``fro``) and left the PARAMETER beside it + # unchecked. The generalisation is the fix: every tolerance this layer + # compares against is validated the way ``operator_zero_tolerance`` and + # ``spectral_zero_tolerance`` already validate theirs. if not np.isfinite(rtol) or rtol < 0.0: raise ValueError(f"rtol must be finite and non-negative, got {rtol}") if not np.isfinite(tp_rtol) or tp_rtol < 0.0: @@ -323,6 +562,13 @@ def certified_eigvals( norm2 = float(np.linalg.norm(L_c, 2)) if L_c.size else 0.0 bound = rtol * float(np.finfo(float).eps) * norm2 + # Round-17 review (PR #121). The primary solve must not be able to end the + # ladder. ``zgeev`` raising ``LinAlgError`` IS a nonconvergence -- exactly + # the failure the real-driver / Schur / balanced routes exist to repair -- + # and letting the exception propagate defeated the ladder precisely in the + # case it was built for. The error is carried instead and re-raised only + # if NO route produced a spectrum at all, so the caller still sees the + # original diagnosis when nothing worked. primary: np.ndarray | None primary_error: sla.LinAlgError | None try: @@ -330,11 +576,44 @@ def certified_eigvals( primary_error = None except sla.LinAlgError as exc: primary, primary_error = None, exc - - if not _tp_certificate_applicable(L_c, tp_defect, fro, tp_rtol): + if ( + not np.isfinite(tp_defect) + # ROUND-21: a non-finite REFERENCE scale is as unusable as a + # non-finite defect. ``defect <= tp_rtol * fro`` against an + # infinite ``fro`` admits every operator, so the band that + # follows is drawn from a scale that was never measured; which + # modes it then swallows is decided by the LAPACK build (2 of 4 + # on the CI runners, 4 of 4 locally). Trace preservation + # RELATIVE to the operator's own scale is not a statement that + # can be made about infinity, so it is refused here. + or not np.isfinite(fro) + # ROUND-23 REVIEW (PR #121), same finding one line further down. + # ``max(fro, tiny)`` floors the REFERENCE SCALE at a constant that + # is larger than the operator whenever the operator is subnormal, + # so the relative test stops being relative: measured for + # ``diag([0, -1e-320, -2e-320, -3e-320])``, defect 3e-320 against a + # floor-derived bound of 2.2e-318, hence 'trace preserving' for an + # operator whose relative defect is 0.80. Repairing the norms above + # without this leaves the rate-unit dependence intact below ~1e-317. + # The floor protected nothing: ``defect <= sqrt(d) * fro`` holds, so + # ``fro == 0`` forces ``defect == 0`` and the comparison ``0 > 0`` is + # False either way -- the exactly-zero generator stays applicable. + or tp_defect > tp_rtol * fro + ): + # No certificate applies without trace preservation, so there is + # nothing to repair TOWARDS and the ladder is not run: the primary + # failure is the honest answer. if primary is None: assert primary_error is not None raise primary_error + # ROUND-22 REVIEW (PR #121). ``zero_mode_count`` kept its dataclass + # default of 1 here, which asserts a stationary mode that nothing + # counted and nothing guarantees -- the whole point of + # ``applicable=False`` is that no zero eigenvalue is implied. Both + # certificate APIs reported ``zero_mode_count: 1`` for + # ``diag([1,2,3,4])``, a spectrum containing no zero at all, and the + # dict is persisted as audit metadata. 0 is the honest count: nothing + # was certified as stationary. return primary, ZeroModeCertificate( applicable=False, certified=False, @@ -346,9 +625,31 @@ def certified_eigvals( ) def _candidates() -> Iterator[tuple[str, np.ndarray]]: + # LAZY by design (twelfth-round review): the ladder stops at the first + # certified spectrum, and the healthy path is the primary solve -- an + # eager list would charge every caller three additional cubic + # decompositions that the loop then never looks at. if primary is not None: yield ("zgeev", primary) + # ``dgeev`` genuinely solves a DIFFERENT matrix unless L is EXACTLY + # real: ``np.allclose`` carries an absolute default atol, so a stiff + # complex generator in small rate units read as "real" and the repair + # silently deleted its Hamiltonian part (measured: D3 moved from 1e-7 + # to ~9.7e-6 under a pure L -> 1e-10 L rescale). Exact realness is the + # only scale-invariant criterion under which the route is valid. if not np.any(L_c.imag): + # GUARDED (round-16 review): an exactly real generator whose + # ``dgeev`` solve raises must not end the ladder. Without this + # suppression the exception propagates out of the generator and + # the two remaining repair routes below -- ``zgees-schur`` and + # ``balanced-zgeev`` -- never run, even when either could recover + # a certifiable spectrum. The route is a REPAIR attempt like the + # others; a repair that fails is a route that did not work, not a + # reason to abandon the remaining ones. + # + # The sibling ladder in ``certified_eig`` already guarded exactly + # this step. The omission here was the inconsistency, not the + # guard there. with contextlib.suppress(ValueError, sla.LinAlgError): yield ("dgeev-real", np.linalg.eigvals(L_c.real).astype(complex)) with contextlib.suppress(ValueError, sla.LinAlgError): @@ -357,6 +658,28 @@ def _candidates() -> Iterator[tuple[str, np.ndarray]]: balanced = sla.matrix_balance(L_c, permute=True)[0] yield ("balanced-zgeev", np.linalg.eigvals(balanced)) + # Split the zero-tolerance band (issue #113). Below the split a mode is + # machine-zero; above it -- but still inside the #108 tolerance, which + # carries a factor ``rtol`` on top -- it is neither machine-zero nor + # resolved. + # + # The factor is set from a MEASURED distribution, not from taste. Across + # 83 healthy generators (single/two-qubit families plus random GKSL, + # unique and degenerate steady states alike) the largest in-band + # |lambda| reaches 2.38 * eps*||L||, median 0.39. Unresolved slow modes + # were measured at 4.87 and 4.87e2. The two populations are therefore + # NOT cleanly separated -- the nearest pair is 2.38 against 4.87, about + # 2x -- so the split is placed at 30x, roughly an order of magnitude + # clear of the healthy maximum rather than midway between the two. + # + # That choice is deliberately conservative: a false NaN on a healthy + # generator destroys a correct analysis, while a missed marginal case + # leaves behaviour exactly as it was before this check existed. The cost + # is bounded reach -- above a spectral spread of ~1e14 the slow modes + # sink below the split and the defect is undetectable by any magnitude + # test, because double precision has genuinely lost the information. + # That residual is documented in issue #113 rather than papered over + # here. split = ZERO_MODE_AMBIGUITY_FACTOR * float(np.finfo(float).eps) * norm2 best: tuple[str, np.ndarray, float] | None = None ambiguous_best: tuple[str, np.ndarray, float, int, int, float] | None = None @@ -364,10 +687,21 @@ def _candidates() -> Iterator[tuple[str, np.ndarray]]: magnitudes = np.abs(ev) residual = float(np.min(magnitudes)) if ev.size else float("inf") if residual <= bound: - in_band, zero_tolerance = refine_zero_band(L_c, ev, magnitudes, bound) + # Issue #113 second axis: a mode the a posteriori bound certifies + # as non-stationary leaves the zero set before it is counted, so a + # genuine slow mode no longer masquerades as a second zero mode and + # hands D1 the NEXT eigenvalue as the gap. + in_band, zero_tolerance = refine_zero_band( + L_c, ev, magnitudes, bound + ) zero_count = int(np.count_nonzero(in_band)) ambiguous = int(np.count_nonzero(in_band & (magnitudes > split))) if not band_discriminates(magnitudes, zero_count, bound): + # ROUND-20: the band swallowed the whole spectrum, so this + # route proved nothing. Carried on as an UNcertified + # candidate, exactly like a route whose residual missed the + # band -- a later route on a different deflation path may + # still produce a spectrum the band can separate. if best is None or residual < best[2]: best = (name, ev, residual) continue @@ -382,14 +716,16 @@ def _candidates() -> Iterator[tuple[str, np.ndarray]]: zero_mode_count=zero_count, zero_tolerance=zero_tolerance, ) + # Round-15 review: an AMBIGUOUS candidate must not end the + # ladder -- a later route can resolve the very modes this one + # cannot (measured: zgeev certified a stiff network with one + # ambiguous mode at 6.3e-7 while dgeev-real returns a clean + # 2.1e-15 stationary mode). It is retained only as the + # fail-closed fallback when EVERY route stays ambiguous, keyed + # by the fewest ambiguous modes (ties: ladder order). if ambiguous_best is None or ambiguous < ambiguous_best[4]: ambiguous_best = ( - name, - ev, - residual, - zero_count, - ambiguous, - zero_tolerance, + name, ev, residual, zero_count, ambiguous, zero_tolerance ) continue if best is None or residual < best[2]: @@ -409,6 +745,7 @@ def _candidates() -> Iterator[tuple[str, np.ndarray]]: zero_tolerance=zero_tolerance, ) if best is None: + # Every route raised (only reachable when the primary solve did too). assert primary_error is not None raise primary_error return best[1], ZeroModeCertificate( @@ -428,7 +765,64 @@ def certified_eig( rtol: float = ZERO_MODE_EPS_FACTOR, tp_rtol: float = 1.0e-10, ) -> tuple[EigenDecomposition, ZeroModeCertificate]: - """Like :func:`certified_eigvals`, also returning left/right eigenvectors.""" + """Like :func:`certified_eigvals`, but also returns left/right eigenvectors. + + Issue #112 follow-up. The eigenvalue-only certificate protects D1/D3/D4, but + the layers that consume eigen*vectors* -- D19's slowest-mode overlap above + all -- were still reading the uncertified decomposition. That matters more + than it sounds: on a stiff generator ``zgeev`` can return a *spurious* slow + mode, and D19 then measures the initial state's overlap with an eigenvector + that corresponds to no physical mode. A wrong overlap on that rung fires a + false A11/F4 Mpemba candidate -- the highest-priority branch of the + classifier, and exactly the failure mode #108 was fixed to prevent. + + Repair ladder, narrower than the eigenvalue-only one + ---------------------------------------------------- + Only routes that actually produce eigenvectors in the ORIGINAL basis are + usable here: + + 1. ``zgeev`` -- the incumbent; + 2. ``dgeev`` on the real part, when ``L`` is real-valued. + + The Schur and balanced routes are deliberately absent: ``zgees`` yields no + eigenvectors directly, and the balanced solve returns them in the balanced + basis, where back-transforming re-introduces exactly the scaling that made + the decomposition unreliable. A narrower ladder is the honest choice -- + it means some spectra that ``certified_eigvals`` can repair are NOT + repairable here, and the certificate then reports ``certified=False`` + rather than pretending otherwise. + + Beyond the eigenvalue certificate, a candidate is accepted only when + every consumed mode's left AND right eigen*vector* residuals pass the + per-mode relative gate ``r_j <= max(VECTOR_RESIDUAL_REL_MAX * |lambda_j|, + bound)`` (round-15 review) -- the vectors are what D19 and the Petermann + factors consume, and a small ``|lambda|`` does not vouch for them. When + the vectors fail on every route, the returned certificate carries + ``certified=False`` with ``residual`` set to the offending (vector) + residual, so the downstream warning shows a number that actually + exceeds the printed bound. Ambiguous candidates likewise no longer end + the ladder (same rule as :func:`certified_eigvals`). + + Returns + ------- + tuple + ``(EigenDecomposition with left_vectors set, certificate)``. + """ + # ROUND-22 REVIEW (PR #121). A tolerance that is not a number cannot + # decide anything. With ``tp_rtol = NaN`` the trace-preservation test + # ``defect > tp_rtol * fro`` is False for EVERY operator (all NaN + # comparisons are), and with ``tp_rtol = inf`` its right-hand side is + # infinite -- so a demonstrably non-trace-preserving generator walks + # straight past the applicability gate and can come back + # ``applicable=True, certified=True``. Measured before this guard: + # ``certified_eigvals(diag([0,-1,-2,-3]), tp_rtol=nan)`` reported + # applicable and certified with a trace defect of 3.0. + # + # Round 21 closed exactly this shape one variable further along (a + # non-finite REFERENCE scale ``fro``) and left the PARAMETER beside it + # unchecked. The generalisation is the fix: every tolerance this layer + # compares against is validated the way ``operator_zero_tolerance`` and + # ``spectral_zero_tolerance`` already validate theirs. if not np.isfinite(rtol) or rtol < 0.0: raise ValueError(f"rtol must be finite and non-negative, got {rtol}") if not np.isfinite(tp_rtol) or tp_rtol < 0.0: @@ -439,6 +833,13 @@ def certified_eig( norm2 = float(np.linalg.norm(L_c, 2)) if L_c.size else 0.0 bound = rtol * float(np.finfo(float).eps) * norm2 + # Round-17 review (PR #121). The primary solve must not be able to end the + # ladder. ``zgeev`` raising ``LinAlgError`` IS a nonconvergence -- exactly + # the failure the real-driver / Schur / balanced routes exist to repair -- + # and letting the exception propagate defeated the ladder precisely in the + # case it was built for. The error is carried instead and re-raised only + # if NO route produced a spectrum at all, so the caller still sees the + # original diagnosis when nothing worked. primary: EigenDecomposition | None primary_error: sla.LinAlgError | None try: @@ -446,11 +847,41 @@ def certified_eig( primary_error = None except sla.LinAlgError as exc: primary, primary_error = None, exc - - if not _tp_certificate_applicable(L_c, tp_defect, fro, tp_rtol): + if ( + not np.isfinite(tp_defect) + # ROUND-21: a non-finite REFERENCE scale is as unusable as a + # non-finite defect. ``defect <= tp_rtol * fro`` against an + # infinite ``fro`` admits every operator, so the band that + # follows is drawn from a scale that was never measured; which + # modes it then swallows is decided by the LAPACK build (2 of 4 + # on the CI runners, 4 of 4 locally). Trace preservation + # RELATIVE to the operator's own scale is not a statement that + # can be made about infinity, so it is refused here. + or not np.isfinite(fro) + # ROUND-23 REVIEW (PR #121), same finding one line further down. + # ``max(fro, tiny)`` floors the REFERENCE SCALE at a constant that + # is larger than the operator whenever the operator is subnormal, + # so the relative test stops being relative: measured for + # ``diag([0, -1e-320, -2e-320, -3e-320])``, defect 3e-320 against a + # floor-derived bound of 2.2e-318, hence 'trace preserving' for an + # operator whose relative defect is 0.80. Repairing the norms above + # without this leaves the rate-unit dependence intact below ~1e-317. + # The floor protected nothing: ``defect <= sqrt(d) * fro`` holds, so + # ``fro == 0`` forces ``defect == 0`` and the comparison ``0 > 0`` is + # False either way -- the exactly-zero generator stays applicable. + or tp_defect > tp_rtol * fro + ): if primary is None: assert primary_error is not None raise primary_error + # ROUND-22 REVIEW (PR #121). ``zero_mode_count`` kept its dataclass + # default of 1 here, which asserts a stationary mode that nothing + # counted and nothing guarantees -- the whole point of + # ``applicable=False`` is that no zero eigenvalue is implied. Both + # certificate APIs reported ``zero_mode_count: 1`` for + # ``diag([1,2,3,4])``, a spectrum containing no zero at all, and the + # dict is persisted as audit metadata. 0 is the honest count: nothing + # was certified as stationary. return primary, ZeroModeCertificate( applicable=False, certified=False, @@ -466,6 +897,7 @@ def certified_eig( ) def _candidates() -> Iterator[tuple[str, EigenDecomposition]]: + # Lazy + exact realness, for the same reasons as the eigenvalue ladder. if primary is not None: yield ("zgeev", primary) if not np.any(L_c.imag): @@ -481,12 +913,36 @@ def _candidates() -> Iterator[tuple[str, EigenDecomposition]]: ) def _vector_residual(decomp: EigenDecomposition) -> float: + """Largest OFFENDING per-mode eigenvector residual; 0.0 when all pass. + + Round-15 review: the eigenvalue certificate alone does not guarantee + the eigen*vectors* -- the quantities D19 and the Petermann factors + actually consume. Measured on a valid stiff classical network: + certified and resolved with ``bound = 3.2e-7``, yet the LEFT vector + of the selected slow mode had residual ``3.1e-4`` -- three orders of + magnitude beyond the certificate, i.e. not an eigenvector of ``L`` + in any usable sense (the right vectors were fine). + + The acceptance test is PER MODE and RELATIVE: + ``r_j <= max(VECTOR_RESIDUAL_REL_MAX * |lambda_j|, bound)`` for every + non-zero-band mode, with ``r_j`` the larger of the unit-normalised + left/right residuals. ``r_j / |lambda_j|`` is the first-order + relative error scale of anything computed from the pair, and unlike + a single operator-scale cutoff it separates the measured + populations: a legitimate ``dgeev-real`` repair carries a few + percent on its slow modes, while a corrupt ``zgeev`` decomposition + carries 22% to 2900% (see + :data:`liouscope._consts.VECTOR_RESIDUAL_REL_MAX` for the measured + calibration). Modes inside the zero band are not gated here: their + vectors span the stationary manifold and are not consumed through + this decomposition. + """ ev = decomp.eigenvalues if not ev.size: return 0.0 vr = decomp.right_vectors vl = decomp.left_vectors - assert vl is not None + assert vl is not None # every ladder route computes left vectors tiny = np.finfo(float).tiny res_r = np.linalg.norm(L_c @ vr - vr * ev[None, :], axis=0) / np.maximum( np.linalg.norm(vr, axis=0), tiny @@ -509,11 +965,14 @@ def _vector_residual(decomp: EigenDecomposition) -> float: tuple[str, EigenDecomposition, float, int, int, float] | None ) = None vector_failed: tuple[str, EigenDecomposition, float, float] | None = None - for name, decomp in _candidates(): magnitudes = np.abs(decomp.eigenvalues) residual = float(np.min(magnitudes)) if magnitudes.size else float("inf") if residual <= bound: + # Issue #113 second axis: a mode the a posteriori bound certifies + # as non-stationary leaves the zero set before it is counted, so a + # genuine slow mode no longer masquerades as a second zero mode and + # hands D1 the NEXT eigenvalue as the gap. in_band, zero_tolerance = refine_zero_band( L_c, decomp.eigenvalues, @@ -525,6 +984,9 @@ def _vector_residual(decomp: EigenDecomposition) -> float: zero_count = int(np.count_nonzero(in_band)) ambiguous = int(np.count_nonzero(in_band & (magnitudes > split))) if not band_discriminates(magnitudes, zero_count, bound): + # ROUND-20, same guard as the sibling ladder. The two loops + # having drifted apart is what produced the round-16 finding; + # this one is added to both in the same commit. if best is None or residual < best[2]: best = (name, decomp, residual) continue @@ -541,24 +1003,25 @@ def _vector_residual(decomp: EigenDecomposition) -> float: zero_mode_count=zero_count, zero_tolerance=zero_tolerance, ) + # Eigenvalues fine, eigenvectors demonstrably not: keep for + # the fail-closed record and try the next route. if vector_failed is None or vec_residual < vector_failed[3]: vector_failed = (name, decomp, residual, vec_residual) continue + # Round-15 review: an ambiguous candidate must not end the + # ladder (same rule as certified_eigvals). if ambiguous_best is None or ambiguous < ambiguous_best[4]: ambiguous_best = ( - name, - decomp, - residual, - zero_count, - ambiguous, - zero_tolerance, + name, decomp, residual, zero_count, ambiguous, zero_tolerance ) continue if best is None or residual < best[2]: best = (name, decomp, residual) if ambiguous_best is not None: - name, decomp, residual, zero_count, ambiguous, zero_tolerance = ambiguous_best + name, decomp, residual, zero_count, ambiguous, zero_tolerance = ( + ambiguous_best + ) return decomp, ZeroModeCertificate( applicable=True, certified=True, @@ -571,6 +1034,10 @@ def _vector_residual(decomp: EigenDecomposition) -> float: zero_tolerance=zero_tolerance, ) if vector_failed is not None: + # No route produced consumable eigenvectors. certified=False keeps + # ``resolved`` False, so every consumer withholds -- the eigenvalue + # residual alone must not be allowed to read as success when the + # vectors it vouches for are not eigenvectors. name, decomp, residual, vec_residual = vector_failed return decomp, ZeroModeCertificate( applicable=True, @@ -582,6 +1049,7 @@ def _vector_residual(decomp: EigenDecomposition) -> float: zero_mode_count=0, ) if best is None: + # Every route raised (only reachable when the primary solve did too). assert primary_error is not None raise primary_error return best[1], ZeroModeCertificate( @@ -596,7 +1064,23 @@ def _vector_residual(decomp: EigenDecomposition) -> float: def require_finite_square_2d(A: np.ndarray, *, name: str = "matrix") -> np.ndarray: - """Validate ``A`` as a finite, square, 2-D array at a public boundary.""" + """Validate ``A`` as a finite, square, 2-D array at a public boundary. + + Returns ``np.asarray(A)`` unchanged when valid, else raises a + :class:`ValueError` that names the offending argument *and* the concrete + defect. This is a fail-closed guard for public entry points: without it a + ``NaN``/``inf``-laden Liouvillian flows into ``scipy.linalg.expm`` / ``svd`` + and surfaces as an opaque, location-blind error deep in LAPACK + (``"array must not contain infs or NaNs"`` / ``"SVD did not converge"``), + leaving the caller no clue which input was malformed. + + Parameters + ---------- + A + Candidate array. + name + Human-readable argument name used in error messages. + """ arr = np.asarray(A) if arr.ndim != 2 or arr.shape[0] != arr.shape[1]: raise ValueError(f"{name} must be a square 2-D array, got shape {arr.shape}") @@ -613,7 +1097,14 @@ def require_finite_square_2d(A: np.ndarray, *, name: str = "matrix") -> np.ndarr def hermiticity_defect(A: np.ndarray) -> tuple[float, float]: - """Return ``(max|A - A^H|, max|A|)``.""" + """Return ``(max|A - A^H|, max|A|)`` -- the absolute defect and its scale. + + Both are elementwise sup-norms, so their quotient is the *relative* + Hermiticity defect: a dimensionless number that is unchanged by a pure + change of units ``A -> cA``. Returning the pair (rather than the quotient) + keeps the zero-operator case in the caller's hands, where the right + semantics are known. + """ A = np.asarray(A) defect = float(np.max(np.abs(A - A.conj().T))) if A.size else 0.0 scale = float(np.max(np.abs(A))) if A.size else 0.0 @@ -626,7 +1117,55 @@ def is_hermitian( *, rtol: float = EPS_HERMITICITY, ) -> bool: - """Return True iff ``A`` is Hermitian, scale-relative by default.""" + """Return True iff ``A`` is Hermitian. + + The gate is SCALE-RELATIVE by default (issue #109): ``A`` is accepted iff + ``max|A - A^H| <= rtol * max|A|``. An operator that is not trace-normalised + -- a Hamiltonian, a jump operator -- carries a physical dimension, so an + ABSOLUTE tolerance makes the validation verdict depend on the caller's + choice of units, and it fails in both directions: + + * **fail-open** (the serious one) -- a genuinely non-Hermitian ``H`` with a + relative defect of ``1e-6`` passes a fixed ``1e-9`` gate as soon as + ``||H|| <~ 1e-3``. The resulting generator is not GKSL at all, and every + downstream diagnostic then describes dynamics that are not a quantum + channel; + * **false rejection** -- an exactly Hermitian ``H`` obtained from a + similarity transform carries round-off ``~eps*||H||``, which exceeds a + fixed ``1e-9`` once ``||H|| >~ 1e8``, so valid input is refused. + + ``rtol`` is deliberately ``EPS_HERMITICITY = 1e-9``, i.e. the historical + absolute value re-read as a relative one. At unit scale (``max|A| ~ 1``, + the regime of every existing fixture) the gate is therefore unchanged; only + operators far from unit scale move, which is exactly the defect. Note that + ``1e-9`` is many orders above round-off (``~1e-16``) *by design*: this is a + validation gate against malformed input, not a round-off discriminator. + + ``rtol=0`` semantics inside :func:`numpy.allclose` are not used here. + ``np.allclose``'s own ``rtol`` compares against ``|b|`` elementwise, which + for ``b = A^H`` makes the tolerance largest exactly where ``A`` is largest + and vanish where ``A`` is zero -- not a statement about the operator's + scale. The scale is taken once, from ``max|A|``. + + Parameters + ---------- + A + Candidate matrix. Non-square (or non-2-D) input returns ``False``. + atol + Legacy ABSOLUTE tolerance. When given, ``A`` is accepted iff + ``max|A - A^H| <= atol``, reproducing the pre-#109 behaviour verbatim. + This is the opt-in for callers whose operator is normalised to a known + scale -- density matrices, where the absolute reading is already + relative to ``tr(rho) = 1`` (see :func:`is_density_matrix`). + rtol + Relative tolerance on ``max|A|``; ignored when ``atol`` is given. + + Notes + ----- + Zero-operator semantics: ``max|A| == 0`` gives ``tol = 0``, so the all-zero + matrix (exactly Hermitian) passes and nothing else does. This mirrors the + documented zero-operator behaviour of :func:`liouscope.numerics.scale.rate_scale`. + """ A = np.asarray(A) if A.ndim != 2 or A.shape[0] != A.shape[1]: return False @@ -641,12 +1180,23 @@ def is_density_matrix( atol_trace: float = EPS_TRACE, atol_psd: float = 1.0e-10, ) -> bool: - """Return True iff ``rho`` is a valid density matrix.""" + """Return True iff ``rho`` is a valid density matrix. + + Checks Hermiticity, unit trace and positive semi-definiteness. + + The Hermiticity gate stays ABSOLUTE here (issue #109). A density matrix is + trace-normalised, so ``atol`` is already implicitly relative to a fixed + scale ``tr(rho) = 1``; the rate-dimension argument that makes an absolute + tolerance wrong for a Hamiltonian does not apply. Reading it as relative to + ``max|rho|`` would instead *tighten* the gate by up to a factor ``d`` for + nearly maximally mixed states -- a behaviour change with no defect behind it. + """ rho = np.asarray(rho) if not is_hermitian(rho, atol=1.0e-9): return False if abs(np.trace(rho) - 1.0) > atol_trace: return False + # Hermitian => use eigvalsh for stability. evals = np.linalg.eigvalsh((rho + rho.conj().T) / 2) return bool(evals.min() >= -atol_psd) @@ -657,7 +1207,14 @@ def support_check( *, eps: float = EPS_SUPP, ) -> tuple[bool, np.ndarray]: - """Verify ``supp(rho_initial) subset supp(rho_steady)``.""" + """Verify ``supp(rho_initial) subset supp(rho_steady)``. + + Returns a tuple ``(ok, regularised_steady)`` where ``regularised_steady`` + is ``rho_steady + eps * I / d`` shifted to avoid singularities when + computing the relative entropy ``D(rho || pi)``. + + Anchor J: this guards against NaNs when ``rho_steady`` is rank-deficient. + """ rho_initial = np.asarray(rho_initial) rho_steady = np.asarray(rho_steady) d = rho_steady.shape[0] @@ -668,6 +1225,7 @@ def support_check( return True, rho_steady null_space = evecs_steady[:, null_mask] + # Projector onto null space P_null = null_space @ null_space.conj().T leakage = float(np.real(np.trace(P_null @ rho_initial @ P_null))) ok = leakage < eps From ac2f74c365be802948d6c0025bd65a9e52398309 Mon Sep 17 00:00:00 2001 From: Marco Hostettler Date: Tue, 1 Sep 2026 21:04:35 +0200 Subject: [PATCH 10/22] fix(#130): add componentwise trace-preservation backward-error gate --- src/liouscope/numerics/linalg.py | 55 ++++++++++++++++++++++++++++++-- 1 file changed, 53 insertions(+), 2 deletions(-) diff --git a/src/liouscope/numerics/linalg.py b/src/liouscope/numerics/linalg.py index 72aa1892..01242e49 100644 --- a/src/liouscope/numerics/linalg.py +++ b/src/liouscope/numerics/linalg.py @@ -27,7 +27,7 @@ ZERO_MODE_APOSTERIORI_MARGIN, ZERO_MODE_EPS_FACTOR, ) -from .norms import scaled_euclidean_norm +from .norms import scaled_cancellation_ratio, scaled_euclidean_norm @dataclass(frozen=True, slots=True) @@ -285,6 +285,44 @@ def trace_preservation_defect(L_super: np.ndarray) -> tuple[float, float]: return scaled_euclidean_norm(defect_vector), scaled_euclidean_norm(L_super) +def trace_preservation_componentwise_error(L_super: np.ndarray) -> float: + """Return the worst componentwise backward error of ``vec(I)^H L = 0``. + + The global ratio ``||vec(I)^H L|| / ||L||`` is necessary but insufficient + as an applicability test: entries of ``L`` that do not participate in a + particular trace equation can make the global operator norm enormous and + thereby hide an order-one violation in another column. Trace preservation + is instead a family of scalar equations, one for every superoperator + column ``j``: + + ``sum_i L[(i, i), j] = 0``. + + For each equation we therefore compute + ``abs(sum terms) / sum(abs(terms))`` and return the largest ratio. This is + the componentwise relative backward-error reading: every equation is + normalised only by the represented coefficients that can actually cancel + in that equation. :func:`scaled_cancellation_ratio` evaluates the ratio in + a power-of-two scaled domain, so the verdict is stable for subnormal and + near-overflow rate units. An exact all-zero equation contributes zero. + """ + arr = np.asarray(L_super) + if arr.ndim != 2 or arr.shape[0] != arr.shape[1] or arr.size == 0: + return float("nan") + if not np.all(np.isfinite(arr)): + return float("nan") + n = arr.shape[0] + d = int(round(np.sqrt(n))) + if d * d != n: + return float("nan") + trace_rows = np.arange(0, n, d + 1, dtype=int) + return float( + max( + (scaled_cancellation_ratio(arr[trace_rows, j]) for j in range(n)), + default=0.0, + ) + ) + + def band_discriminates( magnitudes: np.ndarray, zero_count: int, bound: float ) -> bool: @@ -559,6 +597,7 @@ def certified_eigvals( L_super = require_finite_square_2d(L_super, name="L_super") L_c = np.asarray(L_super, dtype=complex) tp_defect, fro = trace_preservation_defect(L_c) + tp_componentwise = trace_preservation_componentwise_error(L_c) norm2 = float(np.linalg.norm(L_c, 2)) if L_c.size else 0.0 bound = rtol * float(np.finfo(float).eps) * norm2 @@ -587,6 +626,13 @@ def certified_eigvals( # RELATIVE to the operator's own scale is not a statement that # can be made about infinity, so it is refused here. or not np.isfinite(fro) + # ISSUE #130. The normwise ratio above can be diluted by enormous + # entries that do not participate in a violated trace equation. The + # componentwise gate normalises each column only by the diagonal-output + # terms that can cancel in that equation, then takes the worst column. + # Both views must pass before the exact zero-mode theorem is applicable. + or not np.isfinite(tp_componentwise) + or tp_componentwise > tp_rtol # ROUND-23 REVIEW (PR #121), same finding one line further down. # ``max(fro, tiny)`` floors the REFERENCE SCALE at a constant that # is larger than the operator whenever the operator is subnormal, @@ -830,6 +876,7 @@ def certified_eig( L_super = require_finite_square_2d(L_super, name="L_super") L_c = np.asarray(L_super, dtype=complex) tp_defect, fro = trace_preservation_defect(L_c) + tp_componentwise = trace_preservation_componentwise_error(L_c) norm2 = float(np.linalg.norm(L_c, 2)) if L_c.size else 0.0 bound = rtol * float(np.finfo(float).eps) * norm2 @@ -858,6 +905,10 @@ def certified_eig( # RELATIVE to the operator's own scale is not a statement that # can be made about infinity, so it is refused here. or not np.isfinite(fro) + # ISSUE #130. Keep the eigenvector ladder on the exact same + # applicability contract as the eigenvalue-only ladder above. + or not np.isfinite(tp_componentwise) + or tp_componentwise > tp_rtol # ROUND-23 REVIEW (PR #121), same finding one line further down. # ``max(fro, tiny)`` floors the REFERENCE SCALE at a constant that # is larger than the operator whenever the operator is subnormal, @@ -1235,4 +1286,4 @@ def support_check( def matrix_2_norm(A: np.ndarray) -> float: """Spectral norm via SVD; works for non-square arrays.""" - return float(sla.svdvals(np.asarray(A))[0]) + return float(sla.svdvals(np.asarray(A))[0]) \ No newline at end of file From 4d32dab67719ca5634fe2c38ec2b98845e1d39bf Mon Sep 17 00:00:00 2001 From: Marco Hostettler Date: Tue, 1 Sep 2026 21:06:06 +0200 Subject: [PATCH 11/22] test(#130): migrate overflow regression to componentwise TP semantics --- tests/test_pr121_review_round20.py | 123 +++++++++++------------------ 1 file changed, 48 insertions(+), 75 deletions(-) diff --git a/tests/test_pr121_review_round20.py b/tests/test_pr121_review_round20.py index 8ed36204..055f4c67 100644 --- a/tests/test_pr121_review_round20.py +++ b/tests/test_pr121_review_round20.py @@ -139,11 +139,12 @@ def test_a_resolved_certificate_still_publishes_d3_d4_and_the_flag(monkeypatch) def _overflowing_norm_operator() -> np.ndarray: """The reviewer's counterexample, entry for entry. - Finite 4x4, cancelling ``+-1e308`` in the ``vec(I)^H L`` combination, so - the absolute trace defect is ``sqrt(17) = 4.12`` while ``||L||_F`` - overflows to ``inf``. Eigenvalues are exactly ``{1, 2, 3, 4}``: column 0 - and column 2 each carry a single diagonal entry, and the remaining 2x2 on - indices ``{1, 3}`` is lower triangular. + Finite 4x4, cancelling ``+-1e308`` in one trace equation. A naive + sum-of-squares Frobenius norm overflows, but the TRUE norm is the finite + representable value ``sqrt(2)*1e308``. The absolute trace defect is + ``sqrt(17) = 4.12`` and the componentwise trace-equation error is 1: other + enormous entries are not allowed to dilute the violated columns. + Eigenvalues are exactly ``{1, 2, 3, 4}``. """ L = np.diag([1.0, 2.0, 3.0, 4.0]).astype(complex) L[0, 1] = 1.0e308 @@ -152,14 +153,13 @@ def _overflowing_norm_operator() -> np.ndarray: def _nilpotent_wide_band_operator() -> np.ndarray: - """The same failure WITHOUT any overflow, to show the class is general. + """The same band failure WITHOUT a trace-preservation failure. Strictly upper triangular with ``1e10`` entries confined to row 1, so ``vec(I)^H L = 0`` EXACTLY (trace preserving, defect 0.0), ``||L||_F`` is a perfectly ordinary ``1.4e10``, and every eigenvalue is exactly zero. The band ``rtol * eps * ||L||_2`` is then ``3.1e-3`` and swallows all four - -- a nilpotent generator certified as a four-dimensional stationary - manifold, on finite, well-scaled arithmetic. + -- a nilpotent generator whose band has discriminated nothing. """ L = np.zeros((4, 4), dtype=complex) L[1, 2] = 1.0e10 @@ -167,34 +167,36 @@ def _nilpotent_wide_band_operator() -> np.ndarray: return L -def test_the_overflowing_operator_is_the_state_the_reviewer_described() -> None: - """Pins the preconditions, so a numpy change cannot hollow out the test.""" - from liouscope.numerics.linalg import trace_preservation_defect +def test_large_operator_measurement_and_local_tp_error_are_both_real() -> None: + """Pin the corrected arithmetic and the independent policy precondition.""" + from liouscope.numerics.linalg import ( + trace_preservation_componentwise_error, + trace_preservation_defect, + ) L = _overflowing_norm_operator() - assert np.all(np.isfinite(L)), "the input must be FINITE for this to be the bug" - with np.errstate(over="ignore"): - defect, fro = trace_preservation_defect(L) + assert np.all(np.isfinite(L)), "the input must be FINITE for this to be the case" + defect, fro = trace_preservation_defect(L) + local_error = trace_preservation_componentwise_error(L) + assert defect == pytest.approx(np.sqrt(17.0)), "a real, order-one trace defect" - assert not np.isfinite(fro), ( - "the whole failure is that ``defect > tp_rtol * fro`` compares " - "against infinity; if ||L||_F stops overflowing the test is hollow" + assert np.isfinite(fro), "a representable norm must not be turned into infinity" + assert fro == pytest.approx(np.sqrt(2.0) * 1.0e308, rel=2.0e-15) + assert local_error == pytest.approx(1.0), ( + "at least one trace equation has no cancellation at all; unrelated " + "1e308 entries must not make that violation look relatively small" ) assert np.allclose(np.sort(np.linalg.eigvals(L).real), [1.0, 2.0, 3.0, 4.0]) @pytest.mark.parametrize( - "make", [_overflowing_norm_operator, _nilpotent_wide_band_operator] + ("make", "expected_applicable"), + [(_overflowing_norm_operator, False), (_nilpotent_wide_band_operator, True)], ) -def test_a_band_that_swallows_the_spectrum_does_not_certify(make) -> None: - """THE regression, on both ladders. - - Before this guard the certificate read ``certified=True, resolved=True`` - with ``zero_mode_count = 4`` -- every eigenvalue declared stationary, - which for a d=2 system asserts that the entire operator space is a steady - state. ``resolved=True`` is the damaging half: it tells the spectral layer - that D1/D3/D4 may be published. - """ +def test_a_non_tp_input_or_a_nondiscriminating_band_does_not_certify( + make, expected_applicable: bool +) -> None: + """Both ladders fail closed, but the two counterexamples fail differently.""" from liouscope.numerics.linalg import certified_eig, certified_eigvals L = make() @@ -203,64 +205,35 @@ def test_a_band_that_swallows_the_spectrum_does_not_certify(make) -> None: _, cert_vals = certified_eigvals(L) _, cert_vec = certified_eig(L) - # ROUND 21. Both operators must be refused, but for DIFFERENT reasons, and - # this test names which -- asserting only the shared verdict would let - # either mechanism rot unnoticed behind the other. - # - # * finite scale (the nilpotent twin): the certificate legitimately - # applies, and the BAND GUARD catches it -- the band swallowed the whole - # spectrum, so it discriminated nothing. - # * non-finite scale (the overflowing operator): there is no measurable - # scale to be trace preserving *relative to*. ``defect <= tp_rtol * fro`` - # against an infinite ``fro`` admits everything, so the band drawn from - # it was never grounded in a measurement. Which modes it then swallows - # is decided by the LAPACK build -- 4 of 4 on Windows, 2 of 4 on the CI - # runners, and the 2-of-4 case slipped past a guard that only asks - # whether the band took *everything*. The refusal now happens one level - # earlier, where it is platform-independent. - with np.errstate(over="ignore"): - finite_scale = bool(np.isfinite(np.linalg.norm(L, ord="fro"))) - for cert in (cert_vals, cert_vec): assert cert.certified is False assert cert.resolved is False - if finite_scale: - assert cert.applicable is True - assert cert.zero_mode_count == 0 - else: - assert cert.applicable is False, ( - "an operator whose own scale overflows cannot be judged " - "trace preserving relative to that scale" - ) - - -def test_a_non_finite_scale_is_refused_before_any_band_is_drawn() -> None: - """ROUND 21: the fail-closed must key on the SCALE, not on the outcome. - - ``trace_preservation_defect`` returns ``(defect, ||L||_F)`` and the - applicability test is ``defect <= tp_rtol * ||L||_F``. A non-finite - ``defect`` was already refused; a non-finite ``||L||_F`` was not, and that - asymmetry is the whole bug -- a comparison against infinity is vacuously - true, so every operator passed and the band that followed was drawn from a - scale nobody measured. - - This pins the invariant directly rather than through a downstream verdict, - so it cannot be satisfied by accident the way the band guard was. + assert cert.applicable is expected_applicable + assert cert.zero_mode_count == 0 + + +def test_componentwise_tp_gate_refuses_local_defect_before_zero_mode_claim() -> None: + """The #130 regression: finite global norm must not create a fail-open. + + The true Frobenius norm is ~1.4e308, so a GLOBAL normwise quotient makes + ``sqrt(17)/||L||`` microscopic. Trace preservation, however, consists of + separate cancellation equations. Several columns here have order-one + equation error, which invalidates the exact-zero-mode theorem regardless + of unrelated entries elsewhere in the operator. """ from liouscope.numerics.linalg import ( certified_eigvals, + trace_preservation_componentwise_error, trace_preservation_defect, ) L = _overflowing_norm_operator() - with np.errstate(over="ignore"): - defect, fro = trace_preservation_defect(L) + defect, fro = trace_preservation_defect(L) + local_error = trace_preservation_componentwise_error(L) - assert np.isfinite(defect), "the defect itself is finite -- that is the trap" - assert not np.isfinite(fro), ( - "if the Frobenius norm stops overflowing this test is hollow and the " - "asymmetry it guards no longer exists" - ) + assert np.isfinite(defect) + assert np.isfinite(fro) + assert local_error == pytest.approx(1.0) with np.errstate(over="ignore"), warnings.catch_warnings(): warnings.simplefilter("ignore") @@ -512,4 +485,4 @@ def test_neff_does_not_launder_a_nan_series_into_the_full_sample_size() -> None: # Positive control: a finite, genuinely uncorrelated series is unaffected. rng = np.random.default_rng(1) n_eff = estimate_neff_geyer(rng.normal(size=64)) - assert np.isfinite(n_eff) and 1.0 <= n_eff <= 64.0 + assert np.isfinite(n_eff) and 1.0 <= n_eff <= 64.0 \ No newline at end of file From ae0c5b81c5dfaa59f30f93936ec22e30496ffa85 Mon Sep 17 00:00:00 2001 From: Marco Hostettler Date: Tue, 1 Sep 2026 21:07:02 +0200 Subject: [PATCH 12/22] test(#130): require local TP refusal at representable large scale --- tests/test_pr121_review_round23.py | 25 ++++++++++++++----------- 1 file changed, 14 insertions(+), 11 deletions(-) diff --git a/tests/test_pr121_review_round23.py b/tests/test_pr121_review_round23.py index 0e8817ee..e8e22672 100644 --- a/tests/test_pr121_review_round23.py +++ b/tests/test_pr121_review_round23.py @@ -21,9 +21,10 @@ Every guard here is paired with a positive control: a guard that refuses everything discriminates nothing and must not be able to pass this file. Issue -#130 extends finding 12 symmetrically: a mathematically representable large norm -must be measured as finite, while the structural zero-mode certificate remains -responsible for refusing evidence that does not discriminate. +#130 extends finding 12 symmetrically: mathematically representable large norms +must stay finite, while trace-preservation applicability also uses a separate +componentwise backward-error check so unrelated huge entries cannot dilute a +violated trace equation. """ from __future__ import annotations @@ -39,6 +40,7 @@ from liouscope.numerics.linalg import ( certified_eig, certified_eigvals, + trace_preservation_componentwise_error, trace_preservation_defect, ) @@ -111,15 +113,14 @@ def test_a_healthy_generator_in_the_same_rate_units_still_certifies() -> None: def test_representable_large_norm_is_measured_then_rejected_structurally() -> None: - """Measurement arithmetic must not carry the certificate's policy. + """Correct measurement must not create a fail-open applicability claim. The round-20 counterexample has two cancelling ``+-1e308`` entries. Its true Frobenius norm is about ``sqrt(2)*1e308`` and is therefore finite in - float64, even though a naive sum of squares overflows. Issue #130 requires - that finite measurement to stay finite. Safety is preserved one layer - later: the zero-mode band contains the whole measured spectrum, so - ``band_discriminates`` leaves the certificate unresolved instead of - treating arithmetic overflow as a surrogate refusal signal. + float64, even though a naive sum of squares overflows. That finite global + scale cannot prove trace preservation: several trace equations have no + cancellation at all, so their componentwise backward error is 1. The exact + zero-mode theorem must therefore be inapplicable before the band is used. """ op = np.diag([1.0, 2.0, 3.0, 4.0]).astype(complex) op[0, 1] = 1.0e308 @@ -127,14 +128,16 @@ def test_representable_large_norm_is_measured_then_rejected_structurally() -> No assert np.all(np.isfinite(op)), "the input must be FINITE for this to be the case" defect, fro = trace_preservation_defect(op) + local_error = trace_preservation_componentwise_error(op) assert defect == pytest.approx(np.sqrt(17.0)), "a real, order-one trace defect" assert np.isfinite(fro), "a representable Frobenius norm must not overflow" assert fro == pytest.approx(np.sqrt(2.0) * 1.0e308, rel=2.0e-15) + assert local_error == pytest.approx(1.0) with np.errstate(over="ignore"), warnings.catch_warnings(): warnings.simplefilter("ignore") cert = certified_eigvals(op)[-1] - assert cert.applicable + assert not cert.applicable assert not cert.certified assert not cert.resolved @@ -376,4 +379,4 @@ def test_the_default_keeps_the_layer_reporting() -> None: n_haar=3, ) assert np.isfinite(lep.lep_proximity) - assert isinstance(lep.lep_candidate_count, int) + assert isinstance(lep.lep_candidate_count, int) \ No newline at end of file From d4d635f34b0a4fa94eae926693dd7abd626d1403 Mon Sep 17 00:00:00 2001 From: Marco Hostettler Date: Tue, 1 Sep 2026 21:07:19 +0200 Subject: [PATCH 13/22] test(#130): pin componentwise TP backward-error invariants --- tests/test_issue130_componentwise_tp.py | 62 +++++++++++++++++++++++++ 1 file changed, 62 insertions(+) create mode 100644 tests/test_issue130_componentwise_tp.py diff --git a/tests/test_issue130_componentwise_tp.py b/tests/test_issue130_componentwise_tp.py new file mode 100644 index 00000000..8df6ff80 --- /dev/null +++ b/tests/test_issue130_componentwise_tp.py @@ -0,0 +1,62 @@ +"""Issue #130: componentwise trace-preservation backward-error regressions.""" + +from __future__ import annotations + +import numpy as np +import pytest + +from liouscope.numerics.linalg import trace_preservation_componentwise_error +from liouscope.numerics.norms import scaled_cancellation_ratio + + +def test_scaled_cancellation_ratio_exact_cancellation_at_extreme_scales() -> None: + for scale in (1.0e-300, 1.0, 1.0e300): + values = np.array([scale, -scale], dtype=complex) + assert scaled_cancellation_ratio(values) == 0.0 + + +def test_scaled_cancellation_ratio_detects_uncancelled_equation() -> None: + for scale in (1.0e-300, 1.0, 1.0e300): + values = np.array([scale, 0.0], dtype=complex) + assert scaled_cancellation_ratio(values) == pytest.approx(1.0) + + +def test_scaled_cancellation_ratio_is_scale_invariant_for_partial_cancellation() -> None: + reference = np.array([2.0 + 3.0j, -2.0 - 2.0j, 0.25 - 0.5j]) + expected = scaled_cancellation_ratio(reference) + assert 0.0 < expected < 1.0 + for scale in (1.0e-250, 1.0e-100, 1.0, 1.0e100, 1.0e250): + got = scaled_cancellation_ratio(reference * scale) + assert got == pytest.approx(expected, rel=2.0e-15, abs=0.0) + + +def test_tp_componentwise_error_ignores_unrelated_huge_entries() -> None: + # d=2 => trace-output rows are vectorised matrix rows 0 and 3. + op = np.zeros((4, 4), dtype=complex) + op[1, 0] = 1.0e300 # huge but irrelevant to the trace equation + op[0, 2] = 1.0 + assert trace_preservation_componentwise_error(op) == pytest.approx(1.0) + + +def test_tp_componentwise_error_accepts_exact_trace_cancellation_at_huge_scale() -> None: + op = np.zeros((4, 4), dtype=complex) + op[0, 1] = 1.0e300 + op[3, 1] = -1.0e300 + assert trace_preservation_componentwise_error(op) == 0.0 + + +def test_tp_componentwise_error_is_rate_unit_invariant() -> None: + op = np.zeros((4, 4), dtype=complex) + op[0, 0] = 3.0 + op[3, 0] = -2.0 + expected = 1.0 / 5.0 + for scale in (1.0e-300, 1.0e-200, 1.0, 1.0e200, 1.0e300): + got = trace_preservation_componentwise_error(op * scale) + assert got == pytest.approx(expected, rel=2.0e-15, abs=0.0) + + +def test_tp_componentwise_error_fails_closed_on_bad_shape_and_nonfinite() -> None: + assert np.isnan(trace_preservation_componentwise_error(np.ones((3, 3)))) + bad = np.zeros((4, 4), dtype=complex) + bad[0, 0] = np.nan + assert np.isnan(trace_preservation_componentwise_error(bad)) From eed9c0ffd47fc2a4a0d45fc55f2fc80d81f4664a Mon Sep 17 00:00:00 2001 From: Marco Hostettler Date: Tue, 1 Sep 2026 22:00:23 +0200 Subject: [PATCH 14/22] ci: add one-shot byte-preserving EOF repair --- .github/workflows/repair-pr136-eof.yml | 81 ++++++++++++++++++++++++++ 1 file changed, 81 insertions(+) create mode 100644 .github/workflows/repair-pr136-eof.yml diff --git a/.github/workflows/repair-pr136-eof.yml b/.github/workflows/repair-pr136-eof.yml new file mode 100644 index 00000000..000ddb34 --- /dev/null +++ b/.github/workflows/repair-pr136-eof.yml @@ -0,0 +1,81 @@ +name: One-shot PR136 EOF repair + +on: + push: + branches: + - "chatgpt/issue-130-scaled-norms" + +permissions: + contents: write + +concurrency: + group: pr136-eof-repair + cancel-in-progress: false + +jobs: + repair: + if: github.repository == 'marcohost33-maker/Liouscope' + runs-on: ubuntu-latest + steps: + - name: Checkout exact branch head + uses: actions/checkout@11bd71901bbe5b1630ceea73d27597364c9af683 # v4.2.2 + with: + persist-credentials: true + fetch-depth: 0 + + - name: Append exactly one missing terminal newline + shell: python + run: | + from pathlib import Path + + targets = [ + Path("src/liouscope/numerics/linalg.py"), + Path("tests/test_pr121_review_round20.py"), + Path("tests/test_pr121_review_round23.py"), + ] + for path in targets: + before = path.read_bytes() + if before.endswith(b"\n"): + raise SystemExit(f"{path}: expected the known missing-EOF-newline defect") + after = before + b"\n" + path.write_bytes(after) + if path.read_bytes() != after or after[:-1] != before: + raise SystemExit(f"{path}: byte-preservation invariant failed") + + workflow = Path(".github/workflows/repair-pr136-eof.yml") + workflow.unlink() + Path(".pr136-eof-repair.marker").write_text( + "one-shot EOF repair completed; delete this marker via the external GitHub connector to trigger final PR CI\n", + encoding="utf-8", + ) + + - name: Verify change surface + shell: bash + run: | + set -euo pipefail + git diff --check + python - <<'PY' + import subprocess + + expected = { + ".github/workflows/repair-pr136-eof.yml", + ".pr136-eof-repair.marker", + "src/liouscope/numerics/linalg.py", + "tests/test_pr121_review_round20.py", + "tests/test_pr121_review_round23.py", + } + observed = set(subprocess.check_output(["git", "status", "--short"], text=True).splitlines()) + names = {line[3:] for line in observed} + if names != expected: + raise SystemExit(f"unexpected repair surface: {sorted(names)}") + PY + + - name: Commit byte-only repair and self-delete + shell: bash + run: | + set -euo pipefail + git config user.name "liouscope-eof-repair" + git config user.email "actions@users.noreply.github.com" + git add -A + git commit -m "style(#130): restore terminal newlines" + git push origin "HEAD:${GITHUB_REF_NAME}" From 350faea72b63ffe0548966b267909605cabd6200 Mon Sep 17 00:00:00 2001 From: liouscope-eof-repair Date: Tue, 1 Sep 2026 20:00:33 +0000 Subject: [PATCH 15/22] style(#130): restore terminal newlines --- .github/workflows/repair-pr136-eof.yml | 81 -------------------------- .pr136-eof-repair.marker | 1 + src/liouscope/numerics/linalg.py | 2 +- tests/test_pr121_review_round20.py | 2 +- tests/test_pr121_review_round23.py | 2 +- 5 files changed, 4 insertions(+), 84 deletions(-) delete mode 100644 .github/workflows/repair-pr136-eof.yml create mode 100644 .pr136-eof-repair.marker diff --git a/.github/workflows/repair-pr136-eof.yml b/.github/workflows/repair-pr136-eof.yml deleted file mode 100644 index 000ddb34..00000000 --- a/.github/workflows/repair-pr136-eof.yml +++ /dev/null @@ -1,81 +0,0 @@ -name: One-shot PR136 EOF repair - -on: - push: - branches: - - "chatgpt/issue-130-scaled-norms" - -permissions: - contents: write - -concurrency: - group: pr136-eof-repair - cancel-in-progress: false - -jobs: - repair: - if: github.repository == 'marcohost33-maker/Liouscope' - runs-on: ubuntu-latest - steps: - - name: Checkout exact branch head - uses: actions/checkout@11bd71901bbe5b1630ceea73d27597364c9af683 # v4.2.2 - with: - persist-credentials: true - fetch-depth: 0 - - - name: Append exactly one missing terminal newline - shell: python - run: | - from pathlib import Path - - targets = [ - Path("src/liouscope/numerics/linalg.py"), - Path("tests/test_pr121_review_round20.py"), - Path("tests/test_pr121_review_round23.py"), - ] - for path in targets: - before = path.read_bytes() - if before.endswith(b"\n"): - raise SystemExit(f"{path}: expected the known missing-EOF-newline defect") - after = before + b"\n" - path.write_bytes(after) - if path.read_bytes() != after or after[:-1] != before: - raise SystemExit(f"{path}: byte-preservation invariant failed") - - workflow = Path(".github/workflows/repair-pr136-eof.yml") - workflow.unlink() - Path(".pr136-eof-repair.marker").write_text( - "one-shot EOF repair completed; delete this marker via the external GitHub connector to trigger final PR CI\n", - encoding="utf-8", - ) - - - name: Verify change surface - shell: bash - run: | - set -euo pipefail - git diff --check - python - <<'PY' - import subprocess - - expected = { - ".github/workflows/repair-pr136-eof.yml", - ".pr136-eof-repair.marker", - "src/liouscope/numerics/linalg.py", - "tests/test_pr121_review_round20.py", - "tests/test_pr121_review_round23.py", - } - observed = set(subprocess.check_output(["git", "status", "--short"], text=True).splitlines()) - names = {line[3:] for line in observed} - if names != expected: - raise SystemExit(f"unexpected repair surface: {sorted(names)}") - PY - - - name: Commit byte-only repair and self-delete - shell: bash - run: | - set -euo pipefail - git config user.name "liouscope-eof-repair" - git config user.email "actions@users.noreply.github.com" - git add -A - git commit -m "style(#130): restore terminal newlines" - git push origin "HEAD:${GITHUB_REF_NAME}" diff --git a/.pr136-eof-repair.marker b/.pr136-eof-repair.marker new file mode 100644 index 00000000..627e204d --- /dev/null +++ b/.pr136-eof-repair.marker @@ -0,0 +1 @@ +one-shot EOF repair completed; delete this marker via the external GitHub connector to trigger final PR CI diff --git a/src/liouscope/numerics/linalg.py b/src/liouscope/numerics/linalg.py index 01242e49..3708d035 100644 --- a/src/liouscope/numerics/linalg.py +++ b/src/liouscope/numerics/linalg.py @@ -1286,4 +1286,4 @@ def support_check( def matrix_2_norm(A: np.ndarray) -> float: """Spectral norm via SVD; works for non-square arrays.""" - return float(sla.svdvals(np.asarray(A))[0]) \ No newline at end of file + return float(sla.svdvals(np.asarray(A))[0]) diff --git a/tests/test_pr121_review_round20.py b/tests/test_pr121_review_round20.py index 055f4c67..d09c3489 100644 --- a/tests/test_pr121_review_round20.py +++ b/tests/test_pr121_review_round20.py @@ -485,4 +485,4 @@ def test_neff_does_not_launder_a_nan_series_into_the_full_sample_size() -> None: # Positive control: a finite, genuinely uncorrelated series is unaffected. rng = np.random.default_rng(1) n_eff = estimate_neff_geyer(rng.normal(size=64)) - assert np.isfinite(n_eff) and 1.0 <= n_eff <= 64.0 \ No newline at end of file + assert np.isfinite(n_eff) and 1.0 <= n_eff <= 64.0 diff --git a/tests/test_pr121_review_round23.py b/tests/test_pr121_review_round23.py index e8e22672..bc48aa28 100644 --- a/tests/test_pr121_review_round23.py +++ b/tests/test_pr121_review_round23.py @@ -379,4 +379,4 @@ def test_the_default_keeps_the_layer_reporting() -> None: n_haar=3, ) assert np.isfinite(lep.lep_proximity) - assert isinstance(lep.lep_candidate_count, int) \ No newline at end of file + assert isinstance(lep.lep_candidate_count, int) From 59f861a571c5ee06d3c5aa2917079fe10eb22641 Mon Sep 17 00:00:00 2001 From: Marco Hostettler Date: Tue, 1 Sep 2026 22:00:53 +0200 Subject: [PATCH 16/22] ci(#130): finalize byte-only EOF repair --- .pr136-eof-repair.marker | 1 - 1 file changed, 1 deletion(-) delete mode 100644 .pr136-eof-repair.marker diff --git a/.pr136-eof-repair.marker b/.pr136-eof-repair.marker deleted file mode 100644 index 627e204d..00000000 --- a/.pr136-eof-repair.marker +++ /dev/null @@ -1 +0,0 @@ -one-shot EOF repair completed; delete this marker via the external GitHub connector to trigger final PR CI From 561c2bb7bf8653e229f552dfea945a3a03e10c6d Mon Sep 17 00:00:00 2001 From: Marco Hostettler Date: Tue, 1 Sep 2026 23:50:40 +0200 Subject: [PATCH 17/22] fix(#135): make Gaussian profile likelihood scale-safe --- src/liouscope/_types.py | 1 + src/liouscope/diagnostics/relaxation.py | 1 + src/liouscope/fitting/aicc.py | 45 ++++++++-- src/liouscope/fitting/bootstrap.py | 2 + src/liouscope/fitting/gls.py | 58 ++++++++++-- src/liouscope/numerics/norms.py | 34 +++++++ tests/test_issue135_likelihood_scale.py | 113 ++++++++++++++++++++++++ 7 files changed, 243 insertions(+), 11 deletions(-) create mode 100644 tests/test_issue135_likelihood_scale.py diff --git a/src/liouscope/_types.py b/src/liouscope/_types.py index f6893937..2005292e 100644 --- a/src/liouscope/_types.py +++ b/src/liouscope/_types.py @@ -187,6 +187,7 @@ class FitResult: n_eff: float residual_ar1_rho: float success: bool + likelihood_degenerate: bool = False @dataclass(frozen=True, slots=True, kw_only=True) diff --git a/src/liouscope/diagnostics/relaxation.py b/src/liouscope/diagnostics/relaxation.py index d8433bc7..29d334fd 100644 --- a/src/liouscope/diagnostics/relaxation.py +++ b/src/liouscope/diagnostics/relaxation.py @@ -215,6 +215,7 @@ def _fit_with_model( n_eff=n_eff, residual_ar1_rho=fit.rho_ar1, success=fit.success, + likelihood_degenerate=fit.likelihood_degenerate, ) return fit_result, fit.params diff --git a/src/liouscope/fitting/aicc.py b/src/liouscope/fitting/aicc.py index 2130a6de..c40fb40b 100644 --- a/src/liouscope/fitting/aicc.py +++ b/src/liouscope/fitting/aicc.py @@ -11,10 +11,13 @@ from __future__ import annotations +import math from collections.abc import Mapping import numpy as np +from ..numerics.norms import scaled_log_sum_squares + def aicc(log_likelihood: float, k: int, n_eff: float) -> float: """Return ``AICc`` for given log-likelihood, parameter count, and N_eff.""" @@ -29,18 +32,48 @@ def gaussian_log_likelihood( *, sigma: float | None = None, ) -> float: - """Gaussian log-likelihood for ``y - y_hat``. + """Gaussian log-likelihood for ``y - y_hat`` without absolute RSS floors. - Uses MLE sigma if ``sigma`` is omitted. + When ``sigma`` is omitted, evaluate the profile likelihood at the Gaussian + MLE ``sigma_hat**2 = RSS / n`` directly in log-RSS space. Exact zero RSS has + no finite interior MLE for the positive scale parameter and therefore + returns NaN (model-selection likelihood unavailable) instead of inventing + an absolute epsilon variance. With an explicitly supplied finite positive + ``sigma``, zero residuals remain a valid finite likelihood. """ residuals = np.asarray(residuals, dtype=float) n = residuals.size - rss = float(np.dot(residuals, residuals)) + if n == 0: + return float("nan") + + log_rss = scaled_log_sum_squares(residuals) + if math.isnan(log_rss): + return float("nan") + + log_2pi = math.log(2.0 * math.pi) if sigma is None: - sigma_sq = max(rss / n, 1.0e-30) + if log_rss == float("-inf"): + return float("nan") + if log_rss == float("inf"): + return float("-inf") + return float(-0.5 * n * (log_2pi + 1.0 + log_rss - math.log(n))) + + sigma = float(sigma) + if not math.isfinite(sigma) or sigma <= 0.0: + return float("nan") + log_sigma = math.log(sigma) + if log_rss == float("-inf"): + standardised_rss = 0.0 + elif log_rss == float("inf"): + return float("-inf") else: - sigma_sq = sigma * sigma - return float(-0.5 * n * (np.log(2.0 * np.pi * sigma_sq) + rss / (n * sigma_sq))) + log_standardised_rss = log_rss - 2.0 * log_sigma + if log_standardised_rss > math.log(np.finfo(float).max): + return float("-inf") + standardised_rss = math.exp(log_standardised_rss) + return float( + -0.5 * n * log_2pi - n * log_sigma - 0.5 * standardised_rss + ) def choose_model(aiccs: Mapping[str, float]) -> str: diff --git a/src/liouscope/fitting/bootstrap.py b/src/liouscope/fitting/bootstrap.py index b936ee56..f7698a36 100644 --- a/src/liouscope/fitting/bootstrap.py +++ b/src/liouscope/fitting/bootstrap.py @@ -65,6 +65,8 @@ def parametric_bootstrap( + (f" (saturated: {', '.join(base.saturated)})" if base.saturated else "") + (" (the curve carries no resolvable variation, issue #123)" if base.degenerate else "") + + (" (the residual likelihood scale is degenerate, issue #135)" + if base.likelihood_degenerate else "") + "; a bootstrap around a non-estimate has no meaning" ) theta_hat = base.params diff --git a/src/liouscope/fitting/gls.py b/src/liouscope/fitting/gls.py index 73b54f49..07b04446 100644 --- a/src/liouscope/fitting/gls.py +++ b/src/liouscope/fitting/gls.py @@ -12,6 +12,7 @@ from __future__ import annotations +import math import warnings from collections.abc import Callable from dataclasses import dataclass @@ -19,6 +20,7 @@ import numpy as np from scipy.optimize import least_squares +from ..numerics.norms import scaled_log_sum_squares from .aicc import gaussian_log_likelihood from .models import saturation_watch from .neff import _AR1_SMALL_N, ar1_correlation_corrected @@ -41,6 +43,11 @@ class GLSFitOutput: #: reports a fit that ran and ended on a magnitude plateau; here there was #: nothing to fit. Implies ``success`` is False and ``params`` is NaN. degenerate: bool = False + #: True when the residual Gaussian scale has no finite usable MLE for model + #: selection (issue #135). Distinct from ``degenerate`` above: the curve + #: may carry variation and the optimiser may have run, but exact-zero RSS + #: (or an unrepresentable positive MLE scale) cannot support AICc/CI claims. + likelihood_degenerate: bool = False def _whiten(y: np.ndarray, rho: float) -> np.ndarray: @@ -192,14 +199,55 @@ def residual(params: np.ndarray, rho_local: float = rho) -> np.ndarray: ) whitened = _whiten(residuals_raw, rho) n = whitened.size - sigma = float(np.sqrt(max(np.dot(whitened, whitened) / max(n, 1), 1.0e-30))) + log_rss = scaled_log_sum_squares(whitened) + if log_rss == float("-inf") or not math.isfinite(log_rss): + warnings.warn( + "fit_gls_ar1: residual Gaussian scale has no finite positive MLE " + "for model selection; likelihood/AICc/CI evidence is unavailable " + "(issue #135).", + RuntimeWarning, + stacklevel=2, + ) + return GLSFitOutput( + params=p, + residuals=residuals_raw, + rho_ar1=rho, + sigma=float("nan"), + log_likelihood=float("nan"), + success=False, + saturated=tuple(sorted(fired)), + likelihood_degenerate=True, + ) + + log_sigma = 0.5 * (log_rss - math.log(n)) + try: + sigma = float(math.exp(log_sigma)) + except OverflowError: + sigma = float("inf") + if not math.isfinite(sigma) or sigma <= 0.0: + warnings.warn( + "fit_gls_ar1: positive residual MLE scale is not representable as " + "float64; likelihood/AICc/CI evidence is unavailable (issue #135).", + RuntimeWarning, + stacklevel=2, + ) + return GLSFitOutput( + params=p, + residuals=residuals_raw, + rho_ar1=rho, + sigma=float("nan"), + log_likelihood=float("nan"), + success=False, + saturated=tuple(sorted(fired)), + likelihood_degenerate=True, + ) + # Prais-Winsten exact AR(1) likelihood: _whiten keeps observation 0 # (scaled by sqrt(1-rho^2)), so the transform has log-Jacobian - # 0.5*log(1-rho^2). Omitting it makes the reported value a hybrid of the - # exact and conditional likelihoods and biases cross-model AICc (each - # model fits its own rho) toward under-fitting high-rho models. + # 0.5*log(1-rho^2). The profile likelihood is evaluated directly from + # log(RSS), not by squaring ``sigma`` or materialising RSS. jac = 0.5 * float(np.log(max(1.0 - rho * rho, 1.0e-12))) - log_lik = gaussian_log_likelihood(whitened, sigma=sigma) + jac + log_lik = gaussian_log_likelihood(whitened) + jac return GLSFitOutput( params=p, residuals=residuals_raw, diff --git a/src/liouscope/numerics/norms.py b/src/liouscope/numerics/norms.py index 067d8012..712973f9 100644 --- a/src/liouscope/numerics/norms.py +++ b/src/liouscope/numerics/norms.py @@ -73,6 +73,40 @@ def scaled_euclidean_norm(values: np.ndarray) -> float: return float(np.ldexp(scaled_norm, exponent)) +def scaled_log_sum_squares(values: np.ndarray) -> float: + """Return ``log(sum(abs(values)**2))`` without spurious under/overflow. + + The return value is ``-inf`` for an exact all-zero input, ``nan`` when any + component is NaN, and ``inf`` when any component is infinite. For finite, + non-zero float64 input the logarithm remains finite even when the true sum + of squares (or its square root) is outside the representable float64 range. + + This is the likelihood-facing companion to :func:`scaled_euclidean_norm`: + both use the same exact power-of-two scaling contract, but this function + never reconstructs RSS in ordinary floating-point units. + """ + arr = np.asarray(values) + if arr.size == 0: + return float("-inf") + + real = np.asarray(np.real(arr), dtype=float) + imag = np.asarray(np.imag(arr), dtype=float) + if np.any(np.isnan(real)) or np.any(np.isnan(imag)): + return float("nan") + if np.any(np.isinf(real)) or np.any(np.isinf(imag)): + return float("inf") + + scaled = _finite_component_scale(arr) + if scaled is None: + return float("-inf") + scaled_real, scaled_imag, exponent = scaled + sumsq = float( + np.sum(scaled_real * scaled_real, dtype=float) + + np.sum(scaled_imag * scaled_imag, dtype=float) + ) + return float(math.log(sumsq) + 2.0 * exponent * math.log(2.0)) + + def scaled_cancellation_ratio(values: np.ndarray) -> float: """Return ``abs(sum(values)) / sum(abs(values))`` scale-safely. diff --git a/tests/test_issue135_likelihood_scale.py b/tests/test_issue135_likelihood_scale.py new file mode 100644 index 00000000..4e9af428 --- /dev/null +++ b/tests/test_issue135_likelihood_scale.py @@ -0,0 +1,113 @@ +"""Issue #135: scale-safe Gaussian likelihood and degenerate RSS semantics.""" + +from __future__ import annotations + +import math + +import numpy as np +import pytest + +from liouscope.diagnostics import relaxation as relaxation_mod +from liouscope.fitting.aicc import aicc, choose_model, gaussian_log_likelihood +from liouscope.fitting.bootstrap import parametric_bootstrap +from liouscope.fitting.gls import GLSFitOutput, fit_gls_ar1 +from liouscope.fitting.models import M0 +from liouscope.numerics.norms import scaled_log_sum_squares + + +def test_scaled_log_rss_spans_underflow_and_overflow_regimes(): + tiny = np.array([1.0e-320, -2.0e-320]) + huge = np.array([1.0e308, -1.0e308]) + log_tiny = scaled_log_sum_squares(tiny) + log_huge = scaled_log_sum_squares(huge) + assert np.isfinite(log_tiny) + assert np.isfinite(log_huge) + assert log_tiny == pytest.approx(math.log(5.0) + 2.0 * math.log(1.0e-320), rel=2e-5) + assert log_huge == pytest.approx(math.log(2.0) + 2.0 * math.log(1.0e308), rel=1e-14) + assert scaled_log_sum_squares(np.zeros(4)) == float("-inf") + + +def test_profile_loglikelihood_is_scale_covariant_and_delta_invariant(): + x = np.linspace(0.2, 1.2, 64) + r0 = 0.7 + x + 0.03 * np.sin(5.0 * x) + r1 = 0.9 + 1.1 * x - 0.02 * np.cos(3.0 * x) + reference_delta = None + reference_winner = None + for scale in (1.0e-150, 1.0e-40, 1.0, 1.0e40, 1.0e150): + ll0 = gaussian_log_likelihood(scale * r0) + ll1 = gaussian_log_likelihood(scale * r1) + assert np.isfinite(ll0) and np.isfinite(ll1) + delta = ll1 - ll0 + scores = { + "M0": aicc(ll0, k=2, n_eff=64.0), + "M1": aicc(ll1, k=3, n_eff=64.0), + } + winner = choose_model(scores) + if reference_delta is None: + reference_delta = delta + reference_winner = winner + else: + assert delta == pytest.approx(reference_delta, rel=1e-11, abs=1e-10) + assert winner == reference_winner + + +def test_profile_loglikelihood_handles_true_rss_above_float_range(): + ll = gaussian_log_likelihood(np.array([1.0e308, -1.0e308])) + assert np.isfinite(ll) + + +def test_zero_rss_unknown_sigma_is_unavailable_but_known_sigma_is_valid(): + residuals = np.zeros(8) + assert np.isnan(gaussian_log_likelihood(residuals)) + ll = gaussian_log_likelihood(residuals, sigma=2.0) + expected = -0.5 * residuals.size * math.log(2.0 * math.pi * 4.0) + assert ll == pytest.approx(expected, rel=1e-14, abs=1e-14) + + +def test_gls_exact_fit_marks_likelihood_degenerate_and_bootstrap_refuses(): + t = np.linspace(0.0, 4.0, 48) + p0 = np.array([1.25, 0.6]) + y = M0(t, p0) + with pytest.warns(RuntimeWarning, match="likelihood/AICc/CI evidence is unavailable"): + fit = fit_gls_ar1(M0, t, y, p0, n_iters=1) + assert not fit.success + assert fit.likelihood_degenerate + assert not fit.degenerate + assert np.isnan(fit.sigma) + assert np.isnan(fit.log_likelihood) + # The repository treats unexpected warnings as errors. Here the warning is + # part of the intended public contract: the base fit first reports why its + # likelihood evidence is unusable, then bootstrap refuses the non-estimate. + with pytest.warns(RuntimeWarning, match="likelihood/AICc/CI evidence is unavailable"): + with pytest.raises(RuntimeError, match="likelihood scale is degenerate"): + parametric_bootstrap(M0, t, y, p0, B=4) + + +def test_ordinary_noisy_gls_positive_control_is_not_likelihood_degenerate(rng): + t = np.linspace(0.0, 4.0, 80) + p = np.array([1.25, 0.6]) + y = M0(t, p) + 1.0e-3 * rng.standard_normal(t.size) + fit = fit_gls_ar1(M0, t, y, p, n_iters=1) + assert fit.success + assert not fit.likelihood_degenerate + assert np.isfinite(fit.sigma) and fit.sigma > 0.0 + assert np.isfinite(fit.log_likelihood) + + +def test_likelihood_degenerate_state_reaches_fitresult_and_is_nonselectable(monkeypatch): + t = np.linspace(0.0, 1.0, 16) + y = np.exp(-t) + fake = GLSFitOutput( + params=np.array([1.0, 1.0]), + residuals=np.zeros_like(y), + rho_ar1=0.0, + sigma=float("nan"), + log_likelihood=float("nan"), + success=False, + likelihood_degenerate=True, + ) + monkeypatch.setattr(relaxation_mod, "fit_gls_ar1", lambda *args, **kwargs: fake) + fit_result, _ = relaxation_mod._fit_with_model("M0", t, y) + assert not fit_result.success + assert fit_result.likelihood_degenerate + assert np.isinf(fit_result.aicc) From cc2cd2c7cbe0f96c6c6f5391d19243c825fe5121 Mon Sep 17 00:00:00 2001 From: marcohost33-maker Date: Thu, 3 Sep 2026 00:48:50 +0200 Subject: [PATCH 18/22] fix(#135): a fit withheld because a number its likelihood never uses 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) Claude-Session: https://claude.ai/code/session_01FPKDUXBV4qXpWX8ZtpqgAZ --- CHANGELOG.md | 53 ++++++++++++++ CITATION.cff | 17 +++++ src/liouscope/fitting/bootstrap.py | 15 ++++ src/liouscope/fitting/gls.py | 58 +++++++++++---- tests/test_issue135_likelihood_scale.py | 97 +++++++++++++++++++++++++ 5 files changed, 226 insertions(+), 14 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 9968fbca..3375ce73 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,7 +6,60 @@ adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). ## [Unreleased] +### Changed +- **The Gaussian likelihood behind AICc is evaluated in log-RSS space, and an + exact-zero RSS is now an explicit abstention (issue #135).** This is a + METHODOLOGY change with user-visible consequences: it can reorder AICc, + change which of M0..M3b is selected, and therefore change the reported decay + rate for a given run. Three parts: + - `gaussian_log_likelihood` evaluates the profile likelihood at the Gaussian + MLE directly from `log(RSS)` via the new + `numerics.norms.scaled_log_sum_squares`, instead of forming RSS in ordinary + units. Residual series whose true sum of squares lies outside the float64 + range — a rescaled curve is enough — previously produced `0` or `inf` RSS + and an unusable likelihood. The value is now scale-covariant, and + likelihood *differences* between models are invariant under a common + rescaling of the residuals, which is the quantity AICc actually consumes. + - An exact-zero RSS has no finite interior MLE for the positive scale + parameter, so the profile likelihood is `NaN` rather than an invented + absolute epsilon variance. `fit_gls_ar1` reports this as + `likelihood_degenerate=True` with `success=False`, `aicc()` then scores the + model `inf`, and `parametric_bootstrap` refuses: a perfect fit yields no + interval rather than a zero-width one. Passing an explicit positive `sigma` + still gives a finite likelihood for zero residuals, which is the + well-posed case. + - Runs whose reported model or rate changes are those where the previous + absolute-RSS path had over- or underflowed, or where the winner was chosen + against a floored likelihood. Re-run any archived analysis whose selected + model matters; the run manifest does not record per-model likelihoods, so + the change is not visible in `input_hash`. + ### Fixed +- **A fit was withheld entirely because a number its likelihood never uses + could not be materialised (PR #147, round-1 review).** When the whitened + residuals have 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 `fit_gls_ar1` + then returned `success=False`, `log_likelihood=NaN` and + `likelihood_degenerate=True`. `_fit_with_model` scored that model `inf` and + dropped it from selection, which reintroduces exactly the ABSOLUTE SCALE + BOUNDARY into model selection that the change above removed: the same curve + in different rate units either is or is not a candidate. Measured on one + minimum-subnormal residual (`5e-324`) among 64 otherwise-zero points: + `log_rss = -1488.88`, `log_sigma = -746.52`, profile log-likelihood + `+47686.44`, `exp(log_sigma) = 0.0`. + + The fit now stays selectable with its finite log-space likelihood, and only + the scale-dependent evidence is withheld. `GLSFitOutput` gains + `scale_unavailable` (additive, default `False`); `sigma` is `NaN` there, + deliberately not `0.0`, because `_ar1_resample` consumes it as the + innovation standard deviation and `0.0` would generate identical replicates + — a zero-width confidence interval, which is the failure mode of an + uncertainty pipeline, not 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 handler. Three mutations, one per guard line, are + proven to discriminate. - **D16 was published from the spectrum for which D1/D3/D4 had just been withheld (PR #121, round-23 review, finding 13).** Where the zero-mode certificate is applicable but unresolved, the spectral layer reports D1, D3, diff --git a/CITATION.cff b/CITATION.cff index 18b5faca..202b5eb2 100644 --- a/CITATION.cff +++ b/CITATION.cff @@ -84,6 +84,23 @@ abstract: >- # likewise scale-relative rather than absolute (issue #109), which closes a # fail-open path that admitted non-GKSL generators at small ||H||; this changes # which inputs are ACCEPTED, not the numerics of accepted ones. +# Also pending for the next cut: the Gaussian likelihood behind AICc is +# evaluated in log-RSS space rather than by forming the residual sum of squares +# in ordinary units, and an exact-zero RSS is an explicit abstention rather than +# a floored variance (issue #135). This CHANGES NUMERICAL RESULTS and is a +# METHODOLOGY correction, not an addition: AICc ordering, the selected model +# M0..M3b and therefore the reported decay rate can differ for runs whose +# residual scale previously over- or underflowed, and a perfect fit now yields +# no confidence interval instead of a zero-width one. The likelihood value is +# scale-covariant and likelihood DIFFERENCES between models are invariant under +# a common rescaling of the residuals, which is the quantity AICc consumes. It +# must NOT be described as making model selection scale-free in general: only +# the likelihood path was corrected, and the residual-scale parameter itself +# still has a float64 representation limit -- a fit whose MLE scale underflows +# keeps its log-space likelihood and stays selectable, but its bootstrap +# interval is withheld (round-1 review of PR #147). Archived analyses whose +# selected model matters should be re-run; per-model likelihoods are not run- +# manifest fields, so this change is not visible in input_hash. keywords: - open quantum systems - Lindblad diff --git a/src/liouscope/fitting/bootstrap.py b/src/liouscope/fitting/bootstrap.py index f7698a36..1ec84e17 100644 --- a/src/liouscope/fitting/bootstrap.py +++ b/src/liouscope/fitting/bootstrap.py @@ -52,6 +52,21 @@ def parametric_bootstrap( if rng is None: rng = np.random.default_rng(0) base = fit_gls_ar1(model, t, y, p0, bounds=bounds) + if base.scale_unavailable: + # Round-1 review (PR #147). The base fit is a legitimate AICc + # candidate -- its log-space likelihood is finite -- but its innovation + # scale is not representable, so there is nothing to draw replicates + # from. ``sigma`` is NaN there by construction, and ``rng.normal(0, nan)`` + # returns NaN rather than raising, so every replicate would be fitted to + # NaN data and the resulting interval would be an artefact. Refusing + # here routes the case into the caller's existing handler, which reports + # the CI as NaN -- "fit uncertainty UNKNOWN", which is the true state. + raise RuntimeError( + "parametric_bootstrap: the base fit has no representable residual " + "scale (issue #135), so replicates cannot be simulated; the fit " + "itself remains a valid model-selection candidate, only its " + "interval is unavailable" + ) if not base.success: # Round-17 review (PR #121). Every replicate is simulated AROUND # ``theta_hat``; if the base fit ended on the model's magnitude diff --git a/src/liouscope/fitting/gls.py b/src/liouscope/fitting/gls.py index 07b04446..74406439 100644 --- a/src/liouscope/fitting/gls.py +++ b/src/liouscope/fitting/gls.py @@ -48,6 +48,14 @@ class GLSFitOutput: #: may carry variation and the optimiser may have run, but exact-zero RSS #: (or an unrepresentable positive MLE scale) cannot support AICc/CI claims. likelihood_degenerate: bool = False + #: True when the profile likelihood IS computable in log space but the + #: positive MLE scale ``exp(log_sigma)`` is not representable as a float64 + #: (round-1 review of PR #147). ``sigma`` is then NaN while + #: ``log_likelihood`` stays finite: the fit remains a valid AICc candidate, + #: and only the scale-dependent evidence -- the parametric bootstrap, hence + #: the CI -- is withheld. Distinct from ``likelihood_degenerate``, where the + #: likelihood itself has no finite value. + scale_unavailable: bool = False def _whiten(y: np.ndarray, rho: float) -> np.ndarray: @@ -224,23 +232,44 @@ def residual(params: np.ndarray, rho_local: float = rho) -> np.ndarray: sigma = float(math.exp(log_sigma)) except OverflowError: sigma = float("inf") - if not math.isfinite(sigma) or sigma <= 0.0: + scale_unavailable = not math.isfinite(sigma) or sigma <= 0.0 + if scale_unavailable: + # ROUND-1 REVIEW (PR #147). This branch used to return + # ``success=False, log_likelihood=nan, likelihood_degenerate=True``, + # which withheld the whole FIT because a number the likelihood never + # needs could not be materialised. ``log_rss`` is finite here -- the + # gate above already refused the case where it is not -- so the profile + # likelihood is computable directly in log space by + # ``gaussian_log_likelihood``, which never forms ``sigma``. Only + # ``exp(log_sigma)`` left the float64 range. + # + # Measured on the reviewer's construction: one minimum-subnormal + # residual (5e-324) in an otherwise-zero series of 64 points gives + # ``log_rss = -1488.88``, ``log_sigma = -746.52`` and a perfectly + # finite profile log-likelihood of ``+47686.44`` -- while + # ``exp(-746.52)`` underflows to ``0.0`` because the smallest + # subnormal is ``exp(-744.44)``. Marking the fit unsuccessful there + # makes ``aicc()`` return ``inf`` in ``_fit_one`` and the model drops + # out of selection, which reintroduces exactly the ABSOLUTE SCALE + # BOUNDARY into model selection that issue #135 removed: the same + # curve in different rate units is or is not a candidate. + # + # So the fit stays selectable and only the scale-dependent evidence is + # withheld. ``sigma`` is NaN, never 0.0 or inf: it is consumed by + # ``_ar1_resample`` as the standard deviation of the innovation, where + # 0.0 would silently generate a bootstrap of IDENTICAL replicates -- + # a zero-width CI, which is the failure mode of an uncertainty + # pipeline, not a wide one. ``parametric_bootstrap`` refuses on the + # dedicated flag rather than on ``success``, so "no interval" and "no + # estimate" stay distinguishable. warnings.warn( - "fit_gls_ar1: positive residual MLE scale is not representable as " - "float64; likelihood/AICc/CI evidence is unavailable (issue #135).", + "fit_gls_ar1: the positive residual MLE scale is not representable " + f"as float64 (log sigma = {log_sigma:.6g}); the log-space " + "likelihood and AICc remain available, but bootstrap/CI evidence " + "is withheld (issue #135, PR #147 review).", RuntimeWarning, stacklevel=2, ) - return GLSFitOutput( - params=p, - residuals=residuals_raw, - rho_ar1=rho, - sigma=float("nan"), - log_likelihood=float("nan"), - success=False, - saturated=tuple(sorted(fired)), - likelihood_degenerate=True, - ) # Prais-Winsten exact AR(1) likelihood: _whiten keeps observation 0 # (scaled by sqrt(1-rho^2)), so the transform has log-Jacobian @@ -252,8 +281,9 @@ def residual(params: np.ndarray, rho_local: float = rho) -> np.ndarray: params=p, residuals=residuals_raw, rho_ar1=rho, - sigma=sigma, + sigma=float("nan") if scale_unavailable else sigma, log_likelihood=log_lik, success=success, saturated=tuple(sorted(fired)), + scale_unavailable=scale_unavailable, ) diff --git a/tests/test_issue135_likelihood_scale.py b/tests/test_issue135_likelihood_scale.py index 4e9af428..40673b2c 100644 --- a/tests/test_issue135_likelihood_scale.py +++ b/tests/test_issue135_likelihood_scale.py @@ -3,6 +3,7 @@ from __future__ import annotations import math +import warnings import numpy as np import pytest @@ -111,3 +112,99 @@ def test_likelihood_degenerate_state_reaches_fitresult_and_is_nonselectable(monk assert not fit_result.success assert fit_result.likelihood_degenerate assert np.isinf(fit_result.aicc) + + +# -------------------------------------------------------------------------- +# PR #147 round-1 review: an unrepresentable MLE scale withheld the whole FIT, +# not just the scale-dependent evidence. +# -------------------------------------------------------------------------- + +#: Smallest positive float64. One such residual among an otherwise-zero series +#: keeps ``log_rss`` finite while ``exp(0.5 * (log_rss - log n))`` underflows: +#: the profile likelihood is computable, its scale is not. +_MIN_SUBNORMAL = 5.0e-324 + + +def _constant_model(t: np.ndarray, p: np.ndarray) -> np.ndarray: + return np.full(t.shape, float(p[0])) + + +def _underflowing_scale_case() -> tuple[np.ndarray, np.ndarray]: + n = 64 + t = np.linspace(0.0, 5.0, n) + y = np.zeros(n) + y[7] = _MIN_SUBNORMAL + return t, y + + +def test_the_underflowing_case_really_has_a_finite_log_likelihood() -> None: + """Positive control: the fixture must exercise the finding, not a NaN RSS. + + Without this the two assertions below could pass on a case where nothing + was computable in the first place. + """ + _, y = _underflowing_scale_case() + log_rss = scaled_log_sum_squares(y) + assert math.isfinite(log_rss), log_rss + log_sigma = 0.5 * (log_rss - math.log(y.size)) + assert math.isfinite(log_sigma), log_sigma + assert math.exp(log_sigma) == 0.0, ( + "the fixture no longer underflows; the finding cannot be reached" + ) + assert math.isfinite(gaussian_log_likelihood(y)) + + +def test_unrepresentable_scale_keeps_the_fit_selectable() -> None: + """The fit stays an AICc candidate; only ``sigma`` is withheld. + + Measured before the repair: ``success=False``, ``log_likelihood=nan`` and + ``likelihood_degenerate=True`` -- so ``_fit_with_model`` scored the model + ``inf`` and dropped it from selection. That reintroduces an ABSOLUTE scale + boundary into model selection, which is the boundary issue #135 removed: + the same curve in different rate units is or is not a candidate. + """ + t, y = _underflowing_scale_case() + # Recorded rather than ``pytest.warns``: a missing warning must fail this + # test at an ASSERTION, so that removing the guard is attributable. With + # ``pytest.warns`` the test dies inside pytest's own context manager, which + # a mutation run cannot tell apart from an incidental crash. + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + fit = fit_gls_ar1(_constant_model, t, y, np.array([0.0])) + assert any( + issubclass(w.category, RuntimeWarning) + and "scale is not representable" in str(w.message) + for w in caught + ), [str(w.message) for w in caught] + assert fit.success, "the fit was withheld because its SCALE could not be built" + assert math.isfinite(fit.log_likelihood), fit.log_likelihood + assert not fit.likelihood_degenerate + assert fit.scale_unavailable + assert math.isnan(fit.sigma), ( + f"sigma is {fit.sigma!r}; 0.0 would make _ar1_resample draw identical " + "replicates and produce a zero-width interval" + ) + + +def test_unrepresentable_scale_still_withholds_the_bootstrap() -> None: + """Fail-closed control: selectable must not mean an interval was invented.""" + t, y = _underflowing_scale_case() + with pytest.warns(RuntimeWarning): + with pytest.raises(Exception) as excinfo: + parametric_bootstrap( + _constant_model, t, y, np.array([0.0]), B=4, + rng=np.random.default_rng(0), + ) + assert excinfo.type is RuntimeError, excinfo.type + assert "no representable residual scale" in str(excinfo.value) + + +def test_ordinary_fit_is_untouched_by_the_scale_repair() -> None: + """Over-correction control: a normal fit keeps a finite positive sigma.""" + rng = np.random.default_rng(20260903) + t = np.linspace(0.0, 4.0, 64) + p = np.array([1.0, 0.8]) + y = M0(t, p) + 1.0e-3 * rng.standard_normal(t.size) + fit = fit_gls_ar1(M0, t, y, p, n_iters=1) + assert not fit.scale_unavailable + assert np.isfinite(fit.sigma) and fit.sigma > 0.0 From 9808a9a00336a7deeb441412d18c1d021a34385b Mon Sep 17 00:00:00 2001 From: Codie Date: Thu, 3 Sep 2026 02:33:11 +0200 Subject: [PATCH 19/22] fix(#147): a withheld scale the report could not name, and an octave 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) Claude-Session: https://claude.ai/code/session_01FPKDUXBV4qXpWX8ZtpqgAZ --- CHANGELOG.md | 22 ++++++ src/liouscope/_types.py | 9 +++ src/liouscope/diagnostics/relaxation.py | 9 +++ src/liouscope/fitting/aicc.py | 29 ++++++-- tests/test_issue135_likelihood_scale.py | 93 +++++++++++++++++++++++++ 5 files changed, 158 insertions(+), 4 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 3375ce73..9804c9ce 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -35,6 +35,28 @@ adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). the change is not visible in `input_hash`. ### Fixed +- **The reason a confidence interval was withheld did not reach the report (PR + #147, round-2 external review).** When the residual MLE scale is not + representable as float64 the fit stays a valid AICc candidate and only its + interval is withheld — but `_fit_with_model` copied `likelihood_degenerate` + and dropped `scale_unavailable`, so the persisted report showed only + `bca_ci_beta = (nan, nan)`. That is the value ANY bootstrap or jackknife + failure produces, and the distinguishing information lived in a + `RuntimeWarning` an artefact does not keep. `FitResult` carries + `scale_unavailable` now and `_fit_with_model` copies it; `io.export` + serialises `FitResult` field-wise, so the reason travels into dumped + reports. +- **The one-half factor was applied after the overflow guard (PR #147, round-2 + external review).** For the explicit-`sigma` path the standardised RSS never + appears on its own: it enters the log-likelihood as `-0.5 * RSS`, so the + representable range reaches `2 * float64.max`. The guard compared against + `log(float64.max)` and returned `-inf` for the octave above it — measured: + `sigma = 1` with a single residual near `1.4e154` has RSS ~`1.96e308` and a + finite log-likelihood of ~`-9.8e307`, and the model was dropped from + selection on that arithmetic boundary rather than on the data. The bound now + carries `+ log(2)` and the exponential is taken after subtracting `log(2)` + in the octave that needs it; below that octave the arithmetic is unchanged, + so no previously computable likelihood moves by even one ulp. - **A fit was withheld entirely because a number its likelihood never uses could not be materialised (PR #147, round-1 review).** When the whitened residuals have a finite, non-zero RMS below the smallest positive float64, diff --git a/src/liouscope/_types.py b/src/liouscope/_types.py index 2005292e..67ddbf22 100644 --- a/src/liouscope/_types.py +++ b/src/liouscope/_types.py @@ -188,6 +188,15 @@ class FitResult: residual_ar1_rho: float success: bool likelihood_degenerate: bool = False + #: True when the fit is a valid model-selection candidate -- its log-space + #: likelihood and AICc are finite -- but the positive MLE residual scale is + #: not representable as float64, so the parametric bootstrap and therefore + #: the confidence interval were withheld (issue #135, PR #147 round-2 + #: review). Without this field the persisted report showed only + #: ``bca_ci_beta = (nan, nan)``, which is the same value ANY bootstrap or + #: jackknife failure produces: once the warning stream is gone, the reason + #: the interval is missing was unrecoverable from the artefact. + scale_unavailable: bool = False @dataclass(frozen=True, slots=True, kw_only=True) diff --git a/src/liouscope/diagnostics/relaxation.py b/src/liouscope/diagnostics/relaxation.py index 29d334fd..666ddddb 100644 --- a/src/liouscope/diagnostics/relaxation.py +++ b/src/liouscope/diagnostics/relaxation.py @@ -216,6 +216,15 @@ def _fit_with_model( residual_ar1_rho=fit.rho_ar1, success=fit.success, likelihood_degenerate=fit.likelihood_degenerate, + # ROUND-2 REVIEW (PR #147). This conversion carried + # ``likelihood_degenerate`` and dropped ``scale_unavailable``, so the + # newly supported underflow case -- ``success=True`` with a withheld + # ``sigma`` -- reached the report indistinguishable from an ordinary + # fit. Its bootstrap then fails and ``compute_relaxation_layer`` + # records ``bca_ci_beta = (nan, nan)``, exactly what a jackknife + # failure or a non-converged resample produces. The reason lived only + # in a RuntimeWarning, which a persisted artefact does not keep. + scale_unavailable=fit.scale_unavailable, ) return fit_result, fit.params diff --git a/src/liouscope/fitting/aicc.py b/src/liouscope/fitting/aicc.py index c40fb40b..89db3727 100644 --- a/src/liouscope/fitting/aicc.py +++ b/src/liouscope/fitting/aicc.py @@ -63,16 +63,37 @@ def gaussian_log_likelihood( return float("nan") log_sigma = math.log(sigma) if log_rss == float("-inf"): - standardised_rss = 0.0 + half_standardised_rss = 0.0 elif log_rss == float("inf"): return float("-inf") else: log_standardised_rss = log_rss - 2.0 * log_sigma - if log_standardised_rss > math.log(np.finfo(float).max): + # ROUND-2 REVIEW (PR #147). The standardised RSS never appears on its + # own: it enters the result as ``-0.5 * RSS``. Its representable range + # therefore reaches ``2 * float64.max``, and comparing against + # ``log(float64.max)`` refused an entire octave of perfectly finite + # likelihoods. Measured: ``sigma = 1`` with a single residual near + # 1.4e154 has RSS ~1.96e308 and a log-likelihood of ~-9.8e307, and + # this guard returned ``-inf`` -- which drops the model from AICc + # selection on an arithmetic accident rather than on the data, the + # same absolute-scale boundary issue #135 exists to remove. + # + # The one-half factor is now applied BEFORE the decision: the bound + # carries ``+ log(2)`` and, in the octave that only the halved value + # can represent, the exponential is taken after subtracting + # ``log(2)``. Below that octave the arithmetic is left exactly as it + # was (``0.5 * exp(...)``), so no previously computed likelihood + # changes by even one ulp -- widening a range must not perturb the + # values already inside it. + log_max = math.log(np.finfo(float).max) + if log_standardised_rss > log_max + math.log(2.0): return float("-inf") - standardised_rss = math.exp(log_standardised_rss) + if log_standardised_rss > log_max: + half_standardised_rss = math.exp(log_standardised_rss - math.log(2.0)) + else: + half_standardised_rss = 0.5 * math.exp(log_standardised_rss) return float( - -0.5 * n * log_2pi - n * log_sigma - 0.5 * standardised_rss + -0.5 * n * log_2pi - n * log_sigma - half_standardised_rss ) diff --git a/tests/test_issue135_likelihood_scale.py b/tests/test_issue135_likelihood_scale.py index 40673b2c..968bcb89 100644 --- a/tests/test_issue135_likelihood_scale.py +++ b/tests/test_issue135_likelihood_scale.py @@ -13,6 +13,7 @@ from liouscope.fitting.bootstrap import parametric_bootstrap from liouscope.fitting.gls import GLSFitOutput, fit_gls_ar1 from liouscope.fitting.models import M0 +from liouscope.io.export import _to_jsonable from liouscope.numerics.norms import scaled_log_sum_squares @@ -208,3 +209,95 @@ def test_ordinary_fit_is_untouched_by_the_scale_repair() -> None: fit = fit_gls_ar1(M0, t, y, p, n_iters=1) assert not fit.scale_unavailable assert np.isfinite(fit.sigma) and fit.sigma > 0.0 + + +# -------------------------------------------------------------------------- +# PR #147 round-2 review +# -------------------------------------------------------------------------- + + +def test_scale_unavailable_reaches_fitresult() -> None: + """The reason a CI was withheld must survive into the persisted report. + + ``bca_ci_beta = (nan, nan)`` is what EVERY bootstrap or jackknife failure + produces. Without this flag on ``FitResult`` an unrepresentable residual + scale -- a fit that is otherwise successful and selectable -- was + indistinguishable from a non-converged resample once the warning stream + was gone. + """ + t, y = _underflowing_scale_case() + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + fit_result, _ = relaxation_mod._fit_with_model("M0", t, y) + assert any( + "scale is not representable" in str(w.message) for w in caught + ), [str(w.message) for w in caught] + assert fit_result.scale_unavailable + # The fit stays a candidate: the flag records why the INTERVAL is missing, + # it does not withhold the estimate (round-1 review contract). + assert fit_result.success + assert math.isfinite(fit_result.log_likelihood) + # The field is serialised, so an audit artefact carries the reason. + assert _to_jsonable(fit_result)["scale_unavailable"] is True + + +def test_ordinary_fit_reaches_fitresult_without_the_flag() -> None: + """Over-correction control: a normal fit must not be labelled scaleless.""" + rng = np.random.default_rng(20260903) + t = np.linspace(0.0, 4.0, 64) + y = M0(t, np.array([1.0, 0.8])) + 1.0e-3 * rng.standard_normal(t.size) + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + fit_result, _ = relaxation_mod._fit_with_model("M0", t, y) + assert not fit_result.scale_unavailable + + +def test_explicit_sigma_likelihood_uses_the_full_representable_rss_range() -> None: + """``-0.5 * RSS`` is representable up to an RSS of ``2 * float64.max``. + + The guard compared the standardised RSS against ``float64.max`` and + returned ``-inf`` for the octave above it, although the value the formula + actually forms -- the HALVED RSS -- is perfectly finite there. Dropping a + model from selection on that boundary is the absolute-scale dependence + issue #135 exists to remove. + """ + residual = np.array([1.4e154]) + log_lik = gaussian_log_likelihood(residual, sigma=1.0) + assert math.isfinite(log_lik), log_lik + expected = -0.5 * math.log(2.0 * math.pi) - math.exp( + 2.0 * math.log(1.4e154) - math.log(2.0) + ) + assert log_lik == pytest.approx(expected, rel=1e-12) + + # POSITIVE CONTROL that the fixture really sits in the disputed octave: + # the raw standardised RSS is NOT representable, only its half is. + log_rss = scaled_log_sum_squares(residual) + log_max = math.log(np.finfo(float).max) + assert log_max < log_rss <= log_max + math.log(2.0), log_rss + + +def test_a_genuinely_unrepresentable_half_rss_still_returns_minus_inf() -> None: + """Fail-closed control: the bound moved by one octave, it did not vanish.""" + residual = np.array([1.4e308]) + log_rss = scaled_log_sum_squares(residual) + assert log_rss > math.log(np.finfo(float).max) + math.log(2.0) + assert gaussian_log_likelihood(residual, sigma=1.0) == float("-inf") + + +def test_values_inside_the_old_range_are_bit_identical() -> None: + """Widening a range must not perturb the values already inside it. + + The halving stays ``0.5 * exp(...)`` below the disputed octave, so every + likelihood the previous implementation could compute is reproduced exactly + rather than to within an ulp. + """ + for scale in (1.0e-30, 1.0, 3.7, 1.0e30, 1.0e120): + residuals = scale * np.array([0.3, -1.1, 2.0, 0.0]) + log_rss = scaled_log_sum_squares(residuals) + log_standardised = log_rss # sigma = 1 + assert log_standardised <= math.log(np.finfo(float).max) + n = residuals.size + expected = ( + -0.5 * n * math.log(2.0 * math.pi) - 0.5 * math.exp(log_standardised) + ) + assert gaussian_log_likelihood(residuals, sigma=1.0) == expected From cee633f9a48c5ddd4667274850d131448d049686 Mon Sep 17 00:00:00 2001 From: Codie Date: Thu, 3 Sep 2026 02:34:46 +0200 Subject: [PATCH 20/22] test(#147): make the fail-closed control die of an assertion, not an 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) Claude-Session: https://claude.ai/code/session_01FPKDUXBV4qXpWX8ZtpqgAZ --- tests/test_issue135_likelihood_scale.py | 20 ++++++++++++++++++-- 1 file changed, 18 insertions(+), 2 deletions(-) diff --git a/tests/test_issue135_likelihood_scale.py b/tests/test_issue135_likelihood_scale.py index 968bcb89..41a74cb1 100644 --- a/tests/test_issue135_likelihood_scale.py +++ b/tests/test_issue135_likelihood_scale.py @@ -277,11 +277,27 @@ def test_explicit_sigma_likelihood_uses_the_full_representable_rss_range() -> No def test_a_genuinely_unrepresentable_half_rss_still_returns_minus_inf() -> None: - """Fail-closed control: the bound moved by one octave, it did not vanish.""" + """Fail-closed control: the bound moved by one octave, it did not vanish. + + The outcome is captured rather than asserted inline. Without the guard the + call raises ``OverflowError`` from ``math.exp``, and a test that dies of an + exception is indistinguishable from an incidental crash in a mutation run -- + the death must be attributable to an ASSERTION for the proof to count. + """ residual = np.array([1.4e308]) log_rss = scaled_log_sum_squares(residual) assert log_rss > math.log(np.finfo(float).max) + math.log(2.0) - assert gaussian_log_likelihood(residual, sigma=1.0) == float("-inf") + + outcome: object + try: + outcome = gaussian_log_likelihood(residual, sigma=1.0) + except Exception as exc: + outcome = exc + assert not isinstance(outcome, BaseException), ( + f"the likelihood raised {type(outcome).__name__} " + f"({outcome}) instead of returning -inf" + ) + assert outcome == float("-inf") def test_values_inside_the_old_range_are_bit_identical() -> None: From 6507277a53c2fbab1ecd06f1a2d8d70a8ffdffac Mon Sep 17 00:00:00 2001 From: marcohost33-maker Date: Thu, 3 Sep 2026 09:46:12 +0200 Subject: [PATCH 21/22] fix(#147 round-3): decide the halved-RSS boundary on the value, not on 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) Claude-Session: https://claude.ai/code/session_01FPKDUXBV4qXpWX8ZtpqgAZ --- src/liouscope/fitting/aicc.py | 34 ++++++++- tests/test_issue135_likelihood_scale.py | 98 +++++++++++++++++++++++++ 2 files changed, 129 insertions(+), 3 deletions(-) diff --git a/src/liouscope/fitting/aicc.py b/src/liouscope/fitting/aicc.py index 89db3727..6c5d42a9 100644 --- a/src/liouscope/fitting/aicc.py +++ b/src/liouscope/fitting/aicc.py @@ -85,11 +85,39 @@ def gaussian_log_likelihood( # was (``0.5 * exp(...)``), so no previously computed likelihood # changes by even one ulp -- widening a range must not perturb the # values already inside it. + # ROUND-3 REVIEW (PR #147). The round-2 bound was still a decision + # made in LOG SPACE, and at this magnitude log space cannot make it. + # Measured with ``sigma = 1`` and a single residual + # ``sqrt(float_max) * sqrt(2)``: ``log_rss`` comes out at + # 710.475860073944 while ``log_max + log(2)`` rounds to + # 710.4758600739439 -- ONE ulp lower -- so the guard returned + # ``-inf`` although the mathematically relevant half-RSS is + # 1.7976931348623155e+308, a perfectly finite 0.9999999999999999 of + # ``float_max``. Nor can the halved value be recovered by subtracting + # ``log(2)`` and exponentiating: one ulp of a logarithm near 710 is a + # factor of ~1e-16 in the value, which straddles the overflow edge + # exactly where this test is made, and ``math.exp`` raises + # ``OverflowError`` on the same input (measured). + # + # So the halved quantity is MATERIALISED from the residuals instead + # of inferred from their logarithm. Halving each term BEFORE + # accumulating keeps every partial sum inside float64 whenever the + # true half-RSS is inside it, and lets the float64 addition itself + # decide the boundary -- there is no rounded bound left to be one ulp + # wrong about. An overflow here is now a measurement rather than a + # prediction: ``inf`` means the half-RSS genuinely does not fit + # (verified for 4x, 1e3x float_max, and for 64 residuals whose halved + # sum overflows), and only then is the fit dropped. + # + # The lower octave is untouched, so no previously computed likelihood + # changes by even one ulp -- the same promise round 2 made. log_max = math.log(np.finfo(float).max) - if log_standardised_rss > log_max + math.log(2.0): - return float("-inf") if log_standardised_rss > log_max: - half_standardised_rss = math.exp(log_standardised_rss - math.log(2.0)) + scaled = residuals / sigma + with np.errstate(over="ignore", invalid="ignore"): + half_standardised_rss = float(np.sum((0.5 * scaled) * scaled)) + if not math.isfinite(half_standardised_rss): + return float("-inf") else: half_standardised_rss = 0.5 * math.exp(log_standardised_rss) return float( diff --git a/tests/test_issue135_likelihood_scale.py b/tests/test_issue135_likelihood_scale.py index 41a74cb1..8231541a 100644 --- a/tests/test_issue135_likelihood_scale.py +++ b/tests/test_issue135_likelihood_scale.py @@ -317,3 +317,101 @@ def test_values_inside_the_old_range_are_bit_identical() -> None: -0.5 * n * math.log(2.0 * math.pi) - 0.5 * math.exp(log_standardised) ) assert gaussian_log_likelihood(residuals, sigma=1.0) == expected + + +# -------------------------------------------------------------------------- +# PR #147 round-3 review +# -------------------------------------------------------------------------- + + +def test_the_halved_boundary_is_decided_on_the_materialised_value() -> None: + """Round 2 moved the bound; round 3 removes the bound. + + The round-2 guard compared ``log_rss`` against ``log_max + log(2)``, two + separately rounded logarithms. On the reviewer's fixture -- ``sigma = 1`` + with a single residual ``sqrt(float_max) * sqrt(2)`` -- ``log_rss`` lands + exactly ONE ulp above that bound and the fit was dropped, although the + quantity the formula actually forms is 0.9999999999999999 of + ``float_max`` and entirely representable. + + The repair is not a wider bound. It is materialising the halved sum from + the residuals and letting float64 addition decide, so there is no rounded + threshold left to be one ulp wrong about. + """ + residual = np.array([math.sqrt(np.finfo(float).max) * math.sqrt(2.0)]) + + # PREMISE, asserted rather than assumed: the fixture is in the octave the + # round-2 bound rejected, and it misses that bound by a single ulp. + log_rss = scaled_log_sum_squares(residual) + log_max = math.log(np.finfo(float).max) + round2_bound = log_max + math.log(2.0) + assert log_rss > round2_bound, "fixture no longer trips the round-2 bound" + assert (log_rss - round2_bound) == pytest.approx( + math.ulp(log_max), rel=1e-9 + ), "the fixture no longer sits one ulp above the bound" + + log_lik = gaussian_log_likelihood(residual, sigma=1.0) + assert math.isfinite(log_lik), log_lik + + # The value is the one the formula defines, computed the direct way. + scaled = residual / 1.0 + half = float(np.sum((0.5 * scaled) * scaled)) + assert math.isfinite(half) + assert half == pytest.approx(0.9999999999999999 * np.finfo(float).max, rel=1e-15) + assert log_lik == pytest.approx(-0.5 * math.log(2.0 * math.pi) - half, rel=1e-15) + + +def test_the_log_space_route_could_not_have_reached_it() -> None: + """Why the fix is not "subtract log(2) and exponentiate". + + One ulp of a logarithm near 710 is a factor of ~1e-16 in the value, and + that is precisely the width of the decision being made here. Documented + as a live assertion because the obvious smaller repair looks correct and + is not: ``math.exp`` raises on the same input for which the direct + computation returns a finite number. + """ + residual = np.array([math.sqrt(np.finfo(float).max) * math.sqrt(2.0)]) + log_rss = scaled_log_sum_squares(residual) + + raised: BaseException | None = None + try: + math.exp(log_rss - math.log(2.0)) + except Exception as exc: + raised = exc + assert isinstance(raised, OverflowError), ( + "the log-space route no longer overflows; the round-3 rationale " + f"needs rereading (got {raised!r})" + ) + assert math.isfinite(gaussian_log_likelihood(residual, sigma=1.0)) + + +def test_a_half_rss_that_truly_overflows_is_still_dropped() -> None: + """Fail-closed control: removing the bound must not remove the refusal. + + Three shapes of genuine overflow, including one that only overflows in + the SUM -- no single term does -- which a per-term bound would have + missed. Captured and asserted rather than crashed into, so a mutation run + can attribute the death to an assertion. + """ + cases = { + "4x float_max, one residual": np.array( + [math.sqrt(np.finfo(float).max) * 2.0] + ), + "1e3x float_max, one residual": np.array( + [math.sqrt(np.finfo(float).max) * math.sqrt(1000.0)] + ), + "overflows only in the sum": np.full( + 64, math.sqrt(np.finfo(float).max) * math.sqrt(2.0) / 4.0 + ), + } + for label, residuals in cases.items(): + outcome: object + try: + outcome = gaussian_log_likelihood(residuals, sigma=1.0) + except Exception as exc: + outcome = exc + assert not isinstance(outcome, BaseException), ( + f"{label}: raised {type(outcome).__name__} ({outcome}) " + "instead of returning -inf" + ) + assert outcome == float("-inf"), label From 77b3aceaaffe6e47ae27ecfbe0220d41965df298 Mon Sep 17 00:00:00 2001 From: marcohost33-maker Date: Thu, 3 Sep 2026 09:49:03 +0200 Subject: [PATCH 22/22] fix(#147 round-3): test the halved RSS for OVERFLOW, not for non-finiteness 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) Claude-Session: https://claude.ai/code/session_01FPKDUXBV4qXpWX8ZtpqgAZ --- src/liouscope/fitting/aicc.py | 14 +++++++++++++- 1 file changed, 13 insertions(+), 1 deletion(-) diff --git a/src/liouscope/fitting/aicc.py b/src/liouscope/fitting/aicc.py index 6c5d42a9..f0781649 100644 --- a/src/liouscope/fitting/aicc.py +++ b/src/liouscope/fitting/aicc.py @@ -116,7 +116,19 @@ def gaussian_log_likelihood( scaled = residuals / sigma with np.errstate(over="ignore", invalid="ignore"): half_standardised_rss = float(np.sum((0.5 * scaled) * scaled)) - if not math.isfinite(half_standardised_rss): + # Tested for OVERFLOW specifically, not for non-finiteness. + # The mutation run reported a blanket isfinite() check BLIND, and + # measuring why exposed a latent state collapse: for an infinite + # half-RSS the check is provably equivalent (the subtraction + # below already yields -inf), and its only remaining effect would + # be to turn a NaN into -inf -- converting "model-selection + # likelihood unavailable", which this function returns NaN for by + # contract, into "definitively dropped". Those are different + # states. NaN is unreachable here (a NaN residual is caught by + # the isnan(log_rss) gate above, and every term 0.5*s*s is + # non-negative), so this is guarding an invariant rather than a + # path -- which is the reason to state it precisely. + if half_standardised_rss == float("inf"): return float("-inf") else: half_standardised_rss = 0.5 * math.exp(log_standardised_rss)