Improve lineardiff: strong accuracy, but costly to optimize and several rough edges

[pavelkomarov]

Splitting this off from #216, which is about docs and paper wording. This is about lineardiff the method: it looks like it earns its place on accuracy, but it is expensive enough to optimize that it is awkward to include in benchmarks, and it has several rough edges.

It measures up on accuracy

Six simulations, three noise seeds each, every method given its own optimize() run against dxdt_truth so nobody gets hand-picked settings. True RMSE, averaged over seeds:

sim        lineardiff   polydiff   rtsdiff   best
cruise          0.654      0.860     0.854   lineardiff
sine            0.785      1.054     0.815   lineardiff
tri             0.987      0.996     0.943   rtsdiff
pop             0.390      0.444     0.415   lineardiff
linaut          0.713      0.763     0.739   lineardiff
lorenz          0.840      1.030     0.875   lineardiff

Best in 5 of 6, and zero failures across the run. So the method is worth investing in rather than demoting.

But it is expensive to optimize

Mean optimize() wall time on the same run:

lineardiff  1255.5s
polydiff      84.0s
rtsdiff       26.9s

15x slower than polydiff, 47x slower than rtsdiff. Two compounding reasons:

1. Every window is a separate CVXPY solve, so a single lineardiff call is ~0.85-1.2s against 0.002s for savgoldiff and 0.030s for rtsdiff.
2. The search space is the largest of any method. With order searched as a categorical over {1,2,3} and three numerical dimensions (gamma, step_size, window_size), Nelder-Mead gets restarted from 3 x 48 = 144 seed points.

step_size is the dominant cost knob, since it sets how many convex solves happen per call. Worth deciding whether it should be searched at all, or pinned to a fraction of window_size.

Rough edges

No axis parameter. Every other method takes one. This currently blocks the README capability table, which states that every listed method accepts axis.

Not scale equivariant. f(a*x) is not a*f(x); measured relative deviation is 1.00, meaning the output barely scales with the input at all. Cause is that gamma is an absolute regularization weight, so rescaling the data silently changes the balance between fidelity and prior. Tracked generally in #222. An earlier attempt to normalize amplitude inside _lineardiff made things worse (func4 went 5.06 -> 24.4) because it normalized per window, so windows were fit against different scales and the crossfade blended incommensurate pieces. The fix is a single global scale computed once in lineardiff, not per window.

Ill-conditioning at small windows. Each row of the design matrix is one more integration than the row below it, so in raw time the rows differ by powers of the window duration. cond(integral_X) at window_size=11 runs 3.2e3 to 7.2e3 across sims. Working in nondimensional time tau = t/T holds every row at O(1) and brings cond near 30 regardless of window size. That change is written but not yet committed.

Solver failures concentrated at high order. Sweeping 980 parameter combinations, CVXPY fails on 79 without preconditioning, distributed by order as {2:2, 3:15, 4:25, 5:37}, and 90 with it, distributed {3:4, 4:31, 5:55}. Preconditioning removes the order-2 failures entirely but does not help orders 4 and 5, which suggests those orders are intrinsically fragile here rather than merely badly scaled. Searching order over {1,2,3} instead of pinning it at 3 keeps the search away from the worst of it.

Ridge damping is unresolved. A small Tikhonov term markedly helps small windows and hurts large ones:

sim         ws        lam=0  lam=0.0001    lam=0.01       lam=1
cruise      11        3.944       2.051       1.752       1.742
cruise      21        1.747       1.008       0.785       0.779
cruise      41        0.741       0.969       1.408       1.488
lorenz      41        0.751       0.759       1.575       1.982

Since the optimizer picks window_size in the 26-67 range anyway, damping may be solving a problem the search already avoids. Needs a decision rather than a default.

Open question

Can lineardiff's optimization be made speedy enough for the method to go in notebook 4/. On accuracy it deserves to be there; at 1255s per optimize() call it materially lengthens the benchmark.

---

🤖 Written by Claude Code on behalf of @pavelkomarov