Fit with least_squares instead of L-BFGS-B - #45
Merged
Merged
Conversation
Every Bouligand fit now goes through one private method, `_fit`, which is
`scipy.optimize.least_squares` (trust-region reflective) on the existing
vector-valued `residuals` rather than `minimize` on the scalar `min_func`.
The problem is a sum of squares and `residuals` already exposed the vector,
so TRF can use structure L-BFGS-B cannot see. The difference is not
marginal. Measured in CPU time -- wall clock is useless on a shared box --
over synthetic windows, L-BFGS-B spent 3.1 s on the four-parameter fit and
16.5 s on a `dz` profile against 13 ms and 49 ms here. Its cost did not
track the number of function evaluations at all (34x spread at equal
counts), so the time was going into its own machinery, not the forward
model. The objective was 1-3% of a profile.
Answers are unchanged, which is the reason to do it:
- fitted parameters agree to four figures on 8/8 synthetic seeds, and
where they differ least_squares finds the lower misfit
- the four `dz` intervals pinned in the new test are the values L-BFGS-B
produced
- over cached EMAG2 spectra in a downstream project, eight vertices at
2000 km reproduce stock exactly, and 0 of 108 profile intervals across
three window sizes and three targets excluded their point estimate
`_fit` also supplies the two exact Jacobian columns -- `zt` and `C` enter
the forward model linearly, so d r/d zt = -2k/sigma and d r/d C = 1/sigma --
and differences only `beta` and `dz`, with a small cache so the Jacobian
reuses the residual `least_squares` just evaluated. Both are needed: the
columns alone, recomputing the base, gave 1.32x where together they give
1.51x. `beta` enters through the *order* of a Bessel function and has no
closed-form derivative; `dz` has one, but it costs an extra `kv` call,
exactly what the difference it would replace costs.
`_covariance` and `_jacobian` are deliberately left on central differences.
They are 0.16 ms, and the calibrated sigmas depend on them.
One behavioural difference: `least_squares` refuses `lb == ub`, which
L-BFGS-B accepted by pinning the parameter. `_fit` widens a collapsed bound
to the narrowest numerically distinct box, which pins it just as effectively
and stays close enough to the edge for `_warn_on_bounds` to fire.
Two tests drove the fit by monkeypatching `min_func`, which nothing calls
now. The sensitivity one patches `_fit` and counts fits instead, so it still
blows up inside the simulation loop where the prior is copied; the
Metropolis one is unaffected, since `log_posterior` still uses `min_func`,
but its docstring explained a discrepancy that no longer exists.
Suite: 100 -> 107 tests, 293 s -> 12 s.
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01BLGk2RXegPShhDRHcH7FtL
Adversarial review found three things. **The performance numbers were an artefact.** The 3.1 s / 16.5 s figures came from an unpinned BLAS: L-BFGS-B calls a threaded one whose workers spin-wait, and `time.process_time` charges every spinning thread, so on a contended machine the same comparison reads 20x to 400x depending on load. Pinned to one thread it is **3x** -- 12.5 ms against 5.0 ms for a fit, 102 ms against 33 ms for a `dz` profile. The durable number, which no threading can distort, is evaluations of the forward model per vertex: 2012 against 359. CLAUDE.md said to prefer `process_time` over wall clock; that is the advice that produced the wrong number, and it now says to pin the BLAS instead. The change is still worth making. It is 3x and 5.6x fewer model evaluations, not two orders of magnitude. **The Curie chain rule had no test.** Deleting `transform[fixed[0], free.index(1)] = -1.0` or flipping its sign passed the entire suite, because trust-region reflective converges from a wrong Jacobian too -- just by a different route, for a few extra evaluations. The old test asserted `curie.cost == min_func(...)`, which tests `expand()` and holds for any Jacobian at all. Replaced with a test that captures the callable `_fit` hands scipy and differences the residual it hands it alongside, in all three configurations and with priors active. Both mutations now fail it. **An inverted bound was silently accepted.** The widening predicate was `upper <= lower`, so `bounds[2] = (30.0, 10.0)` -- a typo, and `bounds` is documented as reassignable -- was rewritten to `[30, 30+3e-11]` and pinned `dz` at 30 with no error, where L-BFGS-B raised. Now `==` only, and an inverted bound reaches `least_squares` and raises. Also replaced `test_fit_cost_is_min_func`, which was a tautology -- both sides are half the sum of squares of the same vector, and it survived a mutation returning `2 * res.cost` from `_profiled_misfit`, the exact failure its docstring claimed to guard. It now checks that profiling a parameter at its own fitted value recovers the unconstrained misfit, which pins the scale. Suite 107 -> 111. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01BLGk2RXegPShhDRHcH7FtL
The four pinned intervals failed on macOS at rtol=2e-3, by up to 1.3%, with identical code -- the only CI leg that failed. An interval endpoint is where `brentq` crosses the threshold on a deviance curve that is nearly flat there, which is the whole reason dz needs a profile rather than a sigma. A last-ulp difference in `kv`, or in the FFT behind the synthetic, moves the crossing far more than it moves the fit. Pinning to 0.2% was asserting on the platform's libm, not on this package. 5%, which is ~4x the observed cross-platform spread. Still load-bearing: a sign-flipped zt column, zeroed finite-difference columns, and a `_profiled_misfit` off by a factor of two each fail all four cases at the new tolerance. Those move an endpoint by tens of percent. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01BLGk2RXegPShhDRHcH7FtL
Four parallel cleanup reviews. Applied: **The affine map had two implementations.** `transform` said the `dz` row differentiates through `zt`; `expand` separately hand-coded `x[fixed[0]] = value - x[1]`. Same relation, two representations, kept in step by review -- and an inconsistency between them is close to invisible, because the fit converges from a wrong Jacobian anyway, just by a longer route. `expand` is now `transform @ y + offset`, so the Jacobian is the derivative of the expansion by construction. (The matmul stays: measured at 0.98 us against 1.55 us for column indexing, and indexing cannot express the chain-rule term.) **The Jacobian rebuilt its constant part on every call.** The analytic columns depend only on `kh` and `sigma`, and the prior rows only on the prior -- neither on `x`. Hoisting them out drops the lazy `base = None` guard and a four-line comment about not double-counting prior rows, because the differenced columns no longer touch those rows at all. **`needed` probed the matrix it had just built** -- four `ndarray.any()` calls on 3-element rows, measured at 4.4-6.2 us per fit, over half the per-call setup. It is `free` minus the closed-form columns, plus the held one under the Curie constraint, which is also what it means. **The residual cache kept four entries for a one-deep pattern.** Instrumented over six seeds: `residuals` calls are exactly `nfev + 2*njev`, so the base lookup hits every time and nothing else ever does. The comment justified the extra entries by rejected steps being retried, which trust-region reflective does not do -- a shrunken region gives a different point. One entry, same hit rate. (The cache itself earns its keep: ~0.1 ms of hashing against 9-27 ms of forward model per profile.) **`_bound_arrays`** -- `_fit` and `metropolis_hastings` had each written the `None`-means-infinite unpacking, and had already drifted on dtype. **`_ANALYTIC_COLUMNS` keyed by name** rather than by position, so `_PARAMETERS` order is encoded once. Same for `free.index(_ZT)`. Tests: the Jacobian spy no longer runs a fit it discards; the FD reference step uses `_JACOBIAN_STEP` instead of a literal, which would otherwise let a retuned constant pass a test that no longer measures the scheme; parametrised ids; and `assert calls["n"] > 2` is gone, since `exploding` raises on exactly that condition and the assertion could not fail. Mutation-checked after the refactor: chain-rule dropped and sign-flipped, `zt` column sign-flipped, prior rows zeroed, `_profiled_misfit` doubled, and the new `differenced` list missing its Curie entry all still fail. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01BLGk2RXegPShhDRHcH7FtL
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Every Bouligand fit now goes through one private method,
_fit, which isscipy.optimize.least_squares(trust-region reflective) on the existing vector-valuedresidualsrather thanminimizeon the scalarmin_func.Why
The problem is a sum of squares and
residualsalready exposed the vector, so TRF can use structure L-BFGS-B cannot see. Over six synthetic 200 km windows, BLAS pinned to one thread:_fitoptimiseprofile("dz")About 3x, on 5.6x fewer evaluations of the forward model.
Answers
dzintervals pinned in the tests are the values L-BFGS-B producedWhere the likelihood is flat the two can disagree by a lot —
dzof 310 km against 139 km for a misfit difference of 2e-6 on one 100 km synthetic. Both answers are equally good there; that is a property of the window, not of this change.Analytic Jacobian
ztandCenter the forward model linearly, sodr/dzt = -2k/sigmaanddr/dC = 1/sigmaare exact; onlybetaanddzare differenced. A small cache lets the Jacobian reuse the residualleast_squaresjust evaluated — the columns alone gave 1.32x where together they give 1.51x.betaenters through the order of a Bessel function and has no closed-form derivative.dzhas one, but it costs an extrakvcall, exactly what the difference it would replace costs._covarianceand_jacobianstay on central differences: 0.16 ms, and the calibrated sigmas depend on them.Behavioural notes
least_squaresrefuseslb == ub, which L-BFGS-B accepted by pinning the parameter._fitwidens a collapsed bound to the narrowest numerically distinct box. An inverted bound (lb > ub) is left to raise — it is a typo, andboundsis documented as reassignable.min_func, which nothing calls now. The sensitivity one patches_fitand counts fits instead, so it still blows up inside the simulation loop where the prior is copied.optimise's body withminimizewill diverge.Global_CPD'sfit_spectrumdid, and its exactfit_spectrum == optimiseassertion failed until it was routed through_fittoo.Testing
The Jacobian is tested directly, by capturing the callable
_fithands scipy and differencing the residual it hands it alongside — in the unconstrained,dz-held anddz = CPD - ztcases, with priors active. This is not optional rigour: TRF converges from a wrong Jacobian too, so deleting the Curie chain-rule term or flipping its sign passed the entire suite before this test existed.Suite: 100 -> 111 tests, 293 s -> 21 s.
🤖 Generated with Claude Code
https://claude.ai/code/session_01BLGk2RXegPShhDRHcH7FtL