© 2026 IEEE. Personal use of this material is permitted. Permission from IEEE must be obtained for all other uses, in any current or future media, including reprinting/republishing this material for advertising or promotional purposes, creating new collective works, for resale or redistribution to servers or lists, or reuse of any copyrighted component of this work in other works. Accepted to the 2026 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS 2026).
Statistical Learning of Contractive Dynamical Representations for Composite Adaptive Control
arXiv:2609.35758v1 [eess.SY] 28 Sep 2026
Min Kim1,∗ , José Leonardo Brenes1,∗ , Fred Hadaegh1,2 , Soon-Jo Chung1,2 Abstract— We present a representation-learning framework for composite adaptive tracking control under dynamically coupled disturbances. The framework connects classical disturbance-accommodating control (DAC) to recent last-layer adaptive disturbance-rejection methods. Specifically, we introduce a statistically principled hard expectation–maximization (hard-EM) procedure, with a Kalman smoother in the hard E-step, to identify dynamical representations of disturbance whose latent evolution is uniformly contractive. The learned representation evolves a latent disturbance-excitation state from measured plant features and control inputs and decodes that state into the time-varying disturbance acting on the nominal plant, thereby extending prior “fixed-decay” last-layer adaptive methods to a learned, predictive DAC-style formulation. Combined with Bayesian filtering of the learned latent state, this representation yields a composite adaptive tracking controller with predictive capability and provable exponential convergence to a bounded neighborhood. We validate our approach experimentally on a slippery ground vehicle carrying a liquidsloshing tank and a pendulum load, and we further assess its robustness on a system of coupled Duffing oscillators. Across both settings, the method achieves accurate disturbance prediction and improved overall tracking performance relative to fixed-decay representation-learning ablations, LTI disturbanceaccommodating baselines, and model-based PD baselines.
I. INTRODUCTION There is a substantial body of literature on representation learning for last-layer adaptive control and disturbance rejection [1]–[3]; however, many existing applications focus on external disturbances such as aerodynamic effects or terraininduced slip. In such settings, disturbances are captured by a latent variable that is treated as exogenous; its value shifts as the “environment” changes, sometimes with an assumed decay rate, which corresponds to the Ornstein– Uhlenbeck process. In contrast, when the unknown disturbance is dynamically coupled with the controlled plant, such as liquid sloshing in a vehicle’s tank, it is natural to view the latent variable as an augmented dynamical state that evolves according to its own dynamics, in the spirit of disturbance-accommodating control (DAC). In such cases, we expect adaptive controllers with predictive capability to have a comparative advantage over counterparts that do not incorporate prediction. This motivates the learning of dynamical representations tailored to adaptive control; in this paper, we develop a representation that captures the coupled disturbance source’s excitation mechanism and its influence on the main system. * These authors contributed equally. 1 California Institute of Technology (Caltech), Pasadena, CA 91125 USA. 2 Jet Propulsion Laboratory (JPL), Pasadena, CA 91109 USA.
{mink, jbrenes, hadaegh, sjchung}@caltech.edu
Fig. 1. Proposed adaptive control architecture and experimental platform. Left: block diagram of the controller. A latent state is estimated online via a Kalman–Bucy-style filter from the disturbance proxy y, using a latentstate evolution model parameterized by (A, b). The basis functions (Φ, ψ) map latent-relevant features ϕ to an adaptive disturbance-rejection term that augments a nominal tracking controller. The quantities (A, b, Φ, ψ) are parameterized by neural networks and learned using a statistically principled training procedure. Right: tracked mobile robot used for experimental validation.
The contributions of this paper are threefold: (i) we introduce a statistically principled hard expectation–maximization training framework, with Kalman-smoother-based latent inference (Fig. 2), for learning contractive dynamical latent representations of disturbance from data; (ii) we design a neural architecture for the dynamical representation that enforces uniform contraction, is compatible with Bayesian filtering, and yields predictive disturbance estimates with stability guarantees when used in a composite adaptive controller; and (iii) we demonstrate the method on hardware and in simulation, showing accurate disturbance prediction and improved tracking over fixed-decay representation-learning ablations and LTI-based DAC. Taken together, these results show that statistically principled learning of contractive latent disturbance dynamics can bridge representation learning and composite adaptive control for systems with dynamically coupled disturbances. We next position our paper relative to prior work. From the perspective of classical control, our approach is closely related to disturbance-accommodating control (DAC), which represents disturbance using an augmented disturbance state with its own dynamics and estimates this augmented state online for compensation [4]. DAC motivates treating the disturbance as an estimable and predictable dynamical quantity rather than as a memoryless exogenous input. Our perspective is also connected to the robust servomechanism problem
Fig. 2. Statistically principled training via hard-EM with exact latent inference. The procedure alternates between a hard E-step and an M-step: given current parameters (θ, Rd , Qd ), the hard E-step infers the latent trajectory a0:L−1 for each data window, and the M-step updates the noise covariances (Rd , Qd ) from residual statistics and updates the network parameters θ.
[5], where disturbance signals are generated by an exosystem and robust output regulation requires embedding an internal model [6] of the disturbance in the closed loop. From this standpoint, our method can be interpreted as learning a datadriven, contractive analogue of the disturbance generator in the spirit of identification-based internal-model regulation [7]. In contrast to [7], which adapts an internal model by linear regression, we identify an affine-in-latent model via statistical learning and exploit Kalman-based estimation. Learning-based control has been demonstrated on physical robot platforms, including quadrotors under wind disturbance [1], [8], drone racing [9], learning-based MPC for quadrotors [10] and autonomous car racing [11], rapid motor adaptation for bipedal locomotion [12], and robust perceptive locomotion for quadrupedal robots in natural environments [13]. These hardware demonstrations highlight the importance of controllers that remain reliable under modeling errors and disturbances, which are often induced or amplified by robotic maneuvers. Motivated by this need, recent learningbased MPC approaches improve tracking under model mismatch by incorporating learned disturbance or residuals into the MPC framework [14], [15]. Also, [16] learns a Koopman-based disturbance dynamics model and uses the corresponding disturbance observer for compensation. In contrast, we learn a latent disturbance dynamics model by a Kalman-smoother-based latent trajectory inference and hard expectation–maximization (EM). We note that the learning strategy presented in this paper is analogous to previous smoother-based EM methods for system identification [17], [18], although these papers do not involve neural networks. In this paper, we focus on disturbance representations that are uniformly contractive [19], [20] and are affine in the latent variable (see (3)). The contraction property of the dynamical representation is required for insensitivity to the unknown initial latent state (Proposition 1), and is physically plausible for dissipative disturbance dynamics. On the other hand, the affine-in-latent form enables closed-form Bayesian estimation. This motivates us to use a hard-EM procedure to train the neural networks (Fig. 2 and Section III); during training, the latent trajectory is repeatedly estimated by Kalman smoothing. Moreover, using the Kalman–Bucy filter for the latent state, we synthesize a composite adaptive con-
troller [1], [2], [21], whose exponential tracking guarantees are analyzed in Theorems 3 and 4. This analysis relies on a covariance bound for uniformly contractive systems (Lemma 2). In Section IV, we evaluate the learned dynamical representation and the resulting adaptive controller on two different scenarios. First, we study a ground vehicle (“GVR-Bot”) equipped with a liquid-sloshing tank, a pendulum load, and an unknown, inaccessible internal PID controller; vehicle maneuvers excite both the fluid sloshing and the pendulum motion. Then, we consider a system of coupled Duffing oscillators as a numerical testbed, where the robustness of the proposed tracking controller is assessed under increasingly adverse Gaussian-noise conditions. These examples show not only that the method can learn predictive maneuver-todisturbance mappings, but also that the learned mappings can be readily and effectively incorporated into composite adaptive controllers. II. THE CONTROL METHOD AND GUARANTEES The main controller equations are given by the theorems in Section II-B. A. System Description and the Dynamical Representation We consider systems with additive disturbances of the form: M (q)q̈ + C(q, q̇)q̇ + g(q) = u + d (1) where the coefficient functions are continuously differentiable and Ṁ − 2C is skew-symmetric. We also consider the simpler first-order system: v̇ = f (v) + u + d,
(2)
with continuously differentiable f . We aim to learn disturbance estimators of the form (3) for use in composite adaptive control [22]: ȧ = A(ϕ(t), u(t))a + b(ϕ(t), u(t)), dˆ = Φ(ϕ(t))a + ψ(ϕ(t)),
(3a) (3b)
where A : Rdϕ +du → Rda ×da , b : Rdϕ +du → Rda , Φ : Rdϕ → Rdu ×da , ψ : Rdϕ → Rdu are continuous functions to be learned. Intuitively, a ∈ Rda summarizes disturbance-relevant hidden dynamics; A and b define its evolution, while Φ and ψ decode it into a disturbance estimate. We emphasize that ϕ : R≥0 → Rdϕ is a continuous “feature” signal which is determined causally. In most cases ϕ is simply the state: ϕ = (q, q̇), but it could additionally contain measurements from on-board sensors, e.g., a camera. Here, a statistically explains the disturbance d and evolves according to the dynamic model (3a). Unless otherwise stated, the control signals u : R≥0 → Rdu are continuous and ∥ · ∥, when applied to a matrix, denotes the spectral norm. ˆ does not depend on Note that in (3) the estimate d(t) the instantaneous control input u(t). This restriction may be
undesirable in some settings. An alternative is to linearly approximate the map u(t) 7→ d(t); to do so, set dˆk (t) = uaug (t)⊤ Φk a + ψk (4)
Moreover, uniformly contractive linear systems lead to bounded covariance estimates when the Kalman–Bucy filter is used. Specifically, we consider the following covariance evolution equation:
Ṗ = A(t)P + P A(t)⊤ + Q − P Φ(t)⊤ R−1 Φ(t)P (6) for each disturbance coordinate 1 ≤ k ≤ dim(d), where ⊤ ⊤ the augmented input uaug (t) := [u(t) , 1] . This paramewith a fixed initial condition P (0) = P0 ≻ 0 and Q ⪰ 0, terization is convenient for controller synthesis: since u(t) R ≻ 0. Here, A(t) := A(ϕ(t), u(t)) and Φ(t) := Φ(ϕ(t)). ˆ enters (4) affinely, the relation u(t) + d(t) = udesired (t) Lemma 2: Suppose 21 A + A⊤ ⪯ −λI for some concan be easily solved for u(t); in practice, one solves the stant λ > 0. Then the unique solution of (6) is positive Tikhonov-regularized normal equations centered at udesired (t) definite and satisfies for robustness. Finally, we also note that a prototypical class λmax (Q) of examples for (3) consists of multi-component mechanical 1 − e−2λt (7) λmax P (t) ≤ e−2λt λmax (P0 ) + 2λ systems coupled through force elements. Hereafter, we assume access to an online noisy estimate for all t ≥ 0. Moreover, if Φ is bounded, Q ≻ 0, and 21 (A + y(t) = d(t) + ϵ(t), where ϵ(t) is a continuous, additive error. A⊤ ) ⪰ −µI, then λmin (P (t)) ≥ p∗ for some positive p∗ If y(t) were exact and the control loop were sufficiently given by an explicit formula. fast relative to the disturbance timescale, direct cancellation Proof: The positive definiteness claim follows from the would largely suffice; in practice, however, y(t) is often comparison theorem for Hermitian RDE [23, Theorem 4.1.4] obtained by numerically differentiating the measured q or (compare to the same equation but with Q = 0), and v, which amplifies sensor noise, and actuator delays further existence of P on [0, ∞) follows from Theorem 4.1.6. limit the feasibility of instantaneous cancellation for rapidly For the upper bound, consider the linear equation P̃˙ = varying disturbances. AP̃ + P̃ A⊤ + Q with the initial condition P̃ (0) = P0 . The proposed learned dynamical representation instead Theorem 4.1.4 gives P (t) ⪯ P̃ (t) for all t ≥ 0. Since R produces systematically refined estimates of the current dis- P̃ (t) = ΦA (t, 0)P0 ΦA (t, 0)⊤ + t ΦA (t, τ )QΦA (t, τ )⊤ dτ 0 turbance d(t) via the Kalman–Bucy filter, by fusing the full for all t ≥ 0, and ∥ΦA (t, t0 )∥ ≤ e−λ(t−t0 ) for all t ≥ measurement history y(·) and, implicitly, all measurements t0 ≥ 0, we get λmax (P̃ (t)) = ∥P̃ (t)∥ ≤ e−2λt λmax (P0 ) + λmax (Q) available in the training set. When the disturbance source is (1 − e−2λt ). The claim follows from the Weyl 2λ dynamically coupled to the plant, learning an explicit latent monotonicity principle. dynamics is a natural design choice; these learned dynamics For p∗ = √ the lower bound, we have P (t) ⪰ p∗ I with 2 then serve as predictive process models in the Kalman–Bucy −µ+ µ2 +su λmin (Q) Φ ∧ λmin (P0 ), where su := λmin (R) and su filter. Moreover, when constructing disturbance labels for Φ > 0 is an upper bound for Φ; this can be proved by training, one may leverage higher-quality off-board sensing comparing the original equation to a related one: Ṗc = APc + and non-causal signal processing to obtain more reliable ⊤ 2 ⊤ P A + (s p I − p A c u ∗ ∗ − p∗ A ) − Pc (su I)Pc , Pc (0) = p∗ I numerical derivatives than would be available on-board in ⊤ ⊤ −1 2 real time. The refined estimates from the filter are then and noting that Q ⪰ su p∗ I −p∗ A−p∗ A , Φ R Φ ⪯ su I. used for disturbance compensation, thereby yielding composite adaptive tracking controllers with provable convergence B. Composite Adaptive Tracking Controller and Analysis guarantees (see Theorems 3 and 4). While Lemma 2 is of independent interest, its main role Before proving the main Theorems 3 and 4 which give in this paper is to establish exponential tracking results for the controllers ((9) or (12) with (8)), we first formalize various composite adaptive control strategies based on the the intuition that uniformly contractive linear systems are Kalman–Bucy filter. insensitive to initial conditions which will be unknown in First, we consider a tracking problem for (1) with a twice practice. continuously differentiable reference trajectory qd . We use Proposition 1 (Exponential insensitivity to initial condition): Kalman–Bucy-inspired composite adaptation equations of Suppose 12 A + A⊤ ⪯ −λI and ∥Φ∥ ≤ C pointwise for the form: some constants λ > 0 and C > 0. For a fixed pair ϕ and u, â˙ = A(ϕ, u)â + b(ϕ, u) any two trajectories a1 and a2 of (3) satisfy (8a) + P Φ(ϕ)⊤ R−1 (y − Φ(ϕ)â − ψ(ϕ)) + P Φ(ϕ)⊤ s, −λt ˆ ˆ ∥d1 (t) − d2 (t)∥ ≤ Ce ∥a1 (0) − a2 (0)∥ (5) Ṗ = AP + P A⊤ + Q − P Φ⊤ R−1 ΦP (8b) for all t ≥ 0. ˙ where s := q̃+Λq̃, q̃ := q−qd , and P (0) = P0 , Q, R, Λ ≻ 0. 2 d Proof: We have dt ∥a1 (t) − a2 (t)∥ = For the main theorems, we assume d(t) = Φ(ϕ(t))a(t) + ⊤
⊤
2 (a1 (t) − a2 (t)) (ȧ1 (t) − ȧ2 (t)) = (a1 (t) − a2 (t)) (A+ 2 A⊤ ) (a1 (t) − a2 (t)) ≤ −2λ ∥a1 (t) − a2 (t)∥ . Therefore, −λt ∥a1 (t) − a2 (t)∥ ≤ e ∥a1 (0) − a2 (0)∥. Now, the claimed result follows from dˆ = Φa + ψ.
ψ(ϕ(t)) + r(t), so that the disturbance d is realized by continuous signals a(t) and r(t) with a being differentiable. Note that (8) are the standard Kalman–Bucy equations but with P Φ⊤ s, an additional composite adaptation term. Note
√ that when the Kalman covariance P becomes larger, the Differentiating Gc (t) := eαt Wc − c − d/(c∗ α) and tracking error term P Φ(ϕ)⊤ s also becomes larger (along passing to the limit c → 0+ gives: with the innovation term). p p d d √ −αt Remark 1: In (8), A = A(ϕ, u), b = b(ϕ, u), Φ = Φ(ϕ), m∥s(t)∥ ≤ V (t) ≤ e V (0) − ∗ + ∗ . c α c α ψ = ψ(ϕ) are dependent on the measured features and/or (11) control, and therefore are implicitly time-varying. Previous ˙ = −Λq̃ + s in q̃, using By solving the linear system q̃ last-layer adaptive disturbance rejection methods [1], [2] may be interpreted as the case of A = −λI, b = 0, ψ = 0 with ∥ exp(−tΛ)∥ = exp (−tλmin (Λ)) for all t ≥ 0, and choosing −c2 t + hyperparameter λ > 0, although their learning procedures α < λmin (Λ) if necessary, we conclude that ∥q̃∥ ≤ c1 e c for some c > 0 with c , c being independent of â(0), 3 i 2 3 differ from the procedure used here. Remark 2: The Kalman-smoother approach to representa- a(0), q(0), and q̇(0). We also present a simpler analogue of the previous theotion learning (Section III-C) provides data-driven estimates rem in the context of (2). of Q and R that serve as principled starting points for Theorem 4 (Exp. Tracking: 1st-order systems): Consider subsequent gain tuning. the system (2), (8) with the feedback controller Now we present the main theorems demonstrating exponential tracking convergence via the adaptation dynamics (8); u = v̇d − f (v) − Ks − Φâ − ψ, (12) here, exponential convergence to a bounded ball means that there exist constants ci > 0 such that each solution on R≥0 of where K ≻ 0 and now s := v − vd with a continuously dif the closed-loop system satisfies ∥q(t)−qd (t)∥ ≤ c1 e−c2 t +c3 ferentiable reference velocity vd . If −µI ⪯ 1 A + A⊤ ⪯ 2 −c2 t (or ∥v(t) − vd (t)∥ ≤ c1 e + c3 for (2)), where c2 and c3 −λI for some positive constants µ, λ, and if (Aa + b − ȧ), are independent of â(0), a(0), and q(0), q̇(0) (or v(0) for Φ, r, ϵ are bounded by some constants independent of (2)). As expected, the theorems will depend on the bounds â(0), a(0), v(0), then the feedback controller (12) achieves on Aa + b − ȧ and r, i.e., the quality of the dynamical exponential convergence to a bounded ball. representation (3). Proof: We have ṡ = −Ks − Φã + r where ã := Theorem 3 (Exp. Tracking: 2nd-order mechanical systems): â − a. Considering V := s⊤ s + ã⊤ P −1 ã, we again have Consider the system (1), (8) with the feedback controller V̇ = −2s⊤ Ks − ã⊤ P −1 QP −1 + Φ⊤ R−1 Φ ã + 2s⊤ r + 2ã⊤ P −1 (Aa + b − ȧ) + 2ã⊤ Φ⊤ R−1 (r + ϵ). The rest of the u = M v̇r + Cvr + g − Ks − Φâ − ψ, (9) proof is similar to the proof of Theorem 3 once one notes where K ≻ 0 and vr := q̇d − Λq̃. If −µI ⪯ 21 A + A⊤ ⪯ that we have K ⪰ αI and Q ⪰ 2αP for all small α > 0, −λI, mI ⪯ M ⪯ mI for some positive constants µ, λ, which leads to m, m, and if (Aa + b − ȧ), Φ, r, ϵ are bounded by some p d d ∥s(t)∥ ≤ ∗ + V (0) − ∗ e−αt , (13) constants independent of â(0), a(0), q(0), q̇(0), then the c α c α feedback controller (9) achieves exponential convergence to 1/2 a bounded ball. −1 ∗ and d is dewhere c = 1 ∧ (∥P ∥ + ∥Q∥/(2λ)) 0 Proof: First, note that ṡ = q̈ − v̇r and therefore M ṡ = fined as in the proof of Theorem 3. −(K + C)s − Φã + r, where ã := â − a. Moreover, we have P −1 ã˙ = Φ⊤ s + (P −1 A − Φ⊤ R−1 Φ)ã + P −1 (Aa + b − ȧ) + III. NEURAL NETWORK ARCHITECTURE AND Φ⊤ R−1 (r + ϵ). Defining V := s⊤ M s + ã⊤ P −1 ã, we have ITS TRAINING ⊤ ⊤ −1 −1 ⊤ −1 ⊤ V̇ = −2s Ks − ã P QP + Φ R Φ ã + 2s r + ⊤ −1 ⊤ ⊤ −1 2ã P (Aa + b − ȧ) + 2ã Φ R (r + ϵ). In this section, we present a statistically principled training By Lemma 2, for a sufficiently small α > 0, we have procedure for a uniformly contractive disturbance predictor K ⪰ αM and Q ⪰ 2αP for all t ≥ 0, which implies −2K ⪯ whose latent dynamics are governed by feature-dependent −2αM and −P −1 QP −1 − Φ⊤ R−1 Φ ⪯ −2αP −1 . There- matrix coefficients. fore, we conclude that V̇ ≤ −2αV + 2s⊤ r + 2ã⊤ P −1 (Aa + b − ȧ) + 2ã⊤ Φ⊤ R−1 (r + ϵ). A. Contractive Network Parametrization Applying Cauchy–Schwarz gives All maps are parameterized by GELU MLPs, with spectral s normalization applied to the linear layers. To ensure contracV V̇ ≤ − 2αV + 2 tion of the latent dynamics (3a), we parameterize A(ϕ, u) = λmax (P )−1 ∧ m S(ϕ, u) + K(ϕ, u) with S(ϕ, u) = −(β + µ(ϕ, u))I − (10) r L(ϕ, u)L(ϕ, u)⊤ and K(ϕ, u) = 12 M (ϕ, u) − M (ϕ, u)⊤ , × . P −1 (Aa + b − ȧ) + Φ⊤ R−1 (r + ϵ) with β > 0 fixed. We note that this contractive-byconstruction parameterization is inspired by [24]. The scalar Lemma 2 with the assumptions of this theorem shows that the µ(ϕ, u) ≥ 0 and the matrices L(ϕ, u), M (ϕ, u) are produced norm in (10) is bounded by √ a constant, say d ≥ 0. For each by MLPs and bounded element-wise via σ(·) and tanh(·), small c > 0, define Wc = V + c. We have Ẇc ≤ −αW 1/2c + respectively; L(ϕ, u) is lower-triangular and M (ϕ, u) is √ −1 d ∗ := α c + c∗ where c m ∧ (∥P0 ∥ + ∥Q∥/(2λ)) . strictly lower-triangular. Since K is skew-symmetric, it does
not affect the symmetric part of A, and we have the following uniform contraction condition: A(ϕ, u) + A(ϕ, u)⊤ /2 = S(ϕ, u) ⪯ −βI, (14) therefore the latent dynamics are contracting. b(ϕ, u), Φ(ϕ), and ψ(ϕ) are parameterized by MLPs; Φ is bounded via tanh(·). One may optionally enforce anchoring constraints such as b(0, 0) = 0 or ψ(0) = 0 by subtracting the value at the origin.
λa ∥a0 ∥2 regularizes the inferred initial condition and could be interpreted as mitigating the identifiability ambiguity. For fixed θ, the latent rollout {ak } depends affinely on a0 ; consequently, dˆk is affine in a0 and the Gaussian NLL induces a positive definite quadratic objective in a0 . Thus, the a0 block can be optimized exactly, and this optimization can be viewed as a special case of Kalman smoothing; see Section III-C. In view of (15), it is natural to use P0 = λ−1 a I, â(0) = 0, R = (∆t) Rd and a small Q = qI ≻ 0 in (8). C. Process Noise Modeling and the Kalman Smoother
Algorithm 1 Alternating training (hard-EM) Require: initial (θ, Rd , Qd ), choice of λa , data windows 1: while stopping criterion is not met do 2: sample mini-batch B 3: for all w ∈ B do 4: instantiate the model for window w from θ 5: apply Kalman smoothing to infer a0:L−1 6: end for 7: update Rd (from output residuals; sparse/EMA) 8: update Qd (from latent residuals; sparse/EMA)∗ 9: update θ by AdamW with the inferred latents fixed 10: end while ∗ optional
B. Training via Hard Expectation–Maximization We train the architecture from data windows {(ϕk , uk , yk )}L−1 k=0 sampled at time step ∆t; windows are extracted from trajectories with a fixed stride and may overlap. A key challenge in fitting (3) from data windows is that the latent initial condition a0 = a(0) is unobserved, windowspecific, and has no direct physical meaning. As a result, a0 is a missing value that must be inferred jointly with the shared parameters, which can create an identifiability ambiguity between a0 and θ. To regularize this inference problem in a principled way, we place a zero-mean Gaussian prior on the initial latent state, a0 ∼ N (0, λ−1 a I), which induces the quadratic penalty λa ∥a0 ∥2 in the training objective. Given parameters θ, we roll out (3a) by RK2 to obtain ak , form predictions dˆk = Φθ (ϕk )ak +ψθ (ϕk ), and assume a Gaussian observation model yk | a0 , θ ∼ N (dˆk , Rd ) with Rd ≻ 0. This yields a Gaussian negative log-likelihood (NLL) loss min
a0 ,θ,Rd
L−1 X
(yk − dˆk (a0 , θ))⊤ Rd−1 (yk − dˆk (a0 , θ))
k=0
(15)
+ log |Rd | + λa ∥a0 ∥2 , implicitly averaged over a mini-batch of windows. We optimize (15) by an alternating scheme, which can be interpreted as a (generalized) hard-EM procedure [25]: (hard E-step) fix θ and Rd and perform exact minimization with respect to a0 to obtain a MAP estimate for each window; (M-step) Fix the inferred a0 and update Rd from mini-batch statistics; then use AdamW to update θ. The quadratic prior term
The training scheme in the previous subsection corresponds to a deterministic latent rollout within each window: once (θ, a0 ) are fixed, the entire latent trajectory {ak }L−1 k=0 is fixed. We optionally add discrete-time process noise Qd ≻ 0 by treating the entire latent trajectory as missing data; Section III-B corresponds to the degenerate special case of Qd = 0. Concretely, let fθ,k−1 (ak−1 ) denote the one-step RK2 update of (3a) using (ϕk−1 , uk−1 ) (from time k − 1 to time k). Since (3a) is affine in a, fθ,k−1 (ak−1 ) is affine in ak−1 . We replace the deterministic update with ak | ak−1 , θ ∼ N fθ,k−1 (ak−1 ), Qd , k = 1, . . . , L−1, (16) while keeping the same observation model for yk and the same prior a0 ∼ N (0, λ−1 a I). Compared to (15), we replace dˆk (a0 , θ) by dˆk (ak , θ) := Φθ (ϕk )ak + ψθ (ϕk ), and the window NLL loss acquires the additional process-noise term ⊤ PL−1 (L − 1) log |Qd | + k=1 ak − fθ,k−1 (ak−1 ) Q−1 a k − d fθ,k−1 (ak−1 ) . For fixed (θ, Rd , Qd ), minimization over the latent trajectory a0:L−1 (the hard E-step) is exactly solvable via a Kalman smoother [26], yielding an inferred latent trajectory for each window. The Kalman smoother is implemented by the RTS formulas. In the M-step, we fix this inferred trajectory and update Rd and optionally Qd from mini-batch statistics, and then θ by AdamW (see Algorithm 1). The identified discrete-time covariances lead to the continuous-time gains in (8): Q ≈ Qd /∆t and R ≈ (∆t)Rd , as noted in Remark 2. In our implementations, we used diagonal Rd and isotropic Qd = qd I for simplicity. IV. EXPERIMENTAL RESULTS Implementing the controllers (9) and (12) requires discretizing (8). We use a Joseph-form Kalman filter, with forward Euler discretization for the P Φ⊤ s term. For both examples, we used da = 3. We compare our controller with (i) a model-based PD control loop, (ii) the fixed-decay ablations inspired by [1], [2], and (iii) a DAC controller in which an LTI system is identified by [27, Sections 4.3.1, 4.4.2], a variant of the N4SID algorithm [28]. We denote the last controller as N4SIDv-DAC. The fixed-decay ablations, denoted FixedDecay, were trained by the Kalman-smootherbased hard-EM procedure restricted to the function class A = −λI, b ≡ 0, and ψ ≡ 0. For the adaptive controllers, datasets were split into training, validation, and test sets in a 70:15:15 ratio. For the N4SID variant, the number of block
rows i in the past and future Hankel matrices, following the notation of [27], was selected to minimize the validationset MSE. The LTI model was then re-estimated using the combined training and validation sets. The LTI disturbance augmentations are of the form ȧ = Aa + B1 ϕ + B2 u, dˆ = Ca + D1 ϕ + D2 u,
(17a) (17b)
ˆ For example, substituting (17b) into (18) yields with d ≈ d. v̇ ≈ Av v + (B + D2 )u + D1 ϕ + Ca. Therefore, considering (12), the disturbance-compensating command is obtained by solving the equation (B+D2 )u+D1 ϕ+C â = v̇d −Av v−Ks in u. Although the LTI model is presented in continuoustime notation, system identification was performed in discrete time.
epochs using a learning rate of 0.001 and a weight decay of 0.01. The training took less than 3 minutes on an Apple M4 Pro laptop. b) Predictions for Disturbance d: The predictive capability of the architecture was evaluated by rolling out the dynamics in (3a) and (4) using the input, state, and waterlevel measurements from the test set. Figure 3 illustrates the predicted and measured disturbances on two test-set trajectories.
A. A Slippery Tracked Vehicle with a Partially Filled Liquid Tank and a Swinging Pendulum The presented architecture was evaluated on a tracked ground vehicle (“GVR-Bot”) carrying a pendulum and a liquid tank filled with water to approximately 30% capacity. The experiment also provides a test of generalization to an unseen intermediate fill level because the 30% level was absent from the training data. The choice of a liquid-carrying platform was partly motivated by the growing interest in onorbit refueling [29]. To emulate low-traction, slippery terrain, the tracks were wrapped with low-friction tape. Water-tank level sensor measurements w ∈ R were sampled at 18 Hz. The platform had an NVIDIA Jetson AGX Orin as a companion computer, and used visual–inertial odometry (VIO) for localization. The control loop ran at 50 Hz. The vehicle’s velocity is modeled by the following firstorder lag [2]: v̇ = Av v + Bu + d, (18) where v = [vx , ωz ] denotes the vehicle’s forward velocity and yaw rate, u = [vx,cmd , ωz,cmd ] denotes the commanded values, and B = −Av = diag(1/0.2, 1/0.15), determined by system identification on hardware. The vehicle’s internal PID controller is inaccessible, and only the velocity commands are available, which motivates the model (18). The liquid sloshing and pendulum motion are expected to contribute to time-varying changes in the vehicle’s traction/slip behavior and, consequently, in the observed command-to-velocity response. Such unmodeled effects are lumped into d in (18). a) Data Collection and Training: Data were collected for about 11 minutes for each water-tank level (≈ 60%, 0%) in an open testing area, with the vehicle driven at various speeds and in various directions without following a prescribed trajectory. The vehicle was driven through the full operational range of speeds and angular rates, with maneuvers such as abrupt starts and stops or smooth accelerations. We trained the neural network (≈ 4.9k trainable parameters) using (v, w, u, y) data with ∆t = 0.01. Here, y is the estimated d under the zero-order hold assumption, and we used ϕ = (v, w). For both the fixed-decay ablation and the proposed method, the neural networks were trained for 100
Fig. 3. (Tracked vehicle) Test-set rollouts for the predicted disturbance dˆ ∈ R2 as functions of time, on two different test-set trajectories. Black: ˆ disturbance proxy y(t); red: rollout prediction d(t). The proxy y is visibly noisier. The average prediction MSE was 3.53; it was 3.49 after discarding the initial half second for each trajectory as a burn-in period (vertical dashed line).
c) Tracking Controller Test: We compare the proposed controller with three other controllers: our FixedDecay ablation, N4SIDv-DAC, and the model-based PD loop used in [2, Eqs. 14–18] which contains nonlinearities. The vehicle’s commanded trajectory consisted of two consecutive double lane changes, which excited the water and pendulum dynamics while inducing track slip. The model-based PD baseline already provides strong overall tracking, particularly in the forward-velocity channel. The benefits of incorporating disturbance predictors are most apparent in yaw-rate regulation. As shown in Table I, the proposed controller achieves the lowest angular-rate and velocity RMS errors and ties N4SIDv-DAC for the lowest position RMS error. The competitive angular-rate and position tracking performance of N4SIDv-DAC suggests that an LTI augmentation captures a significant part of the model mismatch over the tested regime. In contrast, the fixed-decay ablation performed less favorably in this experiment. Metric
NPD
Ang. rate error RMS [rad/s] 0.449 Vel. error RMS [m/s] 0.151 Pos. error RMS [m] 0.085
N4SIDv-DAC FixedDecay 0.414 0.160 0.078
0.468 0.156 0.081
Ours 0.407 0.149 0.078
TABLE I T RACKED - VEHICLE RMS TRACKING ERRORS . NPD: nonlinear PD loop; N4SIDv-DAC: LTI-based disturbance-accommodating control; FixedDecay: fixed-decay ablation.
B. Coupled Duffing Oscillators: a Robustness Sweep We consider a pair of Duffing oscillators coupled by a nonlinear spring: m1 ẍ1 = −k1 x1 − α1 x31 + u + (d − c1 ẋ1 ), m2 ẍ2 = −k2 x2 − α2 x32 − c2 ẋ2 − d, d = kc (x2 − x1 ) + αc (x2 − x1 )3 ,
(19a) (19b) (19c)
Fig. 5. (Coupled nonlinear oscillator) Test-set rollouts from Eq. (3) for the predicted disturbance as functions of time. Black: disturbance proxy ˆ y(t); red: rollout prediction d(t). The proxy y is visibly noisier, while dˆ is smooth. Mean squared error is 2.54 overall and 0.88 after discarding the first second of each trajectory as a burn-in period (vertical dashed line).
Fig. 4.
(Tracked vehicle) Tracking runs under four controllers.
where we only observe x1 , ẋ1 , u and are unaware of the existence of the coupled x2 -system. We train the neural networks with a synthetic dataset of (x1 , ẋ1 , u, y) (∆t = 0.01). Here, the samples of x1 and ẋ1 are corrupted by zero-mean Gaussian noise with standard deviation 0.005; y is computed using forward differences of the sampled ẋ1 . The parameters used for the nonlinear oscillators were: m1 = 1.0, c1 = 0.30, k1 = 1.0, α1 = 1.7,
m2 = 0.7, c2 = 0.25, k2 = 0.8, α2 = 0.2,
kc = 2.0, αc = 5.0.
(20a) (20b) (20c) (20d)
Our empirical evaluation considers two complementary aspects: (a) we show that (3) with ϕ = (x1 , ẋ1 ) predicts d (or more precisely, its proxy y) accurately on the test set, and (b) we compare the controller (9) with the baseline model-based PD (u = m1 v̇r + k1 x1 + α1 x31 − Ks + c1 ẋ1 ), N4SIDv-DAC, and the fixed-decay representation ablation. a) Predictions for Disturbance d: We rolled out (3) ˆ and compared starting from a0 = 0 to get the estimates d, them to the test-set labels y: see Fig. 5. Note that the initial transients from the fixed initial condition a0 = 0 decay as time passes, as expected from the contraction property (Proposition 1). This predictive capability arises from generalizing prior “fixed-decay” representation-learning approaches to a more flexible learned dynamical model. b) Tracking Controller Test: To illustrate the benefit of predictive capability, we made y available at only 1 Hz. We also injected different levels of Gaussian noise (its standard deviation denoted by σ) into y to test for robustness (Fig. 6). We used Λ = 2.0 and K = 12.0 for all four controllers. In this experiment, the acceleration ẍ1 was estimated by taking backward finite differences of the observed ẋ1 , and then filtered by an exponential moving average. The disturbance proxy y was then calculated as m1 ẍ1 + c1 ẋ1 + k1 x1 +α1 x31 −u; after this, it was corrupted by noise level σ. Note that σ = 0 already results in a noisy proxy y because the measurements for x1 and ẋ1 are noisy. Each reference
Fig. 6. (Coupled nonlinear oscillator) Robustness sweep for tracking control under additional Gaussian noise injected into the sparsely provided disturbance proxy y (averaged over 30 reference trajectories, with each trajectory being 20 seconds long). Mean tracking MSE is plotted versus the added noise level σ, where σ = 0 denotes the baseline noise setting. The proposed dynamical-representation controller performs the best. As σ increases, the neural adaptive controllers degrade because the disturbance proxy becomes increasingly corrupted. The simpler N4SIDv-DAC baseline resulted in a flat graph because the fitted discrete-time A, B1 , and B2 satisfied ∥A∥ ≈ 0.09 ≪ 1, ∥B1 ∥ ≈ 0.08 ≪ 1, ∥B2 ∥ ≈ 0.005 ≪ 1 and hence the disturbance compensation mostly came from dˆ ≈ D1 ϕ + D2 u.
trajectory was generated as a sum of randomized sinusoidal signals. Figure 6 summarizes the robustness sweep for the proposed dynamical-representation controller, the fixed-decay representation ablation, the N4SIDv-DAC controller, and the model-based PD baseline. In the nominal-noise case (σ = 0), the proposed controller achieves the lowest average tracking MSE; in particular, our controller reduced tracking MSE by 27% compared to the closest competitor. In this highly nonlinear setting, N4SIDv-DAC performed worse than model-based PD, suggesting that the identified LTI disturbance model was inadequate. Overall, the proposed method provides the strongest performance and remains strong throughout the sweep. V. CONCLUSIONS This paper presented a composite adaptive control framework based on a representation-learning approach for systems with unobserved (or partially observed),
dynamically coupled disturbance sources. The presented controller complements prior “fixed-decay” last-layer representation-learning approaches by incorporating the classical disturbance-accommodating control (DAC) viewpoint. By modeling the disturbance as the output of a uniformly contractive latent dynamical system, the proposed method introduces a dynamics prior that supports disturbance prediction while mitigating sensitivity to latent-state initialization errors. The latent dynamics are parameterized by neural networks and trained using a statistically principled hard expectation–maximization procedure. Hardware experiments on a tracked vehicle carrying a partially filled liquid tank and a pendulum illustrate the practical feasibility of the proposed disturbance estimator and the composite adaptive controller. The learned dynamical representation improves overall tracking performance during aggressive maneuvering, with the clearest improvement observed in yaw-rate regulation. In simulation, experiments on a system of coupled Duffing oscillators demonstrate accurate disturbance prediction and improved tracking relative to the baselines, while remaining strong under progressively noisier disturbance proxies y. Here, the sparse availability of y allowed us to assess the role of predictive capability for maintaining tracking performance. Taken together, the results indicate that the proposed method holds substantial promise for real-world deployment, particularly when disturbance proxies are sampled at intervals which are not negligible relative to the characteristic disturbance timescale, making predictive compensation important between observations. Future work may include evaluating the method across a wider range of hardware platforms and experimental conditions, as well as systematically studying the training dynamics of the hard-EM procedure, particularly the evolution of covariance estimates. Future work may also include replacing the present hard E-step with the standard EM E-step, which takes expectations with respect to the posterior distribution over latent trajectories. The proposed dynamical representation could also serve as a disturbance model in stochastic MPC with chance constraints [30], potentially useful for real-time trajectory planning. Acknowledgment: The authors thank G. X. Johnson, J. Alindogan, and S. Sukhatme for their help with hardware troubleshooting. The authors also thank the anonymous reviewers. R EFERENCES [1] M. O’Connell, G. Shi, X. Shi, K. Azizzadenesheli, A. Anandkumar, Y. Yue, and S.-J. Chung, “Neural-fly enables rapid learning for agile flight in strong winds,” Science Robotics, vol. 7, no. 66, p. eabm6597, 2022. [2] E. S. Lupu, F. Xie, J. A. Preiss, J. Alindogan, M. Anderson, and S.-J. Chung, “Magic-vfm: Meta-learning adaptation for ground interaction control with visual foundation models,” IEEE Trans. Robot., vol. 41, p. 180–199, Jan. 2025. [3] F. C. Yilmaz, R. Onler, E. Tatlicioglu, and E. Zergeroglu, “A transfer learning based deep neural network adaptive controller for the furuta pendulum subject to uncertain disturbance signals,” Scientific Reports, vol. 15, p. 24012, 2025. [4] C. D. Johnson, “Disturbance-accommodating control; an overview,” in American Control Conf., 1986, pp. 526–536. [5] E. Davison, “The robust control of a servomechanism problem for linear time-invariant multivariable systems,” IEEE Trans. Autom. Control, vol. 21, no. 1, pp. 25–34, 1976.
[6] B. Francis and W. Wonham, “The internal model principle of control theory,” Automatica, vol. 12, no. 5, pp. 457–465, 1976. [7] M. Bin, L. Marconi, and A. R. Teel, “Adaptive output regulation for linear systems via discrete-time identifiers,” Automatica, vol. 105, pp. 422–432, 2019. [8] G. Joshi, J. Virdi, and G. Chowdhary, “Asynchronous deep model reference adaptive control,” in Proceedings of the 2020 Conference on Robot Learning, vol. 155. PMLR, 2021, pp. 984–1000. [9] Y. Song, M. Steinweg, E. Kaufmann, and D. Scaramuzza, “Autonomous drone racing with deep reinforcement learning,” in 2021 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS 2021). IEEE Press, 2021, p. 1205–1212. [10] K. Y. Chee, T. C. Silva, M. A. Hsieh, and G. J. Pappas, “Enhancing sample efficiency and uncertainty compensation in learning-based model predictive control for aerial robots,” in 2023 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS 2023), 2023. [11] J. Kabzan, L. Hewing, A. Liniger, and M. N. Zeilinger, “Learningbased model predictive control for autonomous racing,” IEEE Robot. Autom. Lett., vol. 4, no. 4, pp. 3363–3370, 2019. [12] A. Kumar, Z. Li, J. Zeng, D. Pathak, K. Sreenath, and J. Malik, “Adapting rapid motor adaptation for bipedal robots,” in IEEE/RSJ Int. Conf. Intell. Robot. Sys., 2022. [13] T. Miki, J. Lee, J. Hwangbo, L. Wellhausen, V. Koltun, and M. Hutter, “Learning robust perceptive locomotion for quadrupedal robots in the wild,” Science Robotics, vol. 7, no. 62, p. eabk2822, 2022. [14] P. Krupa, M. Zanon, and A. Bemporad, “Learning disturbance models for offset-free reference tracking,” IEEE Trans. Autom. Control, vol. 70, no. 11, pp. 7683–7690, 2025. [15] H. Zhang, J. Ge, J. Su, K. Gu, F. Wang, W.-H. Chen, and S. Li, “Model predictive control with residual learning and real-time disturbance rejection: Design and experimentation,” Control Eng. Practice, vol. 165, p. 106587, 2025. [16] J. Jia, W. Zhang, K. Guo, J. Wang, X. Yu, Y. Shi, and L. Guo, “Evolver: Online learning and prediction of disturbances for robot control,” IEEE Trans. Robot., vol. 40, pp. 382–402, 2024. [17] Y. Gu, L. Chen, C. Li, and S. Yin, “Multiple-model state-space system identification with time delay using the em algorithm,” J. Franklin Inst., vol. 361, no. 16, p. 107113, 2024. [18] W. Wang, S. Liu, Y. Jiang, J. Sun, W. Xu, X. Chen, Z. Dong, and W. Jiao, “A robust filter and smoother-based expectation–maximization algorithm for bilinear systems with heavy-tailed noise,” Mech. Syst. Signal Proc., vol. 236, p. 112912, 2025. [19] W. Lohmiller and J.-J. E. Slotine, “On contraction analysis for nonlinear systems,” Automatica, vol. 34, no. 6, pp. 683–696, 1998. [20] F. Bullo, Contraction Theory for Dynamical Systems, 1.3 ed. Kindle Direct Publishing, 2026. [21] J.-J. E. Slotine and W. Li, “Composite adaptive control of robot manipulators,” Automatica, vol. 25, no. 4, pp. 509–519, 1989. [22] ——, Applied Nonlinear Control, ser. Prentice-Hall International Editions. Prentice-Hall, 1991. [23] H. Abou-Kandil, G. Freiling, V. Ionescu, and G. Jank, Matrix Riccati Equations in Control and Systems Theory. Birkhäuser Basel, 2003. [24] H. Beik-Mohammadi, S. Hauberg, G. Arvanitidis, N. Figueroa, G. Neumann, and L. Rozo, “Neural contractive dynamical systems,” 2024. [Online]. Available: https://arxiv.org/abs/2401.09352 [25] G. Celeux and G. Govaert, “A classification em algorithm for clustering and two stochastic versions,” Comp. Stat. Data Anal., vol. 14, no. 3, pp. 315–332, 1992. [26] D. Sanz-Alonso, A. Stuart, and A. Taeb, Inverse Problems and Data Assimilation, ser. London Mathematical Society Student Texts. Cambridge University Press, 2023. [27] P. van Overschee and B. De Moor, Subspace Identification for Linear Systems. Springer New York, NY, 1996. [28] P. Van Overschee and B. De Moor, “N4sid: Subspace algorithms for the identification of combined deterministic-stochastic systems,” Automatica, vol. 30, no. 1, pp. 75–93, 1994, special issue on statistical signal processing and control. [29] M. Ramezani, M. Amin Alandihallaj, B. C. Yalçın, M. A. Olivares Mendez, and A. M. Hein, “Fuel-aware autonomous docking using rlaugmented mpc rewards for on-orbit refueling,” Acta Astronautica, vol. 238, pp. 690–705, 2026. [30] A. Mesbah, “Stochastic model predictive control: An overview and perspectives for future research,” IEEE Control Syst. Mag., vol. 36, no. 6, pp. 30–44, 2016.
[31] S. Särkkä and L. Svensson, Bayesian Filtering and Smoothing, 2nd ed., ser. Institute of Mathematical Statistics Textbooks. Cambridge University Press, 2023.
APPENDIX
From this, we can calculate the soft-EM loss averaged over a mini-batch of B windows (indexed by w = 1, . . . , B) as "L−1 B i 1 X 1 Xh y,(w) r tr(Rd−1 Sk ) + log |Rd | Q (η) = B w=1 2 k=0
A. Hard and Soft Expectation–Maximization We briefly outline the changes required to replace hard EM with soft EM. a) Hard vs. Soft EM: Suppose a window of measured features, control inputs, and disturbance observations {(ϕk , uk , yk )}L−1 k=0 is given. Write y := y0:L−1 and η := (θ, Rd , Qd ). In the general case, the missing value is the entire trajectory Z = a0:L−1 . For a fixed η, the statistical model in Section III-C can be written as ak+1 = Fk ak + gk + ξk ,
ξk ∼ N (0, Qd ),
yk = Hk ak + hk + ϵk ,
ϵk ∼ N (0, Rd ),
Ck = Pk+1 G⊤ k.
a,(w)
tr(Q−1 d Sk
) + log |Qd |
i
k=0
# (w) 2 (w) + λa ∥m0 ∥ +tr P0 +L(da + du ) log(2π)−da log λa . (23) Note that the last line in (23) can be omitted in the loss function for the M-step. Now, our M-step is as follows: compute the new (Rd , Qd ) by B L−1
with a0 ∼ N (0, λ−1 a I), where λa > 0 is a fixed hyperparameter. Here, a0 , ξk , ϵk are independent. Note that (Fk , gk ) comes from discretizing (3a) in time by, say, RK2 or forward Euler; the coefficients (Fk , gk , Hk , hk ) depend on θ. At iteration r, let π r (z) := p(z | y, η r ) and Qr (η) := −Eπr [log p(y, Z | η)]. Standard (soft) EM minimizes Qr (η) with respect to η, holding π r fixed, to obtain η r+1 . Hard EM instead selects ẑ r := arg maxz π r (z) and minimizes − log p(y, ẑ r | η) with respect to η, holding ẑ r fixed. In other words, soft EM averages over the latent uncertainty Z ∼ π r , whereas hard EM uses a single point estimate of Z, replacing the soft E-step (expectation) with a hard E-step (posterior maximization). b) The Soft EM Loss: Denote the current parameters as η r = (θr , Rdr , Qrd ), where we assume Rdr ≻ 0 and Qrd ≻ 0. Compute the smoothed means and covariances using a Kalman smoother under parameters η r , and denote the smoothed moments by mk = Eπr [ak ], Pk = Covπr (ak ), and Ck = Covπr (ak+1 , ak ); here, the cross-covariances are obtained from the smoothing gains and the smoothed covariances as follows [31, Eq. (12.13)]: Gk := Pk|k Fk⊤ (Pk+1|k )−1 ,
+
L−2 Xh
(21)
The posterior moments (mk , Pk , Ck ) remain fixed throughout the M-step. For candidate parameters η, let a0:L−1 ∼ π r and define the random variables rky := yk − Hk ak − hk and rka := ak+1 − Fk ak − gk . We have rky = ek − Hk (ak − mk ), rka = vk + (ak+1 − mk+1 ) − Fk (ak − mk ), where ek := yk − Hk mk − hk and vk := mk+1 − Fk mk − gk . Therefore, the residual second moments are ⊤ Sky := Eπr [rky (rky )⊤ ] = ek e⊤ k + Hk Pk Hk ,
(22a)
Ska := Eπr [rka (rka )⊤ ] = vk vk⊤ + Pk+1 + Fk Pk Fk⊤ − Ck Fk⊤ − Fk Ck⊤ .
(22b)
Rd+ = Q+ d =
1 X X y,(w) Sk , BL w=1 k=0 B L−2 X X
1 B(L − 1) w=1
(24a) a,(w)
Sk
,
(24b)
k=0
and then update θ using an optimizer such as AdamW to minimize the loss in (23), omitting its last line. In practice, eigenvalue clipping might be needed to ensure positive definiteness of (Rd+ , Q+ d ). c) The Qd = 0 Variant: In Section III-B, the latent trajectory is determined by θ and the initial state a0 , so we take the missing value as Z = a0 . Hard EM uses the MAP estimate for a0 ; soft EM instead retains the posterior distribution a0 | y, η r ∼ N (m0 , P0 ). The posterior moments (m0 , P0 ) remain fixed throughout the soft-EM M-step. Write the deterministic rollout-plus-readout map as dˆk = Mk a0 + ck where (Mk , ck ) depend on the candidate θ, and define ek := yk − Mk m0 − ck ,
(25)
⊤ Sky := ek e⊤ k + M k P0 M k .
(26)
Assuming Rd ≻ 0, the soft-EM objective reduces to B L−1 i 1 X 1 Xh y,(w) Qr (η) = tr(Rd−1 Sk ) + log |Rd | (27) B w=1 2 k=0
up to additive constants. B. Regularization in the N4SID-Variant Implementation Relative to the algorithm in [27, Fig. 4.7], our implementation applies ridge-regularized pseudoinverse approximations at three locations: (i) in the oblique projections [27, Eq. (1.7)], when computing Oi and Oi+1 (step 1); (ii) Γ†i and Γ†i−1 in determining the state sequences (step 5); ⊤ and (iii) Z † , with Z := [X̃i⊤ , Ui|i ], in the least-squares solve for (A, B, C, D) (step 6). For each ridge-regularized pseudoinverse approximation, we solve (G + λG I)X = BG , where G ∈ Rn×n is the associated Gram matrix, BG denotes the right-hand side, and λG = 0.1 tr(G)/n. Here, 0.1 is a relative regularization parameter chosen for our implementation.