Skip to content

fix(sysid): fit_lm reads a small jump of the residual from a long rejected step #258

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
NicholasEhsanRoy merged 3 commits into from
Oct 7, 2026
Merged

fix(sysid): fit_lm reads a small jump of the residual from a long rejected step #258

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
2e6f3ac
fix(sysid): fit_lm reads a small jump of the residual from a long rej…
NicholasEhsanRoy Oct 7, 2026
8879b31
test(sysid): the spec strategy draws no logit bound a float32 leaf ca…
NicholasEhsanRoy Oct 7, 2026
09e10a1
test(sysid): the ball pin asserts the verdict, not how many candidate…
NicholasEhsanRoy Oct 7, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions CHANGELOG.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -418,6 +418,7 @@ guidance; the itemized changes follow.
The `[verify]` extra now only pulls `hypothesis`.

### Fixed
- **`fit_lm` no longer reports `converged=True` beside a small jump of the residual read from a long rejected step** (MADD-ANO-216's residual case, never released): the one-sided test divided the jump by the linearisation's change over the whole step, so a late bounce 22 float spacings from the iterate read 341 and 456 against a threshold of `2**10` on the stock ball. A rejected step reading above `2**5` is now read again across the one float spacing where the residual departs most (a few residual evaluations, at the end of the run only); the run ends `converged=False` with the existing warning. Action: none.
- **A coupling group under `acceleration="iqn-ils"` or `"iqn-imvj"` returns from `step()` when its iterate leaves float range** (MADD-ANO-233; 0.3.0 and 0.3.1 under `solver="ift"` with `"iqn-imvj"`): the secant least-squares handed LAPACK's SVD a matrix holding an `inf`, and the step never returned. It now runs to `max_iterations` and reports `converged=False`, `residual=inf`, as the other accelerations do; steps of a finite run are unchanged, bit for bit. Action: none; on 0.3.x use `solver="fori"` or `"aitken"`.
- **`gradient_relative_error_bound` covers the constants one pass resolves** (MADD-ANO-234, never released): with a constant in it whose tangent at the returned iterate is rounding (the centre or curve of a nonlinearity evaluated on its centre) it read 1.09 for a relative error of 1.134, usable. Such a constant, and a gain whose whole value moves the pass by less than the residual's float floor, is no longer in the bound; the docstring says which constants are. Action: none.
- **`GraphManager.step`, `run` and `run_adaptive` no longer store a step that reshapes a state leaf** (`MADD-ANO-220`, in 0.1.0 to 0.3.1): a list given for a scalar constant (`BallNode(initial_velocity=[1.0, 2.0])`), or an external input of another shape at a later step, broadcast a scalar leaf, and the state's checkpoint then did not load after `reset_state()`. The step is refused by name (compared once per trace, host-side: no compiled program changes), as `run_scan` always did; the FMU sidecar and bridge refuse it too.
Expand Down
37 changes: 34 additions & 3 deletions docs/validation/known_anomalies.yaml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14716,6 +14716,32 @@ anomalies:
the 90 random starts converged above the minimum, and a targeted
search of 800 drawn fits of the ball finds no fit converged where the
loss still falls along the piece it is on.

A residual case, found by the slow random hunt on 2026-10-07 (CPU,
jaxlib 0.10.2; the same to the last digit on 0.11.0) and fixed
before any release: a small jump read from a long candidate. The
ratio divides the jump by the linearisation's change over the
candidate, which grows with its length, and a late bounce that moves
by one time step changes few samples. Two fits of the ball (the
hunt's, and one of 4,500 uniformly drawn) ended `converged=True` at
a loss of 0.006 with ratios of 456 and 341 -- a jump of 0.015 in the
residual's norm, 22 float spacings of the elasticity from the
iterate, read from a candidate 35 spacings long -- where the loss
still fell along the piece by 11.8 and 5.6 times what its rounding
explains. A candidate above `2**5` and not above `2**10` is now
read a second time: its step is halved down to the one float spacing
of the coordinates across which the residual departs most from its
linearisation, and it is a jump where that departure is `2**8` times
its mirror image's (the whole candidate's, as before) plus the
linearisation's change over one spacing. Measured: 2.0e3 and 1.8e3
at those two endings and 3.4e4 or more at the other 261 jump endings
of the 4,500 draws; at most 52 on smooth fits stopped by the floor
rule (140 endings, none of which is asked: their first reading is
at most 21). After the fix: no fit converged
where the loss still falls along its piece among the 4,500 draws or
in the hunt on three seeds (800 examples each), on either jaxlib;
of the 250 endings that existing per-push and slow tests put to
the floor rule (140 of smooth fits), none changed.
severity: "major"
safety_relevance: "context_dependent"
safety_relevance_rationale: >
Expand Down Expand Up @@ -14743,9 +14769,11 @@ anomalies:
test passes because it is stationary on its own piece; and a fit
converged at the interior minimum of a smooth piece that is not the
lowest (a local minimum, with the lower piece as near as 1e-3 of a
parameter on the ball). The threshold is measured, not derived: a
jump smaller than `2**10` times the linear change of a step of a few
float spacings reads as rounding. The root cause, the missing
parameter on the ball). The thresholds are measured, not derived:
a jump reads as rounding where it is smaller than `2**5` times the
linear change of the rejected step that crossed it, or than `2**8`
times what the residual departs from its linearisation by at that
step's mirror image. The root cause, the missing
event-time derivative, is MADD-ANO-021 and stays open.
verification:
- "tests/core/test_sysid_non_differentiable_residual.py::test_fit_lm_is_converged_on_the_ball_only_at_the_minimum"
Expand All @@ -14755,6 +14783,9 @@ anomalies:
- "tests/core/test_sysid_non_differentiable_residual.py::test_an_iteration_whose_damped_candidates_never_moved_is_asked_again"
- "tests/property/test_sysid_targeted_search.py::test_no_wrong_fit_among_the_fixed_draws"
- "tests/property/test_sysid_targeted_search.py::test_fit_lm_on_the_ball_is_not_converged_where_the_loss_falls_along_its_piece"
- "tests/property/test_sysid_targeted_search.py::test_fit_lm_reads_a_small_jump_of_the_ball_from_a_long_candidate"
- "tests/core/test_sysid_non_differentiable_residual.py::test_a_small_jump_inside_a_long_candidate_is_read_across_one_spacing"
- "tests/core/test_sysid_non_differentiable_residual.py::test_a_smooth_residual_read_a_second_time_is_still_not_a_jump"
github_issue: null

- anomaly_id: "MADD-ANO-217"
Expand Down
15 changes: 13 additions & 2 deletions docs/validation/sysid_fmu_claims.yaml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2624,7 +2624,10 @@ claims:
linearisation on the candidate's side than on the other side and than
the linearisation's own change, the rejection is a jump and not
rounding ... The run then ends converged=False with a RuntimeWarning
naming a residual that is not differentiable"; an iteration "whose
naming a residual that is not differentiable"; "a candidate at more
than 2**5 is read a second time ... and it is a jump where that
departure is 2**8 times the other side's and the linearisation's
change over one spacing"; an iteration "whose
ladder began above lam0 while the undamped step was rejected ...
Where the undamped step promised, on the linear model, 2**10 times
more than the short candidates lost to rounding -- or every candidate
Expand All @@ -2649,7 +2652,8 @@ claims:
an event, a contact or a valve whose timing depends on a trained
parameter, MADD-ANO-021): converged=False and the warning where a
rejected candidate within 2**-10 of a parameter shows the jump against
its mirror image (the stock TableNode -> BallNode graph, float32; that
its mirror image, over its whole length or across one float spacing
inside it (the stock TableNode -> BallNode graph, float32; that
graph does not run under x64); not detected, and not claimed: a kink
(a continuous residual whose slope jumps), a jump whose mirror image
the bounds do not allow, a jump beside an iterate the proposal test
Expand Down Expand Up @@ -2679,6 +2683,13 @@ claims:
- tests/core/test_sysid_non_differentiable_residual.py::test_a_candidate_with_no_mirror_image_or_no_finite_residual_says_nothing
- tests/core/test_sysid_non_differentiable_residual.py::test_an_iteration_whose_damped_candidates_never_moved_is_asked_again
- tests/core/test_sysid_non_differentiable_residual.py::test_candidates_that_never_move_are_asked_again_once_only
- tests/core/test_sysid_non_differentiable_residual.py::test_a_small_jump_inside_a_long_candidate_is_read_across_one_spacing
- tests/core/test_sysid_non_differentiable_residual.py::test_a_jump_read_the_second_time_is_the_one_reported
- tests/core/test_sysid_non_differentiable_residual.py::test_a_small_jump_met_by_one_on_the_other_side_is_not_one_sided
- tests/core/test_sysid_non_differentiable_residual.py::test_a_smooth_residual_read_a_second_time_is_still_not_a_jump
- tests/core/test_sysid_non_differentiable_residual.py::test_a_candidate_outside_the_two_thresholds_is_read_once
- tests/core/test_sysid_non_differentiable_residual.py::test_a_second_reading_through_a_residual_that_is_not_finite_says_nothing
- tests/property/test_sysid_targeted_search.py::test_fit_lm_reads_a_small_jump_of_the_ball_from_a_long_candidate
- tests/property/test_sysid_targeted_search.py::test_no_wrong_fit_among_the_fixed_draws
- tests/property/test_sysid_targeted_search.py::test_fit_lm_on_the_ball_is_not_converged_where_the_loss_falls_along_its_piece
status: verified
Expand Down
138 changes: 125 additions & 13 deletions src/maddening/sysid.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3923,9 +3923,45 @@ def _check_adam_hyper(n_iter, lr, tol, betas, eps, notify_every) -> tuple[int, i
#: coordinate's float spacing moves its value by 6e-4, at 21 -- and on the
#: stock bouncing ball, whose bounce moves by one time step, 2.2e4 to
#: 4.4e5 over 55 endings. ``2**10`` sits between, a factor of 48 above
#: the one and 21 below the other.
#: the one and 21 below the other. A wider draw (4,500 fits of the ball,
#: jaxlib 0.10.2 and 0.11.0) reads 1.9e3 to 5.8e5 at 261 jump endings and
#: 341 at one more: a small jump reads low over a whole candidate, so a
#: ratio above ``2**5`` is read a second time (:data:`_JUMP_SUSPECT`,
#: :data:`_JUMP_ACROSS`).
_JUMP_EXCESS = 2.0 ** 10

#: The ratio above which :func:`_one_sided_excess` reads a rejected
#: candidate a second time, across one spacing of the coordinates
#: (:func:`_excess_across_a_spacing`). ``2**5`` is above every ratio
#: measured on a smooth ending (:data:`_JUMP_EXCESS`: at most 21), so
#: those are read exactly as before.
_JUMP_SUSPECT = 2.0 ** 5

#: The second reading at which a candidate between the two thresholds
#: above is a jump. The ratio of a whole candidate understates a small
#: jump: it divides the jump by the model's change over the candidate's
#: length, and a late bounce that moves by one time step changes few
#: samples. The second reading divides by the model's change over one
#: spacing instead; the measure of the rounding is the same. Measured
#: (CPU, jaxlib 0.10.2 and 0.11.0, the same on both) on the stock ball:
#: two fits ended ``converged=True`` at ratios of 341 and 456 -- a jump of
#: 0.015 in the residual's norm, 22 float spacings of the elasticity from
#: the iterate, read from a candidate 35 spacings long -- where the second
#: reading is 1.8e3 and 2.0e3; it is 3.4e4 or more at the other 261 jump
#: endings of 4,500 drawn fits, and at most 5.2 at the six that ended on a
#: piece's own minimum. On smooth residuals stopped by the floor rule
#: (140 endings of the per-push and slow sysid tests, float32 and x64,
#: none of them asked: their first reading is under ``2**5``) it would be
#: at most 13, but for the wide ``logit`` fit of :data:`_JUMP_EXCESS` at
#: 52. ``2**8`` sits between, a factor of 4.9 above the one and 6.9 below
#: the other.
_JUMP_ACROSS = 2.0 ** 8

#: Halvings :func:`_excess_across_a_spacing` may take to bring a
#: candidate's step down to one spacing of the coordinates: more than the
#: bits of a float64 significand, so the cap is never what stops it.
_JUMP_BISECTIONS = 64


def _norm64(x) -> float:
"""``||x||`` in float64 without squaring a value outside its range."""
Expand Down Expand Up @@ -3961,9 +3997,58 @@ def _gain_in_the_gap(r, J, undamped, theta, rises) -> bool:
and promised > 0.0)


def _one_sided_excess(theta, r, J, rejected, mirror) -> tuple[float, float]:
"""``(excess, move)``: how one-sided the residual is around ``theta``,
read from candidates a Levenberg-Marquardt iteration rejected.
def _excess_across_a_spacing(th, r64, J64, cand, cand_r, behind, mirror, dtype) -> float:
"""The one-sided ratio of :func:`_one_sided_excess` for the rejected
``cand``, with the residual's departure and the model's change taken
across one spacing of the coordinates inside the step from ``th``
instead of over its whole length.

The step is halved, each time keeping the half over which the residual
departs further from the linear model, until its two ends are
neighbours on the coordinates' grid. The departure between those two
is then compared with ``behind`` -- the departure at the mirror image
of the whole candidate, which :func:`_one_sided_excess` measured --
plus the model's change over the one spacing.

A jump the candidate crossed is between the two ends, whole: the
numerator is what it was for the whole candidate, and the model's
change, which diluted it, is a spacing's worth. Of a differentiable
residual one spacing holds only a part of the candidate's departure,
and its mirror image's departure is of the same size: the reading
is of order one however the residual curves. ``behind`` is kept from
the whole candidate, not measured again across one spacing, because
one pair of neighbours is no measure of the rounding: their residuals
often agree exactly (measured: readings of 4e3 to 1e6 at smooth
floors with the mirror image taken across one spacing too).

``0.0`` where a residual is not finite.
"""
lo, lo_r = th, r64
hi, hi_r = np.asarray(cand, dtype=np.float64), np.asarray(cand_r, dtype=np.float64)
for _ in range(_JUMP_BISECTIONS):
mid, mid_r = mirror(jnp.asarray(0.5 * (lo + hi), dtype=dtype))
mid = np.asarray(mid, dtype=np.float64)
mid_r = np.asarray(mid_r, dtype=np.float64)
if np.array_equal(mid, lo) or np.array_equal(mid, hi):
break
near = _norm64(mid_r - lo_r - J64 @ (mid - lo))
far = _norm64(hi_r - mid_r - J64 @ (hi - mid))
if far > near:
lo, lo_r = mid, mid_r
else:
hi, hi_r = mid, mid_r
predicted = J64 @ (hi - lo)
ahead = _norm64(hi_r - lo_r - predicted)
scale = behind + _norm64(predicted)
if not (np.isfinite(ahead) and np.isfinite(scale)) or ahead == 0.0:
return 0.0
return float(ahead / scale) if scale > 0.0 else float(np.inf)


def _one_sided_excess(theta, r, J, rejected, mirror) -> tuple[float, float, bool]:
"""``(excess, move, jumped)``: how one-sided the residual is around
``theta``, read from candidates a Levenberg-Marquardt iteration
rejected, and whether that is a jump.

For a rejected candidate ``theta + d`` the residual's departure from
the linear model, ``||r(theta + d) - r - J d||``, is compared with the
Expand All @@ -3983,6 +4068,17 @@ def _one_sided_excess(theta, r, J, rejected, mirror) -> tuple[float, float]:
jump over those. No scale has to be assumed for the rounding: the
mirror image measures it on the problem itself.

A ratio above :data:`_JUMP_EXCESS` is a jump. The model's change
grows with the candidate and a jump does not, so a small jump read
from a long candidate gives a ratio between the two kinds: a candidate
above :data:`_JUMP_SUSPECT` and not above :data:`_JUMP_EXCESS` is read
a second time, across one spacing of the coordinates
(:func:`_excess_across_a_spacing`), and is a jump where that reading
is above :data:`_JUMP_ACROSS`. ``excess`` and ``move`` are those of
the candidate with the largest reading among the jumps, or among all
of them where none is one; the second reading is reported only for a
candidate it made a jump.

``rejected`` holds ``(candidate, residual, move)``; ``mirror(th)``
returns the coordinates ``th`` is evaluated at (projected onto the
bounds, on the leaves' grid) and the residual there. A candidate whose
Expand All @@ -3992,7 +4088,7 @@ def _one_sided_excess(theta, r, J, rejected, mirror) -> tuple[float, float]:
th = np.asarray(theta, dtype=np.float64)
r64 = np.asarray(r, dtype=np.float64)
J64 = np.asarray(J, dtype=np.float64)
excess, at_move = 0.0, 0.0
excess, at_move, jumped = 0.0, 0.0, False
for cand, cand_r, move in rejected:
d = np.asarray(cand, dtype=np.float64) - th
if not d.any():
Expand All @@ -4008,9 +4104,15 @@ def _one_sided_excess(theta, r, J, rejected, mirror) -> tuple[float, float]:
if not (np.isfinite(ahead) and np.isfinite(scale)) or ahead == 0.0:
continue
ratio = ahead / scale if scale > 0.0 else np.inf
if ratio > excess:
excess, at_move = float(ratio), float(move)
return excess, at_move
is_jump = ratio > _JUMP_EXCESS
if _JUMP_SUSPECT < ratio <= _JUMP_EXCESS:
across = _excess_across_a_spacing(
th, r64, J64, cand, cand_r, behind, mirror, theta.dtype)
if across > _JUMP_ACROSS:
ratio, is_jump = across, True
if (is_jump, ratio) > (jumped, excess):
excess, at_move, jumped = float(ratio), float(move), bool(is_jump)
return excess, at_move, jumped



Expand Down Expand Up @@ -5802,11 +5904,21 @@ def fit_lm(
linearisation's own change, the rejection is a jump and not rounding,
which is alike on both sides. The run then ends ``converged=False``
with a :class:`RuntimeWarning` naming a residual that is not
differentiable; ``params`` is the lowest loss it found. Measured, in
differentiable; ``params`` is the lowest loss it found. A small jump
read from a long candidate falls short of that (the linearisation's
change grows with the step, the jump does not), so a candidate at more
than ``2**5`` is read a second time: its step is halved down to the
one float spacing of the coordinates across which the residual departs
most (a residual evaluation a halving), and it is a jump where that
departure is ``2**8`` times the other side's and the linearisation's
change over one spacing. Measured, in
float32 and under x64: at most 21 on fits of smooth residuals stopped
by the floor rule (the spring over its transforms, units and noise
levels; ``HeartPumpNode`` at noise 0.02 and 1), and 2.2e4 to 4.4e5 on
the ball (float32). Not detected: a *kink* -- a continuous residual whose slope
the ball (float32); the second reading, at most 52 on the former and
1.8e3 or more on the ball's jumps, two of which read 341 and 456 the
first time.
Not detected: a *kink* -- a continuous residual whose slope
jumps -- in whose crease a run ends ``converged=True``, at a point no
step lowers the loss from but whose ``J`` is one-sided; a jump whose
mirror image the bounds do not allow; and a jump beside an iterate the
Expand Down Expand Up @@ -6312,7 +6424,7 @@ def _gauss_newton_stationary(th, before, r, J, loss_at_th=None) -> bool:
# shorter only because it is damped more, so it is not tested.
accepted = False
stationary = False
jump = (0.0, 0.0)
jump = (0.0, 0.0, False)
verdict["representable"] = True
# A loss of exactly 0.0 from a residual that is not: the framed
# sum's unframing underflowed float64 (a float64 residual below
Expand Down Expand Up @@ -6467,10 +6579,10 @@ def _gauss_newton_stationary(th, before, r, J, loss_at_th=None) -> bool:
# loss that can be far above the minimum. The rejected
# candidates tell the two apart (:func:`_one_sided_excess`).
jump = _one_sided_excess(theta, r, J, rejected, _mirrored)
stationary = not jump[0] > _JUMP_EXCESS
stationary = not jump[2]
if again:
continue
if jump[0] > _JUMP_EXCESS:
if jump[2]:
warnings.warn(_JUMP_WARNING.format(move=jump[1], excess=jump[0]),
RuntimeWarning, stacklevel=2)
break
Expand Down
Loading
Loading