System identification for complex dynamical systems: a survey

2.1 Fundamental knowledge

Before introducing classical algorithms for single-system identification, it is necessary to clarify several fundamental concepts, including the main classification schemes, essential analytical tools, and the historical evolution of the field. These elements provide the theoretical foundation and "global picture'' on which the more detailed methods in the next subsection are built.

Classification of single-system identification. Classical single-system identification is most often formulated through prediction-error or likelihood-based estimation for parametric input–output and state-space models. In many widely used settings (e.g., linear regression models, ARX/ARMAX with fixed regressors, or locally approximated nonlinear models), the resulting estimators reduce to least-squares-type problems, while more general nonlinear and closed-loop cases require iterative optimization and dedicated treatments. With modern sensing and high-dimensional measurements, the number of candidate regressors (or basis functions) can be large, leading to redundant representations and ill-conditioned estimation. This motivates sparsity-promoting identification, which seeks parsimonious model structures by enforcing sparsity on parameters (or on a chosen function library), retaining only a small set of active terms that explain the data within an acceptable error level.

According to their methodological principles, the single-system identification methods reviewed in Section 2.2 can be grouped into three families that organize the methods emphasized in this review:

1. Regularized regression methods: These operate in a static regression setting and control model complexity via sparsity-promoting or shrinkage penalties. Greedy algorithms such as OMP^[12] seek sparse supports, whereas penalized estimators such as LASSO and ridge/elastic-net regression^[13,14,15] yield efficient solutions with explicit bias–variance control.

2. Sparse Bayesian methods: These introduce hierarchical probabilistic priors to induce sparsity and quantify uncertainty. Representative examples include the RVM^[16,17], classical SBL^[18,19], and block-structured extensions^[20,21,22].

3. Recursive state-space and filtering methods: These work in a dynamical state-space framework and perform online state and parameter estimation via Bayesian filtering. The KF and its nonlinear extensions (EKF, UKF)^{[23,24,25,26]} are canonical examples, and they form the basis of augmented-state joint state–parameter identification reviewed in Section 2.2^[27,28].

These three families provide an organizing map for the single-system methods emphasized in this review; the specific algorithms in the next subsection can be seen as concrete instances or refinements within this organizing framework. Beyond these three families, recent literature has rapidly expanded a fourth line of development, namely structured recursive identification methods, which are especially important for multivariable, nonlinear, and colored-noise systems. This line includes auxiliary-model identification, multi-innovation identification, hierarchical identification, filtering-based identification, and coupled identification, and will be discussed separately at the end of this section.

Classical backbone and scope of this review. Single-system identification also has a well-established classical backbone, including PEM/ML estimation for parametric input–output and state-space models^[29], subspace identification for state-space realization from input–output data^[30], and closed-loop identification techniques that explicitly account for feedback-induced bias and reduced excitation^[31,32]. These paradigms provide rigorous foundations for consistency, identifiability, and statistical efficiency, and remain central in industrial practice. In this review, we treat them as essential background and focus our technical exposition on sparsity-driven regression, SBL, and recursive filtering.

Analytical foundations. Three analytical tools underpin most identification algorithms: regression analysis, gradient-based optimization, and Bayesian inference.

For a linear model

where $$ y(k)$$ denotes the measured output at sample index $$ k$$ , $$ \boldsymbol{\varphi}(k)$$ is the regressor vector constructed from input, output, or state-related signals, and $$ \boldsymbol{\theta}$$ is the unknown parameter vector to be estimated. Parameter estimation is typically achieved by minimizing a loss function, such as the sum of squared residuals or the negative log-likelihood. Nonlinear cases often adopt iterative optimization schemes such as gradient descent:

where $$ \boldsymbol{\theta}_{k}$$ denotes the parameter estimate at iteration $$ k$$ , $$ \boldsymbol{\theta}_{k+1}$$ is the updated estimate, $$ \alpha$$ is the learning rate controlling the update step, $$ J(\boldsymbol{\theta}_{k})$$ is the objective function evaluated at $$ \boldsymbol{\theta}_{k}$$ , and $$ \nabla_{\theta}J(\boldsymbol{\theta}_{k})$$ denotes the gradient of the objective with respect to the parameters.

From a probabilistic perspective, sparse modeling can be expressed in Bayesian form:

where $$ \hat{\mathbf{s}}$$ is the estimated sparse coefficient vector, $$ \mathbf{y}$$ denotes the observation vector, $$ \mathbf{A}$$ is the sensing or dictionary matrix, $$ \|\mathbf{y}-\mathbf{A}\mathbf{s}\|_2^2$$ measures the data-fitting error, $$ \|\mathbf{s}\|_0$$ counts the number of nonzero entries in $$ \mathbf{s}$$ , and $$ \lambda \ge 0$$ controls the strength of sparsity regularization. This formulation corresponds to a maximum a posteriori (MAP) estimation problem under suitable priors. Different priors yield different regularization effects—for instance, Laplace priors correspond to LASSO-type sparsity, while Gaussian–Gamma priors lead to Bayesian RVM formulations^[33,17,34].

Research progress. The development of single-system identification is often described along three intertwined threads. The classical parametric era (1960s–1990s) established PEM/ML formulations for input–output and state-space models, together with frequency-domain methods and subspace identification, forming the backbone of industrial practice. The regularization and Bayesian era (late 1990s–2010s) emphasized robustness and uncertainty quantification, including regularized least-squares, kernel and Gaussian process viewpoints, and Bayesian identification. The sparsity-aware era (mid-2000s onward) introduced compressed-sensing-inspired methods and SBL (e.g., LASSO/OMP and RVM/SBL) to enable parsimonious structure selection in high-dimensional settings, and to connect naturally to equation discovery and modern learning pipelines.

Today, single-system identification remains a cornerstone for adaptive control, signal processing, and intelligent cyber–physical systems, providing data-driven models for prediction, diagnosis, and decision-making under noise, constraints, and limited excitation.

2.2 Classic algorithms

This subsection reviews representative classical algorithmic families for single-system identification, structured according to the three families outlined above. We first consider sparse regression methods (greedy and convex) for static regression, then sparse Bayesian formulations that provide probabilistic sparsity and uncertainty quantification, and finally recursive state-space and filtering approaches for dynamic, online estimation. These methods form the classical backbone of single-system identification. More recent structured recursive identification methods for multivariable, nonlinear, and colored-noise systems—including auxiliary-model, multi-innovation, hierarchical, filtering-based, and coupled identification schemes—will be discussed separately in the next subsection.

Orthogonal matching pursuit identification. OMP is a prototypical greedy method for sparse approximation. It was initially proposed as a recursive function approximation scheme with applications to wavelet decomposition^[12] and later became a cornerstone algorithm in sparse representation and compressed sensing. In each iteration, OMP selects the dictionary atom most correlated with the current residual and orthogonally projects the data onto the span of the selected atoms. In this way, a sparse representation is constructed through a sequence of locally optimal decisions. Although global optimality is not guaranteed, the method is computationally efficient and has proven effective in high-dimensional regimes.

The modern understanding of OMP is closely tied to the theory of compressed sensing. Sparse recovery problems are typically formulated as underdetermined linear systems with a sparsity prior, $$ A\in\mathbb{R}^{m\times n}$$ , $$ x\in\mathbb{R}^{n}$$ , $$ b\in\mathbb{R}^{m}$$ , where $$ m <n$$ and infinitely many solutions satisfy $$ Ax=b$$ . Imposing that most entries of $$ x$$ are zero defines a well-posed recovery objective—one seeks the sparsest coefficient vector consistent with the data^[35]. Within this setting, a series of theoretical results established that greedy schemes such as OMP can recover the true sparse solution under suitable conditions on the sensing matrix. Tropp^[36] analyzed OMP under mutual-coherence assumptions, while subsequent work refined recovery guarantees in terms of the Restricted Isometry Property (RIP) and related conditions^[20,37]. These developments positioned OMP as a practical alternative to convex $$ \ell_1$$ minimization, offering similar recovery performance with lower per-iteration cost.

Beyond classical signal reconstruction, OMP has been widely adopted for variable selection in linear regression and for model structure discovery in identification tasks. For instance, in high-dimensional settings with $$ p\gg m$$ , OMP iteratively builds a parsimonious subset of predictors by adding the most informative columns of the design matrix and re-estimating coefficients via least squares. In dynamical systems and system identification, the same mechanism is used to select a compact set of basis functions from a predefined library to explain observed dynamics. A representative example is the SINDy framework^[38], where OMP (or closely related sparse regressors) is employed to identify the active terms of governing equations from time-series data, thereby bridging sparse approximation and data-driven equation discovery.

Formally, treat the columns of $$ A=[\varphi_1, \dots, \varphi_n]$$ as atoms. At each iteration, OMP selects the atom with maximal absolute correlation with the current residual, augments the active set, recomputes the least-squares coefficients over the active set, and updates the residual by orthogonal projection. The process terminates when the residual norm falls below a prescribed tolerance or a sparsity budget is reached. A standard version of OMP is summarized in Algorithm 1 (we use $$ b$$ for the observation vector, $$ A$$ for the dictionary, and $$ x$$ for the coefficient vector; in spatio–temporal identification, $$ A$$ corresponds to the library $$ \Theta(\cdot)$$ and $$ x$$ to the sparse coefficient vector).

LASSO identification. LASSO, introduced by Tibshirani^[13], is a foundational method for sparse regression and has become a standard tool for system identification. Compared with ordinary least squares (OLS), which may suffer from high variance and poor interpretability in the presence of multicollinearity, LASSO introduces an $$ l_1$$ penalty that simultaneously shrinks coefficients and performs variable selection. In its Lagrangian form, the objective reads

where $$ Y$$ denotes the response vector, $$ X$$ the design matrix, $$ \beta$$ the regression coefficients, and $$ \lambda \ge 0$$ is a regularization parameter. Equivalently, LASSO can be posed as a constrained least-squares problem

where $$ \alpha$$ is an intercept term, $$ \widehat{\beta}$$ is the estimated coefficient vector, $$ N$$ is the number of samples, $$ p$$ is the number of predictors, $$ y_i$$ denotes the $$ i$$ -th response, $$ \mathbf{x}_i$$ is the corresponding predictor vector, and $$ t \ge 0$$ is a tuning parameter controlling the degree of shrinkage. This formulation makes explicit the "budget" on the $$ l_1$$ -norm of the coefficients.

Conceptually, LASSO may be viewed as a convex relaxation of $$ l_0$$ -penalized subset selection. Whereas ridge regression^[14] employs an $$ l_2$$ penalty to stabilize estimates without enforcing exact sparsity, the $$ l_1$$ penalty of LASSO yields sparse solutions while remaining computationally tractable. Subsequent work developed efficient algorithms and regularization paths, most notably least-angle regression (LARS)^[39] and coordinate-descent implementations such as $$\texttt{glmnet}$$ ^[40], which enable large-scale applications and extensions to generalized linear models.

On the theoretical side, a rich literature has established the statistical properties of LASSO in high-dimensional settings. Elastic Net^[15] combines $$ l_1$$ and $$ l_2$$ penalties to handle strongly correlated predictors, while the adaptive LASSO^[41] achieves oracle properties under suitable conditions. High-dimensional analyses, such as Bickel et al.^[42], provide recovery guarantees and error bounds, clarifying how sparsity, sample size, and noise level jointly determine estimation accuracy.

These developments have made LASSO-based sparse regression a natural building block for system identification. In many physical systems, only a few terms dominate the underlying dynamics, so that the governing equations are sparse with respect to a high-dimensional function library. This principle is operationalized in SINDy^[38], where the response $$ Y$$ corresponds to estimated temporal derivatives ($$ \dot{\mathbf{X}}$$ ), the design matrix $$ X$$ is instantiated as a dictionary of candidate functions ($$ \Theta(\mathbf{X})$$ ), and the coefficient vector $$ \beta$$ encodes a sparse set of active terms defining the governing equations. The resulting optimization problem takes the form

where $$ \dot{\mathbf{X}}$$ denotes the matrix of estimated temporal derivatives, $$ \Theta(\mathbf{X})$$ is the candidate-function library constructed from the state data $$ \mathbf{X}$$ , $$ \Xi$$ is the unknown coefficient matrix, $$ \hat{\Xi}$$ is the identified sparse coefficient matrix, $$ \|\dot{\mathbf{X}}-\Theta(\mathbf{X})\Xi\|_2^2$$ measures the regression error, $$ \|\Xi\|_1$$ promotes sparsity, and $$ \lambda \ge 0$$ controls the sparsity level. In the ODE setting, this formulation recovers the original SINDy model. For systems governed by partial differential equations, the same sparse-regression idea underlies PDE-FIND and its extensions^[43], which augment the candidate library with spatial derivatives and employ tailored regularization strategies. These PDE-oriented variants, together with robustness enhancements such as weak and ensemble formulations, are discussed in more detail in the PDE identification discussion. In system identification practice, LASSO thus plays two complementary roles. First, it offers an interpretable mechanism for model structure selection and complexity control, bridging classical regression and modern compressed sensing. Second, through frameworks such as SINDy, PDE-FIND, and their variants, LASSO-based sparse regression enables the discovery of parsimonious, physically meaningful dynamical models from noisy measurements, which can subsequently support prediction, control design, and model-based diagnostics.

Ridge regression identification. Ridge regression arises from the broader idea of Tikhonov regularization for ill-posed inverse problems, where a quadratic penalty is added to stabilize solutions of noisy or ill-conditioned linear systems^[44]. In the statistical literature, Hoerl and Kennard^[14] introduced ridge regression as a practical remedy for multicollinearity in multiple regression, showing that shrinking coefficients towards zero can substantially reduce estimation variance and improve prediction accuracy when predictors are strongly correlated.

Formally, ridge regression augments the least-squares objective with an $$ L_2$$ -norm penalty on the coefficients,

where $$ Y$$ denotes the response vector, $$ X$$ the design matrix, $$ \beta$$ the regression coefficients, and $$ \lambda \ge 0$$ the regularization parameter controlling the bias–variance trade-off. A larger $$ \lambda$$ enforces stronger shrinkage and improves numerical stability when $$ X^\top X$$ is ill-conditioned or nearly singular, whereas a smaller $$ \lambda$$ recovers the ordinary least-squares solution at the expense of increased variance. Because the problem remains strictly convex, ridge admits a unique closed-form solution and can be efficiently implemented via standard linear algebra routines.

In system identification, ridge-type regularization provides a natural baseline for regularized parameter estimation. Early work adopted it mainly as a pragmatic tool to handle noisy, collinear regressors and limited data records. More recently, ridge has been interpreted through the lens of reproducing-kernel Hilbert spaces and Gaussian priors, where the $$ L_2$$ penalty corresponds to enforcing smoothness or energy constraints on impulse responses or nonlinear kernels^[45,46]. Within this framework, regularized least-squares estimators become central components of kernel-based and Gaussian-process identification schemes for both linear and nonlinear dynamic systems, with provable robustness and uncertainty quantification.

An important extension of ridge regression is the Elastic Net (EN), proposed by Zou and Hastie^[15]. The Elastic Net combines the sparsity-inducing $$ L_1$$ penalty of LASSO with the grouping and stabilization properties of the $$ L_2$$ penalty:

where $$ \lambda_{1}\ge 0$$ controls the $$ L_1$$ sparsity-inducing penalty, $$ \lambda_{2}\ge 0$$ controls the $$ L_2$$ shrinkage penalty, and the other symbols follow the ridge and LASSO formulations above. This hybrid formulation addresses limitations of pure $$ L_1$$ or $$ L_2$$ regularization in high-dimensional settings, encouraging groups of correlated predictors to enter or leave the model together while retaining exact sparsity. For system identification problems with many candidate regressors (e.g., large function libraries, long memory structures, or multiple input channels), Elastic Net–type penalties offer a flexible compromise between variable selection and coefficient shrinkage.

Taken together, ridge regression and its Elastic Net extension form a key branch of regularization-based identification methods. They provide numerically stable and conceptually simple building blocks and integrate seamlessly with more advanced frameworks, such as kernel-based regularization and Bayesian system identification. In practice, the choice between LASSO, ridge, and Elastic Net is driven by data characteristics (e.g., degree of collinearity, signal-to-noise ratio, and dimensionality) and modeling goals (sparsity, grouping, or smoothness), and these methods collectively underpin many modern sparse and regularized system identification pipelines.

Relevance vector machine identification. SBL provides a probabilistic framework for sparse modeling and variable selection, and offers a principled alternative to convex regularization methods such as LASSO and ridge regression. A central milestone in this line of work is the RVM proposed by Tipping^[16], which recasts sparse regression and classification within a Bayesian kernel–machine formulation. In contrast to support vector machines (SVMs), which rely on margin maximization and typically yield dense sets of support vectors, the RVM enforces sparsity through hierarchical priors while retaining fully probabilistic predictions.

In the standard RVM formulation, the system output is modeled as

where $$ \Phi$$ is a design matrix composed of basis or kernel functions, $$ \mathbf{w}$$ is the weight vector, and $$ \mathbf{z}$$ denotes Gaussian noise. SBL imposes independent Gaussian priors with individual precision parameters on each weight, which are themselves governed by hyperpriors. Maximization of the marginal likelihood then drives many precisions to infinity, effectively pruning the corresponding weights and leaving only a small subset of "relevance vectors" active in the final model^[16]. Subsequent work developed efficient training algorithms for this framework, notably fast marginal-likelihood maximization schemes for large-scale problems^[17] and a unifying analysis of SBL/automatic relevance determination (ARD) that clarified its relationship to $$ l_1$$ -type regularization and basis selection^[18].

From a system identification perspective, the RVM/SBL family occupies an intermediate position between purely deterministic regularization and fully Bayesian state-space modeling. On the one hand, it inherits many advantages of convex sparse regression—automatic model-order selection, mitigation of multicollinearity, and robustness in high-dimensional settings. On the other hand, its Bayesian formulation provides uncertainty quantification on both weights and predictions, which is crucial for risk-aware control and experiment design. Empirically, RVM models tend to be substantially sparser than their SVM counterparts for comparable predictive accuracy, reducing evaluation cost and facilitating online adaptation.

These properties have motivated a growing body of applications in identification and diagnostics of dynamical systems. In fault diagnosis and condition monitoring, RVM-based classifiers have been employed for bearing fault detection and fault classification in rotating machinery, often outperforming SVMs under small-sample and noisy conditions^[47]. In environmental and hydrological modeling, multivariate RVMs have been used for spatio–temporal prediction of root-zone soil moisture, demonstrating the ability of SBL to capture complex spatial dependencies with compact models^[48]. More recently, RVM variants have been embedded in equation-discovery and multi-task frameworks, where Bayesian sparsity is exploited to identify governing terms shared across multiple operating regimes while accounting for task-specific noise^[49]. These developments highlight the suitability of RVM/SBL as a bridge between interpretable sparse models and probabilistic system identification.

In this review, we therefore regard RVM-based sparse Bayesian identification as a representative of the broader SBL paradigm. It complements greedy and convex-regularization approaches by providing (ⅰ) automatic relevance determination without manual tuning of sparsity levels, (ⅱ) probabilistic predictions with quantified uncertainty, and (ⅲ) typically very compact model representations. These characteristics make RVM/SBL particularly attractive for single-system identification in data-limited regimes and for multi-system scenarios where uncertainty and model parsimony are central concerns. Figure 2 illustrates the conceptual structure of RVM-based sparse identification.

Figure 2. Conceptual structure of relevance vector machine (RVM) identification.

Classic sparse Bayesian identification. SBL extends sparse modeling into a fully probabilistic framework by introducing hierarchical priors over model parameters. Building on the automatic relevance determination (ARD) idea underlying the RVM^[16], SBL generalizes Bayesian sparse regression to broad classes of linear inverse problems and basis–selection tasks^[18,19]. In contrast to explicit $$ L_{1}$$ -penalized formulations such as LASSO, sparsity here emerges from inference over latent hyperparameters that control the variance of individual coefficients.

A generic observation model can be written as

where $$ \Phi(t)$$ denotes a basis (dictionary) matrix, $$ \mathbf{W}_{ij}$$ is the unknown coefficient vector, and $$ \boldsymbol{\varepsilon}$$ is Gaussian noise. SBL assumes Gaussian priors over each coefficient with individual precision (inverse-variance) hyperparameters, which themselves follow suitably chosen hyperpriors (e.g., Normal–Jeffreys, Normal–Gamma, or Normal–Inverse Gaussian). Under Gaussian likelihood and prior assumptions, the posterior over $$ \mathbf{W}_{ij}$$ remains Gaussian with analytically tractable mean and covariance updates. Different hyperprior choices are equivalent to different implicit penalty structures on $$ \mathbf{W}_{ij}$$ , as summarized in Table 4. In this way, sparse Bayesian identification unifies classical sparse regression and probabilistic inference, achieving both automatic model selection and uncertainty quantification—two properties that are particularly valuable in noisy or data-limited identification settings.

In Table 4, $$ \mu_k$$ denotes a coefficient- or basis-specific variance-related hyperparameter, $$ q_k$$ denotes the corresponding scale variable in the Normal–Inverse Gaussian formulation, $$ N$$ is the number of observations, $$ K$$ is the number of candidate basis terms or coefficient groups, and $$ \alpha$$ and $$ \gamma$$ are hyperparameters controlling the shape and scale of the corresponding sparse Bayesian priors.

Early SBL work focused on establishing these hierarchical formulations and clarifying their relationship to convex and nonconvex sparsity penalties^[18,19]. Subsequent studies exploited SBL's robustness and automatic relevance determination to tackle increasingly challenging inverse problems, including compressive sensing, deconvolution, and structured sparsity. This progression naturally connects SBL to system identification, where ill-conditioned regressors, correlated excitations, and limited data are pervasive.

On the application side, SBL has been deployed across a wide range of domains relevant to identification and modeling. In block-sparse recovery, Sant et al.^[22] proposed a total variation (TV)–regularized SBL framework to handle unknown and nonuniform block structures, combining linear TV and region-aware logarithmic TV penalties to enhance sensitivity to block patterns and improve robustness under mixed sparsity. In high-mobility wireless communications, Wei et al.^[50] developed an off-grid SBL-based channel estimation scheme for orthogonal time–frequency space (OTFS) systems, mitigating grid mismatch via first-order linear approximations and achieving high-precision estimation under low-resolution sampling.

SBL has also emerged as a powerful engine for large-scale data-driven modeling. In radio environment mapping, Wang et al.^[51] combined SBL with environment-aware dictionaries and information-theoretic sensor placement, enabling efficient three-dimensional spectrum reconstruction. For biomedical signal decoding, Wang et al.^[52] proposed the SBLEST framework, which integrates spatial–temporal filtering and classification within a single SBL objective, automatically learning neurophysiologically interpretable EEG features and outperforming deep networks in small-sample regimes. In structural health monitoring, Wang et al.^[53] introduced an improved SBL scheme (iSBL) with Bayesian regularization to suppress overfitting and enhance generalization, delivering accurate modeling of nonstationary vibration signals for infrastructure safety assessment.

Despite these advances, classic SBL still faces challenges in high-dimensional and real-time identification scenarios. Inference can be computationally demanding, and hyperparameter learning often relies on hand-tuned update rules or empirical heuristics^[22,50]. Current research therefore explores faster variational or approximate-inference schemes, more expressive yet tractable hyperprior families, and self-adaptive hyperparameter learning, as well as extensions to multimodal fusion and online intelligent systems. Overall, SBL provides a flexible probabilistic backbone for sparse identification, sitting at the interface between traditional regularization, Bayesian learning, and modern large-scale inverse problems.

Block Bayesian identification. Within the SBL family, block sparse Bayesian learning (BSBL) generalizes element-wise sparsity to group-structured signals. Instead of treating coefficients as independent, BSBL assumes that coefficients within the same block are jointly relevant while sparsity is enforced across blocks. This idea was systematized by Zhang and Rao, who showed that many practical signals (e.g., multichannel or temporally correlated sources) are naturally block-sparse and that explicitly modeling intra-block correlation can substantially improve recovery performance^[20]. Their subsequent extension of SBL algorithms for block-sparse signals with intra-block correlation^[21] formalized the BSBL framework: each block is associated with its own covariance structure and hyperparameters, so that both the block-level sparsity pattern and the within-block correlation are inferred from data.

Following these foundational works, a series of studies focused on algorithmic efficiency and on relaxing assumptions about the block structure. Liu et al. proposed a fast marginalized BSBL algorithm that marginalizes out the coefficients and operates directly on hyperparameters, yielding significant computational savings for large-scale problems while preserving the advantages of BSBL in correlated settings^[54]. Fang et al.^[55] further introduced pattern-coupled SBL for recovery of block-sparse signals with unknown block patterns. By coupling the hyperparameters of neighboring coefficients, their hierarchical prior encourages clustered sparsity without requiring explicit knowledge of block boundaries, thus broadening the applicability of block-type models to scenarios where only loose prior information on structure is available.

More recent developments integrate richer structural priors into the BSBL framework. For example, total-variation–regularized SBL formulations incorporate linear and region-aware TV penalties to better capture piecewise smooth or spatially contiguous blocks, improving robustness to mixed sparsity patterns and nonuniform block sizes^[22]. Collectively, these advances indicate a clear research trajectory: from early BSBL models assuming known block partitions and simple intra-block correlations, toward flexible, data-driven formulations that can learn complex, overlapping, or partially specified group structures while retaining the probabilistic interpretability of SBL.

For system identification, block Bayesian approaches are particularly attractive when the underlying dynamics exhibit modular or subsystem-level sparsity. Basis functions corresponding to the same physical mechanism, spatial neighborhood, or input–output channel can be grouped into blocks, so that BSBL simultaneously identifies which modules are active and how their parameters co-vary. The resulting models provide not only sparse representations but also block-level uncertainty quantification, which is valuable for robust control design, fault diagnosis, and adaptive decision-making under noisy or incomplete observations. Overall, block Bayesian identification represents an important evolution of sparse Bayesian methods toward structured, interpretable, and scalable modeling of high-dimensional dynamical systems.

KF and its extensions. Beyond static sparse regression, Kalman-type filters provide a complementary, model-based route to dynamic system identification and belong to the third family in our organizing framework. Since Kalman's seminal work on linear filtering and prediction^[23], the KF has become the standard recursive estimator for linear Gaussian state–space models, with numerous expositions and engineering-oriented tutorials consolidating its theory and implementation^[24]. For system identification, this line of work is particularly attractive because it naturally accommodates process and measurement noise, supports online updating, and yields uncertainty estimates together with state and parameter trajectories.

In the classical linear setting, the KF addresses the discrete-time state–space model

where $$ x_k$$ denotes the latent system state at time step $$ k$$ , $$ u_{k-1}$$ is the control input, $$ z_k$$ is the measurement, $$ A$$ is the state-transition matrix, $$ B$$ is the input matrix, $$ H$$ is the measurement matrix, and $$ w_{k-1}$$ and $$ v_k$$ denote zero-mean process and measurement noises with covariance matrices $$ Q$$ and $$ R$$ , respectively. The algorithm alternates between prediction and correction steps to update the posterior mean and covariance of $$ x_k$$ , thereby minimizing the mean-square estimation error. Within system identification, this framework is routinely used in an augmented-state formulation, where unknown (possibly time-varying) parameters are appended to the physical state, $$ x'_k = [x_k^\top, ,~\theta_k^\top]^\top$$ , and modeled as a slow random walk. This strategy, systematized in Ljung's monograph on system identification^[29], turns parameter estimation into a recursive filtering problem and lays the groundwork for nonlinear extensions.

To handle nonlinear dynamics and/or measurement relations, the EKF linearizes the system around the current estimate via first-order Taylor expansion. Early analyses by Ljung showed that, under mild regularity and excitation conditions, the EKF can exhibit favorable asymptotic behavior when used as a parameter estimator for (mildly) nonlinear systems^[25]. On the application side, EKF-based identification has been widely adopted in structural health monitoring and vibration analysis, where slowly varying stiffness or damping parameters must be tracked in real time from noisy sensor data^[27]. These studies collectively highlight both the appeal and the limitations of EKF: it is computationally efficient and easy to integrate with existing state–space models, but its performance can degrade when linearization errors and Jacobian sensitivities become pronounced.

Motivated by these shortcomings, Julier and Uhlmann introduced the UKF, which replaces local linearization with the Unscented Transform to propagate a deterministic set of sigma points through the nonlinear dynamics^[26]. The UKF captures posterior means and covariances up to second order without requiring analytical Jacobians, offering improved accuracy and numerical robustness for strongly nonlinear or high-dimensional systems. In the identification literature, UKF-type schemes have therefore become a standard choice for joint state–parameter estimation in control, navigation, and robotics, especially when the underlying process models are only moderately well-known and must be refined from data.

More recently, research has begun to fuse Kalman filtering with deep learning in order to handle complex, nonstationary environments. A representative example is the "deep–Kalman" architecture, where recurrent neural networks—often long short-term memory (LSTM) models—are used to learn nonlinear process dynamics or noise structures, while the Kalman update equations enforce probabilistic consistency and provide uncertainty estimates^[28]. Such hybrid schemes illustrate a broader trend: combining physically motivated state–space models with data-driven components, so that the filter can inherit interpretability and stability from classical control theory while leveraging neural networks to model residual dynamics more expressively. This line of work points toward next-generation identification methods that unify model-based and learning-based paradigms, and will be revisited in the multi-system and learning-augmented sections of this review.

At the same time, another important recent development within single-system identification is the rapid emergence of structured recursive parameter-estimation methods for multivariable, nonlinear, and colored-noise systems. In contrast to the classical sparse-regression and Kalman-filtering paradigms reviewed above, these methods emphasize decomposition, auxiliary modeling, innovation expansion, filtering transformations, and explicit exploitation of coupling structures. Since these directions now constitute an important branch of the contemporary system identification literature, we review them separately below.

2.3 Modern structured recursive identification methods

Beyond regularization-oriented regression, SBL, and Kalman-type state-space filtering, recent single-system identification has developed a rich family of structured recursive and iterative methods. A common goal of these methods is to improve estimation accuracy, computational efficiency, and robustness in multivariable, nonlinear, and colored-noise settings. Auxiliary-model identification introduces surrogate models to handle unmeasurable intermediate variables or noise-corrupted quantities; multi-innovation identification extends the innovation from a single instant to a stacked innovation vector so as to improve data utilization; hierarchical identification decomposes a high-dimensional estimation problem into lower-dimensional subproblems; filtering-based identification transforms the data or model to mitigate the bias induced by colored noise; and coupled identification explicitly exploits shared or partially-coupled information structures in multivariable systems. Recent work increasingly combines these ideas rather than treating them as isolated techniques.

Auxiliary-model identification. Auxiliary-model identification is particularly useful when the identification model contains unmeasurable intermediate variables, internal signals, or nuisance quantities that cannot be handled directly by standard least-squares-type estimators. The basic idea is to introduce an auxiliary model and use its generated signals to facilitate parameter estimation in the original nonlinear or stochastic system. Recent studies show that this idea has become an important component of modern recursive identification. For example, Wang et al.^[5] proposed an auxiliary model-based maximum-likelihood recursive least-squares framework for MIMO systems and further extended it to a multi-innovation version with improved estimation accuracy. For Wiener nonlinear output-error systems, Tian et al.^[56] developed an auxiliary-model gradient-based iterative identification method with a moving data window, illustrating how auxiliary models can be combined with iterative optimization in nonlinear settings. For Hammerstein nonlinear systems with non-uniform sampling, Qin and Ji^[10] proposed two-stage gradient-based iterative algorithms based on key-term separation, further demonstrating the continuing development of decomposition- and auxiliary-model-oriented ideas in nonlinear identification. Overall, auxiliary-model identification provides an effective way to transform difficult estimation problems into more tractable recursive or iterative forms, especially when internal variables are unavailable, or the original estimation equations are difficult to solve directly.

Multi-innovation identification. The central idea of multi-innovation identification is to replace the conventional one-step innovation with a stacked innovation vector over multiple sampling instants, thereby improving information reuse and parameter-estimation accuracy. This strategy is particularly effective when the data are noisy, the model is multivariable, or the recursive estimator suffers from slow convergence under single-innovation updates. In recent literature, multi-innovation methods are frequently integrated with other structured approaches rather than used alone. In a neural-network modeling context, Zhou et al.^[9] developed a dynamic-factor and multi-innovation-based output-input feedback Elman-network modeling approach from measurements, indicating that the multi-innovation principle also extends beyond linear parametric models to recurrent nonlinear structures. For multivariable Hammerstein-input-nonlinear systems, Wang and Liu^[8] proposed hierarchical maximum-likelihood multi-innovation identification methods, showing that multi-innovation updates can be combined with hierarchical decomposition to improve both efficiency and estimation quality. In addition, Ding et al.^[57] proposed hierarchical stochastic-gradient and hierarchical multi-innovation stochastic-gradient identification methods for multivariable ARX models, further confirming that multi-innovation has become an important enhancement mechanism in modern recursive identification.

Hierarchical identification. Hierarchical identification aims to reduce the computational burden of high-dimensional estimation by decomposing a complex identification problem into several lower-dimensional subproblems or virtual subsystems. This principle is especially relevant for multivariable systems, where direct estimation often leads to large parameter vectors, high matrix-inversion costs, and reduced numerical efficiency. Recent studies show that hierarchical ideas are now widely embedded into recursive least-squares, stochastic-gradient, and maximum-likelihood schemes. Xing et al.^[6] developed highly efficient filtered hierarchical identification algorithms for MIMO systems with colored noises by combining filtering transformations, auxiliary models, and subsystem decomposition. For multi-input ARX systems with partially-coupled information vectors, Ding et al.^[58] proposed hierarchical recursive-gradient identification methods that explicitly exploit model structure to reduce redundant computation. Wang and Liu^[8] further showed that hierarchical decomposition can also be integrated with maximum-likelihood and multi-innovation ideas in Hammerstein-input-nonlinear systems. Overall, hierarchical identification has evolved from a computational device for large-parameter problems into a broader structural principle for modern recursive system identification.

Filtering-based identification. Filtering-based identification is motivated primarily by the difficulty of identifying systems in the presence of colored disturbances. When the noise is correlated rather than white, standard recursive estimators may become biased or inefficient. The filtering identification idea addresses this issue by transforming the original model or data into a filtered identification form in which the disturbance structure is easier to handle. Xu et al.^[59] provided a recent formulation of this idea by converting systems with colored noise into identification models with white-noise terms and then deriving a corresponding parameter-estimation algorithm. Building on this line, Xing et al.^[6] proposed filtered auxiliary-model recursive least-squares and four-stage filtered hierarchical recursive least-squares algorithms for MIMO systems with colored noises. Ding et al.^[60] further extended the filtering idea to multivariable equation-error ARMA-like systems and developed filtered hierarchical generalized-extended stochastic-gradient, recursive-gradient, and least-squares methods, together with their multi-innovation counterparts. These studies show that filtering-based identification has become a key tool for combining colored-noise robustness with recursive efficiency, especially in multivariable and noise-correlated settings.

Coupled identification. Coupled identification exploits the fact that, in many multivariable systems, different subsystem identification equations are not fully independent but share common or partially shared information vectors. Instead of estimating each subsystem separately, coupled identification explicitly uses this shared structure to reduce redundancy and improve computational efficiency. Liu et al.^[7] proposed a distributed identification method based on a partially-coupled recursive generalized extended least-squares algorithm for multivariate input-output-error systems with colored noises, showing that partially shared information structures can be embedded directly into recursive estimation. Ding et al.^[58] further developed hierarchical recursive-gradient identification methods for multi-input ARX systems with partially-coupled information vectors. More recently, Ding et al.^[57] derived coupled identification models for multivariable ARX systems and combined them with hierarchical stochastic-gradient and hierarchical multi-innovation stochastic-gradient algorithms. Overall, coupled identification should be viewed as a structural extension of recursive system identification: it is not merely concerned with "multiple systems, " but with how shared information can be explicitly encoded within the estimator for a single multivariable plant.

Taken together, auxiliary-model, multi-innovation, hierarchical, filtering-based, and coupled identification methods represent an important modern branch of single-system identification. They are especially relevant for multivariable, nonlinear, and colored-noise systems, where classical one-shot least-squares formulations are often insufficient. More broadly, these methods illustrate a notable recent trend: modern identification algorithms increasingly arise from structured combinations of modeling, filtering, decomposition, and recursive optimization ideas rather than from a single estimation principle alone.

3.1 Fundamental knowledge

Multi-system identification typically involves spatio–temporal data that characterize coupled dynamic processes across multiple subsystems. Before introducing representative algorithms, it is useful to outline a generic formulation and a taxonomy that jointly cover the methods in Section 3.2.

Preliminary formulation. A common abstraction is to view each spatial location $$ (i, j)$$ as a local subsystem with an associated time series $$ T_{ij} \in \mathbb{R}^M$$ . A shared dictionary of temporal basis functions $$ \phi(t)\in \mathbb{R}^{M\times P}$$ is constructed, and local dynamics are encoded by coefficient vectors $$ \mathbf{W}_{ij}\in \mathbb{R}^P$$ :

where $$ T_{ij}$$ denotes the time series associated with spatial location or subsystem $$ (i, j)$$ , $$ \phi(t)$$ is the shared temporal basis or dictionary evaluated at time $$ t$$ , and $$ \mathbf{W}_{ij}$$ is the local coefficient vector encoding the dynamics of subsystem $$ (i, j)$$ . Sparse identification then seeks the most compact representation of $$ \mathbf{W}_{ij}$$ that accurately reconstructs $$ T_{ij}$$ , thereby capturing the dominant structural dependencies of the multi-system.

To account for model mismatch and measurement noise, the identification problem is often relaxed with an error tolerance:

where $$ \|\mathbf{W}_{ij}\|_0$$ counts the number of nonzero coefficients, $$ \|T_{ij}-\phi(t)\mathbf{W}_{ij}\|_p^p$$ measures the reconstruction error under the $$ p$$ -norm, and $$ \delta>0$$ is the admissible error tolerance. or equivalently expressed as a stochastic model with Gaussian perturbations:

where $$ \xi$$ denotes the stochastic perturbation or measurement noise, $$ \mathcal{N}(0, \sigma^2 I)$$ is a zero-mean Gaussian distribution, $$ \sigma^2$$ is the noise variance, and $$ I$$ is the identity matrix with compatible dimension. The objective is to identify the sparsest $$ \mathbf{W}_{ij}$$ that best explains the observed spatio–temporal evolution across all $$ (i, j)$$ , thereby uncovering the underlying mechanisms that govern the multi-system dynamics.

Research progress and methodological taxonomy. From the perspective of this review, the development of multi-system identification can be organized into three complementary methodological lines that jointly cover all algorithms discussed in Section 3.2:

1. Deep-learning-based multi-system modeling: The entire multi-system is treated as a high-dimensional input–output map, and deep architectures automatically extract spatial, temporal, and nonlinear coupling structures. Representative models include CNNs, RNNs/LSTMs, and their hybrids^{[28,61,62,63,64,65,66,67,68,69]}. These methods emphasize expressive power and scalability, and are particularly suitable for complex nonlinear multivariate systems.

2. Sparse equation discovery and PDE-based frameworks: Here, multi-system or multi-field data are used to identify "shared structure + local parameters'' in governing laws. SINDy^[38] and its AIC/weak-form/ensemble extensions^[70,71,72] are prototypical for ODE-level dynamics, while PDE-FIND and its deep-learning-augmented variants^{[43,73,74,75,76,77,78,79,80]} generalize the same sparse-regression philosophy to spatio–temporal fields. This family focuses on representing multi-system behavior through compact, physically interpretable equations within a unified sparse-regression framework.

3. Hybrid and topology-aware multi-system modeling, including fast adaptation: These methods explicitly characterize network structure, piecewise dynamics, and switching logic. Examples include network reconstruction and two-step topological sampling based on sparse regression^[81,82], as well as the IHYDE framework for hybrid dynamical systems^[83]. Later parts of this review extend this line with meta-learning–based fast identification and LLM-assisted workflows, which operate at the level of system families while retaining structural interpretability and compatibility with control and diagnosis.

This taxonomy provides the "overall'' organizing logic for multi-system identification: Section 3.2 revisits concrete algorithms from each of these three methodological lines and analyzes their respective strengths and limitations.

3.2 Classic algorithms

In line with the taxonomy above, this subsection reviews classic algorithms for multi-system identification by grouping them into three methodological families: (ⅰ) deep-learning-based models for high-dimensional nonlinear dynamics; (ⅱ) sparse equation discovery for ODE and PDE models; and (ⅲ) hybrid, topology-aware, and fast-adaptation methods including meta-learning. For each family, we highlight representative formulations and application domains rather than exhaustively cataloging variants.

Deep-learning-based multi-system identification. The use of deep neural networks (DNNs) for dynamic system identification has gradually evolved from proof-of-concept demonstrations to more structured frameworks for large-scale multi-system settings. Early work showed that deep architectures can exploit hierarchical feature representations to capture highly nonlinear mappings between inputs and outputs, outperforming shallow models such as SVMs or single-layer perceptrons in prediction accuracy and noise robustness.

Within this broader trend, Forgione and Piga^[68] proposed a unified neural modeling framework that couples network structure design with identification criteria, targeting nonlinear dynamical systems measured under realistic, noisy conditions. Their results highlight how carefully designed DNN architectures can approximate continuous-time dynamics while respecting stability and regularization requirements. To mitigate the computational burden of large models, Sahu et al.^[69] introduced a compression-oriented scheme in which learned nonlinear mappings are approximated by piecewise linear surrogates, significantly reducing inference time and memory footprint with only modest loss in identification accuracy. These contributions illustrate a first line of research: making deep models practically deployable for system identification through architecture and complexity control.

In parallel, domain-specific studies have demonstrated that DNN-based identification can be effectively adapted to heterogeneous multi-system environments. Alshathri et al.^[67] applied deep networks to cellular communication systems (GSM, UMTS, LTE), showing that data-driven models can learn complex signal patterns without explicit feature engineering and generalize across operating conditions. In power systems, Cui et al.^[28] proposed a multi-layer LSTM (M-LSTM) framework for time-varying parameter identification in composite load models comprising ZIP components and induction motors, validating robustness on IEEE 23-bus and 68-bus benchmarks. These works collectively demonstrate how deep architectures can be tailored to distinct engineering domains while sharing a common identification viewpoint.

More recent studies have moved toward hybrid deep learning structures that explicitly address challenges such as parameter sensitivity, chaotic behavior, and strong stochasticity. For instance, Wang et al.^[66] developed an evolving interval type-2 fuzzy LSTM (eIT2F-LSTM) model that synchronously identifies chaotic systems under noise and random parameter variations, illustrating how uncertainty-aware modules can be embedded into recurrent architectures. Taken together, these contributions trace a clear research trajectory: from generic DNN predictors to domain-tailored and uncertainty-aware deep models that support large-scale, nonlinear multi-system identification.

Despite these advances, several limitations remain. Deep models often exhibit limited interpretability, making it challenging to translate identified structures into physically meaningful models or to certify stability and robustness. Moreover, training and tuning large architectures can be data- and computation-intensive, particularly when labeled data are scarce or heterogeneous across subsystems. These open issues motivate the hybrid and physics-informed approaches reviewed below, which seek to combine the representational power of DNNs with the structural transparency of classical identification.

Sparse equation discovery and PDE-based multi-system identification. Building on the sparse regression ideas outlined in Section 2.2, Brunton et al. introduced the SINDy framework^[38], which recasts model discovery as a sparse regression problem on a library of candidate nonlinear functions. Instead of assuming a priori knowledge of the governing equations, SINDy starts from time-series measurements of the state $$ \mathbf{x}(t)$$ and seeks a parsimonious representation of the vector field

where $$ \mathbf{x}(t)$$ denotes the system state at time $$ t$$ , $$ \frac{d}{dt}\mathbf{x}(t)$$ is its time derivative, and $$ f(\cdot)$$ is the unknown vector field to be identified. Under the assumption that only a small number of terms in a high-dimensional function library are active. Concretely, the dynamics are approximated in the regression form

where $$ \mathbf{X}$$ denotes the sampled state-data matrix, $$ \dot{\mathbf{X}}$$ collects numerical estimates of the corresponding time derivatives, $$ \Theta(\mathbf{X})$$ is a dictionary of candidate functions such as polynomial or trigonometric terms evaluated on $$ \mathbf{X}$$ , and $$ \Xi$$ is a sparse coefficient matrix whose nonzero entries encode the active terms in the governing equations. Sparse regression—typically via sequential thresholded least squares—then identifies a small set of nonzero entries in $$ \Xi$$ , yielding an interpretable set of governing equations.

Since its introduction, SINDy has evolved along several interconnected directions. One line of work integrates SINDy with control and optimization. Kaiser et al.^[84] proposed SINDy-MPC, which couples SINDy-based model discovery with model predictive control, enabling data-driven control design for systems whose dynamics are initially unknown. By learning low-order, interpretable models that are compatible with MPC constraints, this framework highlights the value of sparse models in bridging identification and control.

A second line of research focuses on model selection and robustness. In practice, SINDy often produces a family of candidate models along a sparsity–accuracy Pareto front, and numerical differentiation can amplify measurement noise. To address the first issue, Mangan et al.^[70] proposed SINDy-AIC, which uses SINDy to generate candidate models and then applies the Akaike Information Criterion for principled, automatic model selection. To mitigate noise sensitivity, Schaeffer^[71] introduced integral and weak-form variants that replace pointwise derivatives with integral constraints against smooth test functions, effectively averaging out high-frequency noise. More recently, Fasel et al.^[72] developed Ensemble-SINDy (E-SINDy), which leverages bootstrap ensembles to quantify term inclusion probabilities and improve robustness in low-data, high-noise regimes.

A third, complementary thread extends SINDy to more complex dynamical settings. PDE-FIND and related methods generalize the sparse regression framework to partial differential equations by augmenting the library with spatial derivatives and enforcing appropriate regularization^[43]. These PDE-oriented variants have been combined with Galerkin projections and reduced-order modeling for high-dimensional fluid flows^[85], with physics-informed constraints in plasma dynamics^[86], and with tailored function libraries in nonlinear optics^[87] and computational chemistry^[88]. Across these applications, the core SINDy principle remains the same: exploit the expected sparsity of governing laws in a rich function library to obtain compact, interpretable models from data.

The SINDy framework was extended from ordinary differential equations (ODEs) to partial differential equations (PDEs) in the seminal PDE-FIND method of Rudy et al.^[43]. PDE-FIND (PDE Functional Identification of Nonlinear Dynamics) adapts the core SINDy idea of sparse regression to spatiotemporal fields. The key change lies in the construction of the candidate library $$ \Theta$$ . Whereas the original SINDy for ODEs builds a library $$ \Theta(X)$$ from functions of the state variables, PDE-FIND constructs a richer library $$ \Theta(U)$$ from the field variable $$ U(x, t)$$ and its partial spatial derivatives. Correspondingly, the target equation shifts from the ODE form $$ \tfrac{d}{dt}x(t) = f(x(t))$$ to the PDE form

where $$ u(x, t)$$ denotes the spatio–temporal field variable, $$ u_t$$ is its time derivative, $$ u_x, u_{xx}, u_{xxx}, \ldots$$ denote spatial derivatives of different orders, and $$ N(\cdot)$$ is the unknown nonlinear function that defines the governing PDE.

In practice, $$ u_t$$ and $$ U$$ are sampled at multiple space–time points. The time derivative $$ U_t$$ is reshaped into a single column vector, while the library $$ \Theta$$ is assembled by evaluating candidate terms (e.g., $$ U$$ , $$ U^2$$ , $$ U_x$$ , $$ U_{xx}$$ , $$ U U_x$$ ) at the same sampling locations. This yields a large, typically overdetermined, linear system

Here, $$ U_t$$ is the vectorized time derivative evaluated at sampled space–time points, $$ \Theta(U, U_x, U_{xx}, \dots)$$ is the candidate library containing field variables and spatial derivatives, and $$ \xi$$ is the sparse coefficient vector associated with the candidate PDE terms. Rudy et al.^[43] advocate Sequential Threshold Ridge Regression (STRidge), which is empirically more robust than standard LASSO when the columns of $$ \Theta$$ are highly correlated.

A major limitation of this baseline PDE-FIND pipeline is its sensitivity to noise. Numerical differentiation of spatial data—especially for higher-order derivatives such as $$ u_{xx}$$ or $$ u_{xxxx}$$ —is notoriously ill-conditioned and strongly amplifies measurement noise, which can easily corrupt the identified equation. This challenge has motivated a series of more robust SINDy-based PDE identification methods. A central idea is the Weak SINDy formulation^[73], which replaces pointwise differentiation with an integral (weak) form: the PDE is multiplied by smooth test functions and integrated over the domain. This operation effectively averages out high-frequency perturbations, substantially improving robustness for noisy spatiotemporal data and enabling the discovery of PDEs in weak form, in line with related ensemble–weak formulations^[72].

A second limitation of sparse-regression approaches like PDE-FIND is their reliance on a pre-defined, "closed" candidate library: if a true term is absent, it cannot be recovered. To relax this constraint, genetic-algorithm (GA)–based symbolic discovery has emerged as a complementary direction. These methods explore the space of equation structures via evolutionary operations such as crossover and mutation, allowing discovery even from incomplete or biased libraries. Xu et al.^[74] proposed DLGA-PDE, a hybrid framework that first trains a deep neural network surrogate to generate denoised meta-data and compute stable derivatives via automatic differentiation, and then applies a GA to identify the PDE terms, achieving robustness under noisy and data-scarce regimes. To further enlarge the search space, Chen et al.^[75] introduced SGA-PDE, a symbolic GA that encodes PDE terms as binary trees composed of basic operators (e.g., $$ +$$ , $$ \partial_x$$ ) and operands (e.g., $$ u$$ , $$ x$$ ), enabling gradient-free discovery of complex "open-form" equations, including fractional or composite structures beyond the reach of traditional sparse regression. More recently, Zhang et al.^[76] proposed MORL4PDEs, which combines multi-objective optimization with reinforcement learning to tackle challenging PDE discovery problems involving higher-order derivatives, such as those targeted by SGA-PDE.

A third line of work addresses noise and differentiation issues by integrating deep learning, particularly physics-informed neural networks (PINNs), into the PDE identification pipeline. The central idea is to first train a neural network to approximate the spatiotemporal field, learning a smooth surrogate $$ \hat{u}(x, t)$$ that effectively denoises the raw measurements; automatic differentiation (AD) is then applied to $$ \hat{u}$$ to compute all required temporal and spatial derivatives with machine precision, providing a clean candidate library for subsequent sparse regression. Two seminal contributions established this direction. Both et al.^[77] proposed DeepMoD, which embeds sparse regression directly into the neural network loss so that the network parameters (approximating $$ \hat{u}$$ ) and the PDE coefficients $$ \xi$$ are optimized jointly, yielding strong robustness to noise and small data. Similarly, Chen et al.^[78] introduced the PINN-SR framework, which also uses a PINN and AD for derivative evaluation but couples them with an alternating direction optimization (ADO) strategy that sequentially updates the network weights and sparse PDE coefficients to handle the nontrivial joint optimization. Building on these ideas, subsequent studies have proposed refined neural architectures and training strategies for PDE discovery^[79,80], further tightening the integration between data-driven surrogates, sparse regression, and physical constraints.

Hybrid, topology-aware, and meta-learning-based multi-system identification. Complementary to purely data-driven DNN approaches and sparse equation discovery, a parallel family of methods exploits system topology, piecewise dynamics, and task-relatedness when identifying families of coupled subsystems.

A first line of work focuses on topology and network structure. From a broader perspective, network–from–dynamics methods aim to reconstruct interaction graphs from observed trajectories by casting the inverse problem as a sparse regression or optimization task^[81]. Within this line of research, "two-step'' frameworks first perform a coarse global sparse regression to screen candidate interactions, then refine models locally using the inferred network structure. A prototypical example is Han et al.^[82], who study networked public-goods games and seek to reconstruct both direct and indirect interactions from time-series observations. In the first, global stage, they apply LASSO-based regression to infer an effective interaction matrix from payoff and strategy trajectories, thereby exploiting sparsity to prune implausible links and obtain a coarse-grained description of the network. This step emphasizes scalability and structural discovery: it reveals which nodes appear to influence which others, but the recovered couplings still mix direct and mediated effects. In the subsequent refinement stage, Han et al. introduce a matrix-transformation and optimization procedure to separate direct from indirect influences and to eliminate spurious edges^[82]. By enforcing additional structural constraints, they sharpen the initial coarse network into a compact, interpretable topology that more faithfully reflects the underlying game interactions. Numerical experiments on simulated networks demonstrate that the two-step pipeline can accurately recover the true adjacency structure from noisy, finite-length time series.

A second line of work focuses on hybrid and piecewise dynamics. Yuan et al. (2019) proposed the IHYDE framework (Identification of HYbrid Dynamical systEms)^[83], which aims to infer both the local dynamical laws and the associated switching logic of hybrid systems directly from data. Rather than assuming a single global model, IHYDE explicitly acknowledges that different operating regimes may follow distinct dynamics, and formulates the identification task as an optimization problem in which hybrid models are selected to minimize prediction error while maintaining parsimony. At a high level, IHYDE proceeds in two stages. First, the observed time series are partitioned into segments that are internally coherent in terms of their dynamical behavior. This segmentation can be viewed as assigning each data point (or short time window) to one of several latent subsystems, effectively uncovering a piecewise structure in the state space. Second, for each subsystem, sparse regression is applied to identify a low-complexity dynamical model from a candidate library, yielding distinct local equations for different regimes. In this way, IHYDE generalizes the sparse-identification philosophy (as in SINDy-type methods) from single-system settings to hybrid dynamical systems, jointly learning regime-dependent models and switching structure from data^[83]. Yuan et al. demonstrate that this piecewise identification paradigm is applicable across a wide range of benchmark and engineering systems^[83], including delayed relays, hysteresis loops, phototactic robots, nonlinear hybrid circuits, autonomous vehicles, industrial process monitoring, power-grid fault detection, real-time smart-grid modeling, and electrophysiological models of human atrial action potentials, as well as non-hybrid dynamical systems.

A third line of work targets fast adaptation across families of related systems via meta-learning. Traditional system identification techniques are typically designed for a single system and aim to construct an accurate model from scratch. This often requires substantial amounts of representative data and extensive parameter tuning. In many real-world scenarios, however, engineers encounter not one isolated system but a family of related systems—for example, batches of robots with slightly different dynamics, power converters with manufacturing tolerances, or vehicles exhibiting distinct driving styles. Repeating full identification for each member of such a family is both data-intensive and computationally inefficient, hindering rapid adaptation and large-scale deployment.

Meta-learning–based identification offers a principled way to address this limitation. Rather than learning a single model, meta-learning seeks to acquire a transferable learning procedure or prior over models from a distribution of related identification tasks $$ p(\mathcal{T})$$ ^[89]. By training across multiple historical tasks, the algorithm captures structural regularities shared within a system family and distills them into meta-knowledge $$ \theta$$ . For a new task $$ \mathcal{T}$$ drawn from $$ p(\mathcal{T})$$ , and given a small support dataset $$ D_{\mathcal{T}}^{\text{sup}}$$ , an adaptation operator $$ \text{Adapt}(\cdot)$$ produces task-specific parameters that are evaluated on a query set $$ D_{\mathcal{T}}^{\text{qry}}$$ :

where $$ \theta$$ denotes the meta-parameters shared across tasks, $$ \theta^*$$ is the optimized meta-parameter vector, $$ \mathcal{T}$$ is an identification task sampled from the task distribution $$ p(\mathcal{T})$$ , $$ D_{\mathcal{T}}^{\text{sup}}$$ and $$ D_{\mathcal{T}}^{\text{qry}}$$ are the support and query datasets for task $$ \mathcal{T}$$ , $$ \text{Adapt}(\cdot)$$ denotes the task-specific adaptation operator, and $$ \mathcal{L}(\cdot)$$ is the task loss evaluated on the query set. Once the optimal meta-parameters $$ \theta^*$$ are obtained, the model can be rapidly adapted to a new system using only a small amount of data, enabling efficient few-shot identification^[90].

A canonical optimization-based approach is Model-Agnostic Meta-Learning (MAML)^[89]. MAML aims to learn an initialization $$ \theta$$ such that a few gradient steps on task-specific data yield a performant model. For each identification task $$ \mathcal{T}_i$$ , the inner update adapts the parameters using the support set $$ D_i^{\text{sup}}$$ :

where $$ \theta$$ is the shared initialization, $$ \theta'_i$$ is the task-adapted parameter vector for task $$ \mathcal{T}_i$$ , $$ \alpha$$ denotes the inner-loop learning rate, and $$ \mathcal{L}_{\mathcal{T}_i}(\theta)$$ is the loss of task $$ \mathcal{T}_i$$ evaluated on its support data. The outer update then refines the shared initialization by aggregating query-set losses across tasks:

where $$ \beta$$ is the meta-learning rate, the summation aggregates losses over tasks sampled from $$ p(\mathcal{T})$$ , and $$ \nabla_\theta$$ denotes differentiation with respect to the shared initialization $$ \theta$$ . Because $$ \theta'_i$$ depends on $$ \theta$$ , the exact gradient involves second-order terms; in practice, first-order variants (FOMAML) are often used to reduce computational cost. Viewed from a system identification perspective, this procedure learns an initialization that is "close" to the optimal parameters for all systems in the family, so that only minor task-specific adaptation is needed.

Meta-learning has been combined with several identification paradigms, ranging from black-box sequence models to grey-box physical models and Koopman operator formulations. In black-box identification, neural sequence models such as RNNs, LSTMs, or GRUs directly approximate input–output dynamics,

where $$ \mathbf{h}_t$$ is the hidden state at time $$ t$$ , $$ \mathbf{h}_{t-1}$$ is the previous hidden state, $$ y_{t-1}$$ is the previous output, $$ u_t$$ is the current input, $$ f_{\theta}(\cdot)$$ is the recurrent state-transition function, $$ g_{\theta}(\cdot)$$ is the output mapping, and $$ \hat{y}_t$$ is the predicted output. with task loss typically defined as the mean squared error between predictions and measurements^[90]. In grey-box and physics-informed identification, meta-learning can be wrapped around parametric physical models via hypernetworks $$ g_\theta$$ that map limited support data to estimates of physical parameters, $$ \hat{\phi}_i = g_{\theta}(D_i^{\text{sup}})$$ , which are then inserted into mechanistic simulators. In Koopman operator–based identification, meta-learning helps discover shared nonlinear embeddings $$ \psi_\theta(\cdot)$$ such that

where $$ x_t$$ and $$ x_{t+1}$$ denote consecutive system states, $$ \psi_\theta(\cdot)$$ is the learned nonlinear lifting or embedding map parameterized by $$ \theta$$ , and $$ \mathbf{K}_i$$ denotes the task-specific linear dynamics matrix in the lifted space^[91]. Bayesian meta-learning further interprets $$ \theta$$ as a hierarchical prior over model parameters^[92], providing uncertainty quantification crucial for risk-aware control, adaptive experiment design, and safe deployment in changing environments.

Taken together, hybrid, topology-aware, and meta-learning-based methods provide a third pillar for multi-system identification. They explicitly exploit structural priors and task-relatedness, complementing deep-learning-based and sparse equation discovery approaches and enabling scalable identification of large networks and families of systems.