Reject step sizes that are spuriously self-consistent in Consistency

[dweindl]

Summary

Consistency (the default success_checker for get_derivative) checked whether the requested methods (e.g. forward/backward/central) agreed with each other at each step size, then blended every self-consistent size's mean into the final value, with one final blanket tolerance check across those means.

The bug: this never compares a step size's estimate against the results from other step sizes. A step size can become small enough that all methods sample points within the target function's floating-point noise floor and become correlated (affected by the same rounding/cancellation error) — spuriously self-consistent, yet biased away from the true derivative. That biased size still got blended into the average while success reported True. Symmetrically, a large step size can also be self-consistently wrong due to higher-order/truncation effects.

This is the confirmed root cause of a long-standing intermittent CI failure in AMICI's PEtab benchmark gradient test (test_benchmark_gradient[Weber_BMC2015-*-unscaled]) — see AMICI-dev/AMICI#3078. That test uses fiddy.get_derivative with Consistency to finite-difference-check an analytically computed gradient for a model parameter (a32) several orders of magnitude smaller than the model's other free parameters. Two independent CI runs, 9 days apart, on identical code/commit, produced a bit-identical wrong finite-difference value (860.83 vs. the stable, cross-checked analytic value ~872.68); other CI runs on the same commit passed. The FD computation for that parameter was on a knife-edge where a negligible amount of run-to-run numerical noise was enough to flip a step size's self-consistency check and pull a biased estimate into the blend — a fragility of Consistency's algorithm, not of AMICI's simulation. A minimal, deterministic, AMICI-free repro of the same mechanism is now test_consistency_rejects_rounding_noise_dominated_step_sizes in this PR.

Fix

Add an order-independent iterative outlier-rejection pass (Consistency._reject_outliers) over the self-consistent step sizes' means: repeatedly drop the single worst-deviating candidate — relative to the median and a robust (MAD-based) spread of the rest — until nothing looks anomalous.

Deliberately does not:

• Treat step size magnitude as a proxy for trustworthiness (an earlier design that trusted the largest step first and cascaded down was dropped: a large step can cross a nonlinear feature and be self-consistently wrong too).
• Reuse the user's own rtol as the cross-size threshold: a relative tolerance loose enough to make sense for the per-size check is too loose to catch this kind of bias no matter where it's applied. A spread estimated from the trusted candidates themselves tightens automatically as they agree more closely, regardless of how loose the user's own rtol is.

This changes Consistency()'s default behavior, intentionally — today's default is the bug. Two new parameters, trend_n_sigma (default 5.0) and min_trend_samples (default 3), gate this; below min_trend_samples self-consistent sizes, behavior is unchanged from before.

Documented limitation (see the updated class docstring): like any purely data-driven consensus check, this has a ~50% breakdown point — if close to half (or more) of the self-consistent step sizes are corrupted, it cannot reliably tell which subset is trustworthy. That's a fundamental property of median/MAD-based statistics, not something engineerable away from data alone. Still a strict improvement over today's ~0% effective breakdown point, where a single corrupted-but-self-consistent size already breaks the check.

Test plan

• pytest tests/ — full suite passes (98 tests, including new tests/test_success.py)
• pre-commit run --all-files — ruff lint/format clean
• Regression test reproducing the flaky-check mechanism from a minimal synthetic function (no AMICI dependency)
• Unit tests for _reject_outliers covering: below-threshold no-op, no-outlier no-op, single/multiple outlier removal, order-independence, trend_n_sigma tuning, vector-valued outputs (drop whole candidate on any-element outlier), and a documented NaN-candidate edge case
• Once merged, AMICI can bump its fiddy pin and drop the exclude_params_unscaled workaround for Weber_BMC2015 in tests/benchmark_models/test_petab_benchmark.py

🤖 Generated with Claude Code

[dilpath]

Hm, we could alternatively name this RobustConsistency: essentially Consistency but with the additional check multi-step-size check. Or AutoRobustConsistency that checks for consistency at any of the step sizes and additionally automatically checks for robustness by adjusting the successful step size to be slightly smaller and checking that the same derivative is computed.

Regarding Weber, did you already check that appropriate simulation tolerances are used? I am not against a more robust Consistency, but previously when I have observed such problematic behavior it was rather due to a user-chosen tolerance being too relaxed.

[dweindl]

Hm, we could alternatively name this RobustConsistency: essentially Consistency but with the additional check multi-step-size check. Or AutoRobustConsistency that checks for consistency at any of the step sizes and additionally automatically checks for robustness by adjusting the successful step size to be slightly smaller and checking that the same derivative is computed.

For me it doesn't matter much whether this goes to an extra class or not, or whether it's opt-in in Consistency. I would anyways only use the more robust one in AMICI. Not sure fiddy is used anywhere else at the moment?

Regarding Weber, did you already check that appropriate simulation tolerances are used? I am not against a more robust Consistency, but previously when I have observed such problematic behavior it was rather due to a user-chosen tolerance being too relaxed.

I can't easily reproduce it to test much more. I am almost certain there would be some tolerance where this specific problem does not occur. However, my strong preference would be having a robust solution in fiddy that can deal with a given function instead of tuning fiddy hyperparameters for every single model and parameter vector. Without that, I don't see how for example AMICI-dev/AMICI#2093 could ever be implemented with fiddy.

[dilpath]

Thanks. The current Consistency check takes all FDs from all step sizes that are self-consistent (all FDs at a self-consistent step size are within atol/rtol of each other), then simply averages them. The new check introduces a new way of comparing FDs by multiple voting and median absolute deviation (MAD).

Instead of introducing MAD, could we just take all FDs from self-consistent step sizes, and then take the largest subset of those FDs that are within rtol/atol of each other, then average those? This would be more aligned with the pre-existing check, re-using the same rtol/atol from the user.

Also a warning should be emitted if a step size appears consistent but is then rejected IMO?

[dweindl]

Thanks. The current Consistency check takes all FDs from all step sizes that are self-consistent (all FDs at a self-consistent step size are within atol/rtol of each other), then simply averages them. The new check introduces a new way of comparing FDs by multiple voting and median absolute deviation (MAD).

Instead of introducing MAD, could we just take all FDs from self-consistent step sizes, and then take the largest subset of those FDs that are within rtol/atol of each other, then average those? This would be more aligned with the pre-existing check, re-using the same rtol/atol from the user.

For the issue I am addressing, this wouldn't help. Everything would be accepted. Or at least again require further fiddling around with the proper rtol/atol values. My feeling is that the MAD would be more robust.

I followed your earlier suggestion and added it via a separate class for now.

Also a warning should be emitted if a step size appears consistent but is then rejected IMO?

Added.

[dilpath]

Thanks. The current Consistency check takes all FDs from all step sizes that are self-consistent (all FDs at a self-consistent step size are within atol/rtol of each other), then simply averages them. The new check introduces a new way of comparing FDs by multiple voting and median absolute deviation (MAD).
Instead of introducing MAD, could we just take all FDs from self-consistent step sizes, and then take the largest subset of those FDs that are within rtol/atol of each other, then average those? This would be more aligned with the pre-existing check, re-using the same rtol/atol from the user.

For the issue I am addressing, this wouldn't help. Everything would be accepted. Or at least again require further fiddling around with the proper rtol/atol values. My feeling is that the MAD would be more robust.

I don't see why... rtol/atol are the precision to which the user wants to compute the derivative, right? Hence, collecting and averaging all derivatives that are within this rtol/atol of each other should only keep approximately "true"/"reproducible across step size" gradient values, up to the user-chosen rtol/atol? I don't see why trend_n_sigma solves the problem better... Do you have specific values from Weber for me to consider?

If you just want to get something merged without digging into it too much then fine for me since it's now a separate class, I can simply review the current state if you like.

[dweindl]

Thanks for pushing on this — let me try to answer with actual numbers rather than just intuition.

Concrete data, from the synthetic repro I added in test_consistency_rejects_rounding_noise_dominated_step_sizes (mimics the a32 case: a linear function plus a fixed-magnitude "noise" term that doesn't shrink with step size, at rtol=0.1/atol=1e-5, same as Weber's rtol_consistency/atol_consistency):

size    mean value
0.5     872.68
0.2     872.68
0.1     872.66
0.05    872.64
0.01    872.64
0.001   874.15   <- borderline-biased
0.0001  856.36   <- badly biased
1e-05   854.36   <- badly biased

With rtol=0.1 on values ~872.68, the tolerance band is ±87.3. Every one of these means is within that band of every other — so "largest subset mutually within rtol/atol" is all 8, identical to today's behavior. Reusing the same rtol that's loose enough for the per-size check (where 10% intra-method disagreement is often legitimate truncation error, not corruption) is too loose to separate "good" step sizes from "corrupted" ones — the corruption itself is only ~0.2-2% relative, well inside a 10% band.

I don't think this is really about picking a better number, though — I think it's structural. Consistency/the new checker gets instantiated once per get_derivative call and its rtol/atol apply across every free parameter in that call (direction_ids=parameter_ids). a32 is several orders of magnitude smaller than Weber's other free parameters — there's no way to give it a tighter cross-size tolerance without a per-parameter knob, which doesn't exist today (that's presumably why the workaround was to exclude it, not retune it). Any fixed threshold — reused from rtol_consistency, or a new tighter constant — has the same problem: it can't be right for parameters spanning many orders of magnitude in the same call. What I want here isn't a better fixed number, it's something that calibrates itself per parameter from that parameter's own step-size data — which is what the MAD-based scale does (it's computed independently per DirectionalDerivative, i.e. per parameter).

Worth being precise about scope too: I don't think fiddy's consistency check should try to assert final gradient accuracy — it structurally can't, it has no ground truth. That job is already rtol_check/atol_check in test_petab_benchmark.py, compared against AMICI's own adjoint/forward-sensitivity derivative. What flaked was the self-consistency gate (Consistency/rtol_consistency) letting a biased blend through with success=True, which then correctly failed the separate, already-existing rtol_check assertion downstream. RobustConsistency is only trying to fix that gate, not duplicate rtol_check's job.

Given it's a separate, opt-in class now, happy for you to just review the current state if that's easier — but wanted to give you the actual numbers first since you asked.

🤖 Generated with Claude Code

[dilpath]

Then rather implement some relative_step_size flag that treats step sizes as being relative to the parameter magnitude? 10% rtol seems too relaxed, but that's a different discussion.

[dweindl]

Then rather implement some relative_step_size flag that treats step sizes as being relative to the parameter magnitude? 10% rtol seems too relaxed, but that's a different discussion.

Well, again, my plan was to not tune tolerances for every single model.

If you don't think this is helpful to have in fiddy, I will move it to amici.

[dilpath]

Ah, by tune I thought you meant manually tune. What is the issue with automatic relative step sizes?

I see the benefit of checking for consistency across step sizes, but I don't yet see the logic in having two different concepts of what consistency means: one for same-size, and one for across-sizes. The only reason seems to be that rtol is too relaxed.