TinyUDE: Solver-Free Universal Differential Equations on Microcontrollers via Lie-Taylor Jet Matching

1.1 Motivation & Problem Statement

• The training pipeline requires solver backpropagation tapes, automatic-differentiation libraries, and offline trajectory access, making it ill-suited to edge deployment.A 354-parameter baseline with K = 40 requires 5.9 MB of memory tape and peaks at 9.1 MB of heap usage.Source: §1.1 · pp. 1–2, §1.2 · p. 2

1.2 Edge Constraints: Memory Limitations and Solver Fragility

• Embedding numerical integration in training couples parameter updates to solver stability, while direct collocation reaches relative field error 0.82 under 5% noise.Solver-free training therefore needs filtering to produce reliable derivative targets.Source: §1.2 · p. 2

1.3 Summary of Contributions

• The method fits hybrid vector fields to first- and second-order derivatives estimated by Savitzky–Golay filtering, using closed-form gradients instead of automatic differentiation.The gradients agreed with complex-step differentiation to machine precision across a 32-check suite.Source: §1.3 · p. 2
• Full-rate phase sampling, reservoir buffering, adaptive filter tuning, and quality gating are used to sustain derivative matching under measurement noise.Full-rate phase sampling extracts approximately 50× more statistically independent derivative targets using a 4.4 kB ring buffer.Source: §1.3 · pp. 2–3
• Jet matching achieves a 0.65× geometric-mean field error ratio to the baseline across 0%–5% noise, while using 108 kB of static memory.On a chaotic double pendulum, field error is 0.0116 versus 0.0137, with damping estimates of 0.388–0.399 against the true 0.40.Source: §1.3 · p. 3

2.1 Universal Differential Equations and Neural ODEs

• Neural ODE and UDE training obtains gradients through numerical integration, while jet matching retains the UDE decomposition but removes the solve-inside-the-loss construction.The paper's baseline uses fixed-step Runge–Kutta integration, multiple shooting, and an exact discrete adjoint.Source: §2.1 · p. 3

2.2 On-Device Machine Learning and Memory Constraints

• Existing microcontroller training methods target classification or update only restricted model components, while SPSA gradient variance grows with parameter count.Jet matching is presented as complementary for dense sensor streams and networks with hundreds of parameters.Source: §2.2 · pp. 3–4

2.3 Derivative Estimation, System Identification, and Collocation

• Jet matching adds an exact-gradient second-order loss that uses acceleration estimates to constrain the vector-field Jacobian, beyond fitting field values alone.It also recovers independent samples discarded by decimated Savitzky–Golay filtering and rejects windows where the polynomial model fails.Source: §2.3 · p. 4
• Table 1 compares differential-equation learning frameworks, with adjoint memory figures measured and jet-matching figures reported as static allocations.The supplied caption does not specify table row or column encodings.Source: Table 1 · p. 4

2.4 Gradient Computation, Optimization, and Verification

• The pipeline composes vector–Jacobian and Jacobian–vector products for analytic gradients and verifies them against complex-step differentiation.Its optimization combines Adam, cosine learning-rate scheduling, iterate averaging, and fixed-memory reservoir sampling.Source: §2.4 · p. 4

3.1 Problem Formulation and Dynamical Systems

• Three UDE configurations omit damping, mis-specify gravity, or combine both errors; field error and regression of network outputs against velocities assess model accuracy and damping recovery.The systems are a damped single pendulum and a chaotic damped double pendulum.Source: §3.1 · p. 5

3.2 Lie–Taylor Jet Matching Loss and Analytic Gradients

• Jet matching fits filtered first- and second-time-derivative estimates to the vector field and its Lie–Taylor second-order prediction, with analytic parameter gradients that avoid ODE solvers.The second derivative prediction is Gθ(x) = JFθ(x)Fθ(x); fixed layer buffers support gradient evaluation without forming the network parameter Jacobian or allocating dynamic memory.Source: §3.2 · p. 6

3.3 Streaming Derivative Estimation: Full-Rate Phase Sampling

• Sliding the filter window by stride s while maintaining tap spacing m yields approximately m/s independent derivative-estimation streams, which are collected in a reservoir for minibatch training.The single-pendulum buffer uses 4.4 kB at m = 50; reservoir capacities are 4096 samples for the single pendulum and 2048 for the double pendulum.Source: §3.3 · p. 7, Figure 1 · p. 7

3.4 Noise-Adaptive Filter Tuning and Quality Gating

• Noise adaptation estimates sensor noise, selects a filter time scale, gates second-order matching by SNR, and rejects windows with excessive polynomial-fit residuals.Second-order matching is disabled below the SNR threshold, while rejected polynomial-misfit samples are excluded from the training buffer.Source: §3.4 · pp. 7–8

3.5 Static Memory Allocation

• Static allocations for parameters, optimizer state, buffers, and filtering total 108.1 kB for the single-pendulum setup, remaining constant during online training.The on-device measurement uses a reduced reservoir capacity of 2048 to fit available board DRAM alongside the runtime.Source: §3.5 · p. 8

4.1 Adjoint Baseline Configuration

• The reference baseline uses RK4, multiple shooting, exact discrete-adjoint backpropagation, and full-batch Adam; its K = 40 configuration requires a 5.9 MB tape and 9.1 MB peak heap.The discrete adjoint was validated against complex-step gradients, with differences of 1.8×10−15 to 5.2×10−13 across tested window lengths.Source: §4.1 · pp. 8–9

4.2 Experimental Control and Data Matching

• Both methods share dynamics, noise generation, initial conditions, network architectures, weight initialization, and field-error code, with matched durations of 60 s and 20 s for the single and double pendulums.Experiments use three random seeds for the single pendulum and four for the double pendulum.Source: §4.2 · p. 9

4.3 Execution Environment and Resource Projections

• Accuracy simulations use 64-bit Python/NumPy, while ESP32 resource projections assume 184,960 MACs per update at three cycles per MAC and 60 cycles per tanhf call.The ESP32 target is specified as a 240 MHz dual-core processor with 320 kB SRAM; the paper reports that Section 5.6 replaces projected timing with a direct hardware measurement.Source: §4.3 · p. 9

5.1 Accuracy Across Sensor Noise Levels

• Across 0%–5% noise, jet matching has a 0.65× geometric-mean error ratio versus the reference baseline and lower mean field error in five of six conditions.At 1% noise, relative field error is 0.0022 for the adaptive method versus 0.0028 for the adjoint baseline under matched 60 s windows.Source: §5.1 · p. 9, Table 2 · p. 9

5.2 Evaluation on Chaotic Double Pendulum Dynamics

• At 1% noise on the chaotic double pendulum, jet matching’s mean relative field error is 0.0115 versus 0.0137 for the adjoint baseline.One of four seeds produced an error of 0.025; rapid chaotic flips can exceed the approximation capacity of fixed-window fourth-order filters, yielding derivative errors that pass the quality gate.Source: §5.2 · p. 11

5.3 System Identification and Long-Term Rollout Dynamics

• Over an 8 s rollout, the trained hybrid field has trajectory error 0.0068, 393× below physics-only models and 10× below unadapted jet matching.When true damping is included in the known field, the estimated residual coefficient is −0.001.Source: §5.3 · p. 11

5.4 Ablation Study and Performance Drivers

• Full-rate phase sampling and larger buffers provide the largest ablation accuracy gains, while time-scale adaptation and quality gating add refinements, particularly at low noise.At 1% noise, the second-order loss remains near its noise floor (∼2×10−2), and disabling second-order matching avoids high-variance estimates affecting updates.Source: §5.4 · p. 12, Table 4 · p. 12, Figure 6 · p. 12

5.5 Resource Usage and Memory Requirements

• Jet matching uses 108.1 kB static memory versus 5.9 MB of baseline tape at K = 40, a 55–109× reduction.The cycle-counted ESP32 projection is 2.85 ms per update, or 7.1% duty cycle at 25 Hz.Source: §5.5 · pp. 12–13, Figure 8 · p. 13

5.6 On-Device Validation on ESP32 Hardware

• On an ESP32, the run reaches relative field error 0.0020 and recovers c = 0.400 (true 0.400) using 61.3 kB static memory.Measured updates take 7.24 ms, corresponding to an 18.1% duty cycle at 25 Hz; the timing is consistent across the 8,000 updates.Source: §5.6 · p. 13, Figure 9 · p. 14

6.1 Summary

• 55–109× lower training memory accompanies comparable or superior vector-field accuracy across multiple test systems.The approach operates within 108 kB of statically allocated memory.Source: §6.1 · p. 14
• Full-rate phase sampling, reservoir buffering, adaptive filter tuning, and quality gating support direct derivative matching under measurement noise.Source: §6.1 · p. 14

6.2 Limitations and Future Work

• Physical ESP32 validation covered one operating point: 1% noise, a single pendulum, and B = 2048.Testing the full noise sweep, double pendulum, and power consumption on-device remains future work.Source: §6.2 · p. 14
• At 0.1% noise, traditional adjoint methods retain an accuracy advantage, with 1.8× lower absolute error.Source: §6.2 · p. 14
• Rapid transitions in chaotic systems can exceed polynomial-filter assumptions, while non-Gaussian noise, latency, and quantization remain untested.The latter conditions are outside the reported focus on simulated physical systems under Gaussian noise.Source: §6.2 · p. 14

A.1 Reverse-over-JVP Derivation for Second-Order Jet Gradients

• The predicted second derivative is Gθ(ˆx) = JFθ(ˆx) Fθ(ˆx), and its parameter gradient is derived from the second-order loss.Source: §A.1 · pp. 17–18
• Reverse-mode differentiation over the Jacobian–vector product JNθv computes parameter gradients without explicitly forming higher-order Jacobian matrices.The forward-mode tangent pass supplies JNθv.Source: §A.1 · p. 18

A.2 Noise Estimator Variance Identity

• The noise-variance estimator is obtained by inverting the variance expression for the scalar residual of a polynomial-preserving Savitzky–Golay fit.The derivation models the signal as an exact polynomial plus Gaussian noise.Source: §A.2 · p. 18

• The bob’s Cartesian position is expressed using θ, and differentiating it gives velocity components used to form kinetic energy.The kinetic-energy expression simplifies using cos^2 θ + sin^2 θ = 1.Source: §A.3 · p. 18
• With y = 0 as the reference plane, potential energy is V = −mgL cos θ; viscous damping enters as Qnc = −b ˙θ in the Euler–Lagrange equation.Source: §A.3 · pp. 18–19
• Setting x1 = θ and x2 = ω = ˙θ reduces the damped second-order equation to a first-order autonomous state-space system using normalized damping c = b mL2.Source: §A.3 · p. 19

A.4 Derivation of Damped Double Pendulum Dynamics

• The double-pendulum derivation defines link coordinates and velocities, then uses them to formulate kinetic and potential energy.The kinetic-energy expression includes a coupling term proportional to cos(θ1 − θ2).Source: §A.4 · p. 19
• Rayleigh dissipation forces are introduced for both angular coordinates before deriving the coupled Euler–Lagrange equations.Source: §A.4 · p. 19
• The angular accelerations are represented as a linear matrix system, whose inversion yields explicit state-space equations using denominator parameter D.The derivation defines angular velocities ω1 and ω2 and angular difference δ = θ1 − θ2.Source: §A.4 · p. 19

B Hyperparameter Reference

• Table 5 lists the configuration parameters used across the single- and double-pendulum experiments.Source: B Hyperparameter Reference · p. 20, Table 5 · p. 20

C Analysis of High-Velocity Error Outliers

• At high angular velocities, polynomial-filter approximation errors coincided with localized double-pendulum discrepancies.In failing seeds, training loss was 5.9×10−3 while field errors exceeded 0.11; one example state was x = [−0.73, −0.30, 4.67, −7.80].Source: C Analysis of High-Velocity Error… · p. 20
• Peak field error fell from 0.135 to 0.025 with polynomial-misfit quality gating.The gate detects filter violations directly from raw sensor windows; standard robust losses such as Huber do not automatically reject biased targets that neural networks fit tightly.Source: C Analysis of High-Velocity Error… · p. 20
• Huber-style robust losses do not automatically reject biased targets that neural networks fit tightly.Polynomial-misfit verification instead detects these violations from raw sensor windows.Source: C Analysis of High-Velocity Error… · p. 20

D Extended Per-Seed Results

• Table 6 reports individual per-seed relative field error values across test configurations.The accompanying prose identifies the table as detailed numerical results across individual random seeds.Source: Table 6 · p. 20, D Extended Per-Seed Results · p. 20

E Supplementary Gradient and Estimator Validation Plots

• Figures 10 and 11 validate analytic jet-matching gradients and baseline discrete-adjoint gradients against complex-step differentiation.The validations cover 32 and 14 test conditions, respectively.Source: Figure 10 · p. 21, Figure 11 · p. 21
• Figure 12 presents Savitzky–Golay derivative-estimation error across effective step sizes and noise levels.Source: Figure 12 · p. 21
• Figure 13 examines how observation sampling rate affects baseline adjoint field error.Source: Figure 13 · p. 22