`tvgamma` can be automatically calculated

[pavelkomarov]

Why make users give tvgamma when it's relatively hard to interpret? Shouldn't they give us instead the cutoff frequency, or "I want to keep __% of power in the signal and consider the rest noise", or even a guessed signal:noise ratio? And from there we can calculate a cutoff frequency, check the sampling rate, and weight the terms of the optimization appropriately, automatically.

[pavelkomarov]

Some measurements that bear on this, from comparing tvgamma-tuned against RMSE-tuned hyperparameters across the notebook 4 sweep (4368 cells per method).

The payoff is real and large. Rank correlation between a simulation's bandwidth and how much accuracy tvgamma tuning gives up is −0.90. Logistic Growth is the extreme: 90% of its true derivative's power sits below 0.5 Hz, but the sweep hands it cutoff_frequency=3, so tvgamma comes out 17.6× too small and every method inherits under-regularized parameters. Its median shortfall against RMSE-tuned parameters is 6.5×, where the well-specified signals sit at 1.1–2.0. Deriving the cutoff per signal would recover most of the largest tuning loss in the study.

But "keep __% of power" does not survive contact with noise. Measuring f90, the frequency below which 90% of power lies, on the noisy signal rather than on the true derivative:

	true f90 of dxdt	f90 of noisy x, ns=1	f90 of noisy x, ns=4
Triangles	4.00	1.25	36.25
Logistic Growth	0.50	1.50	34.50
Cruise Control	4.00	2.50	33.25

At high noise most of the spectrum's power is noise, so the percentile runs off toward Nyquist. Worse, even at low noise it gets the ordering backwards, calling Triangles narrower than Logistic Growth when it is 8× wider.

I also tried the obvious alternative, estimating a noise floor from the top decile of frequencies and taking the last bin above 2× that floor. It returns ~48 Hz for every signal regardless of content, because noise crosses any fixed multiple of its own median by chance all the way to Nyquist.

What might actually work. The library already estimates noise σ robustly, via MAD of second differences over √6, in splinediff (polynomial_fit.py:54) and robustdiff (kalman_smooth.py:353). That gives a principled noise floor with no user guess. Combined with a smoothed periodogram (Welch) and a run-length criterion, rather than a single-bin crossing, a defensible knee looks reachable. Two things to watch: the quantity wanted is the bandwidth of the derivative, which differentiation's ω-weighting makes systematically wider than the signal's, and tvgamma goes as roughly f^-1.6, so a 2× error in the cutoff is a 3× error in the weight.

Related refactor. There are six robust-scale estimates in the package, and they are three distinct quantities rather than six variants:

• noise level: polynomial_fit.py:54 and kalman_smooth.py:353 are identical but for nan_policy; basis_fit.py:218 is the same idea in the wavelet basis, uncentered per Donoho-Johnstone
• data scale: linear_model.py:162 and total_variation_regularization.py:88 are identical
• residual scale: utils/evaluate.py:98 and utils/utility.py:59 are identical

Only the first group earns a shared helper. The /√6 is the part worth writing down once, since it comes from the second-difference stencil (1, −2, 1) inflating variance by 1 + 4 + 1. The other two groups are already a single median_abs_deviation(·, scale='normal') call and would gain nothing from a wrapper.

One note on the good first issue label: separating signal from noise in the spectrum is close to the problem the differentiation is trying to solve, so this may be more involved than it looks.

---

🤖 Written by Claude Code on behalf of @pavelkomarov

[pavelkomarov]

I might prefer just keep plain MAD calls in place, but the trouble is the way noise and data scale are calculated across methods is all the same in spirit but differs in details like whether centering happens, whether to use second differences first, whether to account for possible NaNs. It's annoying and ends up taking up more than just one line in most places we do it, so it may be worth trying to unify so the wrinkles can be bottled up and tested in isolation, much like how integration constant estimation turns out to be shockingly complicated. Readers would then also not have to get derailed by repeated noise scale estimation motifs, when they're really not the point of how any method works.

[pavelkomarov]

When we solve this, the top of notebook 2 should be updated to reflect. I like still showing the frequency plot, but where the red line goes can be automatically selected. That would make the barrier to entry a little lower for a lot of users.

[pavelkomarov]

Re-running notebook 4's sweep with dxdt_truth supplied, so hyperparameters are chosen against true RMSE, gives a median blind-to-truth error ratio of 1.508 across 61,152 cells. That is the ceiling on what any better blind objective could recover.

Most of it is not the objective's shape. The sweep deliberately hands every simulation the same bandlimit=3 to keep the experiment a star in parameter space, but the six signals don't share a bandwidth, so where the assumed cutoff is wrong every method inherits slightly suboptimal hyperparameters. The 50% is therefore mostly the prize for per-signal estimation — the benchmark understates what a practitioner gets estimating the cutoff from their own data.

That handicap does not distort the comparison, though: Every method within each of a plot's frames sees the same signal and the same assumed cutoff, so the rank statistics hold even where the absolute errors are inflated. What differs is how much a wrong cutoff costs each method — going from bandlimit 1 to 5, WaveletDiff's error grows 1.04x while RobustDiff's and TVRDiff's grow 2.8x — which is exactly what sweeping that axis is for.

---

🤖 Written by Claude Code on behalf of @pavelkomarov