99 - signal processing
1010authors :
1111 - name : Pavel Komarov
12- orcid : 0000-0000-0000-0000 # TODO: fill in
1312 corresponding : true
1413 affiliation : 1
1514 - name : Floris van Breugel
16- orcid : 0000-0000-0000-0000 # TODO: fill in
1715 affiliation : 2
1816 - name : J. Nathan Kutz
19- orcid : 0000-0003-0944-900X
20- affiliation : 1
17+ affiliation : 3
2118affiliations :
22- - name : Department of Applied Mathematics , University of Washington, USA
19+ - name : Department of Electrical and Computer Engineering , University of Washington, USA
2320 index : 1
2421 - name : Department of Mechanical Engineering, University of Nevada, Reno, USA
2522 index : 2
23+ - name : Autodesk Research, London, UK
24+ index : 3
2625date : 1 April 2026
2726bibliography : paper.bib
2827---
2928
3029# Summary
3130
3231Computing derivatives from measured data is a foundational requirement across science and
33- engineering. Whether fitting governing equations from experimental observations, designing
32+ engineering. Whether identifying governing equations from experimental observations, designing
3433control laws, or processing sensor streams, researchers routinely need derivative estimates
3534from discrete, noisy measurements. The textbook approach — finite differencing — amplifies
3635noise proportionally to $1/\Delta t$, making it unreliable for real data. Smoothing before
3736differencing helps, but the choice of algorithm and its tuning parameters substantially
3837affect the result, and no single choice is universally best.
3938
40- PyNumDiff is an open-source Python package that addresses this by providing a broad suite of
41- numerical differentiation methods for noisy data under a unified API. Seven families of
42- algorithms are implemented : (1) prefiltering followed by finite difference calculation;
43- (2) iterated finite differencing; (3) polynomial fitting; (4) basis function fitting;
44- (5) total variation regularization; (6) Kalman smoothing and its outlier-robust variant;
45- and (7) local approximation with linear models. All methods return a smoothed signal estimate
46- and a derivative estimate as a matched pair ` (x_hat, dxdt_hat) ` . A companion paper
47- [ @ komarov2025 ] provides a comprehensive theoretical taxonomy of these methods, benchmarks
48- their relative performance, and guides method selection for different application scenarios.
39+ PyNumDiff is an open-source Python package that consolidates a broad suite of numerical
40+ differentiation methods for noisy data under a unified API. Seven families of algorithms are
41+ organized : (1) prefiltering followed by finite difference calculation; (2) iterated finite
42+ differencing; (3) polynomial fitting; (4) basis function fitting; (5) total variation
43+ regularization; (6) Kalman smoothing; and (7) local approximation with linear models. All
44+ methods return a smoothed signal estimate and a derivative estimate as a matched pair
45+ ` (x_hat, dxdt_hat) ` . A companion paper [ @ komarov2025 ] provides a comprehensive theoretical
46+ taxonomy of these methods, benchmarks their relative performance, and guides method selection
47+ for different application scenarios.
4948
5049
5150# Statement of Need
5251
5352Estimating derivatives from noisy measurements arises throughout experimental science,
54- data-driven modeling, system identification, and control. Standard finite differences amplify
55- noise, and while smoothing helps, the field has produced a diverse ecosystem of specialized
56- algorithms — each with different strengths regarding outlier robustness, computational cost,
57- handling of irregular time steps, or treatment of missing data. Without a consolidated
58- library, practitioners either default to the nearest available tool or implement bespoke
59- solutions that are hard to compare or reproduce.
60-
61- PyNumDiff fills this gap. Its unified interface lets users compare methods directly on the
62- same data, exploit specialized capabilities, and select hyperparameters without requiring
53+ data-driven modeling, system identification, and control. The field has produced a diverse
54+ ecosystem of specialized algorithms, each with different strengths regarding outlier
55+ robustness, computational cost, handling of irregular time steps, or treatment of missing
56+ data. Without a consolidated library, practitioners either default to the nearest available
57+ tool or implement bespoke solutions that are difficult to compare or reproduce.
58+
59+ PyNumDiff addresses this gap. Its unified interface lets users compare methods on the same
60+ data, exploit specialized capabilities, and select hyperparameters without requiring
6361ground-truth derivatives. The package is particularly valuable in workflows that use
64- derivatives as regression targets: for example, Sparse Identification of Nonlinear Dynamics
65- (SINDy) [ @brunton2016discovering ] learns governing equations by regressing measured
66- derivatives, making clean, reliable derivative estimates a prerequisite for accurate model
67- identification.
62+ derivatives as regression targets: Sparse Identification of Nonlinear Dynamics (SINDy)
63+ [ @brunton2016discovering ] , for example, learns governing equations by regressing measured
64+ derivatives, making reliable derivative estimates a prerequisite for accurate model
65+ identification. More broadly, any system identification or control design pipeline that
66+ begins by estimating rates of change from sensor data is a natural user of PyNumDiff.
6867
6968
7069# State of the Field
7170
72- Relevant Python tools exist, but none covers the full scope of PyNumDiff. ` numpy.gradient `
71+ Relevant Python tools exist but none covers the full scope of PyNumDiff. ` numpy.gradient `
7372and ` scipy.signal.savgol_filter ` [ @virtanen2020scipy ] are widely available but address only
7473a narrow slice of the method space. The ` findiff ` package provides high-order finite
7574difference stencils for clean simulation data rather than noisy measurements. The standalone
7675` TVRegDiff ` code [ @chartrand2011numerical ] implements a single total variation regularization
77- method; PyNumDiff includes and extends this. No existing Python package combines seven method
78- families with a unified API, NaN and variable-step support, multidimensional capability,
79- circular domain handling, and integrated hyperparameter optimization.
76+ method; PyNumDiff includes and extends this. The ` derivative ` package [ @derivative_pkg]
77+ implements several of the same methods, but without multidimensional support, missing data
78+ handling, or a hyperparameter optimization framework. No existing Python package combines
79+ the breadth of PyNumDiff's seven method families with a consistent API and the full set of
80+ capabilities described below.
8081
81- The original PyNumDiff [ @vanBreugel2022 ] introduced four method families with basic parameter
82- optimization. The present version substantially extends that foundation in scope, capability,
83- and software quality.
82+ The original PyNumDiff publication [ @vanBreugel2022 ] introduced the core concept and a
83+ multi-objective optimization framework for hyperparameter selection. The present version
84+ represents a substantial revision of the codebase: methods were reorganized into a cleaner
85+ taxonomy, corrected, and extended; the interface was unified; test coverage was dramatically
86+ improved; and several new capabilities were added throughout.
8487
8588
8689# Software Design
@@ -91,71 +94,91 @@ and software quality.
9194x_hat, dxdt_hat = method(x, dt_or_t, ** params)
9295```
9396
94- where ` x ` is a NumPy array [ @harris2020array ] of measurements, ` dt_or_t ` is either a
95- scalar step size or an array of sample locations (for variable step size), and keyword
96- arguments configure the method. The two return values always match the shape of ` x ` . This
97- uniformity makes method comparison straightforward and enables drop-in substitution.
98-
99- ** Multidimensional support.** An ` axis ` parameter controls which axis is differentiated,
100- allowing all methods to operate on blocks of time series simultaneously. The implementation
101- iterates over all non-time axes with ` np.ndindex ` , applying the algorithm independently to
102- each vector, so a 2D array of shape ` (N, D) ` is handled without reshaping or looping in
103- user code.
97+ where ` x ` is a NumPy array [ @harris2020array ] of measurements, ` dt_or_t ` is either a scalar
98+ step size or an array of sample locations (for variable step size), and keyword arguments
99+ configure the method. The two return values always match the shape of ` x ` . The revision
100+ moved away from positional parameter lists toward explicit keyword arguments throughout,
101+ making calls self-documenting and eliminating a common source of user error. Deprecated
102+ positional interfaces are retained with warnings.
103+
104+ ** Multidimensional support and interface consistency.** An ` axis ` parameter controls which
105+ dimension is differentiated, allowing all non-deprecated methods to operate on blocks of
106+ time series simultaneously without reshaping or looping in user code. The implementation
107+ iterates over all non-time axes using ` np.ndindex ` , applying the algorithm to each vector
108+ independently.
109+
110+ ** Missing data and variable step size.** Several methods handle NaN-valued entries as
111+ missing observations and non-uniform sample spacing: ` splinediff ` , ` polydiff ` , ` rbfdiff ` ,
112+ ` rtsdiff ` , and ` robustdiff ` . This supports realistic experimental scenarios where sensors
113+ drop samples or operate at irregular rates. For the Kalman-based methods, variable step size
114+ is handled correctly by computing the discrete-time transition matrix via matrix exponential
115+ at each actual time increment, rather than using a fixed $\Delta t$ extracted from the first
116+ two samples — a subtle but important distinction for data with irregular spacing.
104117
105- ** Missing data and variable step size.** Several methods handle NaN-valued entries as missing
106- observations and non-uniform sample spacing: ` splinediff ` , ` polydiff ` , ` rbfdiff ` , ` rtsdiff ` ,
107- and ` robustdiff ` . This supports realistic experimental scenarios where sensors skip samples or
108- operate at irregular rates.
118+ ** Outlier robustness.** ` robustdiff ` solves a convex optimization problem that replaces the
119+ least-squares Kalman cost with Huber loss terms on both measurement residuals and process
120+ residuals, following the robust smoothing framework of @aravkin2013 and using ` cvxpy `
121+ [ @diamond2016cvxpy ] as the optimization backend. The problem is formulated as a sparse
122+ system so that it runs in time linear in the number of samples, making it practical for
123+ long time series. ` tvrdiff ` similarly promotes piecewise-smooth solutions by penalizing
124+ the total variation of the derivative, and the data fidelity term can be replaced with a
125+ Huber loss for additional robustness to outliers. Both methods are available via
126+ ` pip install pynumdiff[advanced] ` .
109127
110128** Circular and wrapped domains.** ` rtsdiff ` accepts a ` circular=True ` flag for quantities
111- like angles that live on a periodic domain. Innovation residuals are wrapped to $[ -\pi, \pi] $
112- before each Kalman update, and ` x_hat ` is returned wrapped to the same range. This avoids the
113- erroneous large-magnitude spikes that naive smoothers produce when a signal crosses the
114- $\pm\pi$ boundary. The wrapping is injected via an optional ` innovation_fn ` hook on the
115- underlying ` kalman_filter ` primitive, keeping the filter itself general.
116-
117- ** Outlier robustness.** ` robustdiff ` solves a convex optimization problem that replaces the
118- least-squares Kalman cost with Huber loss on both measurement and process residuals,
119- dramatically reducing outlier influence. It uses ` cvxpy ` [ @diamond2016cvxpy ] as its
120- optimization backend, available via ` pip install pynumdiff[advanced] ` . ` tvrdiff ` similarly
121- minimizes a total variation penalty on the derivative, promoting piecewise-smooth solutions
122- and tolerating abrupt changes; it supports regularization of the velocity, acceleration, or
123- jerk (orders 1–3).
124-
125- ** Hyperparameter optimization.** Because every method has tuning parameters, PyNumDiff
126- provides a multi-objective optimization framework in ` pynumdiff.optimize ` that minimizes a
127- weighted sum of a smoothness penalty and a data fidelity term [ @vanBreugel2020numerical ] .
128- When ground-truth derivatives are available, the optimizer minimizes derivative error
129- directly. The smoothness weight ` tvgamma ` can be initialized from the signal's estimated
130- cutoff frequency $f_c$ via the empirical formula
129+ like angles that live on a periodic domain. Innovation residuals are wrapped to
130+ $[ -\pi, \pi] $ before each Kalman update step via an ` innovation_fn ` hook on the underlying
131+ ` kalman_filter ` primitive, and ` x_hat ` is returned wrapped to the same range. This avoids
132+ the erroneous large-magnitude spikes that naive smoothers produce when a signal crosses the
133+ $\pm\pi$ boundary.
134+
135+ ** Hyperparameter optimization.** Every method has tuning parameters, so PyNumDiff provides
136+ a multi-objective optimization framework in ` pynumdiff.optimize ` that minimizes a weighted
137+ combination of a smoothness penalty and a data fidelity term [ @vanBreugel2020numerical ] .
138+ When ground-truth derivatives are available, derivative error is minimized directly.
139+ The smoothness weight ` tvgamma ` can be initialized from the signal's estimated cutoff
140+ frequency $f_c$ via the empirical formula
131141$$ \texttt{tvgamma} = \exp(-1.6\ln f_c - 0.71\ln \Delta t - 5.1). $$
142+ The optimization was substantially improved in the current version: intermediate evaluations
143+ are cached to avoid redundant function calls, the loss function was made robust to outliers
144+ via Huber penalty, and the Kalman parameter space was reduced from two independent noise
145+ variances to their log-ratio, which is the quantity the result actually depends on
146+ [ @komarov2025 ] .
147+
148+ ** Testing and continuous integration.** The test suite validates all methods against a set
149+ of analytic test functions with known derivatives, checking both noiseless and noisy cases
150+ across the full dynamic range of expected accuracy. Tests are run automatically on every
151+ push and pull request via GitHub Actions, with coverage tracked via Coveralls. Care was
152+ taken to avoid tautological tests in which the implementation directly determines the
153+ expected result.
132154
133155
134156# Research Impact
135157
136- The original PyNumDiff paper [ @vanBreugel2022 ] has accumulated over 60 citations since its
137- 2022 publication, and the package has been used in experimental biology to estimate flight
138- kinematics from motion capture, in control engineering for observer design, and in
139- data-driven dynamics identification. The companion Taxonomy paper [ @komarov2025 ] , submitted
140- to the Journal of Computational Physics, provides the theoretical analysis underpinning the
141- method collection and benchmarks all included methods across diverse test signals. The
142- package is available on PyPI (` pip install pynumdiff ` ), documented at
143- [ pynumdiff.readthedocs.io] ( https://pynumdiff.readthedocs.io/master/ ) , and archived on
144- Zenodo [ @pynumdiff_zenodo] .
158+ PyNumDiff has been applied in experimental biology to estimate flight kinematics from motion
159+ capture, in control engineering for observer design, and in data-driven dynamics
160+ identification [ @vanBreugel2022 ] . The package is available on PyPI
161+ (` pip install pynumdiff ` ) and documented at
162+ [ pynumdiff.readthedocs.io] ( https://pynumdiff.readthedocs.io/master/ ) . The companion Taxonomy
163+ paper [ @komarov2025 ] , submitted to the Journal of Computational Physics, provides the
164+ theoretical analysis motivating the method collection and benchmarks all included methods.
145165
146166
147167# AI Usage Disclosure
148168
149- The draft of this paper was prepared with assistance from Claude Sonnet 4.6 (Anthropic).
150- This included integrating material from the repository, documentation, and related papers
151- into a coherent draft . The authors reviewed and edited all content and take full
152- responsibility for its accuracy.
169+ The draft of this paper was prepared with assistance from Claude Sonnet 4.6 (Anthropic),
170+ integrating material from the repository, documentation, release history, and related
171+ publications . The authors reviewed and edited all content and take full responsibility for
172+ its accuracy.
153173
154174
155175# Acknowledgements
156176
157177The authors thank Yuying Liu and Bingni W. Brunton for their contributions to the original
158- PyNumDiff package [ @vanBreugel2022 ] .
178+ PyNumDiff package [ @vanBreugel2022 ] , and Sasha Aravkin for discussions on convex
179+ optimization techniques that informed the robust differentiation methods. This work was
180+ supported by the National Science Foundation AI Institute in Dynamic Systems (grant
181+ number 2112085).
159182
160183
161184# References
0 commit comments