diff --git a/.github/workflows/draft-pdf.yml b/.github/workflows/draft-pdf.yml new file mode 100644 index 00000000..97bbba60 --- /dev/null +++ b/.github/workflows/draft-pdf.yml @@ -0,0 +1,19 @@ +on: [push] + +jobs: + paper: + runs-on: ubuntu-latest + name: Paper Draft + steps: + - name: Checkout + uses: actions/checkout@v6 + - name: Build draft PDF + uses: openjournals/openjournals-draft-action@master + with: + journal: joss + paper-path: paper/paper.md + - name: Upload + uses: actions/upload-artifact@v7 + with: + name: paper + path: paper/paper.pdf diff --git a/README.md b/README.md index b6b11fad..5af3865f 100644 --- a/README.md +++ b/README.md @@ -38,6 +38,10 @@ All are ultimately smoothing with similar runtime and accuracy, but some have fl All methods have hyperparameters, described in the [Sphinx documentation](https://pynumdiff.readthedocs.io/master/). We take a principled approach and propose a multi-objective optimization framework for choosing settings that minimize a loss function that balances faithfulness to data with smoothness of the derivative estimate. For more details, refer to [this paper](https://doi.org/10.1109/ACCESS.2020.3034077). +![Three simulated signals and their derivatives, estimated by six of the seven method families, with hyperparameters chosen by `pynumdiff.optimize`.](paper/methods_comparison.png) + +Above, three simulated signals are differentiated by two methods each, drawn from six of the seven families, with hyperparameters chosen automatically by `pynumdiff.optimize`. Reproduce with `python paper/make_figure.py`. + ## Installing Dependencies are listed in [pyproject.toml](https://github.com/florisvb/PyNumDiff/blob/master/pyproject.toml). They include the usual suspects like `numpy` and `scipy`, plus `pywavelets` for `waveletdiff`, `tqdm` for the optimizer, and `cvxpy` for `robustdiff` and `tvrdiff`. diff --git a/paper/make_figure.py b/paper/make_figure.py new file mode 100644 index 00000000..8f26d3d7 --- /dev/null +++ b/paper/make_figure.py @@ -0,0 +1,62 @@ +"""Generates paper/methods_comparison.png. Run from the repo root: python paper/make_figure.py""" +import numpy as np, matplotlib +matplotlib.use('Agg') +import matplotlib.pyplot as plt +from pynumdiff.utils import simulate +from pynumdiff import polydiff, tvrdiff, rtsdiff, spectraldiff, butterdiff, lineardiff +from pynumdiff.optimize import optimize + +# One color per method, fixed wherever it appears. Okabe-Ito, legible to colorblind readers. +STYLE = {spectraldiff:("spectraldiff","#CC79A7"), polydiff:("polydiff","#0072B2"), tvrdiff:("tvrdiff","#D55E00"), + butterdiff:("butterdiff","#5D3A9B"), rtsdiff:("rtsdiff","#009E73"), lineardiff:("lineardiff","#E69F00")} + +# Six of the seven families, paired so each panel contrasts two. Bandlimits are optimizer inputs, round numbers a +# user reads off a power spectrum; Lorenz is the faster signal. Comment out FITTED to re-run the search (minutes). +panels = [("Sum of sines", simulate.sine, [spectraldiff, polydiff], 1), + ("Triangle wave", simulate.triangle, [tvrdiff, butterdiff], 1), + ("Lorenz $x$", simulate.lorenz_x, [rtsdiff, lineardiff], 2)] +FITTED = {"spectraldiff": dict(cutoff_freq=0.0625, even_extension=True, pad_to_zero_dxdt=True), + "polydiff": dict(stride=7, degree=2, window_size=41, kernel='gaussian'), + "tvrdiff": dict(gamma=10.6015625, huberM=6.0, order=1), + "butterdiff": dict(cutoff_freq=0.05, num_iterations=1, filter_order=2), + "rtsdiff": dict(log_qr_ratio=3.81328125, order=1, forwardbackward=False), + "lineardiff": dict(gamma=0.00011125, window_size=53, order=2, kernel='friedrichs')} +#FITTED = {} + +def main(): + fig, axes = plt.subplots(2, 3, figsize=(13, 5.4), sharex=True, + gridspec_kw={'height_ratios':[1, 1.6], 'hspace':0.13, 'wspace':0.2}) + for j, (sname, sim, methods, bandlimit) in enumerate(panels): + x, x_truth, dxdt_truth = sim(duration=4, dt=0.01, noise_parameters=(0, 0.2), random_seed=3) + t = np.arange(len(x))*0.01 + top, bot = axes[0, j], axes[1, j] + + top.plot(t, x, '.', color='0.72', markersize=1.6, label="noisy data") + top.plot(t, x_truth, '-', color='black', linewidth=1.1, label="true $x$") + top.set_title(sname, fontsize=11, pad=6) + bot.plot(t, dxdt_truth, '-', color='0.55', linewidth=2.6, label="true $\\dot{x}$", zorder=1) + + for method in methods: + name, color = STYLE[method] + # Order 1 represents the triangle's piecewise-constant derivative exactly; without it tvrdiff goes spiky. + ssu = {'order':{1, 2, 3}} if (method is tvrdiff and sim is simulate.triangle) else {} + params = FITTED.get(name) or optimize(method, x, 0.01, bandlimit=bandlimit, search_space_updates=ssu)[0] + print(f" {sname:<14} {name:<13} {params}", flush=True) + bot.plot(t, method(x, 0.01, **params)[1], '-', color=color, linewidth=1.1, label=name, zorder=2) + bot.set_xlabel("time (s)", fontsize=10) + + for ax, pad in ((top, 0.17 if j == 0 else 0.04), (bot, 0.28)): # only the headroom each legend needs + lo, hi = ax.get_ylim(); ax.set_ylim(lo, hi + pad*(hi - lo)) + ax.tick_params(labelsize=8) + for s in ('top', 'right'): ax.spines[s].set_visible(False) + bot.legend(fontsize=8, frameon=False, loc='upper center', ncol=3, handlelength=1.6, columnspacing=1.1, + borderpad=0.1, title=f"optimized at bandlimit {bandlimit} Hz", title_fontsize=8) + + axes[0, 0].set_ylabel("$x$", fontsize=11) + axes[1, 0].set_ylabel("$dx/dt$", fontsize=11) + axes[0, 0].legend(fontsize=8, frameon=False, loc='upper center', ncol=2, handlelength=1.6, borderpad=0.1) + fig.savefig("paper/methods_comparison.png", dpi=200, bbox_inches='tight') + print("wrote paper/methods_comparison.png") + +if __name__ == "__main__": # spawned workers re-import this module, so nothing heavy may run at import time + main() diff --git a/paper/methods_comparison.png b/paper/methods_comparison.png new file mode 100644 index 00000000..ddc0362b Binary files /dev/null and b/paper/methods_comparison.png differ diff --git a/paper/paper.bib b/paper/paper.bib new file mode 100644 index 00000000..3e1a584d --- /dev/null +++ b/paper/paper.bib @@ -0,0 +1,181 @@ +@article{vanBreugel2022, + doi = {10.21105/joss.04078}, + url = {https://doi.org/10.21105/joss.04078}, + year = {2022}, + publisher = {The Open Journal}, + volume = {7}, + number = {71}, + pages = {4078}, + author = {{van Breugel}, Floris and Liu, Yuying and Brunton, Bingni W. and Kutz, J. Nathan}, + title = {{PyNumDiff}: A {Python} package for numerical differentiation of noisy time-series data}, + journal = {Journal of Open Source Software} +} + +@misc{komarov2025, + title = {A Taxonomy of Numerical Differentiation Methods}, + author = {Komarov, Pavel and {van Breugel}, Floris and Kutz, J. Nathan}, + year = {2025}, + eprint = {2512.09090}, + archivePrefix = {arXiv}, + primaryClass = {math.NA}, + url = {https://arxiv.org/abs/2512.09090} +} + +@article{vanBreugel2020numerical, + doi = {10.1109/ACCESS.2020.3034077}, + year = {2020}, + author = {{van Breugel}, Floris and Kutz, J. Nathan and Brunton, Bingni W.}, + journal = {IEEE Access}, + title = {Numerical differentiation of noisy data: A unifying multi-objective optimization framework}, + volume = {8}, + pages = {196865--196877} +} + +@article{chartrand2011numerical, + author = {Rick Chartrand}, + title = {Numerical differentiation of noisy, nonsmooth data}, + journal = {ISRN Applied Mathematics}, + year = {2011}, + volume = {2011}, + pages = {164564}, + doi = {10.5402/2011/164564} +} + +@article{brunton2016discovering, + author = {Steven L. Brunton and Joshua L. Proctor and J. Nathan Kutz}, + title = {Discovering governing equations from data by sparse identification of nonlinear dynamical systems}, + journal = {Proceedings of the National Academy of Sciences}, + year = {2016}, + volume = {113}, + number = {15}, + pages = {3932--3937}, + doi = {10.1073/pnas.1517384113} +} + +@article{virtanen2020scipy, + author = {Pauli Virtanen and Ralf Gommers and Travis E. Oliphant and others}, + title = {{SciPy} 1.0: Fundamental algorithms for scientific computing in {Python}}, + journal = {Nature Methods}, + year = {2020}, + volume = {17}, + pages = {261--272}, + doi = {10.1038/s41592-019-0686-2} +} + +@article{harris2020array, + author = {Charles R. Harris and K. Jarrod Millman and St{\'{e}}fan J. {van der Walt} and others}, + title = {Array programming with {NumPy}}, + journal = {Nature}, + year = {2020}, + volume = {585}, + pages = {357--362}, + doi = {10.1038/s41586-020-2649-2} +} + +@article{diamond2016cvxpy, + author = {Steven Diamond and Stephen Boyd}, + title = {{CVXPY}: A {Python}-embedded modeling language for convex optimization}, + journal = {Journal of Machine Learning Research}, + year = {2016}, + volume = {17}, + number = {83}, + pages = {1--5} +} + +@article{savitzky1964, + author = {Abraham Savitzky and Marcel J. E. Golay}, + title = {Smoothing and Differentiation of Data by Simplified Least Squares Procedures}, + journal = {Analytical Chemistry}, + year = {1964}, + volume = {36}, + number = {8}, + pages = {1627--1639}, + doi = {10.1021/ac60214a047} +} + +@article{huber1964, + author = {Peter J. Huber}, + title = {Robust Estimation of a Location Parameter}, + journal = {The Annals of Mathematical Statistics}, + year = {1964}, + volume = {35}, + number = {1}, + pages = {73--101}, + doi = {10.1214/aoms/1177703732} +} + +@article{rauch1965, + author = {Herbert E. Rauch and F. Tung and Charlotte T. Striebel}, + title = {Maximum likelihood estimates of linear dynamic systems}, + journal = {AIAA Journal}, + year = {1965}, + volume = {3}, + number = {8}, + pages = {1445--1450}, + doi = {10.2514/3.3166} +} + +@article{aravkin2013, + author = {Aleksandr Y. Aravkin and James V. Burke and Gianluigi Pillonetto}, + title = {Optimization viewpoint on {Kalman} smoothing with applications to robust and sparse estimation}, + journal = {Journal of Machine Learning Research}, + year = {2013}, + volume = {14}, + pages = {2513--2558}, + url = {https://jmlr.org/papers/volume14/aravkin13a/aravkin13a.pdf}, + doi = {10.1007/978-3-642-38398-4_8} +} + +@software{derivative_pkg, + author = {Andy Goldschmidt}, + title = {derivative: Numerical differentiation in {Python}}, + year = {2021}, + url = {https://github.com/andgoldschmidt/derivative} +} + +@software{findiff, + author = {Matthias Baer}, + title = {findiff: Finite difference derivatives in {Python}}, + year = {2018}, + url = {https://github.com/maroba/findiff} +} + +@software{pykalman, + author = {Daniel Duckworth and {The pykalman developers}}, + title = {pykalman: {Kalman} filters and smoothers for {Python}}, + year = {2012}, + url = {https://github.com/pykalman/pykalman} +} + +@article{kalman1960, + author = {Rudolf E. Kalman}, + title = {A new approach to linear filtering and prediction problems}, + journal = {Journal of Basic Engineering}, + year = {1960}, + volume = {82}, + number = {1}, + pages = {35--45}, + doi = {10.1115/1.3662552} +} + +@article{pysindy, + author = {Brian M. de Silva and Kathleen Champion and Markus Quade and Jean-Christophe Loiseau and J. Nathan Kutz and Steven L. Brunton}, + title = {{PySINDy}: A {Python} package for the sparse identification of nonlinear dynamics from data}, + journal = {Journal of Open Source Software}, + year = {2020}, + volume = {5}, + number = {49}, + pages = {2104}, + doi = {10.21105/joss.02104} +} + +@article{lee2019pywavelets, + author = {Gregory R. Lee and Ralf Gommers and Filip Wasilewski and Kai Wohlfahrt and Aaron O'Leary}, + title = {{PyWavelets}: A {Python} package for wavelet analysis}, + journal = {Journal of Open Source Software}, + year = {2019}, + volume = {4}, + number = {36}, + pages = {1237}, + doi = {10.21105/joss.01237} +} diff --git a/paper/paper.md b/paper/paper.md new file mode 100644 index 00000000..84d67fd4 --- /dev/null +++ b/paper/paper.md @@ -0,0 +1,119 @@ +--- +title: 'PyNumDiff: Practical Numerical Differentiation for Noisy Data' +tags: + - Python + - numerical differentiation + - time series + - denoising + - dynamics + - signal processing +authors: + - name: Pavel Komarov + orcid: 0009-0007-7482-2807 + corresponding: true + affiliation: 1 + - name: Floris van Breugel + orcid: 0000-0001-6538-7179 + affiliation: 2 + - name: Maria Protogerou + affiliation: 3 + - name: J. Nathan Kutz + orcid: 0000-0002-6004-2275 + affiliation: 4 +affiliations: + - name: Department of Electrical and Computer Engineering, University of Washington, USA + index: 1 + - name: Department of Mechanical Engineering, University of Nevada, Reno, USA + index: 2 + - name: Department of Applied Mathematics, University of Washington, USA + index: 3 + - name: Autodesk Research, London, UK + index: 4 +date: 1 April 2026 +bibliography: paper.bib +--- + +# Summary + +Derivatives of measured data are vital across science and engineering: identifying governing equations, designing controllers, and processing sensor streams alike. The textbook remedy, finite differencing, amplifies noise as $1/\Delta t$ and deteriorates rapidly as data grows noisier. Smoothing before differencing helps, but algorithm choice and settings substantially affect the result, and no single approach wins universally. + +PyNumDiff is an open-source Python package consolidating a broad suite of numerical differentiation methods under a unified API, along with an optimization scheme to select hyperparameters even in the absence of ground truth. It implements seven algorithm families: (1) prefiltering followed by finite difference calculation, (2) iterated finite differencing, (3) polynomial fitting [@savitzky1964], (4) global (spectral) and local (wavelet, radial) basis function fitting, (5) total variation regularization [@chartrand2011numerical], (6) Kalman smoothing [@kalman1960; @rauch1965], and (7) local approximation with linear models. All PyNumDiff methods return estimates of the true signal and its derivative as a matched pair, `(x_hat, dxdt_hat)`. A companion survey [@komarov2025, Section 7] situates these methods in the literature, covers the theory behind each method, benchmarks all methods across test signals, and guides selection for different application scenarios. This paper describes the package's second generation, through the v0.3 release series. + + +# Statement of Need + +Numerical differentiation has a diverse ecosystem of algorithms, each with different theoretical foundations and putative strengths. But there has been no consolidated home where methods can be easily swapped or tested head-to-head, and practitioners are left assembling solutions piecemeal from disparate packages [@vanBreugel2022]. + +PyNumDiff addresses this gap: Its unified interface lets users compare methods on the same data, exploit specialized capabilities (e.g. outlier robustness), and tune methods' hyperparameters against known ground truth or to a given signal bandlimit. Derivative estimation trades data fidelity against smoothness, so PyNumDiff uses dual-objective optimization, with a single scalar weight to steer the balance [@vanBreugel2020numerical]. A natural and growing application is SINDy [@brunton2016discovering], which discovers governing equations by regressing measured derivatives, making quality of the derivative estimates a direct determinant of model accuracy. + + +![Three simulated signals with additive noise (top) and their derivatives (bottom). Each panel contrasts two methods from different families, six of the seven in total, with a fixed color per method. Hyperparameters come from `pynumdiff.optimize` rather than hand-tuning, at the bandlimit noted with each derivative panel: 1 Hz for the sines and triangle, 2 Hz for the faster Lorenz signal. That bandlimit is an input to the optimizer, the sort of round figure a user reads off a power spectrum, rather than a claim about the signals, whose frequency content runs past it. The triangle's piecewise-constant derivative separates the methods most clearly, `tvrdiff` reproducing the steps where `butterdiff` rounds and rings at each corner.\label{fig:comparison}](methods_comparison.png) + + +# State of the Field + +Relevant Python tools exist, but each covers only part of the space. `numpy.gradient` and `scipy.signal.savgol_filter` [@virtanen2020scipy] cover only a couple specific cases; `findiff` [@findiff] offers high-order finite difference stencils suited to clean simulation data, but not to noisy measurements; `pykalman` [@pykalman] does not provide state-space models to turn the filter into a differentiator; and total variation regularization [@chartrand2011numerical] was implemented with lagged-diffusivity iteration that approaches the convex minimizer rather than attaining it. Historically, practitioners have had to patch together packages, standalone scripts, and academic mathematics, with no shared API or principled way to compare results. The `derivative` package [@derivative_pkg] attempts to bring more algorithms under one roof but lacks `NaN` handling and hyperparameter optimization. No existing package spans PyNumDiff's seven method families with a consistent interface. + +The original PyNumDiff publication [@vanBreugel2022] established the core method set and optimization framework. This second generation rewrites and consolidates it from the ground up. `kerneldiff`, `rtsdiff`, and `tvrdiff` each cover by parameter what had been separate functions, so the package presents a simplified interface. `rbfdiff`, `waveletdiff`, and `robustdiff` are new, bringing the set to thirteen, all organized into seven modules and callable with consistent keyword-argument signatures. A single optimizer serves every one of them. This generation also adds support for multidimensional data, as well as irregular sample spacing, missing observations, and circular domains where the underlying mathematics allows. + + +# Software Design + +**Package design.** All differentiation methods share the call signature + +```python +x_hat, dxdt_hat = method(x, dt_or_t, **hyperparams) +``` + +where `x` is a NumPy array [@harris2020array] of measurements, `dt_or_t` is either a scalar step size or an array of sample locations, and keyword arguments configure the method, making calls self-documenting. The v0.2 series preserves prior positional signatures behind deprecation warnings; v0.3 removes them, leaving only the keyword-argument interface. + +**Architecture.** PyNumDiff is organized into seven differentiation modules plus shared `utils` and `optimize` in a flat structure. Where strong alternatives exist, PyNumDiff delegates rather than reimplements: SciPy [@virtanen2020scipy] provides spline fitting, Savitzky-Golay filtering, and signal processing routines; NumPy [@harris2020array] provides the FFT; PyWavelets [@lee2019pywavelets] provides the discrete wavelet transform for `waveletdiff`; CVXPY [@diamond2016cvxpy] handles convex optimization for `robustdiff`, `tvrdiff`, and `lineardiff`. The public `kalman_filter` and `rts_smooth` primitives let advanced users with known dynamics bypass `rtsdiff`'s constant-derivative model to build their own bespoke differentiator. `utils` houses common subroutines like estimation of integration constants and data scales, calculation of evaluation metrics, and generation of simulated noisy data. `optimize` includes not only cost minimization but a whole framework: per-method default search spaces, bounds with rounding rules, categorical handling, a caching layer, and a configurable objective that can work even in the absence of ground-truth. + +**Method capabilities.** Every method handles multidimensional data via `axis`. Table 1 groups methods by their additional talents. + +\begin{table}[!ht] +\centering +\begin{tabular}{@{}lp{0.70\linewidth}@{}} +\hline +\textbf{Capability} & \textbf{Methods} \tabularnewline +\hline +Variable step size & \raggedright \texttt{polydiff}, \texttt{splinediff}, \texttt{rbfdiff}, \texttt{rtsdiff}, \texttt{robustdiff}, \texttt{lineardiff} \tabularnewline +Missing data & \raggedright \texttt{polydiff}, \texttt{splinediff}, \texttt{rtsdiff}, \texttt{robustdiff}, \texttt{lineardiff} \tabularnewline +Outlier robustness & \raggedright \texttt{robustdiff}, \texttt{tvrdiff} \tabularnewline +Circular domain & \raggedright \texttt{rtsdiff} \tabularnewline +\hline +\end{tabular} +\caption{Methods by auxiliary capability.} +\end{table} + +**Irregular and incomplete sampling.** Whether and how a method supports variable step size depends on the underlying model. For example, basis spline fits are indifferent to spacing, able to place knots with equal ease anywhere in the domain, but Kalman-based methods must compute a discrete-time transition by matrix exponential at each actual interval, while Butterworth filters inflexibly assume a constant sampling rate. Likewise, missing observations, given as NaN entries, can be excluded from fitting and imputed from the model in select cases, useful for real data where sensors may drop samples. + +**Outlier robustness.** `robustdiff` replaces the quadratic Kalman cost with Huber loss [@huber1964] terms on both measurement and process residuals, following @aravkin2013, with CVXPY [@diamond2016cvxpy] as the optimization backend, operating over a sparse problem formulation to scale linearly with signal length. `tvrdiff` similarly applies Huber loss on data fidelity, with its total variation penalty on the derivative additionally promoting piecewise-smooth solutions. + +**Circular and wrapped domains.** `rtsdiff` accepts `circular=True` for quantities like angles on a periodic domain. Innovation residuals are wrapped to $[-\pi, \pi]$ before each Kalman update via an `innovation_fn` hook, and `x_hat` is returned in the same range, avoiding the large spurious spikes naive smoothers produce when a signal crosses the $\pm\pi$ boundary. + +**Hyperparameter optimization.** `pynumdiff.optimize` selects hyperparameters $\Phi$ by minimizing +$$L(\Phi) = \text{RMSE}\big(\textstyle\int\hat{\dot{x}}(\Phi) + c,\; x\big) + \gamma\,\text{TV}\big(\hat{\dot{x}}(\Phi)\big),$$ +which requires no ground truth, because fidelity is measured by integrating the estimated derivative back against the measured signal [@vanBreugel2020numerical]. Root mean squared error in the first term can be (and is by default since v0.2.1) replaced with a robust variant, so it is less prone to bias in the presence of outliers. The smoothness weight $\gamma$ is derived internally from the data's sampling rate and user-provided `bandlimit`, the highest frequency of meaningful signal. Default search spaces are defined at the top of the `optimize` module, with opportunities to collapse to fewer search dimensions thoroughly explored, for speed. Minimization requires a continuous space, so discrete numbers like polynomial order and window size are handled by rounding, while truly categorical hyperparameters like kernel type are supported by pooling results across runs employing separate choices. Repeated evaluations are cached to avoid redundant work. A single call to `suggest_method` runs this search across nearly every method in the package and reports the best fit, the practical payoff of putting all behind one interface. + +**Testing and continuous integration.** The test suite validates all methods against analytic functions with known derivatives, covering noiseless and noisy cases, scale equivariance, multidimensional application, and missing-data handling. Utilities and optimization are also tested for correct behavior, to verify functionality like robust estimation of data scale and the equivalence of parallel and serial search. The tests avoid tautology, never letting an implementation define its own expected result: Accuracy is asserted against empirically established error bounds, with the suite reporting any method that beats its bound by an order of magnitude so the bound can be ratcheted down. An optional `--plot` flag renders every case for visual inspection. GitHub Actions runs the suite on every push and pull request, and Coveralls tracks line coverage, currently 94%. + + +# Research Impact Statement + +The first generation [@vanBreugel2022] of PyNumDiff is cited in over twenty publications and has been applied in experimental biology (flight kinematics from motion capture) and control engineering (observer design). The PySINDy project [@pysindy] maintains its own differentiation submodule overlapping substantially with PyNumDiff's capabilities, and integration discussions are ongoing. + +The companion taxonomy paper [@komarov2025, Section 7] provides accuracy and bias results for every included method and derives recommendations, based on performance against simulated signals, using loss, search spaces, and other choices now encoded in PyNumDiff v0.3. The evidence is reproducible, with experiments runnable from Jupyter notebooks in PyNumDiff's repo, spanning sampling rate, noise distribution, noise scale, outlier contamination, and signal bandlimit. Additional tutorial notebooks cover basic usage, hyperparameter optimization, automatic method suggestion, and demonstrations of the outlier, multidimensional, and circular-domain capabilities described above. Full API documentation is published at [pynumdiff.readthedocs.io](https://pynumdiff.readthedocs.io/master/). + + +# AI Usage Disclosure + +This paper was drafted with assistance from Claude Opus 5 (Anthropic), which also helped draft code revisions before and during review. All outputs were reviewed and heavily remolded by hand, and the authors take full responsibility for accuracy. + + +# Acknowledgements + +The authors thank Yuying Liu and Bingni W. Brunton for their contributions to the original PyNumDiff package [@vanBreugel2022], and Sasha Aravkin for discussions on convex optimization techniques that informed robust differentiation methods. This work was supported by the NSF AI Institute in Dynamic Systems (grant number 2112085). + + +# References