Learning long range dependencies through time reversal symmetry breaking
Authors: Guillaume Pourcel
Maxence Ernoult
Venue: NeurIPS 2025 oral ()
Keywords: physical learning
continuous adjoint method
temporal credit assignment
Hamiltonian systems
analog computing
training algorithms
state space models
Primary area: machine_learning_for_sciences
TL;DR: We propose a backward-mode AD proxy using only forward passes applying to Hamiltonian recurrent units and stacks thereof (namely, SSMs) with theoretical guarantees and experimental evidence

Abstract: Deep State Space Models (SSMs) reignite physics-grounded compute paradigms, as RNNs could natively be embodied into dynamical systems. This calls for dedicated learning algorithms obeying to core physical principles,  with efficient techniques to simulate these systems and guide their design. We propose \emph{Recurrent Hamiltonian Echo Learning} (RHEL), an algorithm which provably computes loss gradients as finite differences of physical trajectories of non-dissipative, \emph{Hamiltonian systems}. In ML terms, RHEL only requires three ``forward passes'' irrespective of model size, without explicit Jacobian computation, nor incurring any variance in the gradient estimation. Motivated by the potential to implement our algorithm in non-digital physical systems, we first 
introduce RHEL in \emph{continuous time} and demonstrate its formal equivalence with the continuous adjoint state method.
To facilitate the simulation of Hamiltonian systems trained by RHEL, we propose a \emph{discrete-time} version of RHEL which is equivalent to Backpropagation Through Time (BPTT) when applied to a class of recurrent modules which we call \emph{Hamiltonian Recurrent Units} (HRUs). This setting allows us to demonstrate the scalability of RHEL by generalizing these results to hierarchies of HRUs, which we call \emph{Hamiltonian SSMs} (HSSMs). We apply RHEL to train HSSMs with linear and nonlinear dynamics on a variety of time-series tasks ranging from mid-range to long-range classification and regression with sequence length reaching $\sim 50k$. We show that RHEL consistently matches the performance of BPTT across all models and tasks. This work opens new doors for the design of scalable, energy-efficient physical systems endowed with self-learning capabilities for sequence modelling.

Learning long range dependencies through time reversal symmetry breaking
Guillaume Pourcel University of Groningen ", INRIA g a. pourcelérugg.
Maxence Ernoult Google DeepMind mernoult@google. com
Abstract
Deep State Space Models (SSMs) reignite srspiis-g compute paradigms, RNNs could natively be embodied into dynamical systems. This calls for dedicated learning algorithms obeying to core physical principles, with efficient techniques to simulate these systems and guide their design. We propose Recurrent Hamiltonian Echo Leurning (RHEL), an algorithm which provably computes loss gradients finite differences of physical trajectories of non-dissipative, Hamiltonian systems. In MI. terms, RHEL only requires three "forward passes" irrespective of model size, without explicit Jacobian computation, nor incurring any variance in the gradient estimation. Motivated by the potential to implement our algorithm in non- digiial physical systems, we first introduce RHEL in continuous time and demonstrate its formal equivalence with the continuous adjoint state method. To facilitate the simulation of Hamiltonian systems trained by RHEL, we propose discrete-tine version of RHEL which is equivalent Backpropagation Through Time (BPTT) when applied to a class of recurrent modules which we call Hamiltonian Recurrent Unity (HRUs). This setting allows generalizing these results to hierarchies of HRUs, which we call Hamiltonian SSMs (HSSMs). We apply RHEL to train HSSMs with linear and nonlinear dynamics variety of time-series tasks ranging from mid-range to long-range classification and regression with sequence length reaching 50k: We show that RHEL consistently matches the performance of BPTT across all models and tasks 3) This work opens new doors for the design of scalable, energy-efficien physical systems endowed with self-learning capabilities for sequence modelling.
demonstrate the scalability of RHEL by
Introduction
The resurgence of Recurrent Neural Networks (RNNs) for sequence modeling 11. particularly when integrated with State Space Models (SSMs) 1 reopens the debate about the "hardware lottery".
This raises critical question with high stakes 7 should future hardware development continue to
optimize for the long- long-dominant GPU/IPU-Transformer backkrop paradigm {8) or should it explore alternative computational approaches alongside novel training algorithms?
This paper embraces distinctive standpoint in the realm of AI hardware by distinguishing algorithms from hardware [9] and posits that some SSMs could be physically realized in dynamical physical systems and, when endowed with bespoke credit assignment mechanisms, be turned into "self- searniner machines [10] I|I see Section 5 for detailed discussion. To this end, such algorithms should fulfill at least two requirements: use "forward passes" only (i.e. no backward
"Bemoulli Institute & CogniGrron University of I. Lille, Inria, CNRS, Centrale Lille, UMR CRIStAL Ouu code is available on:
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

passes), ii) do not use explicit state-Jacobian (i.e. Jacobian of the system's 's dynamics). The vast majority of existing algorithms endowed with these features revolve around forward-mode automatic differentiation (AD) [2]. The standard forward-mods mode AD algorithm for temporal processing is
Recurrent Real-Time Learning [33 (RTRL.. Due to its cubic memory complexity with respect to the number of neurons, RTRL can only be exactly implemented on architectures of limited size [14]
and either requires low-rank approximations [15]. or hardcoded [16] or metalearned [17 heuristics.
RTRL fails to satisfy criteria ii) because it explic- itly uses the Jacobian to compute directional deriva- tives of the loss L along each direction v2 of the canonical basis of parameters as a supplemen- tary computation during the forward pass. How- ever, it can be reformulated as a zeroth-order pro- cedure by approximating the directional derivatives via Vol. = lim 6 (L + cv) I((O). This recipe theoretically aligns with our three algorithmic requirements but incurs a prohibitive complexity cost by requiring separate forward passes for each pertur- bation along every parameter direction. Alternative methods attempt to circumvent this by sampling ran- dom directions from the parameter space [18 21] rather than exhaustively computing gradients along the entire canonical basis. However, these approaches suffer from either excessive variance [22] or high bias (23. limiting their effectiveness to small-scale appli- cations 24 or fine-tuning tasks (25). To avoid the pitfalls of RTRL and other forward-mode AD prox- ies, a better approach would be to instead emulate backward-mode AD. forward in time [26). Yet, such techniques have only developed to emulate backward- mode implicit differentiation
y
40
Φo
TTo
-To
P-T
2 d
Φ TTg VVgL)
Figure 1: НEB core mechanies [10]. Af- ter following free forward trajectory (blue curve) and undergoing "rebound", the neu- rons exactly travel backward (grey curve). In- stead, HEB prescribes nudging the "echo" " tra- jectory (red curve) closer to y before rebound, with the resulting position gap encoding the error gradient with respect to momentum (the left green arrow).
energy- -bsed models on static inputs [27] and have not been extended out of equilibrium to sequential data.
In this work, we propose an alternative approach that emulates backward-mods AD forward time, inspired by the Hamiltonian Echo Backprop (HEB) algorithm [10), illustrated in Fig. D for a single neuron. Consider a system of neurons described by their position di and momentum ग. which do not dissipate energy. This means, for instance for coupled oscillators (Eq. (8). that the initial energy provided to this system is preserved and transferred across oscillators from kinetic energy to elastic potential or vice versa. Such systems can be modelled by the Hamiltonian formalism (section 2.1 Let us say that we want to learn the initial conditions on this system such that reaches some target y after some time. HEB proceeds in three steps. For the first step, the system evolves freely, but does not reach y. For the second step. HEB perturbs the trajectory for brief duration to drive the system slightly closer to y per some metric L. Then in the last step, the neurons are "bounced"" backwards, causing the system to evolve backward. the perturbation was omitted, the system would exactly symmetry.
retrace its previous trajectory backward in time, a property which is called Because the perturbation breaks time-reversal symmetry, the system no longer retraces its initial trajectory exactly-the resulting deviation encodes the gradient of L with respect to its initial state.
In spite of its significant implications for physical learning. HEB has multiple features which may hinder its broader investigation within the ML community: i) HEB is difficult to compare to standard ML algorithmic basclines as it was derived with dedicated theoretical physics tools and is entangled with fine- aine-grined physics of the systems being trained; ii) HEB theory assumes that model parameters are also dynamical variables and that only their initial state is learnable, which is in stark contrast with standard RNN parametrization; iii) is unclear how HEB would extend to SSMs, f.e. in discrete time, on sequential data with losses defined at all timesteps, on hierarchical recurrent units; iv) finally, HEB was only evaluated on small static problems (e.g. XOR, MNIST) and not proved at larger scale.
Taking inspiration from HEB and using standard ML tools, we propose Recurrent Hamiltonian Echo Learning (RHEL) simple, general and scalable forward-only only proxy of backward-mon mode AD applying to a broad class of dissipative-Fi free Hamiltonian models with the following key contributions:
2

---

With the primary goal in mind to inform the design of self-learning physical systems and facilitate its comparison HEB, we first introduce RHEL continuous time. We show that RHEL generalizes HEB to a broader class of models and problems and demonstrate its equivalence, at every timestep, with the continuous adjoint state method [28] in the limit of small trajectory perturbations (Theorem B.I.. We numerically highlight this property
toy physical model (Eq, ( Fig.
To efficiently simulate this algorithm, extend RHEL discrete time and demonstrate its cquivalence with Backpropagation Through Time (BPTT) (Theorem (.2) on Hamiltonian Recurrent Units (HRUs) which are discrete-time, symplectic integrators of separable Hamiltonian dynamics. We further generalize this result to learning a hierarchy of HRUs, which we call Hamiltonian State Space Models (HSSMs, Fig. 2p. and propose RHL chaining procedure accordingly (Theorem 3.3).
With this simulation toolbox in hand, finally demonstrate the effectiveness of the proposed approach on HSSMs with linear [29] and nonlinear [30] recurrent units. We show that: i) gradients estimates produced by RHEL near perfectly match gradients computed by end-to-end AD (Fig. 4. ii) RHEL remains on par with AD in terms of resulting model
performance across classification and regression long-range tasks (Tables 2
2 Problem statement
2.1 Notations & model description.
Notations. Given a differentiable mapping H RRd - IR, we denote V H E Rd and V Ve, E IRd, xdoz the gradient and Hessian of H respectively. When necessary, we use the Leibniz notation to denote as V,H the gradient of И with respect to the 7th variable. When considering a differentiable mapping M I R", we denote its Jacobian matrix as θ. E gaxd We also use the notation de, to denote total derivative with respect some variable v1 which accounts for both direct and indirect effects of the variable in on the function being differentiated.
Continuous Hamiltonian model. Following 10] we model neurons as a vector & E Re
(Φ T). comprising a position vector & E R de and momentum vector R " such that
We define Σ. R xda such that Σ. + (Φ) -T) and Σg Rdaxxx such that Σ, 4 (T, Denoting 0 θ de and
U Rdd.. the model parameters and inputs, define the Hamiltonian associated to the model as a mapping H S x θ x + R. Taking -T as the origin of time and starting from (--T) x, say that & follows Hamiltonian dynamics under sequence of inputs t + (t) and with some parameters 0 if satifies the ordinary differential equation (ODE):
I
Vt E : (t) V (t). O, (t)), J == -1 0
(1)
Assumptions. The most fundamental requirement the whole proposed approach for the Hamilto- nian models at is time-reversal symmetry. Heuristically, if we were recording the dynamics of conversative system described by Eq. (I) and playing the resulting recording forward and backward, we could not distinguish them both are physically feusible (Lemma A.3) A direct consequence of this fact is that if neurons are let to evolve for some time and then bounced back, i.e. their momenta reversed, they will exactly travel back their initial state (Corollary A.3) 3 Namely (see Fig. T:
+ (--))
(())
x*(0)= = Eg 400
4- (-1'))
(2)
Eq.
Eq. T
Learning as constrained optimization
Problem formulation. Our goal can be framed as constrained optimization problem where we aim find a set model parameters such that, under some input sequence - u(l), ()) the model trajectory approaches much as possible some target trajectory, as measured by a loss function L, under the constraint that the model trajectory is physically feasible in the sense of satifying Eq. (T)
min = (0) dil del (1), 0
s.t. VI E 00 (1)
V (1).
(3)
3

---

(A)
B
Φ (N)
00
0000
Φ (N-1)
(N-1)
(2(()
@
yit
(C)
X20
MICL
I
KK
K
15
TE (
T
Figure (A): Learning Hamiltonian SSMs by RHEL. The forward pass of HSSM (left, Eq. (15 reads composition of Hamiltonian Recurrent Units (HRUs, Eq. (9). During the backward pass (right), the momenta the neurons of the top-most HRU are flipped and nudged by the first error signal, yielding perturbed trajectory (red curve). The contrast with the time-reversed trajectory (which would be obtained without nudging, dotted blue curve) yields an error signal that is passed backward (purple) to previous HRUs (Alg. [T. (B)-(C): checking Theorem 3.1 on toy model. Six coupled harmonic oscillators with learnable parameters under input u and target (Eq. (8) Plots show gradients for spring (dark green) and mass parameter (light green) comparing C-ASM (L + go (t), dotted) and RHEL (1 - (!),), solid), alongside sensitivities for oscillator positions di (orange) and of (blue) under both methods (t Xl), dotted; As (4), solid).
where T. reads as the sum of cost functions t -> E(., t). The most common approach solve Eq. (3) by gradient descent or variants thereof such that the problem boils down to computing Le Note that the problem defined by Eq. 33 is more general than the one solved by the seminal HEB work [[O] as: i) the loss function L is defined over the whole trajectory, ii) the Hamiltonian is time-dependent (through t + ut)) and parametrized by 0 which is shared across the whole computational graph.
Algorithmic baseline. One standard approach to compute dg J. is the continuous state method (ASM) [28]. which can be simply regarded as the continuous counterpart of backward-mooe AD. Given some forward trajectory (+(t))ee()roo) spanned by Eq. (I) this method prescribes solving the following backward ODE:
x0)=00 0,A(t) V 0,
0,u u(-)) JT. A(t) + v t]
(4)
with the gradient of the loss with respect to the initial state of the neurons æ and the model parameters 0 given by (Theorem A.I
deL A(T), doL = S'' (t)dt, with (t ==
JT A(t) (5)
Note that in theory, the forward trajectory (+())){(TT.o) can either be stored recomputed back- wards alongside X yielding ((1) memory cost property which holds beyond Eq. [ specifically BI. However in practice, recomputing variables backward through discrete computational graph is generally inexact and induces bias in the gradient estimation [32]. unless reversible discretization scheme is used. Fortunately, symplectic integrators associated with Hamiltonian flows are reversible 33 a property which has been leveraged yield memory savings in neural networks [Q Therefore, the continuous ASM as well as our proposed algorithm (Alg. [T] also naturally inherit this memory efficiency as a model feature rather than feature of the training algorithms themselves.

---

3 Recurrent Hamiltonian Echo Learning (RHEL)
3.1 Definition in continuous time & equivalence with continuous ASM
We are now equipped to introduce Recurrent Hamiltonian Echo Learning (RHEL.. In comparison to the continuous ASM, RHEL does not require solving a separate adjoint ODE akin to Eq. (4) Instead, RHEL simply prescribes running multiple additional times the forward trajectory Eq. 1 with three modifications: i) the neurons first conjugated, i.e. 1 + x=; ii) the inputs t - u(-t) are processed backwards; iii) the trajectory of the neurons is slightly nudged, with a strength c, towards direction of decreasing loss values with cost functions t - f[ -f processed backwards too. More
precisely, let us assume a "free" forward trajectory Eq. ①T) has been executed between -7 and 0 yielding some neuron state (0). Then, define the echo dynamics of the neurons of through: (0) = 4 (0), Ö, (t, = JV a H [ (tt (),c), u(-t))-ccJV (1, c),
Vt e [0. 7
crucial observation is that when e = O, the echo trajectory simply matches the forward trajectory in reverse: & (1, 0) (-() for i E [0, (Lemma A.3)) reflecting time-reversal symmetry. When e + 0, this symmetry is broken: the resulting trajectory difference ((, - (t, = 0)
implicitly encodes for the error signals carried by the continuous ASM. We introduce this result more formally in Theorem [3I with the help of the following quantities capturing the essence of RHEL as a forward-only, siteeennce-bssee gradient estimator:
AB ARHEL(t. (t, = 2c (V (t, (V2Hl0
c). ct.
O. ((t())]
- -t, -),, 8 u(t)))
ARHEL (1,c) == 2c Eg (L,) $ (1, (1,-c))
(7)
Theorem 3.1 (Informal). Under mild assumptions on the Hamiltonian function and with t -> A(t)
and (t) defined by Eq. (4) 5. the continuous ASM and RHEL are equivalent at all times in the sense that:
X() = lim A WE [0, T),
RHEL
(1,
gl(()
lit A <
ARHEL (1,e)
c0
Proof sketch. Defining A ARHEL(()) := lim A RHEIL (t, e) and noticing that 4 (BHEL((t) Σg ( c))le=o. we can show by differentiating Eq. (6) around = O and with some algebra that t + 4 (t) satisfies the same ODE as + X(t) (Eq. <4 and have same initial conditions. therefore are equal at all times. The other two equalities of Theorem 3.1 can then be easily deduced. The full formal statement and proof are provided Appendix A.1.3
Example. We numerically verify Theorem 3.1 on a toy model of six coupled harmonic oscillators i (Fig, 2) with masses m, and spring parameters k;, ky. The system is driven by an external input force u(() applied to oscillator I and produces an output on oscillator 4, which is nudged toward a target trajectory y(t) during the echo phase. The corresponding Hamiltonian reads:
T
11 L 0, u]
kis 0.)² + udi,
(8)
2m
with parameters 0 = {m,, ki, Kig hig. This yields the following dynamics and associated RHEL gradient estimators (see App. A.4.I)
Vi E 6], 0.0; = mi,
OCT; = k; +
Σ kin - Φi) + da + 8.4 de c(os - v),
Ike (1) lim (0) (1, (0) (c) - (1, (1)
c
2c
where u(t) denotes the external driving force applied oscillator 1, y(t) the target trajectory associated with the output oscillator 4, d.j the Kronecker delta, and ó, an indicator equal to during the echo phase and otherwise.
3.2 Extension to discrete time & equivalence with BPTT
In this section, we propose a discrete-time version of RHEL by first defining family of discrete models as integrators of the Hamiltonian flow which exactly preserve time-reversal symmetry, defining an associated surrogate problem and then solving it. This disscretzze-thenoptiiie approach ensures a better gradient estimation than directly discretizing RHEL theory in continuous time [32].
5

---

Hamiltonian Recurrent Units (HRUs). We
QK
@k+2/3
@k+1
now assume that the neurons & obey the following discrete-time dynamics, starting from $ - and:
V
V V
T
T x++1 MAH 0 Uk
4k K
1.
TCk
T3-1333
TCK+2/3
Tlk+1
where M denotes the integrator associated with the Hamiltonian function H with time dis- cretization 6. We classically pick MH. as the Leapfrog integrator [33] and restrict ourselves to separable Hamiltonians of the form H [ T[]] + V[Ø] yield reversible and explicit integration scheme. One possible parametrization of three explicit Euler integrators, which alternate between updates of the position & (first and third steps) and of the momentum (second step) see Def. AA33 detailed description and Fig. 3 for the associated computational graph. As these two intermediate integrator time steps needed to accurately define the RHEL learning rule in discrete time, we explicitly denote them as fractional time units:
Figure 3: Computational graph of M.
the Leapfrog integrator in this case is as the composition
MH $ - k+1/3 Pk+2/3
4
(10)
We call such models Hamiltonian Recurrent Units (HRUs). Therefore by design, the time-reversal symmetry property defined in Eq. (2) extends in discrete time to HRUs (Corollary A.6).
Surrogate learning problem. Given this modelling choice, the continuous-time problem introduced in Eq. 3) naturally translates here in discrete time as:
min L =
S.t.
4x+1
MH.s 0, U. Vk -
(11)
0
k== K+1
Algorithmic baseline. Eq. defines a classical RNN learning problem where Backpropagation Through Time (BPTT), iee. the instantiation of backward-mode AD in this context, natural algorithmic baseline. BPTT simply amounts to apply the "chain rule" backward through the compu- tational graph spanned by the forward pass of a HRU (Eq. (D). Alternatively, it also be regarded as the discrete counterpart of the continuous ASM previously introduced (Theorem A.3). .. For this reason, we re-use the same notations above and define the following quantities associated BPTT:
K-
L
Σ go (k)) 0=
(k( = do I., (k)) = I ÀE= := do L
(12)
where (k), 9n (k) and A denote the "sensitivity" of the loss L to at time step k; u and $ respectively see App. A.2.2 for more detailed definition and derivation of BPTT.
RHEL discrete time. on HRUs as:
Finally, we extend RHEL in discrete-time by defining the echo dynamics
FC(c) = Tj E. V O,, Φ (+)((c= M
2 (c). O, u (++1)) - cJ Vo- -(k +1] Vk = 0, ... K K-1
(13)
Denoting:
11/ (E) u u- (H 1/3(e), 0, ((+1) 1 ( (E), u (atil). (14)
we can define the discrete-time counterpart of Eq. ( as:
RHEEL c) = V2H (c), 0 u- (k+1) \/H1/2 u u-(++)))
A RHEL ( =
( 0, u-(k))] 1/21 (-e), u-(++)))
A RHEI. (k)
2t E.
(e((c) $ (-c))
and consequently extend Theorem 3.1 in discrete time well.

---

Theorem 3.2 (Informal). Under mild assumptions the Hamiltonian function and with (A*)Rejo,K<11) (ge (go(k)) (k))) and (gu (e)) defined by Eq. 12 BPTT and RHEL are equivalent at all times in the sense that wk = 0, K:
lim A A
RHEL
(k,
(k(k)
lim A ( c-
(k, e), lu(k) lim A RHEL (k. e).
The full formal statement and proof provided in Appendix A.2.3
3.3 Learning stacks of HRUs via RHL chaining
Vectorized notations. Given discrete vector field (Vi)kel- we denote V x''x the forward and backward VTT Ixxv and V := VV
trajectories associated with V. We define the vectorized operator M such that the dynamics of
single HRU unit defined in Eq. (9) and the the nudged dynamies prescribed by RHEL as given in
Eq. 13 rewrite more compactly as:
M
4 u.
T
MH H.5
0, EJ V
L T
Learning Hamiltonian State Space Models (HSSMs) as multi-level optimization. consider hierarchicul models reading composition of N HRU units of the form:
We now
(0) = u, vl [0, N - 1] T (++1)
M
(0)
(++1)
0(000
(15)
H
where the trajectory of the (f)-th HRU unit is fed input sequence into the (E + 1)-th HRU unit (Fig. 2). We call such hierarchies of HRUs Hamiltonian State Space Models (HSSMs). The corresponding optimization problem reads:
min L
E (N)
R:]
(++1)
0000
+
v [O,N 1]
K+1
(16)
Alg. D prescribes an intuitive receipe to compute L L for the above optimization problem (Eq. (6)).
chaining RHEL backward through HRUs. Namely, the echo dynamics of the top-most HRU read as Eq. 13 using the initial learning signal a to nudge its trajectory. On top of estimating its
parameter gradients, also estimate its inpul gradients which used to nudge the echo dynamics of the preceding HRU. This procedure is repeated until reaching the first HRU see Fig 2
Algorithm Recurrent Hamiltonian Echo Learning (RHEL) on a single HRU
Inputs:
o (final state of the forward trajectory), A (incoming gradient), (nudging strength), (timestep)
Outputs: AO (parameter gradient estimate), An (input gradient estimate)
1: A0 Ode
2: A, = (An x." Au,) < Ox
3: T
/6 - =EE-
A
Compute $ in parallel
4:
for k in ..... K do
5:
di. + MHH De
ter
O, u (1++()+)) EEE As.- (k--(
6: AO + AO-
2 (V 1/2 & O, u-(kt) - 12H1/2 c 0, u -((++1))))
Eq.
14
7: Du t (V ch O u- ((+1)) V3H 1/2 . O. u (-(++1)))
8: end for
9:
return AO, A.
Theorem 3.3 (Informal). Given u input sequence and HSSM with layer-wise Hamiltonians (H applying Alg. recursively backward from the top-most HRU solves the multilevel optimization problem defined in Eq. (16) in the sense that:
44 0. N dgin lim ( c++
The full formal statement and proof are provided Appendix
A.4
7

---

4 Experiments
Foreword. We first introduce the two types HRUs at use, empirically assess the validity of Theo- rem 3.3 by statically comparing gradients computed by RHEL and BPTT on HSSMs made up of these two types of HRUs, then perform training experiments across classification and regression sequence tasks. Importantly, the goal of the proposed experiments is not to improve the SOTA performance on these tasks. Rather we want show, on a given model satisfying the requirements of RHEL (i.e. being a HRU or stack thereof), that RHEL maintains training performance with respect to BPTT.
Linear HRU block. Following recent work 22 we introduce the linear HRU block with the parametrization 0 {A(O, ((9) W(O and Hamiltonian:
{ H1 1((0((1 000) q = (012 + 1000" . -0000 B u
(17)
u Flø w('),
where A E R ((da/2)x((ee2)) /2) is assumed to be diagonal and to have positive entries, B R (dde x2., a E R, and F is the nonlinear spatial building block parametrized by W (see
App. A.4.3 Under these assumptions, it can be shown that such HRUs have a bounded spectrum and yield HSSMs which are universal approximators of continuous and causal operators between time-series (29).
Nonlinear HRU block. To demonstrate the validity of our approach beyond linear HRUs, we also introduce a nonlinear HRU block [30] with (() {A W((), B(t), b(t), a and:
H (() ( T ((-1)) = || (0) || + 2
" 111
+ (A log (cosh (1(0) (6 + B(d) (0) + b()).
(18)
a Fiøl-1) W
where A E (dde/2)x(de) /2) is also assumed to be diagonal, ( E 2 ( E R (ds /2).
E R. and F is the nonlinear spatial building block parametrized by ( (see App. A.4.5 It has been shown that these HRUs mitigate by design vanishing and exploding gradients 30 we provide empirical validation in App. A.5.5
4.1 Static gradient comparison
As sanity check of Theorem 3.3 and preamble to training experiments, we first check that RHEL gradients are computed correctly. Given HSSM model, we randomly sample tuple (u: y:) D from the SCP1 dataset (see App. A.5.I for details on this dataset) and pass through the model. Given y. we then run BPTT and RHEL through the model and compare them in terms of cosine similarity and norm ratio of the resulting layer-wise parameter gradients. We run this experiment on a HSSM made up of six linear HRU blocks (which we call "linear HSSM") and another HSSM comprising six nonlinear HRU blocks (resp. "nonlinear HSSM") and obtain Fig. + We observe that in terms of these two comparison metries, RHEL and BPTT parameter gradients are near-perfern sear-perfectly aligned for all parameters and across all HRU blocks, for both the linear and nonlinear HSSMs.
4.2 Training experiments
Classification. To further check our theoretical guarantees, we now perform training experiments with BPTT and RHEL across six multivariate sequence classification datasets with various sequence lengths (from 400 to 18k) and number of classes, on linear and nonlinear HSSMs (see App. A.5. for details on these datasets) and display our results in Table T This series of tasks was recently introduced (34 as a subset of the University of East Anglia (UEA) datasets (35) with the longest sequences for increased difficulty and recently used to benchmark the linear HSSM previously
introduced (29). We observe that models trained by RHEL almost match, on average across the datasets, those trained by BPTT in terms of resulting performance on the test dataset.
Regression & scalability. Finally, to assess the applicability of RHEL beyond classification and scalability to longer sequences, we run training experiments on the PPG-DaLiA dataset, multivariate
8

---

linear HSSM
nonlinear HSSM
Coaloe Simllaity
1.0
1.0
0.5
0.5
1.0
1.0
Ratio
0.5
0.5
A
B
A B
Figure 4: Static comparison between RHEL and BPTT. Given some (u, y:) ~ D, we perform BPTT and RHEL on six deep linear and nonlinear HSSMs. We measure, layer-wise (first layer in purple), the cosine similarity (top panels) and norm ratio (bottom pannels) between RHEL and BPTT parameter gradients of linear HSSM (Eq. (17) and nonlinear HSSM, Eq. (8))
Tasks
Worms
SCP1
SCP2
Ethanol
Hearbeat
Motor
Seq. length
17,984
896
1,152
1,751
405
3,000
classes
5
2
2
2
2
Avg
BPTT
78.3#7.6
86.4
63.9+7.
29.9
73.9
48.1
63. 63.4
Lin
RHEL
75.0
86. 1:2.9
41.4
29.9
73.5
51.625.0
62.9
BPTT 51.1#7.2
86.8+3.2
54.0+4.9
29.9
74.5
56.5#7.6
58.8+44.8
Nonlin
RHEL
50.6
85.6
54.0+2.0
29.9
73.9
33.0
57.8 57.8
Table I: Test mean accuracy (%. higher is better) across five different seeds (# indicates standard deviation) using RHEL and BPTT, nonlinear and linear HSSMs, on six UEA time series classification datasets with various sequence length and number of classes.
time series regression dataset designed for heart rate prediction using data collected from wrist-worn device (36). With sequence length of 50k, this task considered to be a difficult "long-range" benchmark (29]. We display in Table 2 the results obtained when training linear and nonlinear HSSMs with RHEL and BPTT on this dataset. Here again, we observe that on average and with both models, the performance of RHEL matches that of BPTT.
5 Discussion
Limitations. We observe on Tables 2 that RHEL slightly underperforms BPTT, though this gap is statis- tically significant only for nonlinear HSSMs on regression tasks. This performance difference can be attributed i) the approximation bias introduced by finite nudging in our method, a known issue in related works (37), and numerical precision that requires some careful selection of the nudging strength (see App. A.5.4)
PPG-DaLiA
Seq. length
49,920
Input/output dim.
6/1
BPTT RHEL
9.1
Lin
9.5+1.0
BPTT 7.8
"Physical computing"" Every logically irreversible computation, be it on a digital accelerator any alter- native substrate, is sustained by a physical process carried by wires and transistors which fundamentally generates heat (38). Therefore the very notion of "physical comput- ing" may sound pleonastic. However, digital hardware designs abstract core physics away to enable idealized boolean logic relying (among many other things) upon statelessness, unidirectionality, determinism and synchro-
Nonlin
RHEL
8.4
Table 2: Test average mean-square ror (x 10-2, lower is better) across five different seeds (+: indicates standard de- viation) applying RHEL and BPTT to train nonlinear and linear HSSMs on the PPG-DaLiA dataset.

---

nization [39. Compilation pipelines are then typically designed to be as hardware-agnose and as general- as possible [40]. In this generalistic computing paradigm, one of the most popular codesign is to write bespoke low-level kemels for specific workloads to mitigate their compute, memory costs and off-chip memory accesses [6] 41 42], yet remaining largely physics-agnostie Conversely, physical computing firstly targets Applicatiion-Seccffe Inte- grated Circuits (ASICs) rather than general-purpose machines [39]. Secondly, the aforementioned digital design constraints may be relaxed and the resulting analog physics leveraged. Supply voltages may be decreased and the resulting circuits become non-deterministic and be leveraged for stochastic Continuous valued currents may be used for vector-matrix multiplication 44].
practices of "hardware-softwaat
algorithms 44
matrix inversion 45 nearest neighbor search [46] within resistive systems. Constrained [47 or
combinatorial 48
optimization problems may be embodied into a physical system and subsequently solved as the system settles to equilibrium- see more examples in recent review works 39
Self-learning machines and Hamiltonian systems in ML. As a particular instantiation of physical computing, self-learning machines are physical embodiments of neural networks whose inference and gradient computation are, for instance, carried out by relaxing to equilibrium 49 50]. In this regard, RHEL extends the design of scalable self-learning machines beyond dissipative systems emulating implicit difserentiation (51)55) In comparison to these algorithms, RHEL crucially avoids prohibitively long simulation times due to the use of lengthy root-finding algorithms, extends their core mechanics to the broader domain of sequence modelling and possibly larger-scale tasks. Finally, on top of its direct connection with HEB [10], RHEL belongs to a large body Hamiltonian- based ML algorithms dedicated to physics-informed neural networks (56458), memory-eflicient models (30). generative models [59], sampling 660-62] and optimization 63 algorithms.
Forward-only learning. In mainstream ML, zeroth-order optimization (ZO) techniques are pursued for their memory efficiency compared to vanilla backprop with gradient checkpointing [64], which can be further enhanced with appropriate quantization schemes [65]. The relatively recent exploration of ZO techniques on Large Language Models (LLMs) finetuning [25) has spurred a revival of "forward- only" techniques for memory-cheap gradient computation in LLMs on a variety of applications. However, the quest for forward-only algorithms in physical computing is not primarily motivated by memory efficiency. In our context, fforward-only' means that both the forward and backward passes obey the same physical principle, for instance obeying to the same Hamiltonian dynamics as required by RHEL.
Future work. The restriction of RHEL to non-dissipative systems limits model expressivity [29]. While inducing artificial dissipation through the use dissipative integrators [29] or by embedding dissipative system within larger conservative one [10] are convenient workarounds, there is ultimately a need for a temporal credit assignment algorithm tailored for truly dissipative systems. In Appendix I we show that even when assuming that the time-reversed trajectory of dissipative system is physically feasible, perturbations of this trajectory do not encode ASM error signals. Additionally and in the spirit of HEB, idcalized version of RHEL should be entirely "black-box as it still requires explicit Hamiltonian parameter gradients to implement its learning rule (Alg. [ This could possibly be achieved for instance using homeostatic control via control knobs acting on the model parameters [66]- we sketch inside Appendix (aanndiix/A3.22 what such a procedure would look like. Another exciting direction of work is have RHEL operate online it currently requires one forward pass and subsequent ones revisiting the inputs in reverse order. Finally, while HRU hidden units can be recomputed from their final state [30], it still requires storing the input sequence of the HRU so that memory gains are not evident within a stack of HRUs. Looking beyond analog physical systems, endowing RHEL with better memory efficiency and online mechanics could potentially yield compelling alternative to BPTT on digital accelerators like GPUs [1]..
Conclusion. While we simulated Hamiltonian dynamics alongside RHEL on GPUs, the greatest potential of RHEL lies within real Hamiltonian systems, namely on-chip photonic circuits [67]. superconducting circuits [68], or spintronic systems [69]. possible physical implementation of a HSSM would map its feedforward components onto digital circuits (running backprop) and its recurrent components onto analog circuits (running RHEL) [54] see Remark [8 Appendix. We hope this work will incentivize the search for alternative analog hardware capable of training sequence models with much higher energy efficiency, as well as alternative algorithms for digital hardware.
10

---

Acknowledgments and Disclosure of Funding
The authors warmly thank (in alphabetic order) Aditya Gilra, Albert Cohen, Debabrota Basu and James Laudon for their support of the project. ME thanks Vincent Roulet for his careful review which greatly helped improve the quality of this manuseript, as well Amaury Hayat and Alexandre Krajenbrink for useful discussions. We also thank the anonymous reviewers for the very interesting research questions they raised and for greatly helping improve the clarity of the manuscript. Lastly, we thank the donors of the UEA data who indirectly contributed to facilitate the assessment of our algorithm on long-range tasks. GP acknowledges the CHIST-ERA grant for the CausalXRL. project AHIIST-EAA-1--AI-AA ANR-21-CHR4-0007) and the ANR JCJC for the REPUBLIC project (ANR-2-CH23-0003-000
by I.' Agence Nationale de la Recherche, France (grant reference
Author contributions
GP and ME collaborated closely all aspects of this work. Both authors contributed to study design, with GP focusing primarily on experiments and code development, and ME focusing primarily on theoretical development and writing.
References
[1] Antonio Orvieto, Samuel I. Smith, Albert Gu, Anushan Fernando, Caglar Gulcehre, Razvan Pascanu, and Soham De. Resurrecting recurrent neural networks for long sequences. In International Conference on Machine Learning, pages
PMLR, 2023.
[2] Albert Gu, Isys Johnson, Karan Goel, Khaled Saab, Tri Dao, Atri Rudra, and Christopher Ré. Combining recurrent, convolutional, and continuous-time models with linear state space layers. Advances in neural information processing systems, 34:572-585, 2021.
[3] Albert Gu, Karan Goel, and Christopher Ré. Efficiently modeling long sequences with structured state spaces. arxiv 2021. arXiv preprint arXiv:2111.00396, 2021.
[4] Ankit Gupta, Albert Gu, and Jonathan Berant. Diagonal state spaces are as effective as structured state spaces. Advances in Neural Information Processing Systems, 35:22982 22222922-22994, 2022.
[5] Jimmy TH Smith, Andrew Warrington, and Scott W Linderman. Simplified state space layers for sequence modeling. arXiv preprint 2208 04433,, 2022.
[6] Albert Gu and Tri Dao. Mamba: Sinear-time sequence modeling with selective state spaces. arXiv preprini arXiv:2312.00755: 2023.
[7] Sara Hooker. The hardware lottery. Communications of the ACM, 64(12):58 65, 2021.
[8] Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Eukasz Kaiser, and Illia Polosukhin. Attention all you need. Advances in neural information processing systems, 30, 2017.
[9] Jérémie Laydevant, Logan G Wright, Tianyu Wang, and Peter L McMahon. The hardware is the software. Neuron, 11((2):800-88)) 2024.
[10] Victor Lopez-Pastor and Florian Marquardt. Self-learnin machines based on hamiltonian echo backpropagation. Physical Review X, 13(3):031020, 2023.
[11] Ali Momeni, Babak Rahmani, Benjamin Scellier, Logan G Wright, Peter L McMahon, Clara Wanjura, Yuhang Li, Anas Skalli, Natalia G Berloff, Tatsuhiro Onodera, et al. Training of physical neural networks. arXiv preprint arXiv:2406.03.4 030-40060337,,. 2024.
[12] Atilim Gunes Baydin, Barak A Pearlmutter, Alexey Andreyevich Radul, and Jeffrey Mark Siskind. Automatic differentiation in machine learning: survey. Journal of machine learning research, 18(153):1 1-3, 2018.
[13] Ronald Williams and David Zipser. learning algorithm for continually running fully recurrent neural networks. Neural computation, 1(2):270-280 1989.

---

[14] Nicolas Zucchet, Robert Meier, Simon Schug, Asier Mujika, and João Sacramento. Online learning of long-range dependencies. Advances in Neural Information Processing Systems, 36: 10477-10493, 2023.
[15] Corentin Tallec and Yann Ollivier. Unbiased online recurrent optimization. arXiv preprint arXiv:1702.05043, 2017.
[16] Guillaume Bellec, Franz Scherr, Anand Subramoney, Elias Hajek, Darjan Salaj, Robert Leg- enstein, and Wolfgang Maass. A solution to the learning dilemma for recurrent networks of spiking neurons. Nature communications, 11(1):3625, 2020.
[17] Jamie Lohoff, Anil Kaya, Florian Assmuth, and Emre Neftci. A truly sparse and general implementation of gradient-based synaptic plasticity. arXiv preprint arXiv:2501.11407 2025.
[18] Rie Johnson and Tong Zhang. Accelerating stochastic gradient descent using predictive variance reduction. Advances in neural information processing systems, 26, 2013.
[19] David Kozak, Stephen Becker, Alireza Doostan, and Luis Tenorio. Stochastic subspace descent. arXiv preprint arXiv: 19040009004011455 2019.
[20] James C. Spall. Multivariate stochastic approximation using a simultaneous perturbation gradient approximation. IEEE Transactions on Automatic Control, 37(3):332-341 1992. doi: 10.1 109/9. 19632. Describes the simultaneous perturbation stochastic approximation algorthancoottnnkkkeeeennenoccce Jindexx1.
[21] lla Fiete, Michale Fee, and Sebastian Seung. Model of birdsong learning based on gradient estimation by dynamic perturbation of neural conductances. Journal of neurophysioloy, 98(4): 2038-2057 2007.
[22] Mengye Ren, Simon Kornblith, Renjie Liao, and Geoffrey Hinton. Scaling forward gradient with local losses. arXiv preprint arXiv:2210.033) 2022.
[23] David Silver, Anirudh Goyal, Ivo Danihelka, Matteo Hessel, and Hado van Hasselt. Learning by directional gradient descent. In International Conference on Learning Representations, 2021.
[24] Louis Fournier, Stéphane Rivaud, Eugene Belilovsky, Michael Fickenberg, and Edouard Oyallon. Can forward gradient match backpropagation? In International Conference on Machine Learning, pages 10249-10264. PMLR, 2023.
[25] Sadhika Malladi, Tianyu Gao, Eshaan Nichani, Alex Damian, Jason D Lee, Danqi Chen, and Sanjeev Arora. Fine-tuning language models with just forward passes. Advances in Neural Information Processing Systems, 66:53038-53077, 2023.
[26] Maxence Ernoult, Rasmus Høier, and Jack Kendall. A cookbook for hardware-frieedll implicit learning on static data. In NeurIPS 2024 Workshop Machine Learning with new Compute Paradigms, 2024.
[27] Benjamin Scellier and Yoshua Bengio. Equilibrium propagation: Bridging the gap between energy-based models and backpropagation. Frontiers in computational neuroscience, 11:24, 2017.
[28] Lev Semenovich Pontryagin. Mathematical theory of optimal processes. Routledge, 1985.
[29] Konstantin Rusch and Daniela Rus. Oscillatory state-space models. arXiv preprint arXiv:2410.03944, 2024.
[30] T Konstantin Rusch and Siddhartha Mishra. Unicornn: A recurrent model for learning very long time dependencies. In International Conference on Machine Learning, pages 9168-9178. PMLR, 2021.
[31] Ricky TQ Chen, Yulia Rubanova, Jesse Bettencourt, and David K Duvenaud. Neural ordinary differential equations. Advances in neural information processing systems, 31, 2018.
[32] Mathieu Blondel and Vincent Roulet. The elements of differentiable programming. arXiv preprint 2403:22003.14606, 2024.
12

---

[33] Habib Ammari, W Wei, and Y Sanghyeon. Numerical methods for ordinary differential equations. Notes of Course ETH Zürich, 2018.
[34] Benjamin Walker, Andrew D McLeod, Tiexin Qin, Yichuan Cheng, Haoliang Li. and Terry Lyons. Log neural controlled differential equations: The lie brackets make a difference. arXiv preprint wrXiv:2402.185/11 2024.
[35] Anthony Bagnall, Hoang Anh Dau, Jason Lines, Michael Flynn, James Large, Aaron Bostrom, Paul Southam, and Eamonn Keogh. The uea multivariate time series classification archive, 2018. arXiv preprint arXiv:/811.0000 2018.
[36] Attila Reiss, Ina Indlekofer, Philip Schmidt, and Kristof Van Laerhoven. Deep ppg: Large-scale heart rate estimation with convolutional neural networks. Sensors, 19(14):3079) 2019.
[37] Axel Laborieux, Maxence Ernoult, Benjamin Scellier, Yoshua Bengio, Julie Grollier, and Damien Querlioz. Scaling equilibrium propagation deep convnets by drastically reducing its gradient estimator bias. Frontiers in neuroscience, 15:633674, 2021.
[38] David H. Wolpert. The stochastic thermodynamics of computation. Journal of Physics A: Mathematical and Theoretical, 52(19))193001, 2019. doi: 10.1088/1551-81211ab000 Invited review on stochastic thermodynamics and comprtiincoomeelleeeeeeeoooooeer
[39] Maxwell Aifer, Zach Belateche, Suraj Bramhavar, Kerem Y. Camsari, Patrick Coles, Gavin Crooks, Douglas J. Durian, Andrea J. Liu, Anastasia Marchenkova, Antonio J. Martinez, Peter L. McMahon, Faris Sbahi, Benjamin Weiner, and Logan G. Wright. Solving the compute crisis with physics-based ASICs. arXiv preprint arXiv:2507. /0463, July 2025. URL. [httttts //arxiv. wttp:/og///abb/eeeeeeeeee 10463 Discusses relaxing digital constraints and leveraging analog plyssisccomeemkkeeeeeee///ooeess hwwes//eeeeeeeeeeeeeeeeeeeeeeeee (nn.aiiit:ooiiiee6....l.
[40] Chris Lattner, Mehdi Amini, Uday Bondhugula, Albert Cohen, Andy Davis, Jacques Pienaar, River Riddle, Tatiana Shpeisman, Nicolas Vasilache, and Oleksandr Zinenko. Mlir: Scaling compiler infrastructure for domain specific computation. In 2021 /EEE/ACM International Symposium on Code Generation and Optimization (CGO), pages 2-14. IEEE, 2021.
[41] Noam Shazeer. Fast transformer decoding: One write-head all you need. arXiv preprint 1:X1v:191..02150, 2019.
[42] Tri Dao, Dan Fu, Stefano Ermon, Atri Rudra, and Christopher Ré. Flashattention: Fast and memory--fficient exact attention with io-awareness. Advances in neural information processing systems, 35:16344-16359, 2022.
[43] Patrick J. Coles, Collin Szezepanski. Denis Melanson, Kaelan Donaleta, Anto- nio J. Martinez, and Faris Sbahi. Thermodynamic Al and the fluctuation fron- tier. arXiv preprint arXiv:2302.0658/ 2023. URL htttp:///axivv wwwaaii..o///bba..... org/aba/ 2302.06584 Proposes a thermodynamic framework for Al hardware and algo- 10jndexx00000
[44] Geoffrey W. Burr, Robert M. Shelby, Abu Sehastian, Sangbum Kim, Seyoung Kim, Severin Sidler, Kumar Virwani, Masatoshi Ishii, Pritish Narayanan, Alessandro Fumarola. Lucas Sanches, Irem Boybat, Manuel Le Gallo, Kibong Moon, Jiyoo Woo, Hyunsang Hwang, and Yusuf Leblebici. Neuromorphic computing using non-volatile memory. Advances in Physics: X, 2(1):89-124, 2017. doi: 10. 11000/3746199.120 49 1259585. Review article on neuromorphic lindex=11.
computing using non-volatile geemooresconteeneeeeeeeeeeoooooo
[45] Ju-Seog Jang, Soo- Sou-Young Lee, and Sang-Yune Yung Shin. An optimization network for matrix inversion. In Neural information processing systems, 1987.
[46] Geethan Karunaratne, Manuel Le Gallo, Giovanni Cherubini, Luca Benini, Abbas Rahimi, and Abu Sebastian. In-memory hyperdimensional computing. Nature Electronics, 3(6):327-337, 2020.
13

---

[47] Sri Krishna Vadlamani, Tianyao Patrick Xiao, and Eli Yablonovitch. Physics successfully implements lagrange multiplier optimization. Proceedings of the Na- tional Academy of Sciences, 117(43):6633-26665)) 2020. doi: 10.
2015192117. Demonstrates that physical systems implement I. Lagrange multiplier optimiza- 13]ddexx=13.
[48] Kirill P Kalinin, Jannes Gladrow, Jiaqi Chu, James H Clegg, Daniel Cletheroe, Douglas J Kelly, Babak Rahmani, Grace Brenman, Burcu Canakei, Fabian Falck, et al. Analog optical computer for ai inference and combinatorial optimization. Nature, pages 1-8, 2025.
[49] Menachem Stern, Daniel Hexner, Jason W. Rocks, and Andrea J. Liu. Supervised learning in physical networks: From machine learning to learning machines. Physical Review X, 11(2): 02 021045, 2021. doi: 10.1 103/ 11.0210455 Develops local learning rules in physical networrssconteenerfeeeeeeeeenood
15j
[50] Jack Kendall, Ross Pantone, Kalpana Manickavasagam, Yoshua Bengio, and Ben- jamin Scellier. Training end-to-end analog neural networks with equilibrium prop- agation. arXiv preprint arXiv:2006. 01981, 2020. URL https://arxii org/abs/ 20066. 2006.01981 Introduces equilibrium propagation for training analog neural net- (outee/eeeoetttteeeeeeee/ottttee
Yorks:contennReefrceceoooo
[51] Menachem Stern, Daniel Hexner, Jason W Rocks, and Andrea J Liu. Supervised learning in physical networks: From machine learning to learning machines. Physical Review X, 11(2): 02 021045,22 2021.
[52] Axel Laborieux and Friedemann Zenke. Holomorphic equilibrium propagation computes exact gradients through finite size oscillations. Advances in neural information processing systems, 15:12950-129633 2022.
[53] Benjamin Scellier, Maxence Ernoult, Jack Kendall, and Suhas Kumar. Energy-based learning algorithms for analog computing: a comparative study. Advances in Neural Information Processing Systems, 36:527055 52731, 2023.
[54] Timothy Nest and Maxence Emoult. Towards training digitally-tied analog blocks via hybrid gradient computation. Advances in Neural Information Processing Systems, 37:83877-83914 2024.
[55] Benjamin Scellier. A fast algorithm simulate nonlinear resistive networks. arXiv preprint 17114 2024.
[56] Zhengdao Chen, Jianyu Zhang, Martin Arjovsky, and Léon Bottou. Symplectic recurrent neural networks. arXiv preprint arXiv: /909. 13334, 2019.
[57] Samuel Greydanus, Misko Dzamba, and Jason Yosinski. Hamiltonian neural networks. Advances in neural information processing systems, 32, 2019.
[58] Miles Cranmer, Sam Greydanus, Stephan Hoyer, Peter Battaglia, David Spergel, and Shirley Ho. Lagrangian neural networks. arXiv preprint arXiv:2003.0-6533 2020.
[59] Peter Toth, Danilo Jimenez Rezende, Andrew Jaegle, Sébastien Racanière, Aleksandar Botev, and Irina Higgins. Hamiltonian generative networks. arXiv preprint arXiv: 1909. /3799, 2019.
[60] Simon Duane, Anthony D Kennedy, Brian J Pendleton, and Duncan Roweth. Hybrid monte carlo. Physics letters B, 195(2):216-222, 1987.
[61] Radford M Neal et al. Meme using hamiltonian dynamics. Handbook of markov chain monte carlo, 2(11):2, 2011.
[62] Nilesh Tripuraneni, Mark Rowland, Zoubin Ghahramani, and Richard Turner. Magnetic hamiltonian monte carlo. In International Conference on Machine Learning, pages PMLR, 2017.
[63] Son Nguyen, Lizhang Chen, Bo Liu, and Qiang Liu. Memory-effiiiet optimization with factorized hamiltonian descent. arXiv preprint arXiv:2406.09955 2024.
14

---

[64] Tianqi Chen, Bing Xu, Chiyuan Zhang, and Carlos Guestrin. Training deep nets with sublinear memory cost. arXiv preprint arXiv:1604.06174 2016.
[65] Jiajun Zhou, Yifan Yang, Kai Zhen, Ziyue Liu, Yequan Zhao, Ershad Banijamali, Athanasios Mouchtaris, Ngai Wong, and Zheng Zhang. Quzo: Quantized zeroth-order fine-tuning for large language models. arXiv preprint vrxi:22522.12346 2025.
[66] Benjamin Scellier, Siddhartha Mishra, Yoshua Bengio, and Yann Ollivier. Agnostic physics- driven deep learning. arXiv preprint ar arXiv:2205./5021, 2022.
[67] S. Abreu, Boikov, M. Goldmann, T. Jonuzi, A. Lupo. Masaad, L. Nguyen, E. Picco, G. Pourcel, A. Skalli, L. Talandier, B. Vettelschoss, E. A. Vlieg, A. Argyris, P. Bienstman, Daudet, D. Domenech, I. Fischer, F. Horst, S. Massar, C. R. Mirasso, B. J. Offrein, A. Rossi, M. C. Soriano, S. Sygletos, and S. K. Turitsyn. A photonics perspective on computing with physical substrates. Reviews Physics, 12:100093, December 2024. ISSN 2405-4283. doi: 10. r1011/..evip.2222.1100 100093.
D. Brunner, Dambre,
[68] Michael Schneider, Emily Toomey, Graham Rowlands, Jeff Shainline, Paul Tschirhart, and Ken Segall. SuperMind: A survey of the potential of superconducting electronics for neuromorphic computing. Superconductor Science and Technology, 35(5):053001, March 2022. ISSN 2053-2048. doi: 10. 361-6666///4622..
[69] J. Grollier, D. Querlioz, K. Y. Camsari, K. Everschor-Sitte S. Fukami, and M. D. Stiles. Neuromorphic spintronics. Nature Electronics, 3(7):360-370) July 2020. ISSN 2520-1131. doi: 10. (0038//4192--01--06666
[70] Torbjørn Smith and Olav Egeland. Learning dissipative hamiltonian dynamics with reproducing kernel hilbert spaces and random fourier features. IFAC-PPpprsOnLiine, 115-100 2024.
[71] Marc Berneman and Daniel Hexner. Equilibrium propagation for periodic dynamics. arXiv preprint arXiv:2506.20400, 2025.
[72] Dan Hendrycks and Kevin Gimpel. Gaussian Error Linear Units (GELUs), June 2023.
[73] Yann N. Dauphin, Angela Fan, Michael Auli, and David Grangier. Language Modeling with Gated Convolutional Networks. In Proceedings of the 34th International Conference on Machine Learning, pages 933-941, PMLR, July 2017.
[74] James Bradbury, Roy Frostig, Peter Hawkins, Matthew James Johnson, Chris Leary, Dougal Maclaurin, George Necula, Adam Paszke, Jake VanderPlas, Skye Wanderman-Miin and Qiao Zhang. JAX: composable transformations of Python+NumPP programs, 2018. URL
http://githin
15

---

A
Supplementary material
Contents
A.1
Theoretical results in continuous time
17
A.1.1
Definitions assumptions
17
A.1.2 Proof of the continuous adjoint state method (ASM)
18
A 1.3
Proof of Theorem 33.]
21
1.4
Connection to Hamiltonian Echo Backpropagation (HEB)
23
Theoretical results in discrete time
26
A.2.1
Definitions & assumptions
26
A.2.2
Backpropagation Through Time (BPTT)
28
A.2.3
Proof of Theorem 3.2
30
A,244
Proof of Theorem 3.3
34
A.3
Looking ahead
36
A.3.1
Does RHEL break with dissipative dynamics? ................................
36
A.3.2
A black-box version of RHEL?
37
A.4
Models and algorithms details
39
A.4.1
Toy model
39
A.4.2
HSSM architecture
40
A.4.3
Linear HRU Block
41
A.4.4
Relation between Linear HRU and the of LinOsS
43
A.5
Nonlinear HRU
43
A.4.6
Differentiating the timediiscreizzation
44
A.4.7
Echo passes with automatic differentiation
44
A.5
Experimental details
47
A.5.I
Datasets
47
A.52
Simulation details and ressource consumption
47
A.5.3
Hyperparameters)
48
A.5.4
Gradient re-scaling for dealing with numerical instabilities
49
A.5.5
Vanishing and Exploding gradients of the Nonlinear HRU
50
16

---

A.1
Theoretical results in continuous time
Summary. In this section, we present all the results derived in continuous time. More precisely:
We formally define model (Def. A.1 and constrained optimization problem (Def. A.2)
in continuous time.
We state and prove our algorithm baseline, the continuous adjoint state method (ASM) for this constrained optimization problem (Theorem A.II) To ease the comparison of the continuous ASM with Recurrent Hamiltonian Echo Learning (RHEL) in continuous time, state re- peraameetrized version of the continuous ASM where time is indexed backwards (Corollary A.I). Finally, we also state variant of the continuous ASM when the loss function is only defined at the final timestep (Corollary A.2 ease the comparison between the continuous ASM and Hamiltonian Echo Backprop (HEB, [10]).
We introduce two technical results (Lemma A. reversal invariance property of our model (Lemma A.3) We then show that the direct consequence this property is the time-reversibily of our model upon momentum flipping (Corollary
A.2) A.2 which enable us to prove the time-
A.3), the key mechanies which fundamentally underpins algorithm. All these intermediate results allow us to finally introduce RHEL in continuous time and prove its equivalence with the continuous ASM (Theorem[A2
Lastly, we connect HEB [10 with the continuous ASM when the loss is defined only the final time step (Theorem A.4
We highlight some key differences between HEB and RHEL.
A.1.1 Definitions & assumptions
Definition A.1 (Continuous Hamiltonian model). Given O E Rd T E IRC and input sequence
u(t) E ( the continuous Hamiltonian model prediction (t) (R by definition, implicitly given as the solution of the following ODE:
TT,0
0
---T)=x; e [-T,0] (, = (1)
-I 0
We assume that:
1. H fime-reversal invariant:
I
V$ E Rdo,
vO E Rdo, Yu Rdd -|
Σ,
-I
42 + H is twice continuously differentiable,
3. + +H is differentiable,
4. Vi exists and is continuous with respect to
5. u + H|.,,u. u) continuous,
6. $ + V. and + H2.HH are Lipschitz continuous.
Remark 1. While some of these assumptions will be used explicitly in our derivations, they are all needed to guarantee the existence of partial derivatives of s as an implicit function of x and through Eq. 1 and refer (28) for such claims given these assumptions.
17

---

Definition A.2 (Continuous constrained optimization optimization problem). Given continuous
Hamiltonian model (Def. A.I. we consider the following constrained optimization problem:
min = tt & (t), 0]
s.t. vt [-T.0] (,
where we assume that:
1. is time-reversal invariant:
v$ Rdd vO E Rd0 VI [-T,0] L1, e(1 Σ. 0.,
2. + f[t, is continuous,
3. 0 + fl., 0 is differentiable
4. Volt, " - is continuous.
5. ④ + f 4; is twice difserentiable,
Remark 2. Note that while we did not assume in the main part of this manuscript that depended on 0, we assume it in the appendix for the generality of our derivations.
A.1.2 Proof of the continuous adjoint state method (ASM)
Theorem A.1 (Continuous adjoint state method ((28) ) Given assumptions 4.7 A.2 with respect to and x are given by:
the gradients of
doL = 1 go(t)dt, dzL = A(0)
with t + go E (125) ( defined as:
(t) = Vollt, (tt), 0)
JT. A(t) W [-7,0]
and X solving for the adjoint ODE:
{
x(0) 0,A(1) V ((0()))),)))))) JT A(1) V+l(t, ((),,
Proof of Theorem A.7 With slight adaptations, our proof mostly follows that of and also assume the existence of the partial derivatives of + (t, 0) as an implicit function of f, æ and we defer to (28) for the proof of this claim.
We start off defining the Lagrangian associated with the constraint optimization problem:
L( (( = ("]
0
dt (Eft, (t.,),))
where + A(t) E (Rde) (-T.0) denotes the Lagrangian multiplier associated to the constraint.
Derivation of L.. I.. For readability, emphasize the dependence of & / and 0, will leverage the existence of its partial derivatives with respect to these variables. Then, given Assump- tions A.I the total derivative of the Lagrangian with respect to 0 exists and reads:
doL [ dt (8. T)T (tllt, V 0).
-T
+ dt (t) H u(t) JT + V H (t. u(t)] JT - ((,, A(t).
-T
18

---

We can transform the last term of the integrand with 0 0) by applying Schwartz's theorem and and integration by parts as:
dtdg.: (t, x(t) = - { dtd, 0 (t, A(t) - -
1-0 0) x(1) + {"
.,A(()
(0,0 X(0) + [ (t, -.A(()
where the contribution of the first term at - vanishes because:
$ = x = 06 0) =
Plugging this back into do L yields:
doL
-gg
(0.0 X(0)
+ dtdg (1, T V 2tt, 0) (1,, JT A(1) + 0(A(t))
+ L dt (Vell, (t. 0) 0 + V H 2 v 0) J x(t))
Denoting 4, the solution of the (primal) ODE and by defining A, as the solution of the adjoint ODE:
(t).(t) = - H (t, 0) 0, u(t JT (t) - V (t, (), 0), A. 0
have:
dgL = do { tt (t,,)).
deL( An, 0)
-10 dt (Vollt, 4 (t, V (1 J xx(())
Derivation of dzL. Similarly, we can prove that:
deL = -2 (0, x) [V] -)(0)
+ dtd 2 (t, (t,
JT A(() + 0,A(t))
-T
+ A(0)
Using the same 4. and A,, get d X(0).
Reparametrization of the continuous adjoint method. To ease the comparison of the continuous adjoint state method with algorithm, slightly reparametrize the variables introduced in Theorem A.
Corollary A.1. Under the same assumptions as Theorem A.I the gradients of L. with respect to and are given by:
LgL S go(t)dt, d===
with t + go E (1Rd0)100 [O.T! defined as:
(== Vol[-,
V2 1 JT X(I) -(I We [
and X solving for the adjoint ODE:
{
X(0) 0
0,\(() = v (() JT. 2(1) + 4091-1,4(11),,
19

---

Proof of Corollary A.7 Immediately stems from Theorem A.1
Edge case: loss the final timestep only. The backward ODE can be differently parametrized when a loss is defined only at the last time step. We need this formulation for later convenience.
Corollary A.2. Under the same assumptions Theorem A.1 and assuming:
f(t, = 0 Vt E [0, T):
(T,,, == lr,
the gradients of by with respect to 0 and a are given by:
dølr( 0J VolT( + { go(t)dt, deL= - X(T)
with 1 + go E (Rda) [O.T] defined as:
go(1) = , J A(1) W [O
and À solving for the adjoint ODE:
A(O)
a,A(1)
V2 0( (1)] JT. A(11)
Proof of Corollrr A.2
Starting from the Lagrangian:
[ 0
C( x= lr (0,00),
di ( (.) 0,
(J
the proof reads in the exact same fashion as that of Theorem A.1
20

---

A.1.3 Proof of Theorem 3.1
Technical Lemmas. We first introduce two technical Lemmas which will be needed for the derivation our main result.
Lemma A.l (Block-wise Pauli matrices and associated properties). Defining Σz, En, D3 E IR RRd Xd as:
Id
Ld 0
Σx
Ey = i
Σ3 =
Idaa/2
0
I /
where denotes the imaginary unit, the following equalities hold:
1. J =
2. x = = =E² Ids
3. Eg Ey = =iE,, 5, Ey E, E E, E2 = D.
4. E, Ej = =-E E, for any i + [x,
Proof of Lemma A. easily checked.
Because of the block-wise structure of Eg, Ey, Ez, these equalities can be
Lemma A.2.
Under the assumptions of Deft A.I the following equalities hold for all 4; 0 u.
V V Σ, V H 0.
V 1 E. V 4- 0, E.,
V V e x- 0, ul Eg
Proof of Lemma A.2 The above equalities be simply obtained by differentiating through the time-reversal invariance hypothesis which is possible because the differentiability H with respect to & and 0 and using the chain rule. Namely, given some 4, 0, u.
V H4 ul
Va (H 4 0, u))
= (04
V " 0,
0, u.,.. E V " 0,u) = Σ. V H| 00
since := E. +. The other equalities derived in the same way.
Time-reversal invariance. We highlight here how the assumption : H[= $ given inside Def. entails time-reversal invariance of the dynamics up time-reversal the input sequence.
Lemma A.3 (Time-reversal invariance of the dynamics). Under the assumptions given in Def. A.1 if the solution of the ODE:
(() J V V (t)), 0 u())),
the function : Σ + (-11) ④6 (-t), - (-1)) is solution of the same ODE with the time-reversed input sequence t i1)) u -t).
Proof. Let E [-T,0], We have:
0 (t) -E, 0 ((+)
(Def. of x and chain rule)
=-,, J (g),
(by assumption)
+JJ Σ. V=/[(((t)).0.u(-5))
(Lemma A.1
JJ. -. V.. (-1),000-
(Lemma A.2
= J. V H x (t),0,((tt))
(Lemma A.J
21

---

Time reversibility. A direct consequence of the time-reversal of the dynamics under consideration (Lemmal A.3) A.3) is time reversibility: upon flipping the momentum of & at time t (TT(O) (0) + - (0))
and presenting the input sequence in reversed order, the system evolves backward to its initial state.
Corollary A.3 (Reversibility of the dynamics). Under the same assumptions as Lemma A.3
we define
T as the solution of the ODE:
Φ (0) = 44 (0), 0, (1) = V a- (1),0,
VL [O,]
Then:
vt E [0. T] (t) = (t)
Proof. & and oh satisfy: i) the same initial conditions (Ø*(0) (0) (0)), ii) the same ODE (Lemmal A.3). therefore by unicity of the solution of the ODE, they are equal at all time over the domain of definition of +
Main result. We are now ready to state and demonstrate our main result in continuous time.
Theorem A.2 (Equivalence between RHEL and the continuous ASM). Under the assumptions of Def. 7 and Def. A.2 let be the solution of the ODE t E [-T.0]:
T(-T) x. a (t) J V ( 0, u(t)
Given Φ let th be defined as the solution of the other ODE e [0.
0,4% - JV 4 (()), cJV El-1, (0,,,,,111111 (c,
-
Defining:
1
A RHEL
0, an(((())) (1,e) = Voll- x0 (1,e),, 0] + (Voll 10 (4.),0,,u))) ( 0c)
2c
A AHELL (t, = 2c Σx ($ (t, - (t, -())
we have:
Vt E TO X(t) lim A RHEL (t. c), go(t) lim Ap
RHEE (t, c)
c-0
c-
where X and ge are defined in Corollary A.1
Proof.
Defining AB
ARHEL (1) = lim_-0
4
ARHEL (1) and A
ARHBL(I) lim ARHEL e), note
that:
A RHEL (t) = Eg O. @e, (t,))l== +
A RHEL (t) voll-t, 1t 0), 0) + 0.
Derivation of A(t) lim A AHEEL et. Given t E [0, differentiating the ODE satisfied by Te with respect to at e = 0 yields:
0,00. = (,, (t, c))le=
(Schwartz Theorem)
0)) = J. V². H (1, 0)] 0 $ aeeaas((,,,, (1,e) - V a 111, 45 (11110)1 0) 0).
OJ V H H| D (1,e), 0, u(()) - -JV E(t 15 (1, e), (1,e),
By Lemma A.3
$'t,0)== = (t) Vt E [0,
therefore:
= J. V2 H' (-t) 2 (t, 0) - JV JV Elt, (-(-)), A
= J. Σ. 2 (tt, 0) E. d x (1,e))==0 - J. E. Vollt ((t)) 0]
(Lemma A.2)
22

---

Additionally, note that we have by Lemma A.I
J J. = -22
= Ex;
so that:
0, " = - -1 20 (l,
En Vollt (( () (19)
Left multiplying Eq. 119 on both sides by Σ3 yields:
0, (2,. 0,, = E
J. (Ex 0. (1 (1,)))<)) Σ Vo 4 ((()))
- ( 1. c)l==)) V (
(Lemmaa A.I
(= J (En 0. (t, ())=))= V Et 0
(J =-J)
Finally, note that because (0(0) c does not depend c, we have that:
0, (E (0, = 0,
so that all in all, ARHEL satisfies:
{
& ARHELL(O (0)
0
A (HEL(()) (t) =
ARHEL JT RHEL (t) + ((-(),
Therefore A RHEL and À (as defined in Corollary A.1 satisfy the same initial conditions and the same ODE, therefore they equal all times.
Derivation of go (t) = lim A
RHEL (t, c).
Note that by Lemma A.3 and time-reversal invariance
of L:
A ARHEL (1) V Voll- 4 (-(,),)) (VoH (1,s), ul a(-1))) le-o Voll-t, (),, a (Voll (t, O, u( u(- le=o
As the first term of A (R)EL(t) and (t coincide, the remainder of the derivation focuses on the second term of A RHEL
(t). Given & [0. 7]. we have:
(8o0((((--)),,uu))) J A(() = 2
Φe (t,0),, E. JT X(t)
(Lemma A.2
V2 a. (( 0, E. 1E, A(t)
i>,)
= VV () -E, E. A()((t)
(Lemma A.1
-- , "4..HH(0"(11000 (1,00(,)) E, A(1)
(Lemma A.1
- ol ( (t.))((()))) E 0 x (t, e)le=0
(A A RHEL
= (4.))) 0, (t, e)
(Lemma A.
0. (f u(t)) c=
which finishes to prove go (4) A RREL (1) for E [0,7].
A.1.4 Connection to Hamiltonian Echo Backpropagation (HEB)
Remark 3. The above setup and implementation of RHEL is not exactly that of Hamiltonian Echo Backprop (HEB, [10]). In particular:
the loss function in HEB is only defined at the final time step,
the interaction with does not happen simultaneously with H,
finally, we would like to recover the HEB formula giving the gradient estimate of the loss with respect the initial state of the neurons.
23

---

In the following corollary, we make slight algorithmic adjustments to match the seminal HEB implementation as much as possible while preserving the generality of the sequence modelling setting, i.e. dependence of H with 0 and u. Note that in the seminal HEB work, H does not depend on a static set of parameters nor on an input sequence. u.
Corollary A.4. Under the assumptions of Def. A. and Def. A.2 and assuming additionally:
f(6,,-11 =0 vt E (0,T), 4,, = bT,
let be the solution of the ODE. for [-T,0]:
(-1)) x, (1) J V=(l(((),0.u((()
and the solution of another ODE, for t [0,e]:
0, (t) JV er (t).
Let &e be the solution of the following ODE, for t E [c
T (0,c))=((c)) (, = JV *- H| u(-)),
Defining:
1
A RHEEL (1, () =
ul (t)] - VoH( (())))(
2c
1
A$ ARHEL (1, = E, (1,) $111, -t))
2c
we have:
1I E [0, X(1) = lim A RHEL
(t, e), go(() = lim A
RREL (1,6)
c-+0
where X and go defined in Corollary A.2 In particular:
ied CT +
(T, - (-TT)++0) (
where d = dæ denotes the total Wirtinger derivative with respect to x² and the imaginary unit.
Proof of Corollary A.4
The derivation is almost exactly similar to that of Theorem 3.1 with two key
differences:
the version of the continuous ASM against which this version of RHEL is compared is different (Corollary A.2).
the interaction of + with and H do not happen simultaneously but on disjoint intervals.
We will simply show that the interaction with and conjugation + - 4 yields the correct initial conditions and defer to the proof of Theorem 3.1 for the remainder. We will also use the same
notations and denote A RHELL (t) = lim A (t
(t, A AREEL (t) == lim A RHEEL (t. e).
A
Derivation of A(t) lim A ARHEL c). Integrating the ODE satisfied by between 0 and T yields:
(c) (0) + [ dt.J (+) 0 (0) + eJ ())
Therefore, the initial state of df can be written as:
(e) = 4 (0) + E.
(0), = + (0) e (0) 0 + O
(Lemma A.I.
By differentiating the last equality with respect to e at e = 0, we obtain:
A RHEL (0) v+L
+
24

---

Proceeding exactly as in the proof of Theorem A.2 we obtain that:
{
A RHEL
(.AREL (1) V H[ 0, u (-(] JT ARHEL(U) (1) (0)
-v)
&
Therefore A ARHEL and À (as defined in Corollary A.2 satisfy the same initial conditions and the same ODE, therefore they are equal at all times.
Derivation of go lime-+0 ARHEL (t, See proof of Theorem A.2
Connection to HEB formula. In particular, we have:
d (0)), 0 X(T) Σg 0. 0e
(
Dn (T, (T, 0)) + (())
En (T, (-T))) + +(e)
(Lemma A.3
Ey E. (T ("(-T) + +(())
(Lemma
Dy (Φ (T,)" ((-T)) O(e)
Left multiplying on both sides by i and noticing that id = -iE do. finally obtain:
-ied by (0), 0 $- (T.c)* D (022)
25

---

A.2 Theoretical results in discrete time
Summary. In this section, we introduce all the results derived in discrete time. More precisely:
We first formally define Hamiltonian Recurrent Units (HRUs, Definition A.3 HRUs can be regarded the discrete-time counterpart of the continuous model introduced in the previous section (Definition A.I). namely as an explicit and symplectic integrator of the continuous Hamiltonian model which preserves the time-reversal invariance and time- reversibility properties in discrete time. We also introduce the constrained optimization problem naturally associated with HRUs (Definition/A. AA.4). which is the discrete time counter- part of the constrained continuous optimization problem introduced in the previous section (Definition A.2).
We then formally define Hamiltonian State Space Models (HSSMs) as stacks of HRUs (Definition A.5) and the multilevel constrained optimization problem which is naturally associated to these models (Definition A.66) A.6p.
We state and prove our algorithmic baseline, Backpropagation Through Time (BPTT). through the lens the Lagrangian formalism establish clear connection with the continuous ASM. We first introduce and derive BPTT in its general form (Theorem A.3 and then apply more specifically
HRU as defined in Definition A.3 (Corollary A.5
As we did in continuous time, we introduce series of technical Lemmas needed to extend RHEL in discrete time. We first demonstrate the reversiiiiy of HRUs on a single time step (Lemma AA.4) which then enables us extend the time reversibility property derived in continuous time (Corollary (A.3) to discrete time (Corollary A.6) After introducing one last technical result (Lemmal A.5), S), we then state and prove RHEL in discrete time when applied to HRUs (Corollary A.7) As the algorithm prescribed by Corollary A.7 A includes solving an implicit equation, we finally introduce a slight practical (i.e. fully explicit) variant of RHEL in discrete time (Corollary A.8).
Lastly, we show how to estimate gradients end-to-end in HSSMs by using RHEL-chaining (Theorem A.4). We also highlight that in practice, when using feedforward transformations across HRUs, the algorithm prescribed by Theorem A.4 iiiiiiiiiy requires to chain RHEL. through HRUs and automatic differentiation through these feedforward transformations (Remark 8). This remark fundamentally underpins the actual algorithmic implementation of RHEL which was used throughout experiments.
A.2.1 Definitions & assumptions
Definition A.3 (Hamiltonian Recurrent Unit). Given O E Rdo K Nº and an input sequence the Hamiltonian Recurrent Unit (HRU) prediction is given by:
(us)REl-K.of E (Rdu) K
D MH. up] Vk -K ...-- 1,
with H = T + V and:
MH.5 == MT o Mv.8 o MT.8/2
VoT, Mv.6 = 1 +J V
We assume that:
1. H separable, i.e. V and T only depend on & and π respectively:
v 0, u)
V( 0. ni),
T 0, u) T( 0,
2. T and V are time-reversal invariant:
T T T 4, u), {
V d
40 E do Yu R
V 0 = V $...bb. 0,
3. T and V are twice differentiable with respect to C 0 and u.
Remark 4. Note that M is simply Leapfrog integrator associated with H. We justify each of design choices below:
26

---

steps-parametriaation We write the Leapfrog integrator in three-steps fashion to yield reversible integrator.
Separability of the Hamiltonian. In the case where & and Φ can be separated out in the Hamiltonian, function, the Leapfrog integrator becomes explicit (33].
'I' and V as functions of b. Although T and V only depend on π and Φ respectively, we choose to define them functions of so that the proof of RHEL the continuous case seamlessly translates to the discrete case.
Definition A.4 (Constrained optimization optimization problem in discrete time). Given a continuous Hamiltonian model (Def. A.1, we consider the following constrained optimization problem:
min L
S.T.
Dx+1 MM..s 8, usl
vk
K -..-
kmm K+1
where assume that:
1. is time-reversal invariant:
V④ E Rdo 40 E Rdo Yk ..,...
4 0 = -x (
2. is twice differentiable with respect to and
Definition A.5 (Hamiltonian State Space Models), Given (0(1)... ,(N)) E (I20) K E N and
an input sequence (u=)) k.0) E (Teua K Hamiltonian State Space Model (HSSM) is defined as the composition of HRUS defined in Def.
A.3 as:
U. VL [0, 11,
Vk [-K,0]
4
(E+1)
M
(E+1)
109 (
k+1 410,6
in a vectorized fashion as:
T
u, AA 10 N
1]: T (++1)
(
(++1)
0(0)
m
Definition A.6 (Multilevel optimization problem in discrete time). Given a HSSM (Def. consider the following constrained optimization problem:
A.3 we
min Σ
k
T (0) U
k= K+1
vl E [0.N - 1I: T
(1+1) ((+1)
M
(8)
(++1)
0(0) T
H(0)
where we assume that e satisfies the same assumptions as in Def. A.4
27

---

A.2.2
Backpropagation Through Time (BPTT)
General form. We first state and prove Backpropagation Through Time (BPTT) for any integrator MH
Theorem A.3 (Backpropagation Through Time (BPTT)). Given assumptions in Def. A.4 the gradients of the loss with respect to the parameters and the inputs u are given by:
K Σ (k()), k=0
døL
du- gu(() Ak E [0. K - 11.
with:
go(k) Iu (k)
(k+1), -, -k] 82 (((tt), 0, (++1)) d
(..MHol(*-(eee. 0, u- -((++1) Ap.
and where (Ax) satisfy the following recursion relationship:
A Vi 0,
A++1 MH, (1+1), 0, u- u-(+t))) X + Vi (k+1)) - (A: + 1)] vk 0,
K
Proof of Theorem A.3 BPTT can classically be derived through the application of the "chain rule" backward through the inference computational graph defined in Def. A..] Another useful viewpoint though, which directly connects BPTT the discrete counterpart of the continuous ASM and will be useful later in the appendix, is to derive it through the method of Lagrangian multipliers. Namely, the Lagrangian associated the constrained optimization problem Def. reads as:
C( A, = Σ p 4 h O, -b]++]+ \o (MHH [ ((-((+1)) 0, u (K+1)) x
Extremizing & with respect to + and À yield xx and Ap.
Vk: 0, ... K
0x L( X. O u) = - (k+ (k+1) 0 u-
1 -h+ =0,
04. C( pr A., 0, u. = V 20 0,0 A0,- =0
Vk
K - 1: 0 A., O, u = VIEi O, -k]
O,
A d =0
Finally, the total derivative of L with respect to 1 reads as: dg L do 4. X., 0, u) 8g DoL( A., 0, u) + T. 0. A., A. 1, +00A7 -2L((F,, X, U,
2 -k] - 02 (k+1)) O, u (1+1)) A
The total derivative L with respect to u- is derived in the exact same fashion.
Remark 5. Note that using the vectorized notations introduced in subsection 3.3 the constrained optimization problem defined in Def. A.4 re-writes:
the Lagrangian of
L L= 1 0 Tr [MH,s 0.
3)
x
with Tr denoting the trace matrix operator and:
f|Do, 0,0]
E R Kxl
0 :
E kxl
=
[ (K 1), O.-(K- (K-(K-1) -1]]
MH.a Po. 0,u)
MH MH. 0. ū
E IR Kxdx
MH. $ -(K-1)); 0,
28

---

Detailed BPTT. For the needs of our derivation of RHEL in discrete time, we now introduce finer-grained version BPTT given model assumptions given Def. A..3
Corollary A.5 (Detailed BPTT), Given assumptions in Def. A.3 A.4 the gradients of the loss with respect the parameters and the inputs u-k are given by:
K-1
dgL
Σ [k),,,
d (k),
km0
with:
go(k) = V2 k. 0, kk] Vi.277 - 1/--1 u- ( J Ak
8
+oV 1,2 V (x+2233))) 0, u (k+1)) J Ak+1/3
0, a (K+11))) J AR+2/33
2
(k:) = v1.3 T (k+1/3), O, u- (k+11) J A
2
8
+ + 1 (k+2/3) 0, u- (k+)) Akt1/3 + 2
2 - 0, u- u-(k1))]
JT. Ak+2/33
and where (Ax) satisfy the following recursion relationship. with do= = Vol(Pol and Ak E [0, K - 1|:
Ax+1/3 >a + 7 (K+1/33); 0, (k+1)) J A
A4+2/3
V + (3++2/3). 0, u (k+1)] J
Ak+1/3
AR+1 A
v
T + (k+1)) O, u (k+11) J
- - 0)
Proof of Corollary A.5
Direct application of Theorem A.3 with the inference computational graph
details defined inside Def. A.3
29

---

A.2.3 Proof of Theorem 3.2
Time reversibility in discrete time. We first derive the discrete counterpart of Lemma technical pre-requisite for the extension of RHEL the discrete-time setting.
A.3 as
Lemma A.4 (Reversibility of MH,s). Under the assumptions of Def. A.3
A.4
v E Re ve E IR,d Vu E Rdu 4k+1 = MH, 4x, 0 2 = MH o Pi
0 u = 4
Proof of Lemma A.4 Let Pe. Pk+I E Rd be such that:
Pk+1 -HH.
which rewrites, given Def. A.3
@k+1/3 +x+ + V 2.
1+1133 Tkk
Øx+2/3
0k+1/3
MH
(20)
Tx42/3
Tk+1/3
oV (
Φt+1
@x+2/33
T 0, + 6 V T 2++33
Tk+2/3
becomes apparent from Eq. 20) that MI.8 is invertible with respect to its first argument and that inverting MH,5 amounts change to d:
0b+11 (V T 0,
PR42/3
Tk+2/3
Tk+1
M
@x+1/3
@k+2/3
(21)
T++1/3
Tk+2/3 + oV Φ V x 0.
Φ
06+1/3 $ " T T 0,
T2+1/3
and therefore:
1 1x41. 0, u] = MH.- ( 0, u). Φ -M
Denoting π = -t, note that by time-reversal invariance hypothesis in Def. A.3 and Lemma A.2 T "T' 0, u.. Therefore, Eq. (2T) rewrites as:
we have that V VT 0, u]
v
Øk+2/3
Φx+1 + 8 V
T 2
T 0, u]
244/3
x
Φk+1/3
@x+2/3
T
VVVV(oc+2/3 u.u -8
(22)
Φk
-b411/3
+ N.
where equations bearing π have been multiplied on both sides by -1. Finally note that Eq. simply rewrites as:
22
MH.al9i. 0, u] x
Corollary A.6. Under the assumptions of Def. A.3 if T satisfies the following recursive equations:
$- = x, -K
Pk+1 - P O, un]
and Φh is subsequently defined as:
4 Vt 0,... ......K K - 1:
& k+1 MH. 1 0, 0(4)))))))
Then:
== 0, K $ Φ
Proof of Corollary A.6
This result immediately obtained by iterating Lemma A.4 over the whole
trajectory.
D
30

---

A technical pre-requisite. Finally, we need one last technical Lemma to handle subtleties pertaining Jacobian evaluation which only occur in discrete time.
Lemma A.5. Under the ussumptions of Def. A.3 if we have, for some O E IRdo and u E Rdd:
4 x++1
MH Pk. 0, ul,
then:
T T k+1/3: 0. + 0, u/)
V v 22/33 0.ul =V 0 u] T T 0 x+ 0 u
Proof. This can be seen by simply writing M explicitly for & and n:
Φx+1/3
Φb
+ 2 vn TTTA, 0, u.
Tk+1/3
TTK
MH
ΦR+2/3
1/+1133
Tk+2/3
*K++1/3 8V (
Φ1+1
-0x+2/3" + 2 V T 0,u)
TT++1
=
Discrete-time RHEL. We are now ready to state the main result of this section.
Corollary A.7. Under the assumptions of Def. A.3 A.4 ler (#x); satisfy the recursive equation:
① Φ-
Dk+1 MM
D O, un vk ----K -..... -1,
and let th satisfy:
Pf(E) = 1 + E. VollPo. 0,0).
Vk 0.... K 1:
$ ( = -M..//2 $ (c). 0 u (I+1)]
Φ ( = =My; $ 1/3 (t), A u (k-1)
d (e) = Mr.5/2 2 42/33 (c), 0, u (b+1)] - eJ 1
0, - + 1)],
Then defining:
1122 (1)2 (c), 0, u (1+ =
: f1/l//(c)) 0. u u H ( u(u(())))
2
A RHEL (k [2( 0, -k|
uu (x (vgH 1/ 1 (c). 0 u ( V u ((++1))))))))
1/2
O, u -2
2c
A RHEL
(k, = 2c
A (kG))
RHEL
Σ, ( 2x(-c))
2c
we have:
4k 0.... K 1: A lim A (k, c), (o(k) lim A RHEL (k. c),
gu (k)) = lim A c+0
RHEL (k, c),
where (Au) (90(k))) and (gn, (k)) x are defined inside Corollary A.5
31

---

Proof of Corollary A.7
Let & E [0, K 1]. We define:
A $
ARHEL (k) == lim A +
RHELL
(k; c) = Σg 0. ((: c)le=o
c-+
A ABHEL(K) == lim A
AHHEL. (k, E) (2( 0. -k] -80, (22111/2 O, a uu(+t)))) le-o c-0
ARHEL(K) (k) := lim ARHEET (k, -60. (33l 1/2 ( (c) 0, u..uaa++)))))))))))))))))))))u le=0 c-0
Derivation of A lim,- Ap ARHEL (k: e). We proceed exactly in Theorem A.2 with some subtle
adaptations which we highlight. Differentiating the dynamics of d' between and A + 1/3 and
proceeding as in the proof of Theorem A.2 using Lemma A.1 Lemma A.2 and Corollary A.6 (as the discrete counterpart of Lemma A.3 which was used for Theorem
we obtain:
8
ARHE +
(l: 1/3)
A (k) V k. 0 (K+1)))) J A ARHE (k)
However note that this does not correctly match the dynamics satisfied by inside Corollary A.5 between & and k + 1/3: the Hessian V should instead be evaluated at -Kk+1//3) Fortunately,
using Lemma A.5
T k. ((+1)) - -k++1/3)
V 0, u((+++1)]
Therefore we get:
8
Ab
RHELL
(k- - 1/3)
RHEEL (k)
T (k+1/3), 0, (0--(+11)--- JJ A. (HEE((K)
2
Proceeding the same way T
k+2/3 and 1+11 get altogether:
A RHEL (k: 1/3) A ARHEL(N) - O, 0,u-(8+1) JT - ARHEL(K)
$
A RHEL (k 2/3) A AHEEL (k (1+1
+8V V 2k+2/33) (K+ JT A RHEL 1/3)
RHEEL
(k +)) A ARHEL (k 2/3) + V (k+1). 0, u (k+1)) J (
+ 0, - -(k + 1)]
A RHEL +
satisfying the same equations as (A&)a given by BPTT, together with same initial conditions:
RHEL
(0) V Po. 0.
yields the desired equality.
Derivation of go (k) lim, AB ARHEL (ks,). Proceeding the same way as in the derivation of
Theorem A.2 starting from the expression of go (k) derived Corollary A.5 using A (HHEL((K) X
and paying attention to evaluating Jacobian at the right places using Lemma A.5 we obtain Vk:
0, ... K 1:
go (k) = V2 (0), 0, -k] 8
Vo (e),, 0, ux(1))) 2V a (((1))))) + V 0, u (-a+tnl) ==
2
There again, using LemmalA.5
VoT( O, u [k+1)] = VoT (1//())) u a-tpl
V 0, u (u(()1)).
therefore:
go (k) = v2( 0, kk] 0. {
5
T 0, (.
VoV क k+ u.
VoV/ u (k+ VoT (2/ (uu
-2)(0 -k| - 0.
8
0, u. + 110244/2(()) 0, (1422)))
-V)(00(0) 50. (H 1/ (H12((0(())) (e), 0, (++1)))
lim RHEL
AHEL (k,e)
==
c
32

---

Derivation of 9u (k) lime- Au (ke (k.c). c). Strietly identical to the above paragraph.
ARHEL
Remark 6. Note that the echo dynamics prescribed by Corollary A.7 are implicit:
(i(c) 20 + eEg
o.
4k = 0, K 1:
$ 1+133(() (t) Mr,5/2
$ (c),) 0 u u
Φ 1+2/3(()) = My.s
My.s P (4112(()) 0, uuuctt (k-1)) -
(c) = -T & k42/a (c) 0. (o--((+))l=-- cJ V 1,, 0,- + +]],
(c).
where & 1+ appears, highlighted in red, on both sides of the last step. A straightforward way to make the echo dynamics explicit while preserving the theoretical guarantees of Corollary A.7 to linearize the nulging signal, namely using instead the following set of equations:
D((() D6 + Ex Vol 0,0
Vk 0, 0,... K -1
$ (+1/3(()) = -,r,5/2 D (E), O, u u((+++)
$ (+2/3(() = -M,s5 o
1/3 (c), 0 u (a-(+1))
$ (+](()= = =MT,5/) c k+2/a (c), 0, (. (. V Vo.((0-(441)) O + 1)],
Note that the of the original implicit equation becomes
- in its linearized counterpart.
This remark leads us to a slight variant of Corollary
A.7
Corollary A.8. Under the assumptions of Def. A.
A.4
let (Pk); satisfy the recursive equation:
$ Φ x,
1x+1 Mu
D 0.
un vk
-.....-
and let th satisfy:
( = // + c.. yu.
Vk = 0, ... K 1
Φ k+1 1/3 (c) = MT,8/2 (c), 0, u (k+1)]
242/ (c) = .M..s R+ 1/3 (c), 0 u u((())))
(E) = Mr.8/2 (+2/3((), (c), 0, u- + Ex y-(k+1),
where y E Rxxx does not depend on tb Then the same conclusions as Corollary A.7 hold, with (Ne))E(o,K<1] satisfying Ao Yo and vk [0, K 1]:
ARt1/3 = = + V + (k+1/3) O, C-uukkk (k+1)] JJ. -A
Ax12/3
- + VV (k+22/3), O u (k+1)] J >x+11/3
Akt1
= -A=-2/3 + V 1k+ u (k+11) JT Ax42/33 + y (k+1)
Proof Corollary A.8 A.8
Identical to the proof of Corollary A.7
33

---

A.2.4 Proof of Theorem 3.3
Theorem A.4. Assuming HSSM (Def. A.5 and the optimization problem depicted in Def. A.6 we
have:
VL = 0, 0,.0.,,000.. 1: dgcn L = lim (e)
c-
where (c) can be recursively computed backwards, starting from the top-most block as:
$ To,1N(
(e)
M 11(-1), 12, o
(M) (c), O(N-),
V VoorL
A0(N-1(())
V2l f k:
0(1--1) -ki
K-1
Σ
V2H /2
c. (L) (c). M( $ (N-111) V2H 1/2 1/ c. -c),0(N-))
- (-c))
2c km
#
8
V 22H c,( (c), 11)
V Val 2/ (L) (-e), O(N-1) B(N-D1))
Aa (e)
2e
and subsequently for upstream blocks, i.e. A4 = N -2. O:
T Bel(et1)( (e)
M H
(2++1)
(0) (6) + EE - (c)
K-1
(E( = 2E
6
Φ (c), 000) x (1) 22H2/21 (f(++)) (-E), 000) x (
k=0
Agco (c) = 2c
V /2
VaH / (t+1)
V
(-e), ( T101))
Proof. Re-using the vectorized notations introduced in subsection 3.3 and used in Remark 5 the Lagrangian associated to the optimization problem depicted Def. A.6 reads:
[Magond (1+1)) ,00) u a0)) - ml0+1) (2 7
C T (6())
OlL-11+
+
Tr
0=
K- (L))
Σ 0=
L)) -1] +
Σ
(Attti))
((+1) (Mero,, ( (k+1)): 0 x (2)
-k+1)) - $ (((1)
k
t=
We proceed by induction on the block index starting from f= L.
(N) and X (N) for k E [0. K - 1] be the critical points of L. By Initialization (E = L). Let (
Theorem A.3
K-1
dI =
Vzl
M O(N--), -k] + 8,M
(N) -k
(N) (N)
1k+ oll-1), x
x
(N)
-(++1)
k=
with (2 (N)
[e satisfying the following recursion relationship:
X (N) = v (N)
0 ((N-1)..0) 0].
vk: 0..... K
(N) 82
c (N) 1k+1) 1IN-1), $ (1-1) (k++1)
(N)
+ v
c (N)
1(k+1) 0, - (k + 1)
k+1
Given the definition of the dynamics of (N,(NN) by hypothesis, we can directly apply Corollary A.7 to obtain:
K-1 Σ
dgiN 1) T.
Vzl T
(N) k
((N-1),-K)
-k] 0.0. (22H(N-1)312 (NE (c), ((1-11) क (N--1)
-K--1) c==
k=0
d
(N-1) 0g.
s 1 (N)
1k+1)) 0(N-1) Φ (N-11) IT (M
(--
=-60. (V ( (N) (c) $ (N-1)
9(
(k+1) ==
34

---

Induction (e+1 - E). Let us assume that the desired property is satisfied at layer & + 1. We denote again D and 1 for l E [0, K 1I the critical point of L. We have, for k: E [0, K 1]:
A 03M ,8 20
(e+1)
0 20
Ap
(E+1)
vk = 0.... K A (1) M (1-11)6
k+1
$ (k+1))
Oie-1),
x (e-1)
-(++1))
1T
1/ M
p (e+1) 000, x (h+11)
-(++1))
1 ((t11)
Using the induction hypothesis layer (e (e+1): 1):
1 (6+1) (k+11) 009) & (k+11)
(e+1) -8 0. (V311 x
$
( (E), (((ttt)))) ==
lim A (2) (c) c-+
(++1)
Therefore on the one hand, denoting Agw lim A A (c) RKxds the dynamics A rewrite:
2 = (0
vk 0, ... K 1 A 1 = M.M (k+1))
l(--1) 1), Φ (8-1)
2 -1 +
(k+13)
On the other hand, the dynamics (e.(()) read by hypothesis:
(0) +eE A (e((e)
0
Vk: 0. K-1
(.(
M
1 0(e-1) (t-1)
(k-1) + + A a
k+1
(K-1)
using Corollary A.8 with y A p we conclude that:
L -iim (e
c+0
Remark 7, Theorem A.4 and more generally our definition of HSSMs (Def. A.5) assume that the connectivity patlern of the HRU units is linear chain. Note that while we chose this hypothesis for the sake of clarity of our results and their derivations, Theorem could be seamlessly extended to any Directed Acyelic Graph (DAG) of HRUs. This allows, for instance as simple and realistic case, to use skip connections across HRUs within HSSMs.
Remark 8. Note that RHEL chaining as prescribed by Theorem A.4 i.4iiiiic chains RHEL and automatic differentiation. Indeed, if H explicitly parametrizes feramerrizesseeddorvvar mappings across HRUs as:
(()) To
(())
((0) (1)
-(k+1) = H()) T0. (.(et1) (c), F
② (k+1)) 000] (23)
then, denoting u F
$ (2)
000 we have:
-k++1)*
(e.(e+1) a. V33 (1,1111
(c), (t) (()
(((tt)))
-(K+1) ===
F
(+)
000
(
a.
(1,1(0),122
(e(())) (e), ((0) u4 un))
1(++1)>
==
0 F
-(++1))
(VHH(081122
T ()
(4), 0 u²]
12H001112 (e+(()))) (-e), (00) ul)
2
The red part is done by automatic differentiation and the blue part by RHEL. This underpins the implementation of RHEL chaining we used in our own code.
35

---

A.3 Looking ahead
A.3.1 Does RHEL break with dissipative dynamics?
The applicability of RHEL is fundamentally limited by its restriction to conversative systems. Whenever dissipation is introduced, the system is no longer time-reversal invariant so that RHEL does not readily apply. One interesting question though is whether RHEL could extended to the case where the time-reversed trajectory is physically feasible. Namely: the energy which is lost during the forward trajectory could be exactly pumped back into the system during the echo trajectory. Although the system may no longer be time-reversal invariant, we can ask ourselves whether perturbations of the time-reversed trajectory carry relevant error signals.
Foreword.
Dissipative Hamiltonian model. dynamics as
For this purpose, we introduce dissipation inside Hamiltonian
77
+ (--T) - x. At E [-T,0] (t = J R) V 0, u(t)),
with R symmetric and positive definite. We set R = I for simplicity the loss of generality of this analysis not an issue here as we state negative result below.
Time-reversed dynamics.
With this model, the corresponding time-reversed dynamics read:
0 = =- (-t(
=JV a. H (--)) + E. E. V (() 6t,,0,
= =JV a H (1 u( (-t) + V H "1 (t), u u-t))]
(24)
Therefore as expected: i) time-reversal invariance no long holds, ii) the sign of the dissipation in the forward trajectory switched in the time-reversed trajectory.
ASM with dissipative dynamics. For this model, the ASM gradient estimate at time reads
Vt j0, T): go (4) ==
- ( JT A(1) V (
u( 0, x(L)
with À solving for the adjoint ODE:
A(0) = 0
(1A(1)
V (-))))) 0, u(-)) JT. A(1) - ([))),))) A(1) + V [ ((-)), JT
((()),
(25)
We highlight in blue the extra term in the adjoint dynamics which is induced by the dissipation.
RHEL procedure with dissipative dynamics. Following the intuition conveyed earlier, define the echo dynamics as perturbations the time-reversed trajectory (Eq. (24
(0)
(0)
vt E
(t,c)
JV H (t, 0, u(-( V H (t,c),
(t -t]
Proceeding exactly as we did for the derivation of Theorem A.2 we differentiate the above echo dynamics with respect to c 0, (t, En V2 ( 0)] (J Ex) 0, 45 (t,e))==o) (,,))),,)()))))))))))))),)))))",
O, yielding:
+ E. V2 E, a. (t,
ΣS V Φ (-1)),
Multiplying by Σn on both sides and reusing the notation AHEE (1) E,, 0. we
have that:
A RHELL
0 RHEL (t) v ([-100)) JT. A (RHE v H J ARHEL V V ((-1),
(26)
where we have highlighted in red the extra term induced by dissipation.
Comparing Eq. 25 and Eq. 26 we see that the red term inside the perturbed echo dynamics and the blue term inside the ASM dynamics differ. Therefore, cannot be concluded that À and A RHEL equal at all times, as it is the case for time-reversal invariant systems.
36

---

Conclusion. Using our methodology based on the equivalence with the continuous ASM the echo dynamics cannot be seamlessly adapted to the case where the system subject to dissipation. even when "pumping energy" back into the system". However, this does not mean that developing forward-only, perturbation-basee when systems exhibit dissipation. In study released shortly after our submission [7T), it was shown that a class of linear dissipative systems could be described with a Lagrangian formalism (which is closely related to our Hamiltonian formalism) and gradients estimated as finite differences of two dissipative trajectories, provided that the system is subject to periodic boundary conditions.
methods for temporal credit assignment impossible in general
A.3.2 A black-box version of RHEL?
Foreword. RHEL requires knowledge of the analytical expression of the Hamiltonian to implement its learning rule (Alg. []. Based on the Agnostic Equilibrium Propagation [66] algorithm. we briefly sketch the most straightforward application of weight-centrie homeostatic control to RHEL which, under some hypothesis on the model definition, yields an algorithm which may not require exact knowledge of the Hamiltonian. We caution the reader that this proposal is still brittle and would require much more thinking about the right modeling choices, the theoretical guarantees of the resulting algorithm and experimental validation thereof, which goes beyond the scope of this paper.
Model description.
We propose joint dynamics of the neurons and weights as: (-1T) = x, 6- ()) = J
V 2 (t), (t), u(t),
[
(,) J Ve
H | (t), ( u(())
where e denotes a control variable of the same dimensionality as - namely: one control variable per weight- and @ :- (0 ) denote the concatenation of the position and momentum vector of the parameters. Note that this modelling choice similar to that the seminal HEB work where both the neurons and weights are dynamical variables [TO].
candidate procedure. Based on the above model proposal, a plausible procedure could read as follow:
First phase. Throughout the first forward trajectory, we adjust the control variable c dynamically such that 0 remains stationary at its initial value ((0) = Op:
VI E [-T,0] a (() Jo H (()),,(()
(10(() Jo V (1), ((),((()
(0000)++1e())) [lll((() (()())-0(11)) 2
Vt E [-T,0] c(t) = G, + V ( (t), (t), u(t)],
Vn, H (()) u(t) [[(])]]] 0
We denote (co(t) T.0) the resulting control trajectory.
Second phase. During the echo dynamics, we apply the control co in reverse:
(+)] vJ (t),, Vt E [0, 0. (t) Js V H H (t), O (t), ul
a (t) Jo Ve
n (t), (t), ul (()))] -+)]+ + (lco((-t) 0 -
The echo dynamics of the positions and momenta of θ read:
{ 0,0%(1 Vx; H( (t). @1)
u(-()) (t) (t) = - V.H (t), (t), -( (t))) co(tt)
(27)
Note that because we apply the control in reverse and the parameter dynamics are Hamil- tonian, the parameter dynamics are also time-reversal invariant. This means that if 0, executing Eq. (27) will maintain @ at 6g- Now assuming that > O, note that re-using the expression of cu obtained above, the echo dynamics of Ty wr-tt::e:
01T Voll (t), (t), ul( ((-t)]- Voll (oHl( Oo. u(-()l) + 0. (*())
1
-C (t,c) V 2
" - 0.
(28)
37

---

where AgRnEi = (VoH (t) - voH( O. ue a(-t)))) Note that if @ was "close enough" to Oo then we would have A et, & A (t, (t, e). Additionally, note that playing the control Co in reverse explicitly biases 0r towards 0 through an extra elastic contribution which appears in Eq. (28).
Third phase. At end of the previous procedure, Eq. (28) highlights that could impart momentum kick proportional to the gradient of cost function, namely:
L." T MHALLL (t.
(T) [x))==- - -
dt
but what we really want to do to update the parameters This pattern also appears in the original HEB algorithm. As done in [10), would need a decay step which converts the momentum shift into a position shift in order to yield:
(T(T) ~ (-T) - E S/" d (t,e) R (-()) do(-rT
38

---

A.4 Models and algorithms details
Summary. In this section, we provide details about our models and algorithms. More precisely:
first describe toy model used inside Fig. 2 terms of its Hamilto- nian, resulting continuous-ting dynamics and RHEL gradient estimators as prescribed by
In section 1.4.1
Theorem 3.1
In section A.4.2 provide details about the HSSMs which we used in our experiments. We describe in greater details HSSMs made up of linear HRU blocks (section A.4.3) We describe their Hamiltonian, their resulting dynamics, the associated gradient estimators prescribed by RHEL and explain how parallel scan can be used on these models, especially when applying RHEL.. Similarly, we describe HSSMs made up of nonlinear HRU blocks (section A.4.5).
In section A.4.6, we show how the time discretization 8 can itself be trained by absorbing it into the definition of the Hamiltonian function.
Finally, in the light of Remark 8 highlight in section A.4.7 how our implementation hybridizes Automatic Differentiation (AD) and RHEL using Algorithms
A.4.1
Toy model
In this section, we provide the gradient estimators for the parameters of the toy model. The toy model is simple network of six mechanically coupled oscillators. Each oscillator is described by state " "o', x) where d' E IR is the position of the oscillator and x² E R is its momentum. It has a mass parameter m, and spring parameter k; Any pair of oscillators (i, f) is coupled via the spring parameter kij. The input to the models is time-varying external force utt) E TR coupled to oscillator 1. During the echo passes, the nudging force is modelled by spring coupling with parameter E IR
n² =
to an external force y(t)= E R. The Hamiltonian of the system is given by:
H u
Σ
(T+)² 2m4
1
E (0')
ΣΣ jzi
kigld x²))² uu
(29)
Which gives the following equations of motion:
O' x' m
for all i E (1,6)
# k + Ejjys kig - -²), iE (2,
#1 - kI + Exsei kis - -0))== + u
# kgΦz + Ensea Kay (d x+) 6.(( - y)
where d, is the indicator function of the echo pass, it's equal to during the echo pass and otherwise.
RHEL gradient estimators of the model parameters. For the mass m., we have:
RHEL
(V H ( (1, H, (1, -e), 0, ul)
mi:
26
( (t, c)) (a'(t, (t))
(30)
2e
2m
2m
For the spring parameters ki, we haves:
1
A
RHEL
(V
V 11, 0, 0, ul)
2e
1
(1))) - (0 (1, -())
(31)
2c
(32)
39

---

For the coupling parameters kij, we have:
A RHEL
(Vky (1, VR., (1 -t), 0, ul)
ki
26
1
- -6) - (t, -e))
(33)
2e
RHEL gradient estimators of the state sensitivities
For the state position sensitivities, we have:
1
RHEL
Ex (Ø (t, 4 (t, -e)
2E
1
(I) 7' -t))
(34)
26
For the state momentum sensitivities, we have:
1
A RBEL
Eg ( (t, (t, -c)
m²
2
1 (t) (ttt, -()
(35)
2c
A.4.2 HSSM architecture
In this section, we outline the detailed architecture of full multi-layer HSSM architecture. For the inference, we re-use the same stacking architecture of recurrent blocks and feedforward elements as the LinOSS model 29 The model starts by encoding an input sequence IRd.XK via an affine transformation. The transformed sequence then progresses through multiple HSSM blocks, linear (see Appendix A..33). or nonlinear (see Appendix A..45), directly followed by nonlinear transformations. These transtormations include the Gaussian eтror linear unit (GELU) [72) and the Gated Linear Unit (GLU) 773 defined as GLU(x) sigmoid( W wxx)= o W2x where W1.2 represent trainable weight matrices, accompanied by residual connection. The sequence output from the final block undergoes a second affine transformation produce the model output.
The full linear and nonlinear HSSM is further presented Algorithm 2 For clarity, when applying operations to sequence elements denoted with an overline (e.g., W) these operations are implicitly broadcast across the time dimension. Specifically, for any function f applied to 7, we have f(ū):
f(u)) for all time steps / (1, ..., K).
For the inference of the linear HSSMs, we keep the same recurrent block as LinOSS with slight
change in the integrator (see A.4.3). For the nonlinear HSSM, replace the recurrent block by
UniCORRN recurrent block 30 for which we use the same integrator as for the linear HSSM (see
A.4.5
Algorithm 2 HSSM model
1: Input: Input sequence model type type (linear, nonlinear)
2: Output: HSSM output sequence 7
3: (7)) + W Vennu π + bene
4: for f= = l,..., N do
5: if type linear then
6:
T (0)
+ (INEARHRU(((-1I) ,0)
Via parallel scan
7:
else
8:
- + N
> Sequentially
9:
end
if
+- (T( D
10:
T
11:
I (6) + GELU(Z(()
12:
π - GLU( (1) ग (0-1)
end for
14:
0 - W (.) (N) + bdec
40

---

A.4.3 Linear HRU Block
Hamiltonian of the recurrence. The linear HRU block is the composition of a nonlinear spatial transformation and a linear recurrent transformation (see Eq. [7). Here we provide more details about the linear recurrence that is computed with the RHEL gradient estimator. The linear recurrence defined by the following Hamiltonian:
H 0, uj
T(T, 0, ul VlØ,0,
| |
200 AΦ 0 Bu
2
(36)
Dynamics.
The dynamics of the linear HRU block are defined by the following equations:
{ Φ
(37)
π
AΦ + Bu
Which, after -ineeesscetiiationn with the integrator defined in A.2.1
gives the following equations:
8 V
0k+1/3 = Φb +
[TR, O, un]
2
8
= +b +
2 T
Th+1/3 = -
0642/3 = 0x+1/3
Tk42/3 = Th-1/33 + 8V Φ 1/ 0, unl
(38)
=Tk-1/3 - 5A 8Bus
Φ++1 2242/3 + VTT[TTx222 un
2
6
0442/3
2 kk+//3
1&+1
2/3
with the initial condition k = T क
xIx)= x.
For the echo passes, the initial condition $ follow below, with modifications from equation
$ + A.l 0). The dynamics equations
38 highlighted in blue::
Φ $ + 2T
T41133 = T
22 = 02+11/3
(39)
T442/3 Tx+1/3
ò A + +oBu- (k+1)
8
Φ1 = 21+2/3 + т
k4233
te 7 { (k-1)
2
Tf = T 244 tc V c $ (k+1))
41

---

RHEL gradient estimators.
The gradient estimators of the parameters of the linear HRU are:
A A
ARHEL (k:c) (k;
8 2c
( u (x+1i) VAH 1 ( (u-t)))))))
6
[
0i ( of 2 of (
de
( -f
-
8
(oi ( (0i+1)a(e)) 1
(di+1/a(-e)) ( (06+1)a(-e))
(40)
2e
8
A ARHEL
(K,e)
(
VBl1 0,u
2e
5
(dEt1/a(() + (1+22a(()) (=(k1)
4c
(08+1/0l-00 o 02/2/1( (c)) u ((+11
5
(8+1/s(()
(x+1) $i+113( (()
(41)
2c
The gradient estimator with respect to the input of the recurrent transformation is:
ARHEL (k,
V
H
1/2
2 (c). 0. ((+1)) - 1/H1/2/0/ (-c), (a-(++))))
2e
si (1-1//(()) 22 68+2/3(6) B
-BT
p6 (fit1/5(-e) + 2 (8+223((6)
8
B oi o111))))) -B TT (E+1/a(-E))
(42)
Parallel Scan. Similarly to the LinOSS model [29]. we can compute the recurrence of Linear HRU with parallel scan [5] reduce the computational time. To implement the parallel scan, we need to put the discretized dynamics in a form that can be computed in parallel. For this we vectorize the equation 38
4 k+1
0
I
4
2x42/3
I
0
0
2
I
SA
I
2R+1/3
8B u.k
I
8
852
B Ub
2
$ +
T
0
I
-5A
I
0
I
O uu
which gives us the discrete dynamics in matrix form:
x+1 M +
(43)
with:
M =
62 T
A
811
8 A]
FR
8
7
uk
(44)
5A
I- T
A
8B
Uy
To compute the echo passes T k E (1,..., nudging +e:
K}, we only need to adapt Fj. to get for the positive
g B- u- (k+1) + [ [++1)]
F
2
(45)
8B u (k+1) + v E (k++1)
42

---

A.4.4 Relation between our Linear HRU and the one of LinOSS
In this section, we expand on the relationship between Linear HRU and the used in the LinOSS model [29]. We first recall that we are using the base continuous- dynamics as in the LinoSS model (equations) 17
RHEL requires non-dissipative integrator. RHEL in continuous time requires non-dissipative systems (see ). Hence, to extend RHEL to discrete time, we sought an integrator that preserves this property and had rule out the first integrator used in the LinOSS model (the Implicit time integration (IM), equations and 4 of the (29)). We note that this RHEL requirement can be a weakness for some Machine Learning tasks. As highlighted in the empirical results of [29], performance depends the data distribution. On song-sequence tasks like Worms (17,944 timesteps), dissipative systems have an advantage because dissipation enables forgetting of past information. In contrast, non-dissipative Hamiltonian systems must preserve the entire trajectory history and may require larger models for equivalent performance. On the contrary, data distributions coming from non-dissipative Hamiltonian systems, non-dissipative integrators are better models (29). Appendix C.2).
Our adaptation of the IMEX integrator of LinOSS. In the LinOSS model, a non non-dissipative integrator also proposed, the Implicit-cxplici time integrator (IMEX). Adapted our notation (from equations, the equations of the integrator IMEX are:
{
Th+11 n& + 0 AΦ + BuR+1)
Φ ΦE 7
Re-using the Euler integrators defined in A.3 the IMEX integrator consists of doing an Euler step to change the momentum, followed by an Euler step change the position:
IMEX M H.6
Mr.6 Mv.5
In compact matrix form the IMEX integrator writes as:
D =
M IMEX 4b +
FIMEX K++1
(46)
[ T 2
87
8 Bus+1
where
M IMEX
IMEX K++
(47)
8 A
1
5 Bu
This shows that IMEX is still different from the integrator we used in our linear HRU (Eq, 44). A quick calculation shows that:
(MIMEX) -1 (MT, Mv.s) o # Mr.-s o My,-5
M H.-3
IMEX
which violates the discrete-time reversibility required to prove RHEL (Lemma A.4 On the contrary, the leapfrog integrator is the simple rearrangement of update steps:
MH = MT,6/2 o My. MT,8/2,
which has the desired properties, as proven in Lemma A.4
A.4.5 Nonlinear HRU
Hamiltonian of the recurrence. The nonlinear HRU block is the composition of a nonlinear spatial transformation and nonlinear recurrent transformation (see Eq. 18). The Hamiltonian of the nonlinear HRU block is defined by the following equations:
H d 0. u = T( 0, ul + V( 0, u4]
1
||
|| + (A-T log (cosh (A B u b)))
(48)
43

---

Dynamics.
The dynamics of the nonlinear HRU block are defined by the following equations:
π
(49)
(tanh (AΦ Bu + b) + aΦ)
Which, after time-discretization with the integrator defined in A.2.1 gives the following equations:
Φu+1/3 = ~x + V
TR, 0, un]
2
= = +
2 T
T/+1/3 = Tkk
ΦL+2/3 = 1x+1/3
Tk+2/3 = T/-1/3 + OV (0K1/13 un
(50)
= 1/3 8 (tanh [AΦk+1/3* + Bux + b) + aΦx+1/3)
@k+1 = 0242/33 +
V T Tk42233 uz
2
=01-2/3 +
2-333/3
2
Wh+1 Tk-2/3
A.4.6 Differentiating the time-discretization
In training HRUs, one can also train the multidimensional time-dsccretzzatio 8 E R& that is assumed be diagonal. For this, it suffices reparametrize the integrator and the Hamiltonian, so that the time discretization becomes parameter of the Hamiltonian. Hence the HRU equation ? becomes:
k++ MA DE, O, ug 8] Vk = =
(51)
with H is reparametrization of the Hamiltonian И such that MA 1 O, un 6]
MR,a (2x, un).
For instance, for the linear HRU, we can reparametrize the Hamiltonian as:
H u, 2" 8π (50 (6A) (((B)u
1
(52)
And for the nonlinear HRU, we can reparametrize the Hamiltonian as:
H 1. u. π 8m + Top dΦ + (5A- llog (cosh (A 1b B u b)))
(53)
(54)
Note that to match with previous implementations of the nonlinear HRU 33 we used the following parametrization for the discretization step:
1 1 u. 2" To( (00) 2 (a(8)A -. log (cosh (A x B + b)))
(55)
where (6) is diagonal matrix with (6);: 0.5 + 0.5 (anh(6.:/2)
A.4.7 Echo passes with automatic differentiation
As mentioned Remark learning HSSMs with RHEL involves chaining RHEL gradient estimators with Automatic Differentiation (AD). This design makes RHEL implementation readily compatible with modern automatic differentiation frameworks such as PyTorch or JAX. These libraries provide the capability to create custom functions that, when invoked within computational
44

---

Algorithm 3 Custom Reverse Mode Automatic Differentiation for an arbitrary HRU
1: @custom_autodiff
2: function HRU(O RRU, RU,
3: Φ f- HRU-HELPER(O,0O O 0)
4: return T
5: end function
6: function FORWARDHRU(O//R 2
7: Φ f HRU-HELPER(O,, 0 HHRO u.
8: return PK
9:
end function
function BACKWARDHRU(g.. (oowwrd)
11:
HHRU. u <- input_ forward
12:
á + REVERSE
Reverse the input sequence
13:
4 +-
fo Dr + g.
15:
HRU-HELPER ( o HHRU. u.g
16: DOHRU AT t LEARNINGHRU
17: return D Au
Loss gradient with respect the input of FORWARDHRU(;
18:
end function
19: UiRUL.efiie.cooommmaaaatfffttfon
BackwardHHRU)
Algorithm 4 Helper functions for custom Automatic Differentiation
1: function HRU-HELPER( OHRU, T2) See Dynamics. in App. A.4.3 and A.4.5 for
concrete examples
2: Input: Initial State 2o. Parameters of the recurrence OHnU, Inputs प. Nudging
3: Output: State sequence 4
4: end function
5: function @ de τ. > See RHEL gradient estimator in App. A.4.3 for a concrete
example
6: Input: Echo passes T for input sequence time-reversed ū
7: Output: Loss gradient w.rt to the input of FORWARDHRU(;;) )): Ao Am
8: end function
graph, are automatically differentiated. Algorithm 3 demonstrates how to implement RHEL/AD chaining for the recurrent component of an HRU block.
To implement a custom autodiff function HRU, we must define two additional functions required by AD: ForwardHRU and BackwardHRU, The ForwardHRU function is invoked when the HRU function called within a computational graph to differentiated. The ForwardHRU function returns two elements: its output for computing the remainder of the computational graph & and a rexidual to be stored for the backward pass. For the HRU, the residual consists only of the final state & k from the forward pass. The BackwardHRU function is called during graph backpropagation. It receives as input the loss gradient from the upstream portion of the computational graph g. the input from the forward pass ú, and the residual 4x. The BackwardHRU function computes the gradient of the loss with respect to the HRU block parameters, which is subsequently used to update the HRU parameters. It also computes the gradient of the loss with respect to the forward pass input which is then used to propagate the loss gradient backward through the computational graph.
These three functions are then registered with the autodiff library using the HRU. infinee. custom weffine_cutom__mmmmmeee autodiff function, enabling the library to automatically differenti- ate the HRU function when it is called within computational graph.
45

---

The three functions utilize two helper functions; HRU-HELPER and LearningHRU, The HRU-HELPER function implements the HRU dynamics and is employed in both the forward and backward passes. The LearningHRU function implements the RHEL gradient estimator and is used in the backward pass to compute the gradient of the loss with respect to the HRU block parameters.s
46

---

A.5
Experimental details
Summary.
This last section provides additional experimental details. More precisely:
We provide details about the datasets at use, the devices we used, well as detailed hyperparameters.
Importantly, we detail how RHEL can be subject to numerical underflow and how we mitigated this issue see Fig. 5 for details.
A.5.1 Datasets
This series of tasks recently introduced subset of the University of East Anglia (UEA) datasets with the longest sequences for increased difficulty and recently used to benchmark the linear
HSSM previously introduced 29
The classification datasets are drawn from a recently introduced benchmark [34] that selects a subset of the University of East Anglia (UEA) datasets [35, specifically choosing those with the longest sequences to increase difficulty, and which has been recently employed to evaluate the linear HSSM model [29). These datasets include EigenWorms (17,984 sequence length, 5 classes), SelfRegulationnCCPI (896 length, classes), SeliRegulationSCP22 (1,152 length, classes), EtthnnolConcentrationn (1,751 length, 4 classes), Heartheat (405 length, classes), and MotorImagery (3,000 length, 2 classes).
Additionally, we evaluate our HSSMs on the PPG-DaLiA dataset 36 multivariate time series regression dataset designed for heart rate prediction using data collected from a wrist-worn device, includes recordings from fifteen individuals, each with approximately 150 minutes of data sampled at a maximum rate 128 Hz. The dataset consists of six channels: blood volume pulse, electrodermal activity, body temperature, and three-axis acceleration. After splitting the data, sliding window of length 49920 and step size 4992 is applied, representing a challenging very long-range interaction task.
A.5.2 Simulation details and ressource consumption
Simulation details. The code to run the experiments is implemented using the JAX auto- differentiation framework [74]. All experiments were run on Nvidia V100 GPUs, except for the PPG experiments, which were run on Nvidia Tesla A100 GPUs due to larger memory demands.
Computational complexity. We do not claim any superiority of RHEL over BPTT on conventional digital hardware in terms of runtime memory consumption. As mentioned in the "Future work" paragraph in the main text, turning RHEL into a compelling alternative to BPTT on digital accelerators with respect these metrics is research direction of its own. In its current form, RHEL exhibits the same computational and memory complexity as BPTT.
Single Hamiltonian Recurrent Unit (HRU), For HRU with hidden dimension d and T timesteps, both RHEL and BPTT have a time complexity of OTT)) under naive recurrence, corresponding to d² multiply-accumulat operations per step. For linear HRUs that support the use of purallel scan 0. this cost can be reduced to (log(T)dd)) memory complexity of BPTT and RHEL is ( when storing all activations and (dd) when recomputing them backward from the final state, leveraging the time-reves rive-sivesiibility of the symplectic integrator inside HRUs. This memory efficiency is model feature shared by both algorithms, not a key differentiator of RHEL.
both cases, the
Hamiltonian State-Space Model (HSSM). For a HSSM composed of multiple HRUs, the above computational cost scales linearly with the number of HRUs for both RHEL and BPTT. For given HRU inside a HSSM, there is additional ( memory cost for storing the input sequence fed into that HRU.
Empirical runtime and memory usage. We report below the average runtime and peak GPU memory consumption RHEL and BPTT across representative datasets. As expected, the two algo- rithms exhibit comparable performance profiles. In some cases, RHEL yields moderate speedups for nonlinear models likely stemming from implementation details rather than fundamental algorithmic differences.
47

---

Table 3: Runtime and memory comparison between RHEL and BPTT.
Dataset
Model
RHEL Time (s)
BPTT Time (s)
RHEL Mem (GB)
BPTT Mem (GB)
PPG
Linear
11.50::4.72
10.42:1..8
11.387=0.001
11.387
PPG
Nonlinear
1150.32-:26.78
1609.05-124.74
11.728+0.001
11.731
Rigen Worms
Linear
3.16++5.38
3.21+6624
0.716++00001
0.716++0.001
FigenWorms Worms
Nonlinear
65.06++-2.33
117.66++-244
0.759+ 0.004
0.760++0.003
SCP2
Linear
15.59+15.53
13.02++8.87
0.104+0.000
0.103++0.000
SCP2
Nonlinear
27.96=44..0
21.28
0.107=00000
0.108110000
Observations.
Memory usage is comparable between RHEL and BPTT as expected.
Runtime is generally comparable as expected too, with RHEL showing significant speedup for nonlinear models (e.g., 2x faster for FigenWorms Worms with UniCORNN). We hypothesize that this speed-up is implemeenation--eeeedddn and would require further analysis to be conclusive.
A.5.3 Hyperparameters
We adopted the hyperparameters of (29 without modification, their experiments utilized the IMEX integrator, which, like our approach, derives from the LeapFrog integrator. The hyperparameters are: learning rate (Ir), number of layers (#blocks), number of hidden neurons (hidden dim). state- space dimension (state dim), and whether the time dimension is sent as input (include time). These hyperparameters
found by grid search and presented in Table
We note that the implementation of the linear HSSM models is based on the code of [29] and uses complex numbers for implementation, which corresponds to doubling the state-space dimension mentioned in Table We found that this dimensional doubling was necessary to recover the performance reported in the original paper, same doubling parameters for the nonlinear HSSM.
maintained this approach. We did not apply the
We reused the optimization scheme from 29 using the ADAM optimizer with default parameters and an early-stopping procedure based on the validation loss.
For the RH F.T. algorithm, we have two additional hyperparameters: the nudging strength and the scaling factor 7 (see Appx. A....4). The nudging strength e was 10 without prior tunning. For the scaling factor we did a grid search over the values {10%, 10a 102. 10*3 for the regression task (PPG-DaLiA) and found that the best performing parameter was 10 For the classification tasks, we performed grid search over the values {10%, 10 10 10 and found that the best- best-performing scaling was 10 based on the averaged score.
Table 4: Hyperparameters for the linear and nonlinear HSSM model
Dataset
Ir
hidden dim
state dim
#blocks
include time
Worms
0.0001
64
16
2
False
SCP1
0.0001
64
256
False
SCP2
0,00001
64
256
True
Ethanol
0.00001
16
256
False
Heartbeat
0.00001
64
16
True
Motor
0.0001
16
256
True
PPG
0.0001
64
16
2
True
We employ different initialization schemes for linear and nonlinear HSSM variants. For the linear HSSM model, we follow the same initialization scheme as 29
Uniform(O,
Uniform
hidden, dim hiddlen dim
48

---

C Uniform
state_ dim
state dim
D 10,,
.
Uniform(0) (0,
For the nonlinear HSSM model, we adopt the initialization scheme from 30
Uniform(0.5,
B ~ Uniform
hidden dim hidden dim
. C = 0
. D 1))
b 1).
a Un.f,
8 Uniform(-1,1)
A.5.4 Gradient re-scaling for dealing with numerical instabilities
Analysis. From the theoretical results RHEL, the nudging strength should be chosen as small as possible to accurately estimate the gradient of the loss function. However, naive numerical implementation can encounter underflow issues. RHEL, gradient information is encoded small perturbations to the state 4 that generate the echo trajectory T This presents numerical challenges when the perturbations and state values differ by several orders of magnitude.
In practice, when using finite precision representations (such as floating-poit numbers), the perturba- tions may be lost due discretization error, where small values cannot be accurately represented distinguished from zero in finite precision arithmetic.
To address this numerical challenge, employ simple gradient rescaling method. We multiplica- tively scale the output loss .] by a constant y > 1 during the echo dynamics computation. This amplification ensures that the perturbation-indueee changes in the loss remain within the representable range of peating-poim arithmetic. Subsequently, we divide the RHEL parameter gradient estimate (A (HEL, c)) by the same scaling factor y recover the unbiased gradient.
effectiveness this gradient rescaling method on both Linear and Nonlinear HSSM. To evaluate the effect of gradient scaling, we compare the gradients of RHEL and BPTT. Given a HSSM model, we randomly sample tuple (u,, y:) ~ D from the SCP1 dataset (see App. A. and compute the gradients of the loss function with respect to the parameters of the recurrence of the HRUs in Fig. 5A for the linear HSSM and Fig. A for the nonlinear HSSM.
now demonstrate
To gain more fine-grained understanding, we also compute the gradients with respect to the inputs of the third layer of the HSSM (see Eq. 17 and 18 Fig. 5B and Fig. GB. Compared to the parameter gradients, these input gradients are not time-averages and hence reveal more detail about the underflow issues.
Gradient rescaling recovers the underflow issues. As general pattern across both architectures, we observe that the unscaled (7 = 10") parameter and input gradients tend to be biased: they have norm ratio above 1.0 and cosine similarity below 1.0. We also note that for both the linear and nonlinear cases, there exists a scaling factor (10% for linear and 10 for nonlinear) that recovers nearly perfect match between RHEL and BPTT.
Additionally, for the input gradients that are not time-an -verrged (first column Fig. 5B and Fig. 6B), there is large amount of noise in the unscaled gradients, which is eliminated for the best-performing scaling factor (106 for linear and 10 for nonlinear). We conjecture that this noise is due to the underflow effects described above.
Linear HSSM: scaling improves gradient matching. For the linear analysis, we observe general pattern: more we scale the RHEL gradient, the better matches BPTT, both for input gradients and parameter gradients (Fig. AAA and B).
49

---

Nonlinear HSSM: optimal scaling balances underflow and linearization errors. For the nonlin- ear analysis, we also observe that scaling the RHEL gradient improves the match between RHEL and BPTT gradients. However, for higher values of the scaling factor (7 = 10%), we observe that both input and parameter gradients show drop in matching quality. We conjecture that this is due to linearization errors. If the loss gradient is too large, the nudging it drives will be large, and the finite difference method used for the learning rule will no longer be valid.
Solution Implemented. For the experiments of the main text, we conducted grid search over the scaling parameter y. In our initial implementation, we only applied gradient scaling without the corresponding downscaling step, effectively amplifyyng the gradient magnitude throughout training. This provided valuable insights into the sensitivity of RHEL dynamics to gradient scaling and improved performance. The complete rescaling procedure (including both upscaling and downscaling) used in the current section is also implemented in our code base to provide a more comprehensive analysis of the numerical stabilization approach.
A.5.5 Vanishing and Exploding gradients of the Nonlinear HRU
In this section we validate empirically the claim made in section 4 that the Nonlinear HRU mitigate by design the problem of vanishing and exploding gradients. We emphasize that i) this claim has been demonstrated in a prior work [30, as such is not one of ours, ii) this claim is a property of the model, not of the gradient computation algorithms, therefore it equally applies BPTT and RHEL.
Theoretical claims. The gradient stability properties stem from the underlying Hamiltonian struc- ture with symplectic integration. Prior work 30 rigorously demonstrated that the recurrent compo- nent of our nonlinear HRU block (corresponding to the Hamiltonian in Equation 18) possesses:
Non-exploding gradients (Proposition 3.1 in (30]): the total loss gradient (i.e. aggregating loss "sensitivities" across all timesteps) is upper bounded.
Non-vanishing gradients (Proposition 3.2 in 30 ): when looking far enough back in time, the loss "sensitivities" (i.e the additive time-wise contributions to the loss gradient) are strictly non-zero and independent of the timestep.
Experimental setup. We randomly sampled (u,,y:) ~ D from the SCP1 dataset and computed time-dependent sensitivities for both RHEL (A)) and BPTT (ga(k)) across the sequence length. We then performed linear regression analysis:
Analysis
Correlation r
Slope
Interpretation
A AARHEL (k:]] vs [k)]l
0.99998
1.0007
RHEL R BPTT gradients
A RHELL (A)]l k
0.1 1068
2.82 x 10 15
No vanishing/explodingg gradient
[
vs k
0.
2.98 x 10 15
No vamishing/exploding gradient
Table 5: Linear regression analysis between RHEL and BPTT gradient norms and timestep index k.
Key findings.
* Row 1 confirms RHEL's gradient fidelity (r R 1, slope Rs 1).
* Rows 2-3 demonstrate gradient stability: slopes x 0 indicate neither exponential growth (exploding) nor decay (vanishing) over time.
This empirical evidence supports our theoretical claims and demonstrates that RHEL our proposed nonlinear HRU preserves the favorable gradient properties of the underlying Hamiltonian architecture. We note, however, that this analysis only holds at the level of a single HRU; vanishing gradients can still arise across layers in deeper architectures (i.e., within the feedforward transformations coupling successive HRUs). Both RHEL and BPTT are subject to these issues, they compute mathematically equivalent gradients, and both benefit equally from standard architectural remedies such as skip connections (see Appendix C.5-C.5 in (30).
50

---

(A)
A00
Avg A
B
Avg. B(()
540
Avg. 545
10
10
679
0.1
Contae Sumbarty
833
05-
-175
-3.0
-05
-12
30-
Neem Itatoo
10-
1/
10p
Layss
Layer
12
(B)
1.-
1.00
Closisee Similaity
0.5
0.75
0.50
0.25
0.00
200
300
600
800
1000
100
106
Norm Rttti
202
10a
102
001
10a
200
300
600
800
10P
1nr
10c
Figure 5: (A): Parameters gradients comparison between RHEL and BPTT for HSSM. Given some (u, yi) DD, we perform BPTT and RHEL on six blocks-deep linear HSSMs for different scaling factor 7 of RHEL (different colors). We measure, per layer (line plot) and when averaged across layers (bar-plot), the cosine similarity (top panels) and norm ratio (bottom panels) between RHEL and BPTT parameters gradients of linear HSSM (Eq. (17). (B): Inputs gradient com- parison between RHEL and BPTT for HSSM. Same setting as (A) but we focus on the gradients with respect to the inputs of the third layer (3) see Eq. 17). We measure, per time steps (line plot) and when averaged across time (bar-plot), the cosine similarity (top panels) and norm ratio (bottom panels) between RHEL and BPTT inputs gradients.
51

---

(A)
490
A
B
B14
610
Avg. g
0.5
03-
Sindlritt
06
04
3.4
64
0.2
0.3
60
-12
Norta Ratoo
Layer
Layer
Lager
(B)
1.00
Cosite Simllarity
0.75
0.75
0.50
0.25
0.25
0.00
0.00
200
600
800
102
10
10"
10a
Norm Ratio
10a
10a
6 10-1
10
10x
10
Figure 6: (A): Parameters gradients comparison between RHEL and BPTT for HSSM. Given some (u, 9:) ~ 'D2 we perform BPTT and RHEL on six blocks-deep nonlinear HSSMs for different scaling factor 7 of RHEL (different colors). We measure, per layer (line plot) and when averaged across layers (bar-plot), the cosine similarity (top panels) and norm ratio (bottom panels) between RHEL and BPTT parameters gradients of nonlinear HSSM (Eq. (8)) (B): Inputs gradient comparison between RHEL and BPTT for HSSM. Same setting as (A) but we focus on the
gradients with respect to the inputs of the third layer (3), see Eq. 18 We measure, per time steps (line plot) and when averaged across time (bar-plot), the cosine similarity (top panels) and norm ratio (bottom panels) between RHEL and BPTT inputs gradients. For (A) and (B), ? = 10 was also tested but produced numerical instabilities leading NaN values and is therefore omitted from the results.
52

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification:
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made in the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or
NA answer to this question will not perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how be expected to generalize other settings.
much the results
is fine to include aspirational goals as motivation as long is clear that these goals are not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification: included the discussion section.
Guidelines:
The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those not discussed in the paper.
The authors are encouraged to create separate "Limitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification. asymptotic approximations only holding locally). The authors should reflect how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution low images are taken in low lighting. Or speech-to-text system might be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach address problems of privacy and faimess.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't 't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory assumptions and proofs
Question:
For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
Answer: [Yes]
53

---

Justification:
provided in the appendix.
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated referenced in the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but they appear in the supplemental material, the authors are encouraged to provide short proof sketch provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental result reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that affects the main claims and/or conclusions of the paper (regardless whether the code and data are provided or not)?
Answer: [Yes]
Justification: all pseudo-algorithms and hyperparameters are presented in the appendix.
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution is a specific model and empirical evaluation, it may be necessary to either make possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g.. in the case of large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed.
While NeurIPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily a new algorithm, the paper should make clear how to reproduce that algorithm. (b) the contribution is primarily a new model architecture, the paper should describe the architccture clearly and fully.
(c) If the contribution is new model (e.g., a large language model), then there should either way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome to describe the particular way they provide for reprodueibility. the case of closed-source models, may be that access to the model is limited in some way (e.g., to registered users), but it should be possible for other researchers have some path to reproducing or verifying the results.
5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, as described in supplemental material?
54

---

Answer: [No]
Justification: we plan to release the code in an anonymized fashion by the rebuttal phase upon request of the reviewers.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurIPS code and data submission guidelines (https://nips... cc/
htaia///guiiee/ooeeSuniiiii colloooo for more details.
While we encourage the release of code and data, we understand that this might not be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to to reproduce the results. See the NeurIPS code and data submission guidelines (https
ssiicyolicy) for more details.
The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended to the paper) recommended, but including URLs to data and code permitted.
6. Experimental setting/details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: included the appendix. Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper to level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, or as supplemental material.
7. Experiment statistical significance
Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Question:
Answer: [Yes]
Justification: all experimental results are reported across five different seeds with mean and standard deviation.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by error bars, confi- dence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper.
The factors of variability that the error bars capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g., Normally distributed errors).
55

---

It should be clear whether the error bar is the standard deviation or the standard error of the mean.
is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not to show in tables figures symmetric error bars that would yield results that are out range (e.g. negative error rates).
If error bars are reported in tables or plots, The authors should explain the text how they were calculated and reference the corresponding figures tables in the text.
8. Experiments compute resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?
Question:
Answer: [Yes]
Justification: an estimation of the compute ressources is included in the appendix.
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't 't make it into the paper).
9. Code of ethics
Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics http://awaragcccca//a/////aacaaa
Question:
Answer: [Yes]
Justification:
Guidelines:
The answer NA means that the authors have reviewed the NeurIPS Code of Ethics.
If the authors answer No, they should explain the special circumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g., there a special consid- eration due to laws or regulations in their jurisdiction).
10. Broader impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [Yes]
Justification: discussed in the introduction and discussion section.
Guidelines:
The answer NA means that there no societal impact of the work performed.
If the authors answer NA or No. they should explain why their work has no socictal impact why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
56

---

The conference expects that many papers will be foundational research and not tied particular applications, let alone deployments. However, if there is direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.
If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition attacks, mechanisms for monitoring misuse, mechanisms to monitor how system learns from feedback over time, improving the efficiency and accessibility ML).
11. Safeguards
Question:
Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, scraped datasets)?
Answer: [NA]
Justification: we only use public datasets and our models pose no safety risk as the experi- ments were primarily designed for proof-of-concepts.
Guidelines:
The answer NA means that the paper poses such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is challenging, and many papers do require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Question: Are the creators or original owners of assets (e.g.. code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Answer: [NA]
Justification: only used public datasets.
Guidelines:
The answer NA means that the paper does not use existing assets.
The authors should cite the original paper that produced the code package dataset.
The authors should state which version of the asset is used and, if possible, include URL.
The name of the license (c.g., CC-BY 4.0) should be included for each asset.
For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided.
If assets are released, the license, copyright information, and terms of use in the batasets.papersstt paperswi thcode. cww/ datasets package should be provided. For popular datasets,
has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided.
57

---

If this information is not available online, the authors are encouraged to reach out to the asset's creators.
13. New assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [NA]
Justification:
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details of the dataset/code/mmoel as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent obtained from people whose asset is used.
At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsoureing and research with human subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well details about compensation (if any)?
Answer: [NA]
Justification:
Guidelines:
The answer NA means that the paper does not involve crowdsoureing nor research with human subjects.
Including this information in the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included the main paper.
According the NeurIPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional review board (IRB) approvals or equivalent for research with human subjects
Question:
Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
Justification:
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Depending the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and expect authors to adhere to the NeurlPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such the institution conducting the review.
16. Declaration of LLM usage
58

---

Does the paper describe the usage of LLMs if it is important, original, or non-standard component of the core methods this research? Note that if the LLM used only for writing, editing, or formatting purposes and does not impact the core methodology, scientific rigorousness, or originality of the research, declaration is not required.
Question:
Answer: [NA]
Justification:
Guidelines:
The answer NA means that the core method development in this research does not involve LLMs any important, original, non-standard components.
Please refer to our LLM policy for what should or should not be described.
59


Reviews:
[{"review_id": "AnTV8MIA7M", "reviewer": "Reviewer_AkMK", "summary": "This paper proposes Recurrent Hamiltonian Echo Learning (RHEL), a forward-only gradient-estimation algorithm for non-dissipative Hamiltonian systems. RHEL works by perturbing a time-reversed copy of the forward trajectory with a small ``nudge'', and shows that the resulting finite differences converge to the continuous adjoint-state gradients. A discrete variant is derived, and this paper shows its equivalence with BPTT. Empirical validation occurs on six UEA long-sequence classification datasets and the PPG-DaLiA regression task. HRUs trained with RHEL match BPTT (although mean performance is a little lower, this is typically not statistically significant).", "questions": "- Why is it surprising/need to be posited that \"some SSMs could inherently be mapped onto dynamical systems\"?\n- Theorem 3.1, if you take $\\epsilon$ to 0, wouldn't both methods give the echo trajectory simply be the forward trajectory in reverse for both cases?\n\nMinor:\n- Page 2, Line 48: three -> two? I think I only saw two conditions.\n- Page 3, Line 127: let -> left\n- Page 6, Line 219: rewrite -> can be rewritten\n- Page 7, Line 221: spurious \"reading\"?\n- Page 9, Line 296: Remove \"Yet\" \n- It is a bit misleading to bold font all the averaged numbers in Table 1.\n- Recommend to use hyperref so that the refs and cross-refs are clickable.", "limitations": "Yes, the authors are upfront about the limitations of the method and that it slightly underperforms BPTT.", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 2, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strengths:\n- The method proposed offers a hardware-friendly alternative to BPTT and only uses three forward passes,\n- Theorems are also supported with proofs (some sketched in text, others in Appendix).\n- The diagrams in the paper are good, very useful for illustrating the method. \n\nWeaknesses:\n- There are no results/values of runtime, and this reduces the selling point of the paper being \"forward-only\".\n- The writing is not very clear at times and sometimes seem a bit overdramatic. For example, \"the physical realization of our algorithm\" in the abstract is not very clear. and \"this paper embraces an extreme view\" is quite dramatic. \n- The structure and value of the Supplementary materials need to be made clearer. The informal versions of the theorems are stated in the paper and the formal versions are not signposted.\n- Many systems in reality are not conservative, which may reduce the impact.", "was_revised": true, "final_justification": "The authors clarified some misunderstandings and provided additional interesting analysis. The suggested changes would strengthen the paper and I am raising my score from 4 to 5.", "rebuttal": "We thank the reviewer for their valuable time to give constructive feedback and spot typos. We address below the specific concerns they raised and accordingly provide **four propositions** of modifications of the paper in case of acceptance.\n\n> There are no results/values of runtime [...].\n\n**Answer 1**. \nWe do not claim any superiority of RHEL over BPTT on conventional *digital* hardware in terms of run time or memory consumption. As mentioned in the \"Future work\" paragraph L.306-308 of our submission: turning RHEL into a compelling alternative to BPTT on digital accelerators with respect to these metrics is a direction of research of its own.\n\nIn its current form, **RHEL has the *same* memory and computational cost as BPTT**:\n\n- on a **single Hamiltonian Recurrent Unit** (HRU) using $T$ timesteps and with a hidden layer dimension $d$: when using \"naive\" recurrent computation, the time complexity of BPTT and RHEL is $O(Td^2)$, *i.e.* $d^2$ MACs per step. For *linear* HRUs which allow for the use of parallel scan [1], this cost can be reduced to $O(log(T)d^2)$. In both cases, the memory complexity of BPTT and RHEL is $O(Td)$ when naively storing activations of the forward trajectory, and $O(d)$ when re-computing them backwards from the final state leveraging the time-reversibility of the symplectic integrator at use inside HRUs. As mentioned L.150-152 of our submission, this memory efficiency is a *model* feature that equally applies to BPTT and RHEL, not a key differentiator of RHEL.\n\n- On a **Hamiltonian State Space Model** (HSSM): the above computational cost scales linearly with the number of HRUs for both BPTT and RHEL. For the memory cost on a given HRU taken inside a HSSM, there is an extra memory cost of $O(Td)$ for both BPTT and RHEL to store the input sequence which is fed into a given HRU.\n\nThis being said, we acknowledge the value of providing these metrics in general. Below are the runtime and memory comparisons:\n\n| Dataset | Model | RHEL Time (s) | BPTT Time (s) | RHEL Mem (GB) | BPTT Mem (GB) |\n|---------|-------|---------------|---------------|---------------|---------------|\n| PPG | Lin | 11.50±4.72 | 10.42±6.86 | 11.387±0.001 | 11.387±0.002 |\n| PPG | Nonlin | 1150.32±26.78 | 1609.05±24.74 | 11.728±0.013 | 11.731±0.013 |\n| EigenWorms | Lin | 3.16±5.38 | 3.21±6.24 | 0.716±0.001 | 0.716±0.001 |\n| EigenWorms | Nonlin | 65.06±2.33 | 117.66±2.43 | 0.759±0.004 | 0.760±0.003 |\n| SCP2 | Lin | 15.59±15.53 | 13.02±8.87 | 0.104±0.000 | 0.103±0.000 |\n| SCP2 | Nonlin | 27.96±42.50 | 41.28±29.39 | 0.107±0.000 | 0.108±0.001 |\n\n**Key observations:**\n\n* **Memory usage** is comparable between RHEL and BPTT as expected.\n* **Runtime** is generally comparable as expected too, with RHEL showing significant speedup for nonlinear models (e.g., 2× faster for EigenWorms with UniCORNN). We hypothesize that this speed-up is implementation-dependent and would require further analysis to be conclusive. Yet again, as stated above, the goal of our work is not to outperform BPTT in terms of these metrics.\n\n**Proposition 1**. We propose to include both the table runtimes and memory consumptions, as well as the complexity analysis of RHEL and BPTT.\n\n> [...] and this reduces the selling point of the paper being \"forward-only\".\n\n**Answer 2**.\nWe also want to clarify the position of our algorithm being \"forward-only\" with respect to *zeroth-order optimization* (ZO)[2] which the reviewer may have in mind. In the literature of digital hardware-friendly learning, ZO techniques are pursued for their **memory efficiency** compared to vanilla backprop with gradient checkpointing [3], which can be further enhanced with appropriate quantization schemes [4]. The relatively recent exploration of ZO techniques on Large Language Models (LLMs) finetuning [5] has spurred a revival of \"forward-only\" techniques for memory-cheap gradient computation in LLMs on a variety of applications. \n\nHowever, the motivation for the \"forward-only\" feature of RHEL is *not* memory efficiency. Though we make use of the same terminology, \"forward-only\" gradient computation algorithms are very differently motivated in the *physical computing* literature [6] which RHEL belongs to -- see next section of our rebuttal. In our context, \"forward-only\" means that both the forward *and* backward passes obey the *same* physical principles, for instance obeying to Hamiltonian dynamics as required by RHEL. This means that in principle, any piece of hardware that can sustain inference can also sustain gradient computation, out of the same physics, and therefore learn \"by itself\".\n\n**Proposition 2**. We propose to add an extra related work paragraph dedicated to forward-only algorithms and explain how their motivation contrast with ours. \n\n> The writing is not very clear at times and sometimes seem a bit overdramatic. For example, \"the physical realization of our algorithm\" in the abstract is not very clear. and \"this paper embraces an extreme view\" is quite dramatic.\n\n**Answer 3.**\nWe want to acknowledge the lack of clarity of our wording revolving around the motivation for physical computing.\n\nEvery computation, be it on a digital accelerator or any alternative substrate, is sustained by a physical process carried by wires and transistors which creates entropy [7]. Therefore the very notion of \"physical computing\" may sound pleonastic. However, digital hardware designs abstract core physics away to enable idealized computations relying upon statelessness, unidirectionality, determinism and synchronization [6]. Compilation pipelines are then typically designed to be as hardware-agnostic and as general-purpose as possible. In this paradigm, \"hardware-software\" codesign often amounts for instance to write bespoke low-level kernels for specific workloads to mitigate their compute, memory costs and off-chip memory accesses [8]. Still, these approaches remain largely physics-agnostic.\n\n\"Physics-based\" computing significantly deviate from this framework in two ways. First, the designs may target Application-Specific Integrated Circuits (ASICs) rather than general-purpose machines [6]. Second, the aforementioned digital design constraints may be relaxed and the resulting *analog* physics leveraged. For instance, supply voltages may be decreased and resulting circuits become non-deterministic, so that stochastic algorithms can be inherently mapped to such systems [9]. Uni-directionality may be relaxed such that vector-matrix multiplication can be natively implemented by Kichoff laws of currents [10]. Optimization problems may be embodied into a physical system and subsequently solved as the system settles to equilibrium [11]. As a particular case, physical \"self-learning machines\" are physical embodiments of neural networks whose inference and gradient computation are, *for instance*, solved by relaxing to equilibrium [12, 13]. This is how \"physical realization of our algorithm\" should be understood as.\n\n**Proposition 3**. We will remove the sentences which may sound \"overdramatic\" (L. 9 and L.29 ), replace them with clearer and softened statements based on the comments above and add a related work paragraph around physical computing. \n\n> Why is it surprising/need to be posited that \"some SSMs could inherently be mapped onto dynamical systems\"?\n\n**Answer 4.** \nWe acknowledge this is not surprising, as SSMs are themselves derived as discretization of continuous-time dynamics [14]. However, the primary motivation of SSMs is not the design of \"physics-based\" models but a means to overcome the quadratic cost in sequence length of transformer-based training on digital hardware [8]. \n\n> Many systems in reality are not conservative, which may reduce the impact.\n\n**Answer 5**.\nWe acknowledge this limitation. See our answers to Reviewers  WcdN (Answer 1, Proposition 1) and k16k (Answer 1, Proposition 1) for greater details.\n\n> The informal versions of the theorems are stated in the paper and the formal versions are not signposted.\n\n**Proposition 4.** Acknowledged, this will be fixed in case of acceptance. \n> Theorem 3.1, if you take  to 0, wouldn't both methods give the echo trajectory simply be the forward trajectory in reverse for both cases?\n\n**Answer 6**.\nWe assume that the reviewer implies that both $\\Phi^e(t, \\epsilon)$ and $\\Phi^e(t, -\\epsilon)$ converge to the time-reversed trajectory $\\Phi^\\star(-t)$ as $\\epsilon \\to 0$ and therefore the $\\Phi^e(t, \\epsilon) - \\Phi^e(t, -\\epsilon)$ vanishes. Fortunately, the error signal is carried by $\\epsilon^{-1}(\\Phi^e(t, \\epsilon) - \\Phi^e(t, -\\epsilon)) \\to \\partial_\\epsilon(\\Phi^e(t, \\epsilon))|_{\\epsilon=0}$ and converges to a *finite* value as $\\epsilon\\to 0$. Our proof of Theorem 3.1 heavily relies on this fact.\n\n[1] Guy E Blelloch. “Prefix Sums and Their Applications”. In: (1990).\n\n[2] Spall. \"Multivariate stochastic approximation using a simultaneous perturbation gradient\napproximation\". 1992.\n\n[3] Chen et al. \"Training deep nets with sublinear\nmemory cost\". 2016.\n\n[4] Zhou et al. \"QuZO: Quantized zeroth-order fine-tuning for large language models.\" 2025.\n\n[5] Malladi et al. \"Fine-tuning language models with just forward passes.\" 2023.\n\n[6] Aifer, Maxwell, et al. \"Solving the compute crisis with physics-based ASICs.\" 2025.\n\n[7] Wolpert, David H. \"The stochastic thermodynamics of computation.\" 2019.\n\n[8] Gu & Dao. \"Mamba: Linear-time sequence modeling with selective state spaces.\" 2023.\n\n[9] Coles et al, \"Thermodynamic AI and the fluctuation frontier\". 2023.\n\n[10] Burr et al. \"Neuromorphic computing using non-volatile memory\". 2017.\n\n[11] Vadlamani et al , \"Physics successfully implements lagrange multiplier optimization\". 2020.\n\n[12] Stern et al. \"Supervised learning in physical networks: From machine learning to learning machines.\" 2021.\n\n[13] Kendall et al. \"Training end-to-end analog neural networks with equilibrium propagation.\" 2020.\n\n[14] Gu et al. \"Efficiently modeling long sequences with structured state spaces.\" 2021."}, {"review_id": "9PRHQARVeL", "reviewer": "Reviewer_UVRA", "summary": "This paper proposes a novel algorithm for learning continuous-time sequence models on non-dissipative physical device (Hamiltonian systems). In addition, it introduces a discrete-time version called Hamiltonian Recurrent Units (HRUs), providing a means for simulation. Furthermore, by proposing Hamiltonian State Space Models (HSSMs), which scale and extend HRUs in a hierarchical manner, the paper demonstrates that the proposed model can be expanded into sequence models such as State Space Models (SSMs), which have recently gained attention in ML. It also shows that HSSMs can learn very long sequence tasks through simulation.", "questions": "- In line 173, for the time derivative of $\\pi_1$ and $\\pi_4$ in Equation 8, it would be helpful to add explanations about $u$ (external force) and $y$ (target). As shown in Figure 2B, explicitly specifying the external force ($u$) and target ($y$) for the object would make it easier to understand.\n- In line 48, the text mentions three conditions: i) use “forward passes” only (i.e., no backward passes), ii) do not use explicit state-Jacobian (i.e., Jacobian of the system’s dynamics). However, only two conditions are described earlier. This seems to be a typo.\n- Regarding HSSMs with multiple layers, it is stated that activation patterns need to be stored to compute gradients. Does this mean physical implementation is impossible? Or is this limitation only applicable to the discretized time model, while the continuous-time model does not face such physical constraints?\n- The title, \"Learning long-range dependencies through time reversal symmetry breaking,\" seems misleading as there is no mention of time-reversal symmetry breaking in the main text. This could confuse readers. According to Lopez-Pastor & Marquardt (PRX 2023), \"*The self-learning device obeys time-reversal symmetry.*\" The title appears contradictory. Was the authors' intention to highlight that the system being studied has time-reversal symmetry, and that external driving forces applied when starting the backward trajectory create symmetry breaking?", "limitations": "Yes", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 4, "significance": 4, "originality": 3, "strengths_and_weaknesses": "**Strengths** \n- The paper is very clearly written and rigourous proof for each theorem.\n- The paper also compare their algorithm with the previous algorithm in detail and address the important limitation of the previous algorithm.\n- This paper is an excellent contribution as it develops a real on-chip self-learning system that reduces the high energy cost of GPUs while demonstrating scalability, such as with deep SSMs.", "was_revised": false, "final_justification": "", "rebuttal": "We appreciate the reviewer's time for providing a thorough evaluation of our work, making constructive suggestions and spotting typos. In the following, we respond to their specific points and outline **three proposed modifications** to strengthen the paper in case of acceptance. \n\n> In line 173, for the time derivative of $\\pi_1$\n and $\\pi_4$\n in Equation 8, it would be helpful to add explanations about \n $u$ (external force) and \n $y$(target). As shown in Figure 2B, explicitly specifying the external force ($u$) and target ($y$) for the object would make it easier to understand.\n\n**Answer 1**. Thank you for this valuable suggestion to improve readability and strengthen the connection between the mathematical formulation and Figure 2B. In the main text we indeed only presented the Hamiltonian and the dynamics of the raw physical system without its coupling with the input and nudging output, which makes it difficult to understand how the external force and target relate to the visual representation. We will present a condensed version of the full equations from the appendix (Equations 24-25) to make this connection explicit.\n\n\n**Hamiltonian with input**:\n\n$$\nH[\\Phi, \\theta, u] = \\sum_i^6\\frac{\\pi_i^2}{2m_i} + \\frac{1}{2}\\sum_{i}^6k_i \\phi_i^2 + \\frac{1}{2}\\sum_{i}^6\\sum_{j,j>i}^6k_{ij} (\\phi_j-\\phi_i)^2 + u\\phi^1\n$$\n\n**Complete Dynamics Equations:**\n\n$$\n\\partial_t\\phi_i = \\pi_i /m_i, \\quad \\partial_t \\pi_i = -k_i \\phi_i + \\sum_i^6 \\sum_{j, j\\neq i}^6k_{ij}(\\phi_j-\\phi_i) - \\delta_{i1}\\mathbf{u}\\phi^1 + \\delta_{i4}\\delta_e\\epsilon(\\phi^4-y)\n$$\n\nwhere $\\mathbf{u}$ is the external driving force applied to oscillator 1, $y$ is the target trajectory, $\\delta_e$ is the indicator function of the echo pass (equal to 1 during the echo pass and 0 otherwise), and $\\delta_{ij}$ is the Kronecker delta.\n\n**Proposition 1.** We propose to replace the equations of the toy example by the equation mentioned in answer 1 to mention the coupling with input/output of the system.\n\n> Regarding HSSMs with multiple layers, it is stated that activation patterns need to be stored to compute gradients. Does this mean physical implementation is impossible? Or is this limitation only applicable to the discretized time model, while the continuous-time model does not face such physical constraints?\n\n**Answer 2.**\nThis is an excellent question which raises several points:\n\n- this limitation is *solely* tied to the choice of the **Hamiltonian State Space Model** (HSSM) architecture introduced in Eq.16 of our submission. Whether HSSMs are simulated in discrete time *or* in continuous time, or physically implemented (see next bullet), they would *always* require to store (or somehow recompute) the input sequences that are fed into each HSSM in order to estimate gradients using our RHEL chaining algorithm.\n- Whether HSSMs are amenable to a physical implementation, regardless of the above limitation, is an interesting and *independent* question of its own. We could imagine a physical implementation of a HSSM as follows. As the forward pass of a HSSM reads a composition of recurrent and feedforward units, the corresponding \"RHEL chaining\" algorithm (Alg. 1 inside section 3.3 of our submission) reads as an *implicit* backward chaining of backprop and RHEL respectively. Namely: error signals are backpropagated through feedforward modules and \"RHEL-propagated\" through recurrent modules -- see Remark 8, p. 40 inside our appendix where we highlight this aspect. Therefore, an intuitive implementation of a HSSM would be to map feedforward parts onto *digital circuits* and recurrent parts onto *analog circuits*, with backprop and RHEL running on these digital and analog circuits respectively. This mapping is strongly reminiscent of the recently proposed backprop-Equilibrium Propagation hybridization [1]. Yet, this type of hybrid approach potentially comes at the cost of expensive ADC-DAC conversions across layers [2].\n- While we used HSSMs to provide of scalability proof-of-concept of RHEL and leverage most recent SSM architectures, there exists an alternative architecture which may not exhibit this limitation and may be more amenable to a physical implementation. If HSSMs can be regarded as model of hybrid systems comprising digital and analog parts being composed hierarchically, a model corresponding to a **fully analog** Hamiltonian system could be defined as a sum-separable Hamiltonian of the form:\n\n$$\nH = \\sum_{\\ell = 1}^L \\frac{1}{2} \\|\\pi^{(l)}\\|^2 + \\frac{1}{2} \\phi^{(\\ell)^\\top}\\cdot A^{(\\ell)} \\cdot \\phi^{(\\ell)} - \\phi^{(\\ell)^\\top}\\cdot B^{(\\ell)} \\cdot \\phi^{(\\ell - 1)}, \\quad \\phi^0:=u\n$$\n\n - Applying Hamiltonian dynamics (Eq. 1 in our submission) with this choice of Hamiltonian yields:\n $$\n \\dot{\\phi}^{(\\ell)}= \\pi^{(\\ell)}, \\quad \\dot{\\pi}^{(\\ell)}=-A^{(\\ell)}\\cdot \\phi^{(\\ell)} + B^{(\\ell)}\\cdot\\phi^{(\\ell - 1)} + B^{(\\ell + 1)^\\top}\\cdot \\phi^{(\\ell + 1)}\n $$\n \n - Note that the resulting equations are very different from that of a HSSM: instead executing $L$ trajectories hierarchically, a *single* trajectory updating all layers *simultaneously* at each time step is executed. As such, the only sequence which would need to be stored is just the *data sequence* $\\phi^{(0)}:=u$ itself.\n- Note also we could seamlessly introduce nonlinearities inside the proposed Hamiltonian above in the same way as it is done in the nonlinear HRU introduced in our submission (Eq. 19 of our submission).\n\n**Proposition 2.** We propose to add a detailed discussion around physical implementations of these various Hamiltonian models and RHEL in a dedicated paragraph.\n\n> The title [...] Was the authors' intention to highlight that the system being studied has time-reversal symmetry, and that external driving forces applied when starting the backward trajectory create symmetry breaking?\n\nYes, this was *exactly* our intention. And this is also consistent with the requirement that the system must obey time-reversal symmetry, as highlighted in the paragraph \"Assumptions\" of our submission (starting at L.123), and therefore also consistent with the seminal HEB work [3]. This being said, we acknowledge that the term \"symmetry breaking\" is never mentioned again in the rest of the paper which may hinder clarity. \n\n**Proposition 3.** We propose to make an explicit connection to \"symmetry breaking\" when introducing the RHEL algorithm. \n\n\n[1] Nest et al. \"Towards training digitally-tied analog blocks via hybrid gradient computation.\" 2024.\n\n[2] Li et al. \"Merging the interface: Power, area and accuracy\nco-optimization for rram crossbar-based mixed-signal computing system\". 2015.\n\n[3] Lopez-Pastor & Marquardt. \"Self-learning machines based on Hamiltonian echo backpropagation\". 2023."}, {"review_id": "EJlDyadfz9", "reviewer": "Reviewer_WcdN", "summary": "This paper introduces Recurrent Hamiltonian Echo Learning (RHEL), a novel and elegant algorithm for training a class of recurrent models known as Hamiltonian State Space Models (HSSMs). The core idea is to leverage the physical principle of time-reversal symmetry, which is inherent to conservative Hamiltonian systems. Instead of using backpropagation, RHEL computes gradients through a clever, physics-inspired procedure: it runs the model forward, applies a small \"nudge\" to the final state in the direction of the loss gradient, and then lets the system evolve backward in time (an \"echo\"). The deviation of this echo trajectory from the true time-reversed path provides a direct, unbiased estimate of the loss gradients.\n\nThe authors make several key contributions:\n1.  They formalize RHEL first in continuous time, proving its equivalence to the continuous adjoint state method.\n2.  They derive a practical discrete-time version for a specific class of models they call Hamiltonian Recurrent Units (HRUs), proving its equivalence to Backpropagation Through Time (BPTT).\n3.  They generalize the algorithm to deep, stacked models (HSSMs) via a procedure called RHEL chaining.\n4.  They provide compelling empirical evidence that RHEL perfectly matches the gradients computed by BPTT and achieves equivalent training performance on challenging long-range sequence modeling tasks.\n\nUltimately, the work provides a powerful new perspective on credit assignment, framing it as a physical process and paving the way for the design of novel, energy-efficient \"self-learning\" physical hardware.", "questions": "1.  **Extending to Dissipative Systems:** The restriction to conservative Hamiltonian systems is the most significant limitation. Could the authors elaborate on the fundamental challenges of extending this framework to dissipative systems where time-reversal symmetry is broken? Does the echo mechanism break down completely, or could it be adapted, perhaps by having the echo system model the dissipation in reverse (i.e., by actively pumping energy back into the system)?\n2.  **Sensitivity to the Nudging Parameter `ε`:** The algorithm introduces a new hyperparameter, the nudging strength `ε`. Could you provide more intuition on its role and sensitivity? How was it selected in your experiments? Is there a principled way to set it, perhaps by adapting it based on the magnitude of the loss or the learning rate?\n3.  **Path to a \"True\" Physical Implementation:** The vision of a self-learning physical system is compelling. However, Algorithm 1 still requires explicitly computing the gradient of the Hamiltonian with respect to the parameters (`∇_θ H`), which may be non-trivial for a physical system to compute. Could you expand on the feasibility of achieving a truly \"black-box\" physical implementation, perhaps by elaborating on the \"homeostatic control\" idea mentioned in the discussion? What would this look like in practice?", "limitations": "Yes. The authors have done an excellent job of discussing the limitations of their work in the discussion section. They are upfront about the restriction to non-dissipative systems, the offline nature of the algorithm, the computational cost on digital hardware, and the potential bias from the finite nudging strength. Their discussion is honest, thorough, and points clearly to important directions for future research.", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 4, "significance": 3, "originality": 4, "strengths_and_weaknesses": "**Strengths:**\n\n1.  **Exceptional Originality and Significance:** The paper's central contribution is not an incremental improvement but a paradigm-shifting idea. By elegantly connecting the mathematical formalism of backpropagation to the physical principle of time-reversal symmetry, it opens a genuinely new avenue for research in AI. The ultimate vision of building physical systems that learn *in-situ* using their own dynamics is profound and could have a major long-term impact on the field of energy-efficient AI and neuromorphic computing.\n\n2.  **High Technical Quality and Rigor:** The paper is built on a solid theoretical foundation. The authors provide clear, formal theorems (Theorems 3.1, 3.2, 3.3) that establish the equivalence of RHEL with standard gradient-based methods (Adjoint State Method and BPTT). The derivation is logically sound, progressing from the continuous-time idealization to a practical, discrete-time algorithm and finally to deep hierarchical models.\n3.  **Outstanding Clarity and Presentation:** The paper is exceptionally well-written and easy to follow. The core concept of RHEL is introduced with clear intuition, supported by excellent illustrations (Fig. 1 and 2). The logical flow from problem statement to continuous theory, discrete algorithm, and stacked models makes the complex material highly accessible. The authors do a great job of situating their work with respect to prior art like HEB while clearly articulating their novel contributions.\n\n**Weaknesses:**\n\n(These are minor and mostly represent avenues for future work, which the authors acknowledge.)\n\n1.  **Limited Model Class:** The primary limitation is that RHEL is applicable only to conservative Hamiltonian systems that exhibit time-reversal symmetry. This is a restrictive subset of dynamical systems, and many real-world processes (and potentially more expressive models) are dissipative. While the authors suggest workarounds, extending the core theory to non-conservative systems remains an open challenge.\n\n2.  **Practicality on Conventional Hardware:** On current digital hardware (GPUs/TPUs), RHEL requires three \"forward-like\" passes (one free evolution and two parallel echo passes) to compute gradients. This may not offer a significant speed or memory advantage over a standard BPTT implementation (one forward, one backward pass). The true promise of the algorithm lies in its potential implementation on novel analog or physical hardware, a point the authors rightly emphasize.\n\n3.  **Offline Nature:** The current formulation of RHEL is offline, as it requires processing the entire sequence forward before initiating the backward-in-time echo pass. This makes it unsuitable for online or real-time learning scenarios.", "was_revised": true, "final_justification": "I’m maintaining my Strong Accept (6) rating. The authors have done an outstanding job addressing all of my questions, and this is by far the best rebuttal response I’ve seen.", "rebuttal": "We thank Reviewer WcdN for the time they spent to carefully review our submission, raise important questions and we are very happy that they liked our paper. We make **two propositions** of modifications of our paper in case of acceptance along the lines suggested by Reviewer WcdN. \n\n> Extending to Dissipative Systems: [...] fundamental challenges of extending this framework to dissipative systems where time-reversal symmetry is broken? Does the echo mechanism break down completely, or could it be adapted, perhaps by having the echo system model the dissipation in reverse [...]?\n\n**Answer 1.** \nThis is an excellent question! To figure it out, let us model dissipation in the Hamiltonian formalism and follow the exact same steps as in our submission.\n\n- **Modelling dissipation in the Hamiltonian formalism.** One common way to model dissipation in the Hamiltonian formalism is to modify Eq. 1 of into [1]:\n\n$$\n   \\forall t \\in [-T, 0]: \\ \\partial_t \\Phi(t) = (J -R)\\cdot \\nabla_{\\Phi} H[\\Phi(t), \\theta, u(t)],\n$$\n\nwhere $R$ must be *symmetric* and *positive definite*. For simplicity of the derivation, we choose $R=I_\\Phi$ (the identity matrix).\n\n- **Time-reversed dynamics.** With this new model, the time-reversed dynamics, *i.e.* the dynamics of $\\widetilde{\\Phi}^\\star(t):=\\Sigma_z\\cdot\\Phi(-t)$ read:\n\n$$\n    \\partial_t \\widetilde{\\Phi}^\\star(t) =(J + I)\\cdot  \\nabla_{\\Phi^\\star}H[\\widetilde{\\Phi}^\\star(t), \\theta, u(-t)] \n$$\n\nAs the reviewer rightfully noted: time-reversal symmetry is broken since the dissipation ($-R\\cdot \\nabla_{\\Phi} H$) yields, in *reverse*, \"pumping energy back into the system\"($+R\\cdot \\nabla_{\\Phi^\\star} H$). Yet, one could wonder whether RHEL could be extended such that deviations from this time-reversed trajectory, *even if it does not match the forward-time trajectory*, capture the error signals carried by the continuous Adjoint State Method (ASM).\n\n- **ASM with dissipative dynamics**. To answer this question, we must first derive the ASM on the dissipative dynamics introduced above which yields, with $\\lambda(0)=0$:\n\n$$\n\\partial_t \\lambda(t) = \\nabla^2_{\\Phi}H[\\Phi(-t), \\theta, u(-t)] \\cdot J^\\top \\cdot \\lambda(t) + \\nabla_{\\Phi}\\ell[-t, \\Phi(-t), \\theta]\n-\\nabla^2_{\\Phi}H[\\Phi(-t), \\theta, u(-t)] \\cdot \\lambda(t)\n$$\n\nDissipation introduces the new last term in the above equation. We will therefore pay particular attention to this term in the last step.\n\n- **Echo dynamics with pumping.** As suggested by the reviewer, we now define the echo dynamics as, taking again $\\Phi^e(0) = \\Phi^\\star(0)$:\n\n$$\n    \\partial_t \\Phi^e(t) =(J + I)\\cdot  \\nabla_{\\Phi^e}H[\\Phi^e(t), \\theta, u(-t)] - \\epsilon J \\cdot \\nabla_{\\Phi^e}\\ell[-t, \\Phi^e(t)], \n$$\n\nwhere we actively pump energy back into the system. These echo dynamics are designed such that when $\\epsilon=0$, we recover the time-reversed dynamics derived above. Proceeding exactly as we did in the proof of Theorem 3.1 (p. 27 in the Appendix) the derivative of the above echo dynamics with respect to $\\epsilon$ evaluated at $\\epsilon=0$ yields the following equation, denoting $\\Delta_\\Phi(t) := \\Sigma_x \\cdot \\partial_\\epsilon(\\Phi(t,\\epsilon))|_{\\epsilon=0}$:\n\n$$\n\\partial_t \\Delta_\\Phi(t)= \\nabla^2_{\\Phi}H[\\Phi(-t), \\theta, u(-t)] \\cdot J^\\top \\cdot \\Delta_\\Phi(t) + \\nabla_{\\Phi}\\ell[-t, \\Phi(-t), \\theta]\n$$\n\n$$\n-J\\cdot\\nabla^2_{\\Phi}H[\\Phi(-t), \\theta, u(-t)] \\cdot J \\cdot \\Delta_\\Phi(t)\n$$\n\nUnfortunately, we observe that $\\Delta_\\Phi(t)$ no longer satisfies the adjoint dynamics because the last term of the above equation ($-J\\cdot\\nabla^2_{\\Phi}H \\cdot J \\cdot \\Delta_\\Phi$) does not match the last term of the adjoint dynamics ($-\\nabla^2_{\\Phi}H \\cdot \\lambda$).\n\n**Conclusion.** Using our methodology based on the equivalence with the continuous ASM, the echo dynamics cannot be seamlessly adapted to the case where the system is subject to dissipation, even when \"pumping energy back into the system\". However, this does *not* mean that developing forward-only, perturbation-based methods for temporal credit assignment is impossible in general when systems exhibit dissipation. In a very recent work released shortly after our submission, it was shown that a class of *linear* dissipative systems could be described with a Lagrangian formalism (which is closely related to our Hamiltonian formalism) and gradients estimated as finite differences of two *dissipative* trajectories, provided that the system is subject to periodic boundary conditions [2].\n\n\n**Proposition 1**. We propose to mention this important discussion in the main of the paper, and elaborate in greater details with these equations in Appendix. \n\n\n>Sensitivity to the Nudging Parameter ε: [...] Is there a principled way to set it [...]?\n\n**Answer 2.** As reviewer k16k raised the same question, we refer to our comprehensive response in section k16k (Answer 2).\n\n\n>Path to a \"True\" Physical Implementation [...] Algorithm 1 still requires explicitly computing the gradient of the Hamiltonian with respect to the parameters (∇_θ H) [...] Could you expand on the feasibility of [...] the \"homeostatic control\" idea mentioned in the discussion?\n\n**Answer 3.**\nThis is also a very important question. We elaborate below how an *agnostic* version of RHEL (\"aRHEL\") would look like in practice by closely following the methodology proposed in the *agnostic Equilibrium Propagation* paper [3]. While we sketch these ideas in the rebuttal to provide some intuitions, *we acknowledge upfront that they would require a more rigorous treatment in a separate paper*. \n\n- **Definition of the joint neuron and parameter dynamics**. We define the *joint* dynamics of $\\theta$ and $\\Phi$ over $t\\in[-T, 0]$ as:\n\n$$\n\\partial_t \\Phi(t) = J\\cdot \\nabla_{\\Phi} H[\\Phi(t), \\theta(t), u(t)],\n$$\n$$\n   \\partial_t \\theta(t) = - \\alpha  \\nabla_{\\theta}\\left[H[\\Phi(t), \\theta(t), u(t)] + \\frac{1}{2} \\|c(t) - \\theta(t)\\|^2\\right]\n$$\n\nwhere $c(t)$ is a *control* variable acting on the parameters and $\\alpha$ a learning rate. Note that: i) the dynamics of $\\theta$ are purely dissipative, ii) we do not model the *momentum* of $\\theta$ but only its position. We acknowledge that this modeling choice is far from obvious and could be subject to further discussions, but we adopt this one as this is *the most simple choice* to readily adapt the agnostic mechanism from [3] to RHEL.\n\n- **First phase of the agnostic procedure.** Assuming that $\\theta(-T):=\\theta_\\star$, we adjust the control variable $c$ during the free phase such that $\\theta$ remains at $\\theta_\\star$ - for instance using a simple proportional controller as in [3]. Namely, $\\forall t \\in [-T, 0]$:\n\n$$\n \\partial_t \\Phi(t) = J\\cdot \\nabla_{\\Phi} H[\\Phi(t), \\theta(t), u(t)],\n$$\n\n$$\n   \\partial_t \\theta(t) = - \\alpha  \\nabla_{\\theta}\\left[H[\\Phi(t), \\theta(t), u(t)] + \\frac{1}{2} \\|c(t) - \\theta(t)\\|^2\\right] \\approx 0 \\\\\n$$\n\n$$\n  \\Rightarrow \\forall t \\in [-T, 0]: \\ c(t)\\approx \\theta_\\star + \\nabla_{\\theta}H[\\Phi(t), \\theta_\\star, u(t)]  \n$$\n\nWe denote $\\{c_\\star(t)\\}_{t\\in[-T, 0]}$ the resulting control trajectory.\n\n- **Second phase of the agnostic procedure.** During the echo phase, the control trajectory $c_\\star$ is played in reverse, yield the following dynamics $\\forall t \\in [0, T]$:\n\n$$\n\\partial_t \\Phi^e(t) = J\\cdot \\nabla_{\\Phi^e} H[\\Phi^e(t), \\theta^e(t), u(-t)] - \\epsilon J\\nabla_{\\Phi^e}\\ell[-t, \\Phi^e(t)], \\\\\n$$\n\n$$\n\\partial_t \\theta^e(t) = - \\alpha  \\nabla_{\\theta}\\left[H[\\Phi^e(t), \\theta^e(t), u(-t)] + \\frac{1}{2} \\|c_\\star(-t) - \\theta^e(t)\\|^2\\right] \n$$\n\nTaking a closer look at the dynamics satisfied by $\\theta^e$ and using the definition of $c_0$ above:\n\n$$\n\\partial_t \\theta^e(t) = -\\alpha \\left[\\nabla_{\\theta}H[\\Phi^e(t), \\theta^e(t), u(-t)] - c_\\star(-t) + \\theta^e(t)\\right]\n$$\n\n$$\n\\approx -\\alpha \\left[\\nabla_{\\theta}H[\\Phi^e(t), \\theta^e(t), u(-t)] - \\nabla_{\\theta}H[\\Phi(-t), \\theta_\\star, u(-t)] - \\theta_\\star + \\theta^e(t)\\right]\n$$\n\n$$\n= \\alpha \\epsilon \\Delta_{\\theta}^{\\rm aRHEL}(t, \\epsilon, \\alpha) - \\alpha \\nabla_{\\theta}\\left[\\frac{1}{2} \\|\\theta^e(t) - \\theta_\\star \\|^2\\right]\n$$\n\nwith $\\Delta_{\\theta}^{\\rm aRHEL}(t, \\epsilon, \\alpha):=\\epsilon^{-1}(\\nabla_{\\theta}H[\\Phi^e(t), \\theta^e(t), u(-t)] - \\nabla_{\\theta}H[\\Phi(-t), \\theta_\\star, u(-t)])$. Note that $\\Delta_{\\theta}^{\\rm aRHEL}$ would recover $\\Delta_{\\theta}^{\\rm RHEL}:=\\epsilon^{-1}(\\nabla_{\\theta}H[\\Phi^e(t), \\theta_\\star, u(-t)] - \\nabla_{\\theta}H[\\Phi(-t), \\theta_\\star, u(-t)])$ if $\\theta^e(t)$ was \"close enough\" to $\\theta_\\star$, which the second term of the above dynamics naturally biases towards.\nTherefore, the proposed agnostic mechanism approximates the desired weight update $\\Delta \\theta^{\\rm RHEL} = \\lim_{\\epsilon, \\alpha \\to 0}\\int_0^T \\Delta_{\\theta}^{\\rm aRHEL}(t, \\epsilon, \\alpha)$, without requiring the exact knowledge of $H$ to compute the learning rule.\n\n**Proposition 2**. We propose to include this analysis in the appendix of the paper.\n\n[1] Smith et al. \"Learning dissipative Hamiltonian dynamics with reproducing kernel Hilbert spaces and random Fourier features.\" 2024.\n\n[2] Berneman et al. \"Equilibrium Propagation for Periodic Dynamics.\" 2025.\n\n[3] Scellier et al. \"Agnostic physics-driven deep learning.\" 2022."}, {"review_id": "tlJigq19b9", "reviewer": "Reviewer_k16k", "summary": "This paper extends and generalizes Hamiltonian Echo Backprop (HEB) into the proposed Recurrent Hamiltonian Echo Learning (RHEL). The method is shown to be formally equivalent to the continuous-time adjoint state method (ASM) and Backpropagation Through Time (BPTT) in discrete time. The authors introduce Hamiltonian Recurrent Units (HRUs) and stack them into Hamiltonian State Space Models (HSSMs) to build deep sequence models that are compatible with RHEL.", "questions": "- Can you discuss the computational cost of running RHEL vs BPTT in a conventional setting?\n- Is there a recommended heuristic for setting epsilon in practice?", "limitations": "yes.", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strengths:\n-This paper is well written with precise exposition of the RHEL method and supporting results.\n-The experiments on long-range tasks demonstrate that RHEL matches the performance of BPTT while being useful for continuous learning or alternate hardware implementations, with potential memory and computation savings.\n- The theoretical equivalences to both BPTT and ASM are well developed in the limit of small perturbations and supported with clear empirical results.\n\nWeaknesses:\n- As presented this method is only applicable to reversible Hamiltonian  systems which can be a restrictive class, and relies on symplectic \nintegration.\n- The theoretical equivalence between RHEL and BPTT and ASM is  achieved as the perturbation strength, epsilon, goes to zero. However accurate gradient computations rely on epsilon > 0, particularly in lower precision, which may result in practical implementation being sensitive. It would be useful to see some consideration of continuity / similarities for a non zero perturbation.\n-The claim of eliminating / reducing exploding and vanishing gradients would be strengthened by a supporting result. There is no \ndemonstrated improvement over BPTT (which is expected and fine) but could be the case when BPTT experience gradient maladies.\n-Some of the phrasing in the introduction and conclusion, motivating this method through connections with hardware advancements, is a \nlittle odd and does not seem to fully compliment the method begin presented.\n\nSome of the language, particularly in sentence one of the abstract, could be toned down. Other examples: \"\"This paper takes an extreme view in the world of ai hardware...\". The clarity score would be higher if these were addressed.", "was_revised": false, "final_justification": "", "rebuttal": "We thank Reviewer k16k for their thoughtful evaluation. We appreciate their constructive feedback on several important aspects of RHEL. Below we address their specific concerns and provide detailed empirical analysis to strengthen our claims as well as **two proposals** to incorporate in our paper in case of acceptance. \n\n>As presented this method is only applicable to reversible Hamiltonian systems [...].\n\n**Answer 1.** \nWhile we fully acknowledge that extending RHEL to dissipative systems represents an important avenue for future work, *the restriction to Hamiltonian systems is less limiting than it may seem at first glance*.\n\nFrom a physical viewpoint, strictly reversible dynamics can be relaxed into *a suitable interplay between the time constants involved*. Namely: the timescale of dissipation within the self-learning device must be much larger than the timescale of interaction between the input and self-learning system. This timescale separation allows the system to behave *approximately* as a reversible Hamiltonian system during the computation window. As elaborated in the original HEB work [1], many existing physical systems can satisfy this condition: optical systems with Kerr nonlinearity, cold atom clouds in optical lattices, spin wave systems, and superconducting circuits. The key requirement is that dissipation timescales are much longer than the computation timescales, which is achievable in many well-isolated physical systems.\n\nAlternatively, we also investigate in details the possibility of extending RHEL to dissipative systems in our answer to Reviewer WcdN (Answer 1, Proposition 1).\n\n**Proposition 1.** We propose to incorporate the above discussion on the physical feasibility of reversibility constraints to implement RHEL on real systems.\n\n>The theoretical equivalence between RHEL and BPTT and ASM is achieved as [...] epsilon, goes to zero. However accurate gradient computations rely on epsilon > 0, particularly in lower precision [...] Is there a recommended heuristic for setting epsilon in practice?\n\n**Answer 2.** Both reviewers k16k and WcdN correctly identified that our theoretical results assume infinitesimally small nudging strength $\\epsilon$, which requires careful consideration in finite precision implementations. We conducted extensive empirical analysis in **Appendix A.4.4 to inform our selection of $\\epsilon$**, revealing a fundamental trade-off:\n\n*Numerical precision challenges*. In practice, gradient information in RHEL is encoded in small perturbations to the state $\\Phi$ that generate the echo trajectory $\\Phi^e$ (Theorem 3.1 of our submission). When perturbations and state values differ by several orders of magnitude, finite precision arithmetic suffers from *numerical underflow* issues whereby small perturbation values are flushed to zero. To address this, we employ a *gradient rescaling method* where we multiplicatively scale the output loss $\\ell[\\cdot,\\cdot]$ by a constant $\\gamma > 1$ during echo dynamics computation. Note that: i) it is equivalent to scaling the nudging strength, ii) gradient rescaling has become standard and integrated inside machine learning libraries. However, some subtleties arise with RHEL. \n\n*Accurate gradient estimation trade-offs*: \n- **Linear case**: We observed that increasing the scaling factor consistently improves gradient matching between RHEL and BPTT. This is because linear systems maintain linear response regardless of $\\epsilon$ magnitude, so larger scaling only helps with numerical precision without introducing approximation errors. Note that vanilla backprop, and more generally the adjoint method, fall in this category, as error propagation is *linear* with respect to the adjoint variables.\n- **Nonlinear case**: There exists an optimal scaling factor that balances two competing error sources:\n  - *Underflow errors* ($\\epsilon$ too small): Perturbations lost to finite precision\n  - *Linearization errors* ($\\epsilon$ too large): RHEL uses finite difference approximation $\\nabla f \\approx [f(x + \\epsilon) - f(x - \\epsilon)]/(2\\epsilon)$, which assumes linear response to perturbations. In nonlinear systems, large nudging causes the response to deviate significantly from linear behavior, making the finite difference approximation inaccurate. A solution to this (which we did not explore in our initial submission) would be to use *higher order* estimators, for instance: $\\nabla f \\approx (1/12\\epsilon)[-f(x+2\\epsilon) + 8 f(x + \\epsilon) - 8 f(x - \\epsilon) + f(x - 2\\epsilon)]$ \n\n*Implementation*: Through grid search, we found optimal scaling factors for both the Linear and Nonlinear case. \n\nThis analysis demonstrates that while theoretical equivalence requires $\\epsilon \\to 0$, practical implementation with appropriate scaling maintains gradient accuracy while avoiding both numerical instabilities and linearization errors.\n\n> -The claim of eliminating / reducing exploding and vanishing gradients would be strengthened by a supporting result\n\n**Answer 3**.\n*Clarification of our claims*. We understand the need for such evidence and realize that our wording in the initial submission was not clear enough. When we wrote that \"it has been shown that these HRUs mitigate by design vanishing and exploding gradients\" (L.248-249) we meant that: i) this claim has been demonstrated in a *prior work* [2], as such it is not one of ours, ii) this claim is a property of the *model*, not of the gradient computation algorithms, therefore it *equally* applies to BPTT and RHEL. This being said, we appreciate the value of providing details about the theoretical grounding of these claims as well an empirical validation thereof for self-containedness of our work. \n\n*Theoretical guarantees*. The gradient stability properties stem from the underlying Hamiltonian structure with symplectic integration. Prior work [2] rigorously demonstrated that the recurrent component of our nonlinear HRU block (corresponding to the Hamiltonian in Equation 18) possesses:\n- **Non-exploding gradients** (Proposition 3.1 in [2]): the *total* loss gradient (*i.e.* aggregating loss \"sensitivities\" across all timesteps) is upper bounded. \n- **Non-vanishing gradients** (Proposition 3.2 in [2]): when looking far enough back in time, the loss \"sensitivities\" (*i.e* the additive time-wise contributions to the loss gradient) are strictly non-zero and independent of the timestep. \n\n*Empirical validation*: Since RHEL provides accurate BPTT gradient approximations (Figure 4 of our submission), it must also inherit these stability properties. To demonstrate this empirically, we analyzed gradient magnitudes over time for the recurrence matrix $A$ of the nonlinear HRUs.\n\n*Experimental setup*: We randomly sampled $(u_i, y_i) \\sim \\mathcal{D}$ from the SCP1 dataset and computed time-dependent sensitivities for both RHEL ($\\Delta_A^{\\text{RHEL}}(k)$) and BPTT ($g_A(k)$) across the sequence length. We then performed linear regression analysis:\n\n\n| Analysis | Correlation \\(r\\) | Slope | Interpretation |\n|-----------|----------------|-------|----------------|\n| $\\|\\Delta_A^{\\text{RHEL}}(k)\\|$ vs $\\|g_A(k)\\|$ | 0.99998 | 1.0007 | RHEL ≈ BPTT gradients |\n| $\\|\\Delta_A^{\\text{RHEL}}(k)\\|$ vs $k$ | 0.11068 | 2.82×10⁻¹⁵ | No vanishing/exploding gradient |\n| $\\|g_A(k)\\|$ vs $k$ | 0.11680 | 2.98×10⁻¹⁵ | No vanishing/exploding gradient |\n\n**Key findings**: \n- **Row 1** confirms RHEL's gradient fidelity (r ≈ 1, slope ≈ 1)\n- **Rows 2-3** demonstrate gradient stability: slopes ≈ 0 indicate neither exponential growth (exploding) nor decay (vanishing) over time.\n\nThis empirical evidence supports our theoretical claims, those of [2] and demonstrates that RHEL maintains the beneficial gradient properties of the underlying Hamiltonian architecture.\n\n> There is no demonstrated improvement over BPTT (which is expected and fine) but could be the case when BPTT experience gradient maladies.\n\nThis is an interesting suggestion. However, we do not expect RHEL to improve over BPTT regarding gradient pathologies, as these are inherent properties *of the model architecture* rather than the learning algorithm as mentioned earlier. At the single HRU level, both methods inherit well-behaved gradients (non-exploding/non-vanishing) from the Hamiltonian dynamics. While vanishing gradients can occur *across* layers in deep architectures (i.e. inside feedforward transformations binding HRUs), both RHEL and BPTT face identical challenges since they compute mathematically equivalent gradients. Both methods benefit equally from standard architectural solutions like skip connections -- see Appendix C.4-C.5 in [2].\n\n**Proposion 2.** We propose to add the table presented in answer 2 to validate the non-vanishing/exploding gradient properties of the nonlinear HRU.\n\n> Some of the phrasing in the introduction and conclusion [...] is a little odd [...] Some of the language, particularly in sentence one of the abstract, could be toned down.\n\nWe fully acknowledge this point. As phrased in details in our answer to Reviewer AkMK (Answers 2-3, Propositions 2-3), we will eliminate these sentences and replace them with a more detailed background on physical computing and contrast this literature with standard digital design with a softened tone. We will also include detailed discussions about possible physical implementations of RHEL (Reviewer UVRA, Answer 1, Proposition 2) as well as a more physically realistic version of RHEL which does not require the exact knowledge of the Hamiltonian to perform the weight update (Reviewer WcdN, Answer 3, Proposition 2). \n\n> Can you discuss the computational cost of running RHEL vs BPTT in a conventional setting?\n\nWe carry out this comparison in details in our answer to Reviewer AkMK (Answer 1, Proposition 1).\n\n[1] Lopez-Pastor & Marquardt. \"Self-learning machines based on Hamiltonian echo backpropagation. 2023.\n\n[2] Rusch & Mishra. Unicornn: A recurrent model for learning very long time dependencies. 2021."}]

Author rebuttal:

Meta-review: This paper introduces a new approach for training recurrent models without relying on Backpropagation Through Time (BPTT). Specifically, it proposes Recurrent Hamiltonian Echo Learning (RHEL), an algorithm that computes loss gradients as finite differences of trajectories of non-dissipative Hamiltonian systems.

A key downside of the proposed approach is that RHEL requires the underlying system to be non-dissipative. This is a somewhat restrictive assumption that can limit applicability and hurt performance. As a result, models trained with BPTT currently tend to outperform RHEL in terms of accuracy. Nevertheless, the fact that only three forward passes are required for learning is an exciting property. As mentioned by the authors, its full potential will likely depend on the development of improved analog hardware, and it remains to be seen whether such hardware will become available. Despite these limitations, the work provides a fresh perspective and opens an exciting new research direction motivated by dynamical systems theory.

The reviews are overall very positive and confirm that the underlying mathematical framework is sound. Open questions and concerns have been effectively addressed during the rebuttal and discussion phases. I agree with some of the comments that the work remains somewhat limited in its current form, but this does not diminish its potential impact. I would encourage the authors to more explicitly acknowledge current limitations and to expand the discussion of promising directions for future work. Additionally, I suggest moderating the use of strong adverbs (e.g., “radically,” “tremendous”) to maintain a more neutral academic tone.

In light of these considerations, I strongly recommend this paper for acceptance.
Decision: Accept (oral)