Paper deep dive
Variational meta-learning inference for low dimensional neural system identification
Matteo Rufolo, Dario Piga, Marco Forgione
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 88%
Last extracted: 7/22/2026, 3:19:58 AM
Summary
The paper proposes a variational meta-learning framework for low-dimensional neural system identification to address overfitting and lack of uncertainty quantification in deep learning models under low-data regimes. It extends a deterministic manifold meta-learning approach by introducing amortized Variational Inference to learn a generative prior over the low-dimensional parameter manifold. The method combines Maximum A Posteriori (MAP) estimation with Laplace approximation to provide calibrated uncertainty bounds, evaluated on static regression and the Bouc-Wen dynamical system benchmark.
Entities (7)
Relation Signals (6)
Variational Meta-Learning → uses → Variational Inference
confidence 95% · based on amortized Variational Inference, where a generative prior over the low-dimensional parameter manifold is learned.
Variational Meta-Learning → evaluatedon → Bouc-Wen System
confidence 90% · Evaluated on a static regression task and the Bouc--Wen dynamical system benchmark
Manifold Meta-Learning → extendsto → Variational Meta-Learning
confidence 90% · We propose a fully probabilistic extension of the manifold meta-learning framework, based on amortized Variational Inference
Manifold Meta-Learning → addresses → Overfitting
confidence 85% · The recently developed manifold meta-learning framework addresses the data efficiency problem
Variational Meta-Learning → uses → Laplace Approximation
confidence 85% · we combine Maximum A Posteriori estimation with the Laplace approximation to yield a mathematically grounded posterior approximation.
Neural Network → suffersfrom → Overfitting
confidence 80% · heavily parameterized neural networks are prone to overfitting in low-data regimes
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Deep learning has proven highly effective for nonlinear system identification, but heavily parameterized neural networks are prone to overfitting in low-data regimes and lack reliable uncertainty quantification. The recently developed manifold meta-learning framework addresses the data efficiency problem by restricting the model parameters to a meta-learned low-dimensional manifold. However, that method is purely deterministic. We propose a fully probabilistic extension of the manifold meta-learning framework, based on amortized Variational Inference, where a generative prior over the low-dimensional parameter manifold is learned. During task-specific adaptation, we combine Maximum A Posteriori estimation with the Laplace approximation to yield a mathematically grounded posterior approximation. Evaluated on a static regression task and the Bouc--Wen dynamical system benchmark, the proposed approach achieves predictive accuracy comparable to its deterministic counterpart while successfully providing calibrated uncertainty bounds in severely low-data regimes.
Tags
Links
- Source: https://arxiv.org/abs/2607.18965v1
- Canonical: https://arxiv.org/abs/2607.18965v1
Trouble viewing inline? Open PDF directly →
Full Text
44,248 characters extracted from source content.
Expand or collapse full text
Variational meta-learning inference for low dimensional neural system identification Matteo Rufolo1, Dario Piga1, Marco Forgione1 Corresponding author: matteo.rufolo@supsi.ch.1SUPSI-DTI-IDSIA, Dalle Molle Institute for Artificial Intelligence, Lugano, Switzerland. Abstract Deep learning has proven highly effective for nonlinear system identification, but heavily parameterized neural networks are prone to overfitting in low-data regimes and lack reliable uncertainty quantification. The recently developed manifold meta-learning framework addresses the data efficiency problem by restricting the model parameters to a meta-learned low-dimensional manifold. However, that method is purely deterministic. We propose a fully probabilistic extension of the manifold meta-learning framework, based on amortized Variational Inference, where a generative prior over the low-dimensional parameter manifold is learned. During task-specific adaptation, we combine Maximum A Posteriori estimation with the Laplace approximation to yield a mathematically grounded posterior approximation. Evaluated on a static regression task and the Bouc–Wen dynamical system benchmark, the proposed approach achieves predictive accuracy comparable to its deterministic counterpart while successfully providing calibrated uncertainty bounds in severely low-data regimes. I Introduction In recent years, system identification has increasingly adopted deep learning techniques, leveraging the high expressivity of black-box models such as neural networks to capture complex, nonlinear dynamics when the underlying physics are unknown [16]. However, a primary challenge in deploying these heavily parameterized models is their reliance on large datasets to ensure reliability across diverse operating conditions. When data is scarce, training these architectures becomes prone to overfitting and poor generalization. This contrasts with physics-based models, which are typically data-efficient and capable of generalizing beyond the training regime, though deriving such first-principles models is often unfeasible for complex systems. Meta-learning techniques have recently gained traction within the system identification and control communities to facilitate rapid adaptation across related tasks, spanning distributionally robust task weighting [17], probabilistic transformer system identification [18], transformer-based filtering [6] and controller-tuning transfer [12]. A key architectural question shared across this line of work is how to summarize a variable-length dataset into a fixed-size task representation for conditioning: permutation-invariant set encoders such as Deep Sets [19] and the Set Transformer [13] enforce exchangeability when the ordering of observations carries no information, while sequence encoders such as (bidirectional) recurrent networks are preferred when temporal ordering is itself informative. This amortized-encoding paradigm also underlies the neural process family [9], which learns to map a context set directly to a predictive distribution rather than solving an explicit inner optimization per task. In particular, this paper extends a recently proposed manifold meta-learning framework [8], which leverages an amortized learning strategy, instantiated as a permutation-invariant or sequence encoder depending on whether the task data is order-invariant or sequential, as discussed above, to automatically extract simplified, efficient model architectures. The approach relies on an over-parameterized base architecture capable of capturing the dynamics of all systems within a meta-dataset. Since the physical factors varying across systems correspond to far fewer degrees of freedom than the total parameter count, it is natural to seek a low-dimensional manifold within the base architecture’s parameter space. Restricting the model parameters to this meta-learned manifold yields a reduced-complexity representation that achieves performance comparable to the full model in data-rich scenarios, while substantially outperforming it when training data is scarce. Unlike neural processes, which predict directly from the amortized representation, the manifold framework retains an explicit low-dimensional parameter ϕφ and lifting function, using the encoder’s output only to initialize, rather than replace, a subsequent optimization-based refinement (Section I-B). The primary contribution of this work is the extension of [8] into a probabilistic setting, drawing on recent advances in probabilistic machine learning [18, 9]. Specifically, we adopt a variational learning approach inspired by the Variational Autoencoder (VAE) framework [11], in which an inference network processes each training dataset and outputs an approximate Gaussian posterior over the model’s reduced-complexity parameters in the low-dimensional manifold. Despite the Gaussian restriction on the latent space, the highly nonlinear transformation from the manifold to the full-complexity model output can induce arbitrarily complex, non-Gaussian distributions over the predicted outputs. This probabilistic formulation enables rigorous uncertainty quantification and enhances model interpretability without sacrificing predictive accuracy. The remainder of the paper is organized as follows. Section I reviews the deterministic manifold meta-learning framework. Section I details the mathematical formulation, training, and testing procedures for the proposed probabilistic extension. Finally, Section IV presents two numerical examples: a static regression and a dynamical system identification problem, demonstrating the framework’s capabilities in uncertainty quantification and interpretability while maintaining similar performance to its deterministic counterpart. I Manifold Meta-Learning for Neural System Identification I-A Data Distribution In line with the meta-learning setting for system identification, we assume access to a meta-dataset =D(i)i∈ℕD=\D^(i)\_i , where each dataset D(i)=((i),(i))D^(i)=(u^(i),y^(i)) consists of input-output sequences of length N generated by a specific dynamical system S(i)S^(i). We assume the underlying data-generating mechanism across datasets is governed by a low-dimensional latent variable ∈ℝnzz ^n_z. This variable corresponds to the abstract factors of variation across systems (e.g., varying payload masses in a robotic arm or varying friction coefficients). The inputs (i)u^(i) are independently sampled from a distribution p()p(u), while the outputs depend on both the inputs and the system-specific latent variable via p(∣,)p(y ,z). Overall, the collection of datasets can be formalized as sampling from: p(D)=p()∫p()p(∣,).p(D)=p(u) p(z)p(y ,z)dz. (1) I-B Manifold Meta-Learning and Amortized Optimization The goal of [8] is to learn a model structure of the input-output dependency ↦u across datasets governed by a low-dimensional parameter vector ϕφ, consistent with the low-dimensionality hypothesis in (1). Concretely, this is achieved by restricting the weights θ of an overparameterized neural network, referred to as base architecture: ^=F(,θ),θ∈ℝnθ y=F(u,θ), θ ^n_θ (2) to a low-dimensional manifold via a learned lifting function P:ℝnϕ→ℝnθP:R^n_φ ^n_θ, with nϕ≪nθn_φ n_θ. This yields the reduced-complexity model structure: ^=F(,P(ϕ)),ϕ∈ℝnϕ. y=F\! (u,\,P(φ) ), φ ^n_φ. (3) In [8], the lifting function is parameterized: P=PγP=P_γ, where γ∈ℝnγ ^n_γ are meta-parameters. The latter are meta-learned across datasets to ensure that each input-output dependency in D can be accurately approximated by F(,Pγ(ϕ))F(u,P_γ(φ)), for some dataset-dependent ϕφ and for a shared γ. This is equivalent to meta-learning a low-dimensional manifold ℳ=θ∈ℝnθ:θ=Pγ(ϕ),ϕ∈ℝnϕM=\θ ^n_θ:θ=P_γ(φ),φ ^n_φ\ in base architecture’s parameter space that preserves the base architecture’s approximation capabilities within the meta-dataset. By splitting each dataset D into a training portion Dtr=(tr,tr)D_tr=(u_tr,y_tr) of length NtrN_tr, and a test portion Dte=(te,te)D_te=(u_te,y_te) of length NteN_te (with N=Ntr+NteN=N_tr+N_te), the objective is formalized as a bi-level optimization problem: γ γ =argminγp(D)[ℒ(te,F(te;Pγ(ϕ^))] = _γE_p(D) [L(y_te,F (u_te;\,P_γ( φ) ) ] (4a) ϕ φ =argminϕℒ(tr,F(tr;Pγ(ϕ))), = _φL(y_tr,F(u_tr;\,P_γ(φ))), (4b) where ℒL is the mean squared error loss. In essence, (4) explicitly seeks the manifold that maximizes the average test performance, when fitting of reduced-complexity models is done by optimizing ℒL on each training split. Because exact bi-level optimization is computationally prohibitive, the framework employs an amortized optimization procedure [1]. The inner optimization (4b) is replaced by an encoder network EψE_ψ that directly predicts the low-dimensional parameter ϕφ from the training data. This allows for end-to-end meta-learning according to the criterion: γ^,ψ^=argminγ,ψp(D)[ℒ(te,F(te;Pγ(Eψ(tr,tr)))], γ, ψ= _γ,ψE_p(D) [L(y_te,F(u_te;\,P_γ(E_ψ(u_tr,y_tr))) ], (5) with the expectation approximated by the sample average over datasets in D. At test time, the encoder Eψ^E_ ψ provides a fast, zero-shot estimation of ϕφ. However, because Eψ^E_ ψ is trained to approximate the solution of a nonconvex inner optimization (4b) via a single amortized forward pass, its output is treated throughout this paper as an initializer rather than a final estimate: at meta-test time (Section I-B) it is always followed by gradient-based MAP refinement, so the encoder’s role is to place the optimizer in the right basin of attraction rather than to be trusted in isolation. via gradient descent. Limitations of the Deterministic Framework: While this manifold meta-learning strategy successfully reduces model complexity, it is fundamentally deterministic. The learned encoder Eψ^E_ ψ maps a dataset into a single point estimate ϕφ on the manifold, failing to quantify the epistemic uncertainty inherent in limited-data scenarios. In practical system identification, especially for safety-critical applications, knowing the uncertainty in both the parameter estimates (ϕφ) and the resulting output predictions ( y) is paramount. This fundamental limitation motivates a probabilistic reformulation of the manifold meta-learning discussed so far. I Variational Extension: Formulation I-A Learning the probability of the latent variable To explicitly account for the uncertainty lacking in the deterministic framework, we transition to a generative probabilistic model. To this aim, we introduce learnable prior pω(ϕ)p_ω(φ) and likelihood pγ((i)∣(i),ϕ)p_γ(y^(i) ^(i),φ) with meta-parameters ω and γ, respectively. The likelihood is chosen as: pγ((i)∣(i),ϕ)=(F((i),Pγ(ϕ)),σe2),p_γ(y^(i) ^(i),φ)=N(F(u^(i),P_γ(φ)), _e^2), (6) i.e., a Gaussian whose mean has the same structure of the point prediction in the deterministic manifold framework, and characterized by a noise variance σe2 _e^2. 111In this section, the noise variance σe2 _e^2 is assumed to be known to simplify notation and derivations. In the numerical experiments, σe2 _e^2 is learned together with the other likelihood parameters γ. The prior pω(ϕ)p_ω(φ) is also Gaussian and it is parameterized by its mean vector and a covariance matrix. For given γ,ωγ,ω, the posterior distribution pγ,ω(ϕ∣D(i))=pγ((i)∣(i),ϕ)pω(ϕ)pγ,ω((i)∣(i))p_γ,ω(φ D^(i))= p_γ(y^(i) ^(i),φ)p_ω(φ)p_γ,ω(y^(i) ^(i)) (7) describes the posterior belief of the latent parameters of system (i)(i), given the observation of D(i)=((i),(i))D^(i)=(u^(i),y^(i)). However, (7) is a meaningful characterization of uncertainty only if the prior and likelihood accurately describe the generative process underlying the meta-dataset D. The meta-parameters γ,ωγ,ω must therefore be optimized accordingly. Moreover, because the base architecture F is a neural network, the marginal likelihood, or evidence pγ,ω((i)∣(i))=∫ϕpγ((i)∣(i),ϕ)pω(ϕ)pγ,ωϕp_γ,ω(y^(i) ^(i))= _φp_γ(y^(i) ^(i),φ)p_ω(φ)p_γ,ω\;dφ appearing in the denominator of (7) is intractable, and consequently the posterior cannot be computed analytically. To address both challenges, we rely on Variational Inference [2]. Remark 1 The Gaussian prior assumption is not restrictive, as the nonlinear transformations induced by F and PγP_γ can map this simple prior into arbitrarily complex distributions over the output space. Moreover, one can always choose a fixed isotropic Gaussian prior (,)N(0,I), since the effect of a different mean and covariance can be absorbed into the first layer of PγP_γ. We adopt this isotropic prior in the remainder of the paper, and accordingly omit the prior parameters ω. I-A1 Amortized Variational Inference and the ELBO The core idea of Variational Inference is to approximate the intractable posterior pγ(ϕ∣D(i))p_γ(φ D^(i)) with a parameterized variational distribution qψq_ψ. In this paper, qψq_ψ is chosen as a multivariate Gaussian with diagonal covariance. To avoid learning a separate set of variational parameters for each dataset in our meta-dataset, we employ amortized variational inference. We introduce an inference network (encoder) EψE_ψ, with weights ψ, that takes as input a dataset and outputs the sufficient statistics of the approximate posterior. Specifically, taking into account the train-test split within each dataset D(i)=(Dtr(i),Dte(i))D^(i)=(D_tr^(i),D_te^(i)), the encoder processes the training portion DtrD_tr and provides a mean vector and the diagonal elements of the covariance: Eψ(Dtr(i))=((i),(i))E_ψ(D_tr^(i))=( μ^(i), σ^(i)). The variational distribution qψ(ϕ∣Dtr(i))q_ψ(φ D^(i)_tr) is ((i),diag((i)))N( μ^(i),diag( σ^(i))). In Variational Inference, the Kullback-Leibler (KL) divergence between the approximate and the exact posterior is minimized. Mathematically, this is equivalent to maximizing (w.r.t. the variational parameters ψ) the Evidence Lower Bound (ELBO) [5], that for a single dataset is defined by: ELBO(i)(γ,ψ)=qψ(ϕ∣Dtr(i))[logpγ(te(i)∣te(i),ϕ)]−DKL(qψ(ϕ∣Dtr(i))∥(,)),ELBO^(i)(γ,ψ)=E_q_ψ(φ D_tr^(i)) [ p_γ(y_te^(i) _te^(i),φ) ]\\ -D_KL (q_ψ(φ D_tr^(i)) (0,I) ), (8) Crucially, the ELBO is a lower bound of the log-evidence logpγ((i)∣(i)) p_γ(y^(i) ^(i)), and the gap is the KL divergence between qψ(ϕ∣Dtr(i))q_ψ(φ D_tr^(i)) and pγ(ϕ∣Dtr(i))p_γ(φ D_tr^(i)), namely: logpγ((i)∣(i))=ELBO(i)(γ,ψ)+DKL(qψ(ϕ∣Dtr(i))∥pγ(ϕ∣Dtr(i))) p_γ(y^(i) ^(i))=ELBO^(i)(γ,ψ)+\\ D_KL (q_ψ(φ D_tr^(i)) p_γ(φ D_tr^(i)) ) (9) Therefore, by jointly maximizing the ELBO with respect the variational parameters ψ and the generative model parameters γ, it is possible to recover (i) a probabilistic model that maximizes the data evidence, and thus is optimal in the (marginal) maximum likelihood sense and (i) as a by-product, an approximate expression of the otherwise intractable posterior (7). Remark 2 The ELBO expression (8) balances two competing objectives. The first term is the expected log-likelihood (reconstruction error), which encourages the latent representation to accurately reproduce the output dynamics via the lifting function. The second term acts as a regularizer, preventing the approximate posterior qψ(ϕ∣D(i))q_ψ(φ D^(i)) from deviating too much from the prior. I-A2 Meta-training By maximizing the ELBO across the dataset distribution p(D)p(D), we learn a generative model of p(D)p(D) governed by a low number of latent factors ϕφ, aligning with the modeling assumption (1). Overall, the meta-ELBO objective is: (γ,ψ)=p(D)[qψ(ϕ∣Dtr)[logpγ(te∣te,ϕ)]−βDKL(qψ(ϕ∣Dtr)∥(,)).J(γ,ψ)=E_p(D)[E_q_ψ(φ D_tr) [ p_γ(y_te _te,φ) ]\\ - _KL (q_ψ(φ D_tr) (0,I) ). (10) Here, the hyperparameter β (standard in β-VAEs [10]) scales the relative weight of the KL-divergence regularization. Tuning β proved critical in our experiments to actively prevent posterior collapse [14] and ensure a meaningful latent representation. During training, the outer expectation over the dataset distribution is approximated via Monte Carlo sampling across mini-batches of size b=128b=128. For the inner expectation, we draw a single stochastic sample per dataset from the approximate posterior qψq_ψ. To ensure that gradients can back-propagate through this random sampling step into the encoder parameters ψ, we employ the standard reparameterization trick [11]. This enables end-to-end gradient-based optimization of the joint meta-parameters: γ^,ψ^=argmaxγ,ψ(γ,ψ). γ, ψ= _γ,ψJ(γ,ψ). (11) The full meta-training framework pipeline is illustrated in Figure 1. Figure 1: Overview of the meta-training pipeline. The encoder EψE_ψ maps a training dataset to (μ, log σ); ϕφ is sampled via reparameterization and passed through the lifting function PγP_γ and base architecture F to produce y^te y_te. The KL term (orange, left) regularizes only at training time; the reconstruction term (orange, right) uses ytey_te. I-B Meta-Testing, MAP and Laplace Approximation At meta-test time, the framework is evaluated on a new, unseen system generating a dataset Dnew=Dtrnew,DtenewD new=\D new_tr,D new_te\. While the encoder Eψ^E_ ψ provides an immediate, zero-shot Gaussian approximation of the posterior qψ^(ϕ∣Dtrnew)q_ ψ(φ D new_tr), relying solely on this amortized network can lead to sub-optimal accuracy and miscalibrated uncertainty for out-of-distribution dynamics. To refine the prediction, we first perform Maximum A Posteriori (MAP) estimation, using the encoder’s mean output as an initialization point. Specifically, we obtain the MAP point estimate ϕMAP _MAP as: ϕ^=argminϕℒnlp(ϕ), φ= _φL nlp(φ), (12a) where the loss ℒnlp(ϕ)=‖trnew−F(trnew;Pγ^(ϕ))σe‖2+‖ϕ‖22L nlp(φ)= y_tr new-F(u new_tr;P_ γ(φ)) _e ^2+ φ _2^2 is the negative log-posterior, up to an additive factor. The optimization mirrors the task-specific fine-tuning of the deterministic framework [8], yielding a single point estimate for the low-dimensional parameter. Then, to recover an accurate uncertainty measure, we apply the Laplace approximation around the MAP estimate. By taking the second-order Taylor expansion of the log-posterior at ϕMAP _MAP, we approximate the posterior as a Gaussian (ϕMAP,ΣLaplace)N( _MAP, _Laplace), where the covariance matrix is the inverse of the Hessian of the negative log-posterior loss: ΣLaplace=(∂2ℒnlp(ϕ)∂ϕ2|ϕ=ϕMAP)−1. _Laplace= ( . ∂^2L nlp(φ)∂φ^2 |_φ= _MAP )^-1. (12b) Because the dimensionality of the manifold nϕn_φ is small, computing and inverting this Hessian is computationally feasible. This yields a principled, fast, and mathematically grounded quantification of uncertainty. Remark 3 The Laplace approximation in (12b) is a local approximation: it fits a single Gaussian around the mode ϕMAP _MAP found by (12a), and by construction cannot represent multimodality or heavy-tailed structure in the true posterior pγ(ϕ∣)p_γ(φ ). This is a meaningful restriction whenever distinct latent configurations ϕφ can explain the same observed trajectory equally well, a situation that is not merely theoretical for hysteretic systems such as Bouc–Wen, where several combinations of the shape parameters (γ,β,δ,ν)(γ,β,δ,ν) can produce near-identical force–deformation loops, particularly under short or narrowband excitation. In such regimes the single-mode Gaussian in (12b) may underestimate total predictive uncertainty by ignoring probability mass around competing modes, or misrepresent the local curvature if ϕMAP _MAP lands near a saddle between two comparably good optima. The empirical coverage results in Fig. 5(b), evaluated on the nominal-parameter benchmark trajectory, do not stress-test this failure mode directly, since the adaptation data there is long enough (L≥100L≥ 100) to identify a well-separated mode; shorter or less-informative excitations are where multimodal collapse is most likely. Addressing this is a natural direction for future work, e.g., replacing the single Laplace posterior with a mixture of Laplace approximations fit at multiple MAP restarts, in the spirit of [7], or a full-covariance variational posterior that is not restricted to a diagonal Gaussian at inference time. IV Numerical Examples The proposed variational framework is evaluated on two examples: a synthetic toy example concerning static regression of sinusoidal functions, and a dynamical system identification problem using the Bouc–Wen benchmark [15]. In both cases, the variational approach is compared against its deterministic counterpart [8]. The primary objective is to demonstrate that the probabilistic formulation achieves comparable predictive performance while additionally providing mathematically rigorous, well-calibrated uncertainty estimates for both the output predictions and the low-dimensional latent parameters. The code was implemented in JAX and is publicly available in the Github repository https://github.com/mattrufolo/VAE-sysid-neural-manifold. All computations were executed on a server equipped with an NVIDIA RTX 3090 GPU. The network structures used in the numerical examples are summarized in Table I. TABLE I: Architectures employed in the numerical examples. Benchmark Encoder (EψE_ψ) Lifting Function (PγP_γ) Base Architecture (F) Sines DeepSets: ∙ Shared MLP: 2 hidden layers, ReLU ∙ Aggregation: Mean pooling ∙ Output MLP: 1 hidden layer, ReLU → parallel Dense layers for μ and log σ Affine Map: Dense layer (no activation) MLP: Hidden dimensions: [64, 32, 16, 1], ReLU activations Bouc-Wen RNN + MLP: ∙ RNN: Bi-directional GRU (128 units) ∙ Aggregation: Average pooling ∙ Output MLP: 1 hidden layer (128 units), tanh → parallel Dense layers for μ and log σ (nϕ=20n_φ=20) Affine Map: Dense layer (no activation), expanding nϕ=20n_φ=20 to nθ=244n_θ=244 Neural State-Space: ∙ nx=3,nu=1,ny=1n_x=3,n_u=1,n_y=1 ∙ Nf,NgN_f,N_g: MLPs with 1 hidden layer (16 units), tanh activations IV-A Static regression The meta-training dataset comprises 10001000 sine wave regression tasks of length N=100N=100. For each dataset D(i)D^(i), outputs are generated from the sinusoidal function over the input with an additive white noise, via (i)=A(i)sin((i)+α(i))+(i)y^(i)=A^(i) (u^(i)+α^(i))+e^(i), with inputs (i)∼(−5,5)u^(i) (-5,5), amplitudes A(i)∼(0.5,2.0)A^(i) (0.5,2.0), and phases α(i)∼(0,π)α^(i) (0,π). The term e(i)e^(i) denotes additive white noise with 1%1\% relative amplitude. Each sequence is evenly split into training and test sets D(i)=Dtr(i),Dte(i)D^(i)=\D_tr^(i),D_te^(i)\ (Ntr=Nte=50N_tr=N_te=50). Because the underlying generative process has exactly two degrees of freedom (amplitude and phase), the latent manifold dimensionality is fixed to the theoretical minimum nϕ=2n_φ=2. Architectural details are provided in Table I. We employ a permutation-invariant DeepSets encoder [19], aligning with the structural prior of static regression. Meta-training is performed by minimizing (5) and (11) for the deterministic and probabilistic settings, respectively, using the Adam optimizer with a learning rate of 10−210^-2 over 10,00010,000 iterations. To ensure a fair comparison, all hyperparameters and architectural choices are shared between the two frameworks, with the sole addition of the ELBO regularization weight β=0.1β=0.1 for the probabilistic setting. IV-A1 Unconditional Generation from the Prior We first assess the ability of the learned probabilistic model to generate synthetic datasets that resemble the meta-dataset. To this aim, we draw 3030 latent vectors samples ϕ~(i) φ^(i) from the prior (,)N(0,I) and evaluate the functions F(⋅,Pγ^(ϕ~(i)))F(·,P_ γ( φ^(i))) over a dense grid dsu_ds of 10001000 linearly spaced points in the domain [−5,5][-5,5]. As shown in Figure 2, these unconditionally-sampled parameters generate plausible trajectories that exhibit a clear sinusoidal pattern. For a quantitative verification, we fit a parametric sinusoidal model to these generated curves to extract their empirical amplitudes and phases. We obtained generated amplitudes A^∈[0.332,2.101] A∈[0.332,2.101] and phases φ^∈[0.006,3.165] ∈[0.006,3.165]. These empirical ranges are aligned with the ranges [0.5,2.0][0.5,2.0] and [0,π][0,π] of the uniform distributions of A and α used to generate the meta-training set. This illustrates that the learned prior effectively captures the true generative process of the sines meta-dataset. Figure 2: Unconditional generation of 3030 trajectories from the prior ϕ∼(,)φ (0,I). The outputs exhibit the expected sinusoidal pattern characterizing the meta-training distribution. IV-A2 Performance and Uncertainty in Testing During meta-testing, the shared meta-learned representation is used to find the optimal low-dimensional parameters ϕφ that fit the training portion of an unseen dataset. This adaptation is implemented for the deterministic and probabilistic cases by solving Equations (4b) and (12), respectively. This procedure was executed across 100 new test trajectories. To evaluate the robustness of the adaptation, we assessed the models under three different limited-data regimes by varying the length of the training sequence (Ntr=1N_tr=1, Ntr=2N_tr=2, and Ntr=5N_tr=5). To measure the true functional accuracy of the adapted models, the predictive Root Mean Square Error (RMSE) is computed over the same high-resolution grid dsu_ds utilized to generate Fig. 2. The quantitative results, presented in Table I as the mean and standard deviation across all 100100 test trajectories, demonstrate that the predictive performance of the proposed VAE formulation remains highly competitive with the baseline deterministic approach. TABLE I: Continuous functional RMSE for the Sines benchmark evaluated over a dense 1000-point grid. Results are reported as the mean ± standard deviation across 100 unseen test trajectories for different training sequence lengths (N). Training Length (N) Deterministic RMSE VAE RMSE N=2N=2 0.548 ± 0.723 0.634 ± 0.909 N=4N=4 0.121 ± 0.273 0.086 ± 0.180 N=10N=10 0.062 ± 0.124 0.026 ± 0.044 Figure 3: Laplace uncertainty bounds on a novel test sequence. Top: Adaptation with 5 points (N=10N=10) yields a tight fit. Bottom: Adaptation with 1 point (N=2N=2) results in a poor pointwise fit, yet the ±3σ± 3σ bounds correctly maintain the inherent sinusoidal structure learned from the prior. The uncertainty quantification capability of the probabilistic meta-learning approach is visualized in Figure 3 for Ntr=5N_tr=5 and Ntr=1N_tr=1. In these plots, the discrete blue markers represent the sparse observations used for task adaptation, while the solid lines represent the continuous functional predictions evaluated over the high-resolution grid. For Ntr=5N_tr=5, the training points are sufficient to accurately anchor the phase and amplitude, yielding a tight fit and narrow uncertainty bounds. Conversely, for Ntr=1N_tr=1, adaptation relies on a single data point, leading to a predictably poor pointwise fit. However, the uncertainty bands in this data-starved regime remarkably retain a clear, bounded oscillatory structure. This confirms that the meta-learned probabilistic model successfully captured the fixed-frequency pattern across tasks, even when the data is insufficient to pinpoint the exact phase or amplitude. IV-B Bouc-Wen benchmark We perform a comparison between the probabilistic meta-learning framework introduced in this paper and its deterministic counterpart on the same Bouc-Wen benchmark [15] utilized in [8]. The benchmark models the vibrations of a 1-DoF hysteretic mechanical system whose dynamics are described in state-space form by: [p˙(t)v˙(t)z˙(t)] bmatrix p(t)\\ v(t)\\ z(t) bmatrix =[v(t)1mL(u(t)−kLp(t)−cLv(t)−z(t))αv(t)−β(γ|v(t)||z(t)|ν−1+δv(t)|z(t)|ν)] = bmatrixv(t)\\ 1m_L (u(t)-k_Lp(t)-c_Lv(t)-z(t) )\\ α v(t)-β(γ|v(t)||z(t)|^ν-1+δ v(t)|z(t)|^ν) bmatrix (13a) y(t) y(t) =p(t), =p(t), (13b) where p(t)p(t) (m) is the measured position, v(t)v(t) (m/s) the unmeasured velocity, z(t)z(t) (N) an unmeasured hysteretic force, and u(t)u(t) (N) the known input force. TABLE I: Bouc–Wen coefficients: nominal and uniform min/max bounds. mLm_L cLc_L kLk_L α β γ δ ν nom 2 10 5.0×1045.0× 10^4 5.0×1045.0× 10^4 1000 0.8 -1.1 1 min 1 5 2.5×1042.5× 10^4 2.5×1042.5× 10^4 500500 0.5 -1.5 1 max 3 15 7.5×1047.5× 10^4 7.5×1047.5× 10^4 4500 0.9 -0.5 1 The benchmark provides a noise-free test dataset generated by simulating (13) with the nominal coefficients in Table I, excited by a random-phase multisine with length N=8192N=8192, sampled at frequency fs=750f_s=750 Hz, that excites the frequency range [5,150][5,150] Hz, and has as root mean square of 5050 N. A training dataset (N=40960N=40960) is also provided, where the system is excited by a distinct multisine signal with the same properties, and the output is corrupted by additive Gaussian noise (bandwidth [0−375][0-375] Hz, with amplitude 8⋅10−38· 10^-3). Following [3], we extend the original scope of the benchmark and use (13) to generate a meta-training dataset composed of random multisine trajectories from systems with coefficients uniformly sampled within the min-max bounds in Table I. Crucially, to prevent meta-overfitting, the reduced-complexity architecture is optimized exclusively on this synthetic meta-dataset. The original benchmark dataset (simulated at nominal values) is strictly reserved as meta-test set to evaluate the final adaptation performance. IV-C Meta-Training the System Manifold To enable robust identification in the low-data regime, we meta-learn a low-dimensional manifold (nϕ=20n_φ=20) embedded within the full-complexity parameter space (nθ=244n_θ=244). The specific configurations for the encoder EψE_ψ, lifting function PγP_γ, and base architecture F are detailed in Table I. The lifting network is affine: Pγ(ϕ)=Vϕ+θbiasP_γ(φ)=Vφ+ _ bias, with tunable parameters γ=vec(V,θbias)∈ℝ5124γ=vec(V, _ bias) ^5124. The encoder EψE_ψ, which comprises nψ=136,340n_ψ=136,340 parameters, utilizes a Bi-Gated Recurrent Unit (GRU) [4] followed by an average-pooling mechanism to aggregate input-output sequences into a compact latent representation. We use as base architecture the same linear-plus-residual neural state-space model from [8]: xk+1 x_k+1 =Axk+Bux+Nf(xk,uk;Wf) =Ax_k+Bu_x+N_f(x_k,u_k;W_f) (14a) yk y_k =Cxk+Ng(xk;Wg), =Cx_k+N_g(x_k;W_g), (14b) with nθ=244n_θ=244 parameters (further details in Table I). For the probabilistic setting, we use β=0.1β=0.1 in (10). The meta-dataset D consists of input-output sequences of length N=2000N=2000, generated on-the-fly by sampling Bouc–Wen systems within the bounds of Table I and exciting them with random multisines. Optimization is performed using Adam, running for 200,000200,000 iterations with a learning rate decayed from 1⋅10−41· 10^-4 to 1⋅10−51· 10^-5 via cosine scheduling. Because datasets are generated dynamically, the meta-training is performed over 25.6 million unique trajectories and took ∼ 25.3 hours. IV-D Task-Specific Adaptation and Uncertainty Analysis The primary objective of this numerical evaluation is to compare the proposed variational framework against its deterministic counterpart [8]. For a broader analysis comparing the baseline manifold meta-learning strategy against other state-of-the-art system identification methods, we refer the reader to [8]. Here, we focus exclusively on assessing how the probabilistic extension affects predictive accuracy and enables uncertainty quantification. In line with [8], the performance is evaluated in terms of the FITFIT index. To systematically evaluate the proposed variational framework in the low-data regime, we execute a comprehensive Monte Carlo adaptation study. Following the experimental setup established in [3], we extract nmc=100n_ mc=100 independent subsequences of varying lengths L∈100,200,400,500,600,800,1000,2000,3000,4000,5000L∈\100,200,400,500,600,800,1000,2000,3000,4000,5000\ from the official 40960-sample Bouc–Wen training dataset. For each subsequence, we adapt the meta-learned model to the specific system instance by performing MAP estimation over the low-dimensional latent space to identify the task-specific parameter ϕMAP _MAP. Our evaluation serves two primary objectives. First, we demonstrate that transitioning to a probabilistic generative formulation does not sacrifice the predictive accuracy of purely deterministic methods. Figure 4 reports the test FITFIT indices evaluated on the benchmark dataset. The results confirm that the variational framework remains highly competitive with its deterministic counterpart across all data regimes. To provide comprehensive context, we include as baselines an LTI and a full-complexity model (13) trained on the complete 40960-sample training dataset; detailed formulations for these baselines are in [8]. Figure 4: Deterministic and VAE reduced-order models: test FIT vs. training sequence length L. Second, and most critically, the variational approach equips the identified models with a mathematically rigorous quantification of epistemic uncertainty. By computing the inverse Hessian of the adaptation loss at ϕMAP _MAP via the Laplace approximation, we extract a localized Gaussian posterior for the latent system parameters. To illustrate the practical utility of this formulation, Figure 5(a) displays a benchmark test prediction alongside its ±3σ± 3σ confidence intervals, obtained after adapting the meta-model via MAP estimation with L=100L=100. The generated bounds successfully encapsulate the true system dynamics, actively compensating for regions where the mean prediction slightly deviates from the ground truth. Furthermore, to explicitly validate the statistical consistency of these bounds, Figure 5(b) evaluates the empirical coverage against the expected theoretical probability across various adaptation lengths, confirming that the framework yields properly calibrated uncertainty estimates. (a) Time-domain prediction with ±3σ± 3σ uncertainty bounds. (b) Empirical coverage probability vs. adaptation length L. Figure 5: Performance and uncertainty validation of the variational framework. (a) Benchmark test prediction displaying calibrated confidence intervals after MAP adaptation and Laplace uncertainty computation (L=100L=100). (b) Evaluation of empirical coverage against theoretical probability across different adaptation regimes. Finally, to explicitly analyze how this uncertainty scales with data availability, Figure 6 plots the sequence-averaged predictive standard deviation as a function of the adaptation length L. The mean standard deviation exhibits a clear, monotonic decay as L increases. This behavior perfectly aligns with Bayesian principles: as the model is provided with more task-specific information, the epistemic uncertainty inherently shrinks, yielding progressively more confident and reliable parameter estimates. Figure 6: Evolution of model output uncertainty as a function of the adaptation sequence length. V Conclusion This paper introduced a fully probabilistic framework of manifold meta-learning for nonlinear system identification. By formulating the amortized optimization procedure via Variational Inference, the architecture learns a generative prior over a low-dimensional system manifold. While maintaining the high predictive accuracy of its deterministic counterpart, the proposed variational approach provides calibrated epistemic uncertainty bounds over both the latent physical parameters and the predicted trajectories. The framework was systematically validated on an interpretable synthetic regression task and the Bouc-Wen dynamical benchmark. In both cases, the model demonstrated a consistent, monotonic reduction in predictive variance as adaptation data increased, ensuring robust and reliable identification even in severely low-data regimes. Future work will focus on integrating these data-driven uncertainty bounds into risk-aware control applications, such as stochastic Model Predictive Control, and embedding physics-informed structural priors into the generative latent space to further improve interpretability and performance. References [1] B. Amos (2023) Tutorial on amortized optimization. Foundations and Trends in Machine Learning 16 (5), p. 592–732. Cited by: §I-B. [2] D. M. Blei, A. Kucukelbir, and J. D. McAuliffe (2017) Variational inference: a review for statisticians. Journal of the American statistical Association 112 (518), p. 859–877. Cited by: §I-A. [3] A. Chakrabarty, G. Wichern, V. M. Deshpande, A. P. Vinod, K. Berntorp, and C. R. Laughman (2025) Meta-learning for physically-constrained neural system identification. Neurocomputing, p. 130945. Cited by: §IV-B, §IV-D. [4] J. Chung, C. Gulcehre, K. Cho, and Y. Bengio (2014) Empirical evaluation of gated recurrent neural networks on sequence modeling. arXiv preprint arXiv:1412.3555. Cited by: §IV-C. [5] C. Doersch (2016) Tutorial on variational autoencoders. arXiv preprint arXiv:1606.05908. Cited by: §I-A1. [6] Z. Du, H. Balim, S. Oymak, and N. Ozay (2023) Can transformers learn optimal filtering for unknown systems?. IEEE Control Systems Letters 7, p. 3525–3530. Cited by: §I. [7] R. Eschenhagen, E. Daxberger, P. Hennig, and A. Kristiadi (2021) Mixtures of laplace approximations for improved post-hoc uncertainty in deep learning. Cited by: Remark 3. [8] M. Forgione, A. Chakrabarty, D. Piga, M. Rufolo, and A. Bemporad (2025) Manifold meta-learning for reduced-complexity neural system identification. under revision to Automatica. Cited by: §I, §I, §I-B, §I-B, 12, §IV-B, §IV-C, §IV-D, §IV-D, §IV. [9] M. Garnelo, D. Rosenbaum, C. Maddison, T. Ramalho, D. Saxton, M. Shanahan, Y. W. Teh, D. Rezende, and S. A. Eslami (2018) Conditional neural processes. In International conference on machine learning, p. 1704–1713. Cited by: §I, §I. [10] I. Higgins, L. Matthey, A. Pal, C. Burgess, X. Glorot, M. Botvinick, S. Mohamed, and A. Lerchner (2017) Beta-VAE: learning basic visual concepts with a constrained variational framework. In International Conference on Learning Representations, Cited by: §I-A2. [11] D. P. Kingma and M. Welling (2014) Auto-encoding variational bayes. In International Conference on Learning Representations, Vol. 2014. Cited by: §I, §I-A2. [12] B. Lakshminarayanan, F. Dettú, C. R. Rojas, and S. Formentin (2025) Inverse supervised learning of controller tuning rules. Automatica 178, p. 112356. Cited by: §I. [13] J. Lee, Y. Lee, J. Kim, A. Kosiorek, S. Choi, and Y. W. Teh (2019) Set transformer: a framework for attention-based permutation-invariant neural networks. In International conference on machine learning, p. 3744–3753. Cited by: §I. [14] J. Lucas, G. Tucker, R. B. Grosse, and M. Norouzi (2019) Understanding posterior collapse in generative latent variable models. In Deep Generative Models for Highly Structured Data Workshop, ICLR, Cited by: §I-A2. [15] J. Noel and M. Schoukens (2016) Hysteretic benchmark with a dynamic nonlinearity. In Workshop on nonlinear system identification benchmarks, p. 7–14. Cited by: §IV-B, §IV. [16] G. Pillonetto, A. Aravkin, D. Gedon, L. Ljung, A. H. Ribeiro, and T. B. Schön (2025) Deep networks for system identification: a survey. Automatica 171, p. 111907. Cited by: §I. [17] M. Rufolo, D. Piga, and M. Forgione (2025) Distributionally robust minimization in meta-learning for system identification. IEEE Control Systems Letters 9, p. 1652–1657. Cited by: §I. [18] M. Rufolo, D. Piga, G. Maroni, and M. Forgione (2025) Enhanced transformer architecture for in-context learning of dynamical systems. In 2025 European Control Conference (ECC), p. 819–824. Cited by: §I, §I. [19] M. Zaheer, S. Kottur, S. Ravanbakhsh, B. Poczos, R. R. Salakhutdinov, and A. J. Smola (2017) Deep sets. Advances in neural information processing systems 30. Cited by: §I, §IV-A.