TinyUDE: Solver-Free Universal Differential Equations on Microcontrollers via Lie-Taylor Jet Matching
Authors: Pranavanath Balamurali, Hrishi Kamireddy
Organizations: Department of Electrical and Computer Engineering The University of Texas at Austin · Department of Computer Science The University of Texas at Austin
Training Universal Differential Equations (UDEs) traditionally relies on backpropagating through numerical ODE solvers, creating memory footprints far exceeding the capabilities of edge microcontrollers. We present Lie-Taylor jet matching, a solver-free training framework that fits a hybrid vector field directly to the first and second time-derivatives of observed system states. These derivatives, the truncated Lie-Taylor jet, are estimated online via Savitzky-Golay filtering, yielding fully analytic gradients without automatic differentiation software. We evaluate whether eliminating the solver compromises accuracy against a conventional baseline (fixed-step RK4 integration, multiple shooting, exact discrete adjoints, Adam) sharing identical dynamics, noise models, network architectures, and metrics. While naive derivative matching degrades under sensor noise, our noise-adaptive mechanisms close and reverse this gap: full-rate phase-shifted sampling, a reservoir buffer, cosine-annealed optimization with weight averaging, on-device noise estimation, and polynomial-misfit quality gating. On a damped pendulum and chaotic double pendulum, our method matches or exceeds baseline accuracy at matched data windows and recovers unmodeled damping coefficients. Across noise levels from 0% to 5%, it attains a geometric-mean relative field error of 0.65x that of the baseline within 108 kB of static memory, compared with megabytes of solver tape. On an ESP32 microcontroller, the on-device run reaches a field error of 0.0020 and recovers the damping coefficient to c = 0.400 (true 0.400) within 61.3 kB of static memory and 7.24 ms per update (18.1% duty cycle at 25 Hz), confirming real-time on-device training is feasible without a numerical solver.
Figures & tables
Framework
Optimization Paradigm
Noise Handling
MCU Training
Training Memory
Adjoint UDE [ 1 ]
Solver in loss, reverse-mode
Multiple shooting
Non-viable
5.9 – 11.8 MB (tape)
PINN [ 17 ]
Autodiff on implicit repr.
Loss regularization
Offline only
High (grid-dependent)
SINDy [ 40 ]
Direct numerical diff.
Smoothing + sparsity
Feasible (linear)
Low
Jet Matching (Ours)
Lie–Taylor jet, analytic grads
Full-rate SG + adaptive gates
Microcontroller
108 kB (static)
Table 1: Taxonomy of differential equation learning frameworks. Memory figures for the adjoint baseline reflect measured usage, while jet matching figures represent static memory allocations.
Figure 1: Overview of standard decimated sampling versus full-rate phase sampling. Maintaining tap spacing m while sliding the filter window by stride s yields m/s independent derivative estimation streams.
Noise
Adjoint Baseline
Baseline ( K=80 )
Jet Matching (v2)
Ratio (Ref.)
Ratio ( K=80 )
Jet c^ (True: 0.40)
0%
0.00133±0.00013
0.00133±0.00013
0.00068±0.00022
0.51
0.51
0.399
0.1%
0.00123±0.00012
0.00123±0.00012
0.00223±0.00060
1.81
1.81
0.394
0.5%
0.00158±0.00019
0.00158±0.00019
0.00155±0.00022
0.98
0.98
0.398
1%
0.00283±0.00134
0.00283±0.00134
0.00218±0.00017
0.77
0.77
0.397
2%
0.00695±0.00328
0.00355†
0.00266±0.00031
0.38
0.75
0.395
5%
0.02477±0.00925
0.01180†
0.00642±0.00151
0.26
0.54
0.388
Table 2: Relative field error ( μ±σ over 3 seeds), single pendulum Case A, matched 60 s data window. Bold text indicates lowest error per noise level. † : 1–3 seeds. Error ratios <1.0 indicate lower relative error for jet matching.
Figure 2: The solver-based training loop that jet matching removes. Each gradient requires 4K vector-field evaluations forward and again in reverse, and the checkpoint tape retains every solver stage and activation for all K steps. Jet matching replaces the entire dashed optimization loop and the tape with the streaming pipeline of Algorithm 1 .
Metric
Solver-Based Baseline
Jet Matching (Ours)
Relative field error @ 1% noise
0.0028
0.0022
Relative field error @ 5% noise
0.0248
0.0064
Double pendulum field error @ 1%
0.0137
0.0115
Recovered damping c^ (true: 0.40)
0.398
0.397
Gradient error vs. complex-step
1e−13
1e−13
Numerical ODE solver in loop
Yes (RK4, 4K evaluations)
No
Table 3: Detailed comparison under matched physical plants, neural network structures, noise models, and evaluation routines.
Figure 3: Evaluation of the adjoint baseline: field error versus noise level at K=40 (A), accuracy across various shooting window lengths K (B), and trade-offs between computational cost and field accuracy (C).
Figure 4: System identification across UDE configurations under 1% noise: reconstructed residual functions versus true dynamics (A, B) and estimated damping parameters across random seeds (C).
Figure 5: Autonomous trajectory rollouts over an 8 s horizon (A), phase space trajectories across multiple initial conditions (B), and relative state error progression over time (C).
Configuration
Field Error
Error vs. Baseline ( 0.0028 )
v1: Decimated SG, minibatch size 32, basic SGD
0.0090
3.2× higher
+ Reservoir buffer ( B=2048 ) + Cosine Adam + EMA
0.0070
2.5× higher
+ Full-rate phase sampling
0.0037–0.0045
1.3 – 1.6× higher
+ Extended update steps ( 8,000 ) on matched dataset
0.0023
0.81× (lower)
+ Adaptive filtering, SNR gating, and quality gating (v2)
0.0022
0.77× (lower)
Table 4: Ablation study showing cumulative relative field error improvements on the single pendulum under 1% noise across 3 random seeds.
Figure 6: Evaluation of derivative loss components during training at 1% noise. The second-order loss term remains near a noise floor (A), and disabling it yields equivalent performance (B), supporting the SNR gating strategy.
Figure 7: Performance comparison of second-order loss contributions across noise levels and filter decimation steps: field error ratios ( λ2=0.1 vs. λ2=0 ) relative to calculated SNR thresholds.
Figure 8: Ablation analysis under 1% noise (A), noise-free parameter recovery performance (B), and the ESP32 resource account, projected alongside the physical measurement of Section 5.6 (C).
Figure 9: On-device validation on ESP32 hardware at 1% noise, B=2048 : relative field error and recovered c^ over training (A), measured per-update latency against the projected value and the 25 Hz timing budget, with a zoomed inset on the measured spread (B), and the measured resource account against the simulated projection (C).
Appendix figures & tables6 assets
Supplementary material from the paper’s appendix.
Appendix
Parameter
Single Pendulum
Double Pendulum
State dimension n / Network
n=2 , 2→16→16→2 ( p=354 )
n=4 , 4→16→16→4 ( p=420 )
Sensor rate / Training rate
fs=1 kHz, ftrain=25 Hz, 8,000 updates
SG window / Order / Stride
2M+1=11 , P=4 , s=5
Decimation step m
Automatic ( m≈39 at 1% noise)
Automatic ( m≈17 – 18 at 1% noise)
Reservoir capacity B / Batch
4096 / 64
2048 / 64
Optimizer / Schedule
Adam ( lr=0.02 ), Cosine decay to 2%, Norm clip 5.0
Appendix
Table 5: System parameters and configuration settings.
Test Case
Individual Seed Values
Mean
Jet v2, 0% noise
0.00048, 0.00091, 0.00064
0.00068
Jet v2, 0.1% noise
0.00163, 0.00284, 0.00222
0.00223
Jet v2, 0.5% noise
0.00163, 0.00172, 0.00131
0.00155
Jet v2, 1% noise
0.00205, 0.00212, 0.00236
0.00218
Jet v2, 2% noise
0.00249, 0.00302, 0.00247
0.00266
Jet v2, 5% noise
0.00810, 0.00603, 0.00515
0.00642
Appendix
Table 6: Individual per-seed relative field error values across test configurations.
Figure 10: Validation of analytic jet matching gradients against complex-step differentiation across 32 test conditions.
Figure 11: Validation of baseline discrete adjoint gradients against complex-step differentiation across 14 test conditions.
Figure 12: Savitzky–Golay derivative estimation error across effective step sizes Δteff and noise levels.
Figure 13: Effect of observation sampling rate on baseline adjoint field error.
Recovering continuous-time dynamics from discrete observations is difficult because local supervision (e.g., pointwise regression targets, derivative approximations, or equation residuals) loses fidelity as the observation interval grows. We replace local supervision with a global structural constraint: any flow representing autonomous dynamics must satisfy the semi-group property under time translation. We train a time-conditioned secant velocity field whose deviation from this property, which we call Symmetry Rupture, serves two purposes. As a training regularizer, it confines the hypothesis space to flows that compose consistently across temporal scales. As an inference oracle, it lets the solver select the largest step size that preserves internal consistency, replacing the local truncation error that conventional adaptive solvers depend on. On the diffusion-reaction benchmark under time-informed inference, our method reduces rollout RMSE by 87% while using 5x fewer function evaluations than a Neural ODE baseline. In the more demanding direct auto-regressive setting, where the model must predict distant future frames without intermediate temporal cues, our adaptive solver allocates compute based on local geometric complexity -- maintaining the lowest rollout RMSE on two of three PDE benchmarks while baselines either diverge or require up to an order of magnitude more function evaluations to remain stable.
Yuxiang Luo, Andrew Perrault
Department of Computer Science and Engineering The Ohio State University
We introduce a technique that enables Neural-ODEs to approximate arbitrary velocity fields with a priori planted fixed-points. Specifically, a recipe is given to explicitly accommodate for a finite collection of points in the reference multi-dimensional space of the Neural-ODE where the velocity field is exactly equal to zero. In this way, the gradient-based training is rigorously constrained inside the prescribed hypothesis class while leaving the expressive power of the Neural-ODE unaltered. We rigorously prove the universality of the Neural-ODE under any local constraints in the velocity field and give a computationally convenient way of imposing the fixed points. Our method is then tested on two paradigmatic physical models.
Feliciano Giuseppe Pacifico, Duccio Fanelli, Lorenzo Buffoni +3
Department of Informatics and Computer Science, University of Pisa, Italy · Department of Physics and Astronomy and INFN, University of Florence, Sesto Fiorentino, Italy
Neural surrogates for stiff differential-algebraic equations (DAEs) face two barriers: soft-constraint methods leave algebraic residuals that stiffness amplifies into errors, and hard-constraint methods require trajectory data from stiff integrators. We introduce an extended Newton implicit layer that enforces algebraic constraints exactly and reduces fast dynamics to their quasi-steady-state values in a single differentiable solve. Embedded in a physics-informed DeepONet, the layer recovers all fast and algebraic states exactly from slow-state predictions, removes the per-window stiffness-amplification pathway, and yields a stiffness-scaled Implicit Function Theorem gradient absent from penalty methods. Cascaded implicit layers extend this to multi-component systems with provable convergence. On a grid-forming inverter (stiffness ratio of about 4712), extended Newton attains 1.42% error versus 39.3% (penalty) and 57.0% (standard Newton); augmented Lagrangian and feedback linearization diverged. Two independently trained models compose without retraining (0.72% to 1.16% error, exact constraint satisfaction). Cross-domain validation on the Robertson stiff DAE (stiffness ratio up to 105) confirms generalization. Conformal prediction provides 90% coverage with automatic out-of-distribution detection.
Huy Hoang Le, Haoguang Wang, Christian Moya +2
School of Mechanical Engineering, Purdue University, West Lafayette, IN 47907 USA · Department of Mathematics, Purdue University, West Lafayette, IN 47907 USA · Department of Electrical and Computer Engineering, New Jersey Institute of Technology, Newark, NJ 07102 USA