Deep Time Series in Structural Health Monitoring of Civil Structures: A Review of Architectures, Applications and Challenges

1. Introduction

Many civil infrastructures are in a poor state [1,2,3]. This is mainly due to poor maintenance practices as well as exposure of the infrastructure to harsh environmental conditions [4]. Besides these, structural deterioration has been associated with misuse of infrastructure [5]. Regardless of the cause, structures in such a state usually pose risk to life, property and economic performance. Talking of negative impacts of structural deterioration, South Africa is a good example. To put this into context, South Africa loses almost 50% of the water pumped into its piped networks to leakages, among other factors [6]. This is in part due to deteriorating infrastructure which is rarely maintained. Problems related to deteriorating infrastructure are also encountered in South Africa’s power sector and transportation system [2]. The issue of deteriorating structures is further highlighted in South Africa’s Infrastructure Development scenarios for 2050 [7]. While the problem of deteriorating infrastructure is worse in the developing world, similar challenges are also observed in developed countries. A recent report by MakeUK indicates that approximately 54% of manufacturers in the United Kingdom perceive a decrease in road infrastructure quality over the last decade [8]. These statistics are a clear indication that interventions have to be taken to ensure serviceability of the existing infrastructure. Considering the economic value of civil infrastructure, their exposure to harsh environmental influences, and progressive material deterioration, conventional inspection and maintenance practices, which often rely on manual and periodic interventions, are insufficient to ensure structural integrity and long-term serviceability. As such, there is a need for continuous monitoring to enable early detection of deterioration and timely maintenance decisions [9,10].

Structure Health Monitoring (SHM) appears to be a promising solution to the issues discussed. The field of SHM was developed in the 1990s and has since matured over the years [11]. These days, SHM is extensively studied and applied to civil infrastructure. The aim of SHM is to enable damage identification, localisation, quantification and prognosis. This is accomplished through instrumentation, whereby sensors are deployed to collect data from a structure, followed by data analysis. Fundamentally, SHM is guided by 7 Axioms [12]. Of particular interest to this review is Axiom IV(a), which states that sensors cannot measure damage as such feature extraction is necessary to convert sensor data to damage information. In this regard, statistical and deep learning methods are used.

SHM methods are usually classified into two main categories: global and local [13]. Global methods typically involve acquisition of data across a full structure. The data is then processed to assess damage. Vibration based SHM methods fall under this category. Global methods can be used to detect, locate, and quantify damage on a structure. On the other hand, local methods focus on a particular area of a structure. In most cases, the location is known a priori and the goal is to assess the level of damage. Some methods under non-destructive evaluation are in this category. Considering the wide availability of cheap sensors, some modern SHM systems are hugely global and a common characteristic of such methods is that the generated data is time series. It is therefore important that time series analysis is well understood.

Traditionally, time series analysis uses methods such as Auto Regressive Moving Average(ARMA), AutoRegressive Moving Average with eXogenous excitation (ARMAX), and state space models (SSMs), among others [6]. This is so because time series data exhibit correlations and dependencies which are otherwise absent in static data. Despite the successes of these methods on SHM systems, current data scale in SHM makes them inefficient in some cases. These traditional methods require assumptions which do not fit reality. As is the case with real world systems, noise and missing data are a common issue. Furthermore, modern SHM systems generate huge amounts of highly dimensional data and traditional methods have always found it difficult to analyse such data due to a phenomenon called the curse of dimensionality.

Recently, researchers have turned to deep learning to resolve these issues. Deep learning has proven to work extremely well in high dimensional spaces. With regards to time series, deep learning methods designed for sequential modeling have proven to be effective [14]. However, despite this, the application of deep learning to time series in SHM proves to be an outstanding issues for a number of reasons. Firstly, despite the generation of huge amounts of data by SHM systems, most of this data is inconsistent, contain significant amount of noise and is unlabeled. Secondly, in the real world, the accuracy of models is significantly poor which among other issues is due to domain shift. This is the case because structures are exposed to extremely harsh environments whose nuances might have not been captured in the model. A highlight on these issues is also the lack of integration of existing knowledge in developed models. Another issue with deep learning models is their high computational costs which hinders their development and deployment [15]. Deep learning models are known to be computationally demanding due to the enormous matrix operations which are performed during training and inference. A fourth and less commonly considered reality are regulatory restrictions. The advent of tools like ChatGPT has brought this issue to light [16]. As for structural systems, which by the end of the day can impact life and economies, their deployment is trivial. This is a serious issue as classical deep learning is inherently deterministic and black-box, in which case no rationale is provided for the decision made by the system.

The aim of this paper is therefore to provide a mathematically grounded review of deep learning for time series in an attempt to address the current challenges and scepticisms. The review provides a discussion of neural network architectures which are common in time series, generative deep learning models, and uncertainty quantification in deep learning.To ensure the review is self-contained, we include derivations of key results where appropriate.Although the first part of this review tackles deep learning in general, the primary target audience of this review are SHM practitioners and researchers. Thus, the review provides a discussion of deep learning for time series in SHM, including challenges.

The remainder of the paper is structured as follows: It begins with the discussion of the related works in Section 2. To provide context, Section 3 provides the background to deep learning. This is followed by a discussion of foundational architectures in Section 4 and generative models in Section 5. A review of uncertainty quantification in deep learning is provided in Section 6 .The final discussion in this review deals with applications in Section 7, state of DL in SHM in Section 8 and a concluding synthesis in Section 9

2. Related Works

In recent years, deep learning (DL) has gained significant traction in structural health monitoring (SHM) due to its ability to model complex, high-dimensional, and often nonlinear data [17,18]. Since 2020, several review studies have emerged to synthesize the expanding range of DL applications in SHM. This section analyzes 24 key publications that collectively reflect the current state of research, with a focus on methodological trends and persistent research gaps.

Several reviews have aimed to map the application of DL in SHM. Jia and Li [19], for example, presented a taxonomy based on data types, emphasizing the prevalence of vibration- and vision-based methods. However, their review lacked systematic coverage of sequential modeling or generative approaches. Similarly, Afshar et al. [20] emphasized the predictive power of neural networks for the prediction of structural responses, but provided minimal theoretical grounding.

A focused subset of the current literature addresses vibration-based SHM. Notable among these are works by Spencer et al. [17], Toh and Park [21], Azimi et al. [22], Indhu et al. [23], Wang et al. [24] and Avci et al. [25]. While insightful regarding modal data, these reviews largely overlooked temporal models such as recurrent neural networks (RNNs), temporal convolutional neural networks (TCNNs), and long short-term memory networks (LSTMs), resulting in an application-centric rather than theory-driven perspective.

In vision-based SHM, studies by Hamishebahar et al. [26], Chowdhury and Kaiser [27], Gomez-Cabrera and Escamilla-Ambrosio [28], Sony et al. [29] and Deng et al. [30] reviewed CNN-based methods for surface crack detection, defect classification , among others. Though these works offered comprehensive discussions of image-based SHM and relevant architectures, they paid little attention to temporal modeling, hybrid strategies, or generalization beyond specific materials. Chowdhury and Kaiser [27], for example, focused solely on concrete structures.

A few studies have attempted a more integrative approach. Khan et al. [31] explored hybrid methods combining DL with physics-based models, identifying avenues for methodological convergence, but their review lacked coverage of generative models and uncertainty quantification. Zhang et al. [32] focused on DL-based imputation techniques for incomplete SHM data, though without a rigorous theoretical foundation for the models used.

Despite the rising interest in generative models, only Luleci and Catbas [33] and Luleci et al. [34] have examined their use in SHM, and both lacked mathematical depth. Other reviews by Spencer et al. [17], Cha et al. [35], Abedi et al. [36], Zhang et al. [37], Tapeh and Naser [38] and Xu et al. [39] have made notable contributions, yet often remained narrowly application-focused and did not thoroughly address the theoretical principles underlying DL models.

Collectively, several limitations emerge across the current reviews. Most reviews remain confined to either vision-based or vibration-based SHM, with little integration of multimodal approaches or comparative learning paradigms. While CNNs are extensively covered, sequential models such as LSTMs, gated recurrent units (GRUs), and attention-based architectures receive insufficient attention, despite their common usage in SHM.This analysis finds no prior review dedicated exclusively to these models. Generative models are rarely discussed in depth, and theoretical treatments of uncertainty quantification are largely absent.

CNNs continue to dominate existing literature, while key areas such as uncertainty quantification and generative modeling are rarely covered, and when they are, theoretical treatment is minimal. This review addresses these gaps by providing a unified synthesis of deep learning methods, grounded in theory and explicitly tied to SHM applications. While some theory-focused or theory-practice-focused reviews exist, such as those on RNNs [14], generative models [40], time series [41,42,43], and uncertainty quantification [44,45],this work stands out by offering an integrated perspective across multiple DL paradigms. To the best of our knowledge, it is the first to do so.

3. Background

3.1. Time Series

Since Structural Health Monitoring (SHM) involves continuous monitoring of structures, a significant portion of the generated data is time series [46]. This data may include accelerometers, strain gauges, displacement transducers, and piezoelectric elements that capture the structural response to operational and environmental loading [47] . SHM data is often multivariate, noisy and might contain missing values. As such, modeling this data requires taking these issues into consideration.

Formally, a time series can be represented as $\mathbf{X}=\{x_1,x_2,\dots ,x_{L-1},x_L\}\in {\mathbb{R}}^{\mathrm{D}\times \mathrm{L}},$ where $\mathrm{D}$ denotes the number of dimensions and L represents the length of the time series [42].

Time series are grouped into two major categories, i.e. univariate and multivariate. In a univariate time series, D=1. A typical example would be water level in a dam. For a multivariate time series, D> 1, and thus additional variables are considered. The presence of multiple variables introduces complexity as they may exhibit correlations which need to be taken into account during analysis. In the dam context, water levels are most likely influenced by daily temperature fluctuations and seasonal variations in rainfall.

Time series models typically serve four key purposes, namely; classification [48], forecasting[49], anomaly detection [50,51], and data imputation [52]. In an SHM system, classification involves assigning labels to time series data based on patterns identified in structural behaviour. A common classification task might be the categorization of the state of a structure based on inspection records. On the other hand, forecasting deals with the prediction of future values based on historical trends. In most cases, forecasting deals with continuous values and provides a quantitative outlook on how a system is likely to evolve over time. Anomaly detection (Figure 1) plays a crucial role in identifying deviations from normal behaviour that may indicate a number of issues such as structural damage, shock events e.g. an earthquake, sensor malfunction or simply false alarms. This is important to ensure that events are correctly detected and an appropriate action is taken. Most importantly, anomaly detection can help save a structure and life of its users or occupants.Finally, data imputation focuses on estimating and filling in missing values to ensure the integrity of a dataset. In real-world scenarios, measurements are often incomplete due to factors such as sensor malfunctions, transmission errors, or the high cost of continuous data collection. For robust analysis, it is important to ensure that the missing data is reconstructed with reasonable accuracy. Imputation directly impacts the results of other time series analysis discussed before. Thus, imputation supports the development of robust predictive models and enhances the reliability of decision-making.

3.1.1. Components of a Time Series

Time series data can be decomposed into four primary components, namely: trend, seasonality, cyclic variations, and irregular fluctuations [53].To begin with, the trend represents long-term changes that are often observable as a gradual shift in the system’s behaviour. In structural systems, a trend may manifest as a slow change in the natural frequency of vibrations due to material aging or increasing stress. Such changes could shift the structural dynamics and thus affect structural response to external loads.

On the other hand, seasonality refers to repetitive patterns that happen at fixed time intervals. In a structure, this could appear as periodic vibrations caused by factors, such as daily traffic patterns or seasonal wind fluctuations.

However, not all oscillations are strictly periodic. The cyclic component captures long-term, non-periodic fluctuations. For example, temperature changes over the year can cause materials to expand and contract, altering the vibrational characteristics of a structure over an extended period. This type of fluctuation, while repetitive, does not follow a fixed schedule, making it distinct from seasonal patterns, yet it still reflects periodic behaviour over the long run.

Lastly, irregular fluctuations, represent random, unpredictable disturbances in a time series. In a signal, irregular fluctuations remain after the seasonal, cyclic, and trend components have been removed. Generally, this component is associated with random noise and and cannot be easily explained.

All in all, time series models typically assume an additive $Y_t=T_t+S_t+C_t+R_t$ or multiplicative $Y_t=T_t·S_t·C_t·R_t$ approach to capture the interactions between these components.

3.1.2. Methods of Time Series Analysis

Different methods for time series analysis can be categorised into time domain and frequency-domain [54,55]. Time domain methods consider raw signals relative to the time variable and tend to be useful for forecasting and understanding trends, seasonal patterns, and autocorrelations within the data.However, there are cases where it is important to transform the raw data into the frequency domain to capture other important features. This is done using tools such as such as the Fourier Transform [56] or Wavelet Transform [57]. Frequency-domain analysis is effective at detecting cyclical or periodic patterns, and thus provides insights into the dominant frequencies that influence the behavior of the time series.

3.2. Deep Learning

3.2.1. A Brief History

The foundations of modern neural networks date back to the McCulloch-Pitts model [58] which laid the conceptual groundwork for viewing computation in terms of interconnected logical units. Building on this idea, the perceptron was introduced by Rosenblatt Rosenblatt [59] in the late 1950s. Following that, the development of neural networks had several successes and winters. Despite several challenges, which among others include lack of funding, several researchers persisted on neural network research. The progression of research in this area is in part linked to the the works of Werbos [60], Fukushima [61], Hopfield [62], Rumelhart et al. [63], LeCun et al. [64], Hochreiter and Schmidhuber [65], Bengio et al. [66], Hinton et al. [67] and Glorot and Bengio [68]. Although much progress was made in 1980s, late 1990s and early 2000s, it was not until 2012 that a breakthrough in image classification was achieved using a neural network model [69]. The major reasons for this advancement include data availability, improvements in computational hardware and software, and novel research in machine learning algorithms [70]. Unlike traditional machine learning approaches, neural networks have been empirically shown to perform extremely well on large, highly-dimensional datasets. This is linked to the problem of curse of dimensionality in high dimensional spaces [71,72]. Besides data availability, and hardware and software innovations, he successes of deep learning are attributed to the ability of neural networks to learn a broad class of functions with arbitrary precision. Several results have proven this in the context of multi layer perceptrons [73,74]. The goal of this section is to dissect the inner workings of the neural network model and establish the fundamental principles necessary for their understanding.

3.2.2. Artificial Neuron

A neuron is the basic computational block of an artificial neural network. A conceptual model of a neuron is provided in Figure 2.

Clause 3.4.9 of ISO/IEC 22989:2022 [75] defines a neuron as a primitive processing element which takes one or more input values and produces an output value by combining the input values and applying an activation function, i.e.

$$o_j=\varphi \left(ne{t}_j\right)=\varphi \left(\sum _{i=1}^n{x}_i{w}_{ij}+b_j\right),$$

(1)

The use of an activation function in a neural network is to learn non-linearity in the target function. With a few exceptions, commonly used activation functions are non-linear. An activation function is required to be continuous or piecewise differentiable in order to facilitate the update of network weights, which are usually learned through backpropagation [76]. The choice of an activation function is typically guided by empirical performance, with ReLU [77] being the most commonly used due to its effectiveness in practice. Typical examples of activation functions are shown in Table 1.For a comprehensive treatment of activation functions, the works by Kunc and Kléma [78], Dubey et al. [79] can be consulted.

3.3. Artificial Neural Network

To compose a neural network, a group of neurons is organized into layers. Neurons in one layer are connected to those in the next layer through weights. To build a deep neural network, several layers are arranged in a particular organization called an architecture. While the traditional architectures sequentially stack layers one after the other, it is also common to have architectures with recurrent loops [80] and skip connections [81] which are introduced to solve issues with sequential data and efficient learning in very deep layered networks. For most practical applications, deep architectures are preferred due to their ability to learn rich and expressive feature representations [82].

3.3.1. Learning in Neural Networks

The goal of a neural network is to learn a mapping $f:{\mathbb{R}}^D\to {\mathbb{R}}^Q,$ where D and Q represent the dimensions of the input and output spaces, respectively. The network is designed to handle inputs of arbitrary dimensionality and approximate complex functions. To measure if the network has accurately learned an optimal mapping, its performance is evaluated using a loss function, also known as an objective function, which measures the discrepancy between the neural network’s prediction $f(x_t;\theta )$ and the actual target $y_t$. This is achieved by minimization of the objective function,

$$\widehat{\theta }=arg\underset{\mathit{\theta }}{min}\mathcal{L}\left(\theta \right).$$

(2)

Loss functions differ based on the task at hand. Generally, there are three foundational approaches to machine learning, i.e. supervised learning, unsupervised learning, and reinforcement learning [83]. The discussions in this review will only consider the first two.

In a supervised learning setting, the model is trained to predict an output $y_t\in {\mathbb{R}}^Q$ given an input $x_t\in {\mathbb{R}}^D$. The data is labelled prior to training where each input is set to predict a given output, i.e., $\mathcal{D}={\left\{(x_t,y_t)\right\}}_{t=1}^{\mathrm{T}},$ with $\mathrm{T}$ being the size of the dataset. In this setup, the loss function is designed in such a way that the predicted outputs $f(x_t;\theta )$ closely match the target labels $y_t$. For regression tasks, the mean squared error is a commonly used metric

$$\mathcal{L}\left(\mathit{\theta }\right)={\displaystyle \frac{1}{\mathrm{T}}}\sum _{t=1}^{\mathrm{T}}{\left(y_t-f(x_t;\mathit{\theta })\right)}^2.$$

(3)

For classification tasks, a commonly used loss function is the cross-entropy loss, given by

$$\mathcal{L}\left(\mathit{\theta }\right)=-\sum _{t=1}^{\mathrm{T}}{y}_{t}lnf(x_t;\mathit{\theta })$$

(4)

Cross entropy has its origins in information theory.In ML it is used to measure the deviation of the predicted distribution to the true distribution. The goal is to achieve a low cross entropy, where the predicted probability is close to the true label. In the case where the target output has only two labels, the binary cross entropy is used. For multiple classes, the categorical cross entropy is more appropriate.

It is possible to derive the above loss functions using a probabilistic framework where learning is considered an attempt at modeling the probability distributions over the target variables. With such a framework in place, loss functions for regression are formulated considering a continuous probability distribution while those for a classification task consider a discrete distribution [84]. However, it has to be noted that discrete distributions can also be used to formulate loss functions for regression problems.

In unsupervised learning, the goal is to learn patterns from the data which are then used to cluster, reduce data dimensionality, or generate new data samples. Thus, the objective function for unsupervised learning is crafted with these tasks in mind. These are tackled in detail when discussing generative models.

3.3.2. Regularization

The parameter space is highly multidimensional, encompassing various possible parameter sets. This complexity makes learning challenging, as the learning process must identify parameters that generalize well to unseen data. To address this, a regularization term, $\Omega \left(\theta \right)$, can be incorporated into the loss function, helping steer the learning process toward reasonable parameters

$$\tilde{\mathcal{L}}\left(\mathit{\theta }\right)=\mathcal{L}\left(\mathit{\theta }\right)+\Omega \left(\mathit{\theta }\right).$$

(5)

Model regularization is a way of restricting the model parameters and helps avoid overfitting [85]. With overfitting, the network learns the training data so well but fails to generalize to unseen data during inference. Common explicit regularization terms are the ${\ell }^1$, ${\ell }^2$, and elastic net[86]. The ${\ell }^1$ regularization term is given by, ${\lambda \parallel \mathit{\theta }\parallel }_1$, where $\lambda$ is the regularization parameter, and $\theta$ represents the model weights. In practice, ${\ell }^1$ regularization produces sparse models [87]. Due to this property, ${\ell }^1$ is effective at feature selection and dimensionality reduction in high-dimensional datasets with correlated features. In contrast, ${\ell }^2$ regularization, expressed as, ${\lambda \parallel \mathit{\theta }\parallel }_2^2$, imposes a large penalty on large weights, encouraging the model to select smaller weights. With ${\ell }^2$ regularization, most weights are close to zero. In other settings, the best of the two approaches are combined to produce the elastic net. This creates a more flexible regularization scheme. Besides these explicit regularization techniques, neural networks employ other approaches such as early stopping [88], batch normalization [89], and mixup [90]. The discussion of these approaches is outside the scope of this review.

3.3.3. Update of Network Parameters

During training, neural networks are initialized with weights, which get updated over multiple iterations. The initial weights can be categorized as weak or strong. Strong initial weights impose stricter assumptions compared to their weak counterparts.In this light, various weight initialization techniques exist, including Xavier [68], He [91], LeCun [92], and Orthogonal [93,94].

Typically, weights are updated using the following expression

$${\mathit{\theta }}_{n+1}={\mathit{\theta }}_n-\eta {\nabla }_{\mathit{\theta }}\tilde{\mathcal{L}}\left({\mathit{\theta }}_n\right),\forall n\ge 0.$$

(6)

Here, ${\nabla }_{\mathit{\theta }}\tilde{\mathcal{L}}\left({\theta }_{\mathbf{n}}\right)$ evaluates the sensitivity of the objective function to changes in the parameters.This expression is computed using backpropagation, a specialized case of automatic differentiation [95]. Backpropagation applies the chain rule to compute this term. Following this, parameter updates can occur after a full pass over the dataset (batch gradient descent) or over a smaller subset of data points (mini-batch or stochastic gradient descent(SGD)). Among these two approaches, SGD is preferred as it is more efficient and scalable.

The learning rate $\eta$ governs the step size toward the optimal solution, i.e, ${\nabla }_{\mathit{\theta }}\tilde{\mathcal{L}}\left({\theta }_{\mathbf{n}}\right){|}_{\mathit{\theta }={\mathit{\theta }}^{*}}=0$. Convergence of the model is hugely dependent on the choice of the learning rate. A small learning rate may lead to slow convergence, while a large learning rate risks overshooting minima or diverging altogether. This is so because the objective function encountered in real world problems is a highly complex, non-convex function characterized by multiple local minima and saddle points.

Despite its use, the basic update rule is somewhat naïve in the sense that the update strength remains constant, regardless of the model’s position in the learning trajectory. To overcome this limitation, a momentum term can be introduced to the update rule. Momentum helps accelerate updates in consistent directions and suppresses oscillations in directions with fluctuating gradients. The modified update rule is given by

$$\begin{array}{cc}\hfill v_n& =\gamma v_{n-1}+(1-\gamma ){\nabla }_{\mathit{\theta }}\tilde{\mathcal{L}}\left({\mathit{\theta }}_n\right)\hfill \\ \hfill {\mathit{\theta }}_{n+1}& ={\mathit{\theta }}_n-\eta v_n\hfill \end{array}$$

(7)

Here, $v_n$ is the velocity term, $\gamma \in [0,1]$ is the momentum coefficient (typically around 0.9). This formulation effectively averages gradients over time, enhancing movement in stable directions while dampening updates in volatile regions. However, one issue with traditional momentum is the potential to overshoot minima. A refined version known as Nesterov Accelerated Gradient (NAG) anticipates the future position of parameters by calculating the gradient at a lookahead point,

$$\begin{array}{cc}\hfill v_n& =\gamma v_{n-1}+\eta {\nabla }_{\mathit{\theta }}\tilde{\mathcal{L}}({\mathit{\theta }}_n-\gamma v_{n-1})\hfill \\ \hfill {\mathit{\theta }}_{n+1}& ={\mathit{\theta }}_n-v_n\hfill \end{array}$$

(8)

This subtle change improves convergence in practice and reduces the risk of overshooting. Rosebrock [96] likens Nesterov momentum to a child rolling down a hill who decelerates early upon seeing a brick fence at the bottom. Besides these two approaches, other more advanced optimization algorithms such as Adagrad [97], RMSProp [98], and Adam [99] have also been proposed to further address the challenges in training deep networks.

4. Foundational Architectures

This section introduces the commonly used architectures for time series analysis. Generative models have been reserved for the next section. The aim is to maintain the natural flow of the discussion, since the architectures covered in this section can serve as the backbone for either generative or discriminative models. Our discussion begins with multilayer perceptrons as they are extensively covered in deep learning literature.

4.1. Multilayer Perceptrons

A multilayer perceptron (MLP) is a feedforward neural network whose architecture is defined by an input layer, one or more hidden layers, and an output layer. A typical compact representation of an MLP with an added observation error is shown in Figure 3.

MLPs have found success in modeling highly complex functions for both classification [101,102] and forecasting tasks [103]. Unlike RNNs, as will be seen later, MLPs do not have an inherent memory mechanism to handle sequential data. In order to process such data, the problem is transformed into a supervised learning problem by creating a lagged version of the data. This approach is called the sliding window approach, where to predict the value at time $t+1$, the input consists of a fixed-size window of past observations, typically represented as $(x_t,x_{t-1},x_{t-2},\dots ,x_{t-k})$, and this works for both univariate and multivariate time series.

In addition to the sliding window technique, features inherent to time series can be manually engineered and incorporated into model training. Furthermore, network architectures can be designed to autonomously extract essential features from sequential data. For instance, frequency-domain features can be integrated into the model. While several approaches exist, the focus here is on a novel architecture known as the Fourier Analysis Network (FAN) [104]. FANs activate a certain subset of the neurons in a layer using cosine and sine functions, while the remaining neurons are activated using standard functions such as ReLU. Compared to traditional multilayer perceptrons (MLPs), FANs offer the advantage of reduced trainable parameters and, consequently, lower computational complexity due to their design. In principle, FANs are simply performing a Fourier series analysis, as seen by the compact representation of the layer operation,i.e.

$$\mathbf{B}+{\mathbf{W}}_{\mathbf{c}}cos\left({\mathbf{W}}_{\mathrm{in}}x\right)+{\mathbf{W}}_{\mathbf{s}}sin\left({\mathbf{W}}_{\mathrm{in}}x\right),$$

(9)

which is structurally equivalent to

$$a_0+\sum _{n=1}^{\infty }{a}_{n}cos\left({\displaystyle \frac{2\pi nx}{T}}\right)+\sum _{n=1}^{\infty }{b}_{n}sin\left({\displaystyle \frac{2\pi nx}{T}}\right),$$

(10)

where $\mathbf{B}$ corresponds to the constant term $a_0$, ${\mathbf{W}}_{\mathbf{c}}$ and ${\mathbf{W}}_{\mathbf{s}}$ represent the learned amplitudes analogous to $a_n$ and $b_n$, respectively, and ${\mathbf{W}}_{\mathrm{in}}$ plays the role of the frequency-modulated phase ${\displaystyle \frac{2\pi nt}{T}}$.

Beyond feature engineering and activation design, recent research has demonstrated that carefully structured MLP architectures can achieve competitive performance on time series forecasting tasks. NBEATS [105] introduced a purely MLP-based approach by stacking fully connected layers into forecast and backcast modules. This enables the network to model trend and seasonality components directly from the data without the need for recurrent or convolutional operations. Building on this idea,N-HiTS [106] extended the architecture by introducing hierarchical interpolation strategies that effectively capture multi-scale temporal patterns. N-HiTS was proven to improve accuracy with almost 20% over Transformer architectures for time series. Similarly, TSMixer [107] employ stacked MLPs to separately mix temporal and feature dimensions. This work was inspired by earlier vision models which showed that simple linear projections over time and features can match or surpass more complex sequential models. The last architecture in our discussion is gMLP [108]. This architecture incorporates spatial gating mechanisms into MLP layers to allow the network to model interactions across time steps more effectively without explicit recurrence.

4.2. Recurrent Neural Networks

Over the last 30 years, recurrent neural networks (RNNs) have been the standard deep neural architecture for time series. RNNs are inherently designed to remember past information using memory cells. Thus, RNNs could also be termed memory or stateful networks due to their ability to remember the past. Different variations of RNNs have so far been developed and the following sections tackle these in detail.

4.2.1. Simple RNNs

Simple RNNs use RNN cells which are stacked together to create an RNN. An RNN cell accepts two inputs at each time step; the previous hidden state, which encodes previous information, and the input at the current time step.

Figure 4. A representation of a simple Recurrent Neural Network cell (a) and a single-layered unrolled Recurrent Neural Network (b).

Figure 4. A representation of a simple Recurrent Neural Network cell (a) and a single-layered unrolled Recurrent Neural Network (b).

These inputs are combined using Equation 11 to produce the new hidden state at time step t,

$$h_t=\varphi (\mathbf{W}{h}_{t-1}+\mathbf{W}{x}_t+{\mathbf{b}}_{\mathbf{t}}).$$

(11)

$h_t$ contains a summary vector of all inputs up to and including t[109].Based on the overall structure of the neural network, the hidden state can be used to compute the output, which is given by Equation 12, or it can be passed to the next RNN cell at time step $t+1$,

$$y=\varphi (\mathbf{W}{h}_t+\mathbf{b}),$$

(12)

4.2.2. Deep RNN

The typical RNN discussed so far only considers a single layer. However, as indicated earlier, neural networks perform better with depth. Just as MLPs, this is possible with RNNs. Multiple RNN layers are stacked on top of each other as shown in Figure 5 to compose a Deep RNN. Cells in the top layers get input from the previous cell of the same layer as well as from the bottom layer.

4.2.3. Bidirectional RNN

To help address issues like missing values, RNNs use an innovative architecture called bidirectional recurrent neural network(BiRNN). This architecture was proposed by Schuster and Paliwal [111]. BiRNN combines two RNNs oriented in different directions with each RNN treated separately, i.e., $\overrightarrow{h_t}=f({\mathbf{W}}_{\mathbf{h}}\overrightarrow{h_{t-1}}+{\mathbf{W}}_{\mathbf{x}}{x}_t+{\mathbf{b}}_{\mathbf{h}})$ for the forward RNN and $\overleftarrow{h_t}=f({\mathbf{W}}_{\mathbf{h}}\overleftarrow{h_{t+1}}+{\mathbf{W}}_{\mathbf{x}}{x}_t+{\mathbf{b}}_{\mathbf{h}})$ for the backward RNN. The outputs of these networks are concatenated, i.e., $h_t=\overrightarrow{h_t}\oplus \overleftarrow{h_t}$, to produce the final output.

Figure 6. Architecture of a Bidirectional Recurrent Neural Network (rolled(a) and unrolled(b)). Adapted from Vuong [110].

Figure 6. Architecture of a Bidirectional Recurrent Neural Network (rolled(a) and unrolled(b)). Adapted from Vuong [110].

4.3. Comparisons of Standards RNNs

The three architectures of RNNs discussed thus far have some similarities and differencies. Simple RNNs offer a lightweight structure capable of modeling short-term dependencies; however, their effectiveness diminishes on long sequences due to vanishing gradients [112] and limited memory capacity . On the other hand, deep RNNs are able to learn hierarchical temporal features and tend to perform better than simple RNNs. Despite improving model accuracy, the added depth adds a layer of complexity and there is always a tradeoff between accuracy and model generalization. Thus, careful regularization and initialization are needed in deep RNNs to mitigate training instability. BiRNN’s introduce a layer of complexity with their dual encoding feature. This architecture is useful in applications such as anomaly detection, classification or data imputation tasks where the entire sequence is available a priori . However, bidirectionality is inherently noncausal and thus unsuitable for real-time forecasting, where future observations are not accessible. Despite these differences, all RNNs share a common advantage, weight sharing. Weight sharing reduces the number of trainable parameters in RNNs and thus reducing the computational costs

Table 2. Comparison of Simple RNN, DRNN, and BiRNN in Time Series Modeling.

Aspect	Simple RNN	Deep RNN	Bidirectional RNN
Temporal direction	Forward only	Forward only	Both
Long-term memory	Poor	Improved	Strong (non-causal)
Depth and abstraction	Shallow	Hierarchical	Context-rich
Suitability for forecasting	Online / short horizon	Long horizon forecasting	Not suitable for real-time use
Computation and training	Efficient, stable	Expensive, harder to train	High overhead
Best use case	Basic time series tasks	Multiscale or nonlinear time series	Offline classification or anomaly detection

Table 2. Comparison of Simple RNN, DRNN, and BiRNN in Time Series Modeling.

Aspect	Simple RNN	Deep RNN	Bidirectional RNN
Temporal direction	Forward only	Forward only	Both
Long-term memory	Poor	Improved	Strong (non-causal)
Depth and abstraction	Shallow	Hierarchical	Context-rich
Suitability for forecasting	Online / short horizon	Long horizon forecasting	Not suitable for real-time use
Computation and training	Efficient, stable	Expensive, harder to train	High overhead
Best use case	Basic time series tasks	Multiscale or nonlinear time series	Offline classification or anomaly detection

4.4. Shortfalls of RNNs

There are two main problems with the discussed RNN architectures. Firstly, RNNs are trained through a process called backpropagation through time (BPTT), which involves going back in time n steps to compute the derivative of the loss. Thus, the gradient of the loss function involves repeated multiplication of the Jacobian matrix of the hidden state transition,

$${\displaystyle \frac{\partial {\mathbf{h}}^{\left(\mathrm{T}\right)}}{\partial {\mathbf{h}}^{\left(1\right)}}}=\prod _{t=2}^{\mathrm{T}}{\displaystyle \frac{\partial {\mathbf{h}}^{\left(t\right)}}{\partial {\mathbf{h}}^{(t-1)}}}$$

(13)

If the largest eigenvalue $\lambda$ of this matrix is less than 1, gradients tend to vanish. On the other hand, if $\lambda$ is greater than 1, gradients grow exponentially, leading to the exploding gradient problem (unstable training) [14,113].

There are different ways to resolve these issues including gradient clipping, proper weight initialization, and use of gating mechanisms (LSTMs and GRUs). This issue is also tackled via training technique called teacher forcing, where ground truth at t is fed as input to the cell at $t+1$ [70].

A second shortfall of RNNs, which is also common to their advanced variants, is their lack of parallelizability, as they process data sequentially. This can significantly slow down training times. Attempts to resolve this have led to the development of advanced architectures such as transformers and parallelizable RNNs.

4.5. Long Short-Term Memory

To handle long and short sequences effectively without issues such as exploding or vanishing gradients, Long Short-Term Memory(LSTM) networks were introduced [65]. LSTM neural networks use a cleverly designed cell (LSTM cell) instead of the simple RNN cell. These gates are depicted in Figure 7.

The input to the forget gate is the previous hidden state and current input

$$f_t=\sigma ({\mathbf{W}}_{\mathbf{f}}{x}_t+{\mathbf{U}}_{\mathbf{f}}{h}_{t-1}+{\mathbf{b}}_{\mathbf{f}}).$$

(14)

$f_t$ acts as a mask for the cell memory since a sigmoid function squashes all outputs into values between 0 and 1. Zero entries erase irrelevant previous memory contents at a specific memory location.

To update the cell memory, an LSTM uses the following equations,

$$\widetilde{c_t}=tanh({\mathbf{W}}_{\mathbf{c}}{x}_t+{\mathbf{U}}_{\mathbf{c}}{h}_{t-1}+{\mathbf{b}}_{\mathbf{c}}),$$

(15)

$$i_t=\sigma ({\mathbf{W}}_{\mathbf{i}}{x}_t+{\mathbf{U}}_{\mathbf{i}}{h}_{t-1}+{\mathbf{b}}_{\mathbf{i}}),$$

(16)

$$c_t^{+}=i_t\odot \widetilde{c_t},$$

(17)

$$c_t=f_t\odot c_{t-1}+c_t^{+}.$$

(18)

The input gate $i_t$ controls how much new information, which is represented by the candidate memory $\widetilde{c_t}$, should be stored in the cell state $c_t$. It decides what proportion of the candidate memory will be added to the current memory, effectively controlling the flow of information into the cell state. The candidate memory $\widetilde{c_t}$ represents a potential new memory value that could be added to the cell state. Its value is computed based on the current input and previous hidden state, capturing the new information to be incorporated into the cell state. The input gate modulates this candidate memory and determines how much of it should influence the memory update.

The current hidden state is then calculated based on the updated memory. A copy of the current memory $c_t$ is passed through a tanh activation function to normalize it between -1 and 1. This result is then multiplied by $o_t$, which is computed as

$$o_t=\sigma ({\mathbf{W}}_{\mathbf{o}}{x}_t+{\mathbf{U}}_{\mathbf{o}}{h}_{t-1}+{\mathbf{b}}_{\mathbf{o}}).$$

(19)

The current hidden state, given by

$$h_t=o_t\odot tanh\left(c_t\right),$$

(20)

along with the current cell memory $c_t$, are passed to the next LSTM cell.

4.6. Other LSTM Variants

Different variations of the LSTM cell exist. The focus here is on the peephole LSTM [115] and xLSTM [116] The peephole LSTM extends the standard LSTM architecture by adding connections from the internal cell state to the input, forget, and output gates. The main benefit of this approach is improved gradient flow, which aids in mitigating the vanishing gradient problem. In contrast, xLSTM approaches a similar problem by incorporating additional gates for finer control over information retentions. The architecture also includes an exponential gating mechanism in which exponential activation functions are used in the input and forget gates to better manage data flow

4.7. Gated Recurrent Unit

The Gated Recurrent Unit(GRU) [117] was designed to simplify the LSTM while maintaining comparable accuracy. The GRU removes the need for a separate memory cell and reduces the number of gates to two; the reset gate $r_t$ and the update gate $z_t$.

Figure 8. Internal structure of a Gated Recurrent Unit(GRU) [114].

Figure 8. Internal structure of a Gated Recurrent Unit(GRU) [114].

The reset gate $r_t$ controls how much of the previous hidden state should be considered irrelevant for the current time step. The input to the reset gate is the current input $x_t$ and the previous hidden state $h_{t-1}$

$$r_t=\sigma ({\mathbf{W}}_{\mathbf{r}}{x}_t+{\mathbf{U}}_{\mathbf{r}}{h}_{t-1}+{\mathbf{b}}_{\mathbf{r}}).$$

(21)

After applying the logistic sigmoid activation, $r_t$ determines how much of the previous hidden state should be reset. Next, the candidate current hidden state ${\tilde{h}}_t$ is calculated by multiplying the reset gate output $r_t$ with the previous hidden state $h_{t-1}$, and then adding the current input

$${\tilde{h}}_t=tanh({\mathbf{W}}_{\mathbf{h}}(r_t\odot h_{t-1})+{\mathbf{U}}_{\mathbf{h}}{x}_t+{\mathbf{b}}_{\mathbf{h}}).$$

(22)

The update gate $z_t$ determines the balance between the previous hidden state and the candidate current hidden state. A $z_t$ value of 0 completely ignores the previous hidden state and only considers the current input. Conversely, a $z_t$ value of 1 completely ignores the current input and only considers the previous hidden state. The update gate is computed as

$$z_t=\sigma ({\mathbf{W}}_{\mathbf{z}}{x}_t+{\mathbf{U}}_{\mathbf{z}}{h}_{t-1}+{\mathbf{b}}_{\mathbf{z}}).$$

(23)

Finally, the current hidden state $h_t$ is calculated as a linear interpolation between the previous hidden state $h_{t-1}$ and the candidate hidden state ${\tilde{h}}_t$ based on the value of $z_t$,

$$h_t=z_t\odot h_{t-1}+(1-z_t)\odot {\tilde{h}}_t.$$

(24)

4.8. Comparison Between LSTM and GRU

Both the LSTM and GRU were designed to solve the issue of vanishing gradients and enable long term modeling. The subtle difference between the two architectures is mostly in the cell complexity, where the number of gates in the GRU is reduced to two. While both models perform well on many sequence modeling tasks, GRUs often generalize comparably to LSTMs with less computational cost, making them a preferred choice in scenarios with limited resources or smaller datasets. However, in tasks that require modeling more intricate temporal dynamics, LSTMs may offer better performance due to their more expressive gating mechanisms. For more expressive modeling, in practice, it is also common to use deep or bidirectional RNNs with the LSTM or the GRU as the central computing mechanism.

4.9. Transformer Models

Transformers [118] were introduced to address the limitations of recurrent neural networks (RNNs). The issues referred to here are long-term memory and parallelization. Tranformers rely on self-attention to capture contextual dependencies across input sequences. Mathematically, self-attention is defined as

$$\mathrm{Attention}(\mathbf{Q},\mathbf{K},\mathbf{V})=\mathrm{softmax}\left({\displaystyle \frac{\mathbf{Q}{\mathbf{K}}^{\top }}{\sqrt{\mathbf{d}}}}\right)\mathbf{V},$$

(25)

where $\mathbf{Q}$, $\mathbf{K}$, and $\mathbf{V}$ represent the query, key, and value matrices, respectively, and d denotes the dimensionality of the key vectors. To ensure causality, which is a key issue in time series, masking is used to prevent information from future time steps from influencing the current prediction.

Original transformer models were designed for machine translation tasks. However, over the past few years, transformers have been adapted to time series analysis due to their ability to capture long-range dependencies and enable parallel computation. Notable variants include Temporal Fusion Transformer (TFT) [119], Informer [120], Autoformer [121], FEDformer [122], and PatchTST [123].

4.10. Mamba

Mamba [124] is a recent architecture for efficient and expressive sequence modeling. This architecture was developed to overcome the quadratic complexity $\mathcal{O}\left({\mathrm{L}}^2\right)$ of self-attention in Transformers. It builds on discrete-time state space models (SSMs), represented by

$$h_k=\overline{\mathbf{A}}{h}_{k-1}+{\overline{\mathbf{B}}}_{\mathbf{k}}{x}_k,y_k={\mathbf{C}}_{\mathbf{k}}{h}_k.$$

(26)

Unlike classical SSMs, Mamba introduces input-dependent parameters $\overline{{\mathbf{B}}_{\mathbf{k}}}=B\left(x_k\right)$, ${\mathbf{C}}_{\mathbf{k}}=C\left(x_k\right)$, and a dynamic step size ${\Delta }_k=\Delta \left(x_k\right)$, enabling selective compression of relevant information and improved adaptability. To extend memory and enable parallelization, Mamba uses structured initialization, (HIgh-order Polynomial Projection Operator) [125] for $\overline{\mathbf{A}}$ and a hardware-aware parallel scan algorithm.

4.11. Convolutional Neural Networks

Convolutional neural networks (CNNs) are among the most popular architectures and have significantly influenced the modern era of deep learning. CNNs were initially developed to handle grid-like data such as images. However, researchers have found CNNs to perform well on sequential data, such as time series. It has to be noted that time series can be considered one-dimensional grid data [109]. The fundamental operation in a CNN is the convolution. Mathematically, for the one dimensional discrete case, this is defined as,

$$(G★H)\left(i\right)=\sum _{m}G\left(m\right)H(i-m),$$

(27)

where G is the input signal, H is the filter or kernel, ★ denotes convolution operator, i is the output index, m is the summation index, $G\left(m\right)$ is the value of G at m, and $H(i-m)$ is the flipped and shifted version of H.

Convolutions can also be multidimensional. For 2D convolutions, the operation becomes,

$$(G★H)(i,j)=\sum _m\sum _{n}G(m,n)H(i-m,j-n).$$

(28)

As can be seen, a typical convolution operation involves flipping the function H (kernel).However, in CNNs, it is uncommon to flip function H before applying it, so the operation is more accurately described as cross correlation[126]. Despite this, the term convolution is still commonly used. The purpose of CNNs is to learn a set of functions H, also known as filters or kernels. These filters are typically small tensors (e.g., 1x1, 3x3, 5x5, or 7x7). For multi-dimensional inputs, the filters are extended accordingly. Unlike MLPs, CNNs rely on local connections and share the same filter across the entire input to create the output feature map. Multiple filters can be used in one layer, and their outputs are concatenated to form a complete feature representation. The dimensions of the feature map are dependent on the number of kernels, padding, stride, dilation, pooling, and upsampling. Since the same filter is reused across the input, CNNs are able to learn spatial relationships invariant to transformations. Additionally, CNNs require fewer parameters compared to fully connected networks, making them well-suited for high-dimensional data. Unlike RNNs, CNNs are parallelizable.

4.11.1. 1D CNN

The most basic CNNs used for time series are 1D CNNs [127]. This involves sliding the filter along the sequence length from the beginning. Notably, 1D CNNs can be used for multivariate time series modeling as well, provided that the filter width matches the width of the sequence. Despite their effectiveness, 1D CNNs suffer from future data leakage, violating causality.

4.11.2. Causal CNN

To overcome the issue of future data leakage in 1D CNNs, causal convolutions are used. Mathematically, causal convolution is written as,

$$(G★H)\left[i\right]=\sum _{m=0}^{\mathrm{M}-1}G\left[m\right]H[i-m]$$

(29)

Causal convolutions only use past data to make predictions. To handle long sequences, dilated convolutions, introduced by Yu and Koltun [128] are used to increase the receptive field with fewer layers,

$$(G★H)\left[i\right]=\sum _{m=0}^{\mathrm{M}-1}G\left[m\right]H[i-d·m]$$

(30)

This generalization allows flexible dilation without modifying the filter. The standard convolution is a special case where $d=1$.

4.12. Graph Neural Networks

In several practical applications, time series emerge within networked structures such as road networks, water distribution systems, and railway infrastructures. In such settings, the data exhibit dependencies that are both spatial and temporal in nature. To capture these dependencies accurately requires their explicit modeling. Graph Neural Networks(GNNs) are a natural architecture for this as they provide a more flexible way to model such dependencies. Graph Neural Networks that take into account the aspect of time are called dynamic. In this section, we will firstly discuss the fundamentals of GNNs and will then consider how they fit in time series analysis.

Typically, the input to a GNN is a graph, e.g. a protein structure, social network, transportation network, etc. A graph $\mathcal{G}=(\mathcal{V},\mathcal{E},\Phi ),$ is defined by a set of nodes $\mathcal{V}$, a set of edges $\mathcal{E}$, and associated mappings $\Phi$ that describe the relationships between nodes. Mathematically, a graph is represented either through an adjacency list or adjacency matrix (A). The latter is widely used in Graph Neural Network applications. Each element in A indicates the existence or absence of an edge between a pair of nodes, i.e., $A_{ij}=1$ if there is an edge from node i to node j, and $A_{ij}=0$ otherwise. The specific structure of A depends on the type of graph. Graphs can have single or multiple edges between a pair of nodes. It is also possible to have edges which start and terminate on the same node (self-loops). Further categorisation of a graph depends on edge- directionality and edge-weight. Thus, graphs can be classified as directed or undirected as well as weighted or unweighted. A special case which will be our focus are undirected graphs without self-loops or multiple edges. The analysis of such graphs is simple since the adjacency matrix is symmetric.

Node ordering in a graph is arbitrary as such different formulations of the adjacency matrix are possible [84]. For this reason, a permutation matrix P is used to ensure that the output of the GNN remains invariant to the ordering of nodes. This is achieved through permutation-invariant embeddings $X^{\prime }=XP$ and a transformed adjacency matrix $A^{\prime }=P^{T}AP$.

In general, the goals of a GNN span a variety of tasks, including node classification [129], link prediction [130], graph classification [131], and graph generation [132] . A GNN aims to map the input node or edge features into an embedding space that offers a rich representation for further modeling [133]. These node or edge features are continuously updated by aggregating information from neighbouring nodes or edges through a process known as message passing,

$$h_v^{\left(t\right)}=\varphi \left({\mathbf{W}}^{\left(\mathbf{t}\right)}·\mathrm{AGG}\left(\left\{h_u^{(t-1)}\mid u\in \mathcal{N}\left(v\right)\right\}\right)\right),$$

(31)

where $\varphi$ is a non-linear activation function, ${\mathbf{W}}^{\left(\mathbf{t}\right)}$ is the learnable weight matrix at layer t, AGG is a neighborhood aggregation function such as mean or sum, and $\mathcal{N}\left(v\right)$ denotes the set of neighbors of node v.

Message passing in a Graph Neural Network (GNN) occurs iteratively over several iterations. After a specified number of iterations t, each node encodes information about all other nodes in the graph. The number of iterations is considered optimal to ensure that the GNN has a sufficiently large receptive field while preventing oversmoothing; a scenario where all embeddings become indistinguishable [134]. The basic approach assumes all nodes and edges are of equal importance; however, more advanced methods have been developed to account for varying levels of importance.

4.12.1. Temporal Modeling in GNN

To make predictions, a GNN typically utilizes the final embedding layer. For temporal modeling, it is possible to integrate a GNN with other network architectures, where the final embedding layer is passed as input to subsequent neural network blocks. This approach is commonly used in temporal GNNs [135], where the final embedding is fed into another network block, such as an RNN or CNN. This allows the network to capture both spatial patterns using the GNN and temporal features using the RNN or CNN.

GNNs can also be used for data which is not spatial by nature. This is possible by first transforming the data into a graph structure. The approach seems to work quite well with time series. Liang et al. [136] employed this technique to impute missing values in traffic datasets. For a comprehensive overview of GNN applications in time series tasks, refer to Jin et al. [137].

4.13. Physics-Informed Neural Networks

Traditionally, the approach to training neural networks has been heavily data-centric. However, this approach may have drawbacks, such as noisy data that deviates from actual physical reality.To resolve this, researchers have been experimenting with the idea of incorporating the physics of a problem into a model. Enforcing physics into a model is a well studied field of machine learning under the umbrella of Scientific Machine Learning (SciML). Several studies have so far shown that models in scientific settings perform better on unseen data when the physics of the phenomenon being studied is considered [138] . Although SciML has been achieved through different techniques, our consideration in this review are Physics Informed Neural Networks (PINNs) [139] . PINNs are a household name in SciML hence their discussion. PINNs enforce physics through a physics regularization term which is added to the loss function. This is achieved through the problem’s differential equation and its associated conditions.

Consider the general form of a partial differential equation,

$$\mathcal{F}\left(u\right(x,t\left)\right)=f(x,t),$$

(32)

where $\mathcal{F}$ represents a differential operator, $u(x,t)$ is the solution to the PDE, and $f(x,t)$ denotes the forcing term. The solution to this PDE must satisfy both initial and boundary conditions. Specifically, the initial condition is,

$$u(x,t_0)=g\left(x\right),x\in \Omega ,$$

(33)

where $t_0$ is the initial time, and $g\left(x\right)$ describes the initial state of the system. The boundary condition is expressed as,

$$u(x,t)=h(x,t),t\in [t_0,T],x\in \partial \Omega ,$$

(34)

where $h(x,t)$ prescribes the behavior of the solution on the boundary $\partial \Omega$ of the spatial domain $\Omega$. Here, the spatial variable x spans the domain $\Omega \subset {\mathbb{R}}^n$, and the temporal variable t lies within the interval $[t_0,T]$.

The first step in a PINN is to calculate the residual of the PDE as follows,

$$\mathcal{R}(x,t,u)=\mathcal{F}\left(u\right(x,t\left)\right)-f(x,t).$$

(35)

Residuals are also calculated for initial and boundary conditions, and the total residual is termed the physics loss. Minimizing this loss across the domain and over the specified time interval is crucial for validating the accuracy and consistency of the solution under the given initial and boundary conditions. This ensures that the proposed solution $u(x,t)$ adheres to both the dynamics described by the PDE and the constraints imposed by the initial and boundary conditions. However, in practice, exact satisfaction is rarely achievable; hence, PINNs do not perfectly satisfy the governing physics [140,141].

PINNs are usually trained by minimizing the following composite loss function,with each term adequately weighted,

$$\mathcal{L}\left(\mathit{\theta }\right)={\lambda }_{ic}{\mathcal{L}}_{\mathrm{ic}}\left(\mathit{\theta }\right)+{\lambda }_{bc}{\mathcal{L}}_{\mathrm{bc}}\left(\mathit{\theta }\right)+{\lambda }_r{\mathcal{L}}_{\mathrm{r}}\left(\mathit{\theta }\right)+{\lambda }_{data}{\mathcal{L}}_{data}\left(\mathit{\theta }\right),$$

(36)

where,

$${\mathcal{L}}_{\mathrm{ic}}\left(\mathit{\theta }\right)={\displaystyle \frac{1}{{\mathrm{N}}_{\mathrm{ic}}}}\sum _{i=1}^{{\mathrm{N}}_{\mathrm{ic}}}{\left|u_{\mathit{\theta }}(x_{\mathrm{ic}}^i,t_0)-g\left(x_{\mathrm{ic}}^i\right)\right|}^2,$$

$${\mathcal{L}}_{\mathrm{bc}}\left(\mathit{\theta }\right)={\displaystyle \frac{1}{{\mathrm{N}}_{\mathrm{bc}}}}\sum _{i=1}^{{\mathrm{N}}_{\mathrm{bc}}}{\left|u_{\mathit{\theta }}(x_{\mathrm{bc}}^i,t_{\mathrm{bc}}^i)-h(x_{\mathrm{bc}}^i,t_{\mathrm{bc}}^i)\right|}^2,$$

$${\mathcal{L}}_{\mathrm{r}}\left(\mathit{\theta }\right)={\displaystyle \frac{1}{{\mathrm{N}}_{\mathrm{r}}}}\sum _{i=1}^{{\mathrm{N}}_{\mathrm{r}}}{\left|\mathcal{R}(x_{\mathrm{r}}^i,t_{\mathrm{r}}^i,u_{\mathit{\theta }})\right|}^2.$$

${\lambda }_{ic}$ is the Initial Condition Weight, ${\lambda }_{bc}$ is the Boundary Condition Weight, ${\lambda }_r$ is the Physics Weight, and ${\lambda }_{data}$ is the Data Weight, all balancing their respective influences in the PIN model.

PINNs have been successfully applied to time series problems, such as the estimation of cuffless blood pressure [142]. A recent application of PINNs by Park et al. [143] involved weather data forecasting, where an RNN served as the backbone network and the harmonic oscillator equation was used to enforce physical constraints (see Figure 9).

Due to the nature of the harmonic motion solutions, the neural PDE component effectively modeled the seasonal component of the time series. This approach bears some similarity to the Fourier Analysis Network described earlier.

4.13.1. Issues with PINNs

Similar to other neural architectures, PINNs have shortfalls. A major shortfall lies in the training instability and difficulty in balancing the contributions of data loss and physics loss, especially when dealing with stiff PDEs or multi-scale phenomena [144]. This imbalance often leads to poor convergence or solutions that satisfy the data but violate the governing equations, or vice versa. Moreover, PINNs typically rely on automatic differentiation to compute residuals, which can be computationally intensive . Another notable challenge is the difficulty PINNs face in enforcing complex boundary conditions, especially in domains with discontinuities or sharp gradients [145]. A detailed investigation into the failure modes of PINNs and strategies to address them is provided by Krishnapriyan et al. [140].

5. Generative Modeling

In practical applications, there are scenarios where it is essential to model the input data itself. In such cases, a different class of models, generative models, may be more appropriate. Generative models learn the joint probability distribution $p(x,y)$ and can use this distribution to generate new data points that resemble the training data. The aim of this section is to introduce various classes of generative models that are commonly employed for tasks such as data imputation and synthetic data generation in time series applications.

5.1. Autoencoders

The goal of an autoencoder is to learn a low-dimensional latent space of the data, and then to reconstruct the original data from this compressed space. An authoencoder achieves this using an encoder $z=g(\mathbf{W}x+\mathbf{b})$, for compression, and $x^{*}=f(\mathbf{W}z+\mathbf{b})$, for reconstruction. Learning is done by minimizing the reconstruction loss, which is given by

$$\mathcal{L}(x,x^{*})={\parallel x-x^{*}\parallel }_2^2.$$

(37)

Exact reconstruction is not possible because the latent space has a lower dimension than the input and output. However, this architectural design is necessary to encourage the autoencoder to learn useful features required for reconstruction.Different researchers have successfuly applied autoencoders to data denoising [146], data imputation [147], anomaly detection [148], and more. Despite its wide usage in literature, the standard autoencoder cannot perform meaningful data interpolation, and this arises from the discrete nature of the learned latent space [149]. Among other approaches, this limitation is addressed by variational autoencoders.

5.2. Variational Autoencoders

Unlike traditional autoencoders, variational autoencoders (VAEs) [150] adopt a probabilistic framework to learn a more structured latent space.

Figure 10. Architecture of a Variational Autoencoder (VAE). The input x is passed through an encoder network $g(x;\theta )$ to predict the variational parameters: $\mu$ and $\Sigma$, which parameterize a Gaussian variational distribution $q_{\theta }(z\mid x)$. A latent variable $z^{*}$ is sampled from this distribution using the reparameterization trick: $z^{*}=\mu +{\Sigma }^{1/2}\odot ϵ$, where $ϵ\sim \mathcal{N}(0,\mathrm{I})$. This sample is then fed into the decoder $f(z^{*};\varphi )$, to predict the data x. Adapted from Prince [84].

Figure 10. Architecture of a Variational Autoencoder (VAE). The input x is passed through an encoder network $g(x;\theta )$ to predict the variational parameters: $\mu$ and $\Sigma$, which parameterize a Gaussian variational distribution $q_{\theta }(z\mid x)$. A latent variable $z^{*}$ is sampled from this distribution using the reparameterization trick: $z^{*}=\mu +{\Sigma }^{1/2}\odot ϵ$, where $ϵ\sim \mathcal{N}(0,\mathrm{I})$. This sample is then fed into the decoder $f(z^{*};\varphi )$, to predict the data x. Adapted from Prince [84].

First, the input is mapped to a latent space z using an encoder $z=g(x;\mathit{\theta })$, where $z\sim p_{\mathit{\theta }}(z\mid x)$. Since $p_{\theta }(z\mid x)$ is generally intractable due to the intractability of the marginal likelihood $p\left(x\right)$, VAEs introduce a more tractable approximate posterior $q_{\mathit{\theta }}(z\mid x)$. Then, to reconstruct the data from the latent space, a decoder $x^{*}=f(z;\varphi )$ is used.

VAEs are trained by maximizing the marginal likelihood over the reconstructed sample, $p_{\theta }\left(x\right)$, which is obtained by marginalizing over the latent space,

$$p_{\varphi }\left(x\right)=\int p_{\varphi }(x,z)dz.$$

(38)

This integral is generally intractable and is typically simplified to provide tractable proxy objectives. By introducing a variational distribution $q_{\mathit{\theta }}(z\mid x)$, we can apply Jensen’s inequality to get lower bound of the log marginal likelihood,

$$ln{p}_{\varphi }\left(x\right)=ln\int q_{\mathit{\theta }}(z\mid x){\displaystyle \frac{p_{\varphi }(x,z)}{q_{\mathit{\theta }}(z\mid x)}}dz\ge \int q_{\mathit{\theta }}(z\mid x)ln\left\{{\displaystyle \frac{p_{\varphi }(x,z)}{q_{\mathit{\theta }}(z\mid x)}}\right\}dz.$$

(39)

This gives rise to the evidence lower bound (ELBO),

$$\begin{array}{cc}\hfill \mathcal{L}(\mathit{\theta },\varphi )& =\int q_{\mathit{\theta }}(z\mid x)ln\left\{{\displaystyle \frac{p_{\varphi }(x\mid z)p_{\varphi }\left(z\right)}{q_{\mathit{\theta }}(z\mid x)}}\right\}dz\hfill \\ \hfill & =\int q_{\mathit{\theta }}(z\mid x)ln{p}_{\varphi }(x\mid z)dz-\mathrm{KL}(q_{\mathit{\theta }}(z\mid x)\Vert p_{\varphi }\left(z\right)).\hfill \end{array}$$

(40)

Intuitively, the first part of this expression represents the reconstruction term, which encourages the decoder to accurately reconstruct the input data. The second term encourages the variational distribution to be as close as possible to the prior. It is common not to sample directly from $q_{\theta }(z\mid x)$, but rather to use a reparameterization trick (Equation 41),

$$\mathbf{z}=\mu +{\Sigma }^{1/2}\odot ϵ\mathrm{with}ϵ\sim p\left(ϵ\right).$$

(41)

This allows the expectation over $q_{\mathit{\theta }}(z\mid x)$ to be rewritten as an expectation over a fixed noise distribution $p\left(ϵ\right)$. Expression (40) can therefore be rewritten as a reparameterized objective

$$\mathcal{L}(\mathit{\theta },\varphi )={\mathbb{E}}_{ϵ\sim p\left(ϵ\right)}\left[ln{p}_{\varphi }(x\mid z)+ln{p}_{\varphi }\left(z\right)-ln{q}_{\mathit{\theta }}(z\mid x)\right].$$

(42)

According to Gundersen [151], without the reparameterization trick, backproping would not compute an estimate of the derivative of the ELBO and would thus give no guarantee that sampling large numbers of z will help converge to the right estimate of the gradient. The gradient of the ELBO is computed with respect $\mathit{\theta }$ and $\varphi$ via

$${\nabla }_{\mathit{\theta },\varphi }{\mathbb{E}}_{ϵ\sim p\left(ϵ\right)}\left[ln{p}_{\varphi }(x\mid z)+ln{p}_{\varphi }\left(z\right)-ln{q}_{\mathit{\theta }}(z\mid x)\right]$$

(43)

At this stage, standard stochastic optimization algorithms (e.g., SGD) are used to optimize the variational parameters. Typical VAE architectures for time series include HyVAE [152], TimeVAE [153], and Time-NeighbourVAE [154].

5.3. Generative Adversarial Network

A Generative Adversarial Network (GAN) [126] aims at generating realistic data similar to real-world data present during training by fooling another network into believing that the generated samples are real. This is achieved by simultaneously training two networks: a generator G and a discriminator D.

Figure 11. Schematic of a Generative Adversarial Network (GAN) training loop. A latent variable z is sampled and transformed by the generator $G(z;\theta )$ to produce samples. These, along with real data $x^{*}$, are evaluated by the discriminator $D(x;\varphi )$, which outputs the probability of each sample being real. The discriminator is updated to better distinguish real from generated, while the generator updates its parameters to fool the discriminator.Adapted from Calin [155].

Figure 11. Schematic of a Generative Adversarial Network (GAN) training loop. A latent variable z is sampled and transformed by the generator $G(z;\theta )$ to produce samples. These, along with real data $x^{*}$, are evaluated by the discriminator $D(x;\varphi )$, which outputs the probability of each sample being real. The discriminator is updated to better distinguish real from generated, while the generator updates its parameters to fool the discriminator.Adapted from Calin [155].

The generator creates samples $G(z_j,\theta )$ from a latent variable $z_j\sim \mathcal{N}(·)$. The objective of the generator is to produce data that is indistinguishable from real data. Mathematically, this is equivalent to minimizing the following expression with respect to $\theta$,

$${\mathcal{L}}_G={\mathbb{E}}_{z\sim p_z\left(z\right)}\left[ln(1-D(G(z;\theta )\left)\right)\right].$$

(44)

Technically, if samples are assigned a higher probability, i.e., if the discriminator classifies them as being real, the above expression approaches zero. However, if generated data is classified correctly, i.e., a low probability is assigned to generated data, the expression tends to infinity. This means that the generator loss can only be reduced if the discriminator incorrectly classifies generated data, in which case the generator achieves a minimum loss.

As already highlighted, the discriminator is responsible for classifying data as being real or generated(binary classifier). The input to the discriminator is either a generated or real sample, and it returns a probability that is higher when a sample is real (i.e., close to 1). The objective of the discriminator is usually to maximize the following expression with respect to $\varphi$

$${\mathcal{J}}_D={\mathbb{E}}_{x\sim p_{\mathrm{data}}}\left[logD(x;\varphi )\right]+{\mathbb{E}}_{z\sim p_z\left(z\right)}\left[log\left(1-D\left(G\right(z;\theta \left)\right)\right)\right]$$

(45)

Despite the fact that the underlying principle behind GANs is conceptually straightforward, they are difficult to train in practice. A major challenge is mode collapse, where the generator ignores the variability in the latent space and produces limited or nearly identical output [84,156]. To address this issue, several techniques have been proposed. They include feature matching and mini-batch discrimination [157], unrolled GANs [158], Wasserstein GANs (WGAN) [159], spectral normalization [160], and gradient penalty methods (e.g., WGAN-GP) [161].

Although GANs were initially developed for image data, they have since been extended to handle sequential data such as time series data. Notable adaptations in this space include TimeGAN [162], SeriesGAN [163]; FinGAN [164] to name a few. For a detailed survey of GAN applications in time series, the works of Brophy et al. [165] and Zhang et al. [166] can be consulted.

5.4. Normalizing Flows

Instead of implicitly modeling the data distribution like GANs, normalizing flows explicitly model the probability distribution over the data [84]. Normalizing flows transform a simple distribution over a latent variable $z\sim \mathcal{N}(·)$ into a more complex distribution using differentiable and bijective transformations $f_{\mathit{\theta }}$ (see Figure 12).

The transformation function is usually parameterized by a neural network. Transformations are required to be differentiable and bijective so as to ensure a smooth transition between the latent and data space while enabling efficient learning via gradient-based methods. Several techniques have been developed to design transformations with these desired properties. Among them are coupling flows [168,169], neural ODEs [170], and autoregressive flows [171].

Since the data distribution is defined through a transformation, its density is determined using the change of variables formula,

$$p_X\left(x\right)=p_Z\left(z\right){\left|\left({\displaystyle \frac{\partial f_{\mathit{\theta }}\left(z\right)}{\partial z}}\right)\right|}^{-1},$$

(46)

where $z=f_{\mathit{\theta }}^{-1}\left(x\right)$.

The objective is then to maximize the likelihood of the data under this model, i.e.,

$$arg\underset{\mathit{\theta }}{max}\sum _{i=1}^{\mathrm{N}}ln{p}_X\left(x^{\left(i\right)}\right),$$

(47)

Recent studies have applied normalizing flows to capture complex temporal dependencies and uncertainty in time series. Dai and Chen [172] proposed a graph-augmented normalizing flow (GANF) using Bayesian networks to model inter-series causality.Flow based modeling has also been combined with temporal and attribute-wise attention for high-dimensional, label-scarce data [173]. In forecasting, Rasul et al. [174] integrated autoregressive models with conditional flows to enhance extrapolation and capture temporal correlations. Fan et al. [175] proposed a decoupled formulation called instance normalization flow (IN-Flow), which addresses distribution shifts via an invertible transformation network.

5.5. Diffusion Models

Diffusion models (DM) [176] are trained by sequentially adding noise to a data sample over multiple time steps and then attempting to remove this noise. The noise is added according to a Markovian process

$$\begin{array}{cc}\hfill {\mathbf{z}}_1& =\sqrt{1-{\beta }_1}·\mathbf{x}+\sqrt{{\beta }_1}·{ϵ}_1\hfill \\ \hfill {\mathbf{z}}_t& =\sqrt{1-{\beta }_t}·{\mathbf{z}}_{t-1}+\sqrt{{\beta }_t}·{ϵ}_t\forall t\in \{2,\dots ,\mathrm{T}\},\hfill \end{array}$$

(48)

with its associated distribution,

$$\begin{array}{cc}\hfill q({\mathbf{z}}_1\mid \mathbf{x})& =\mathcal{N}\left(\sqrt{1-{\beta }_1}\mathbf{x},{\beta }_1\mathbf{I}\right)\hfill \\ \hfill q({\mathbf{z}}_t\mid {\mathbf{z}}_{t-1})& =\mathcal{N}\left(\sqrt{1-{\beta }_t}{\mathbf{z}}_{t-1},{\beta }_t\mathbf{I}\right)\forall t\in \{2,\dots ,\mathrm{T}\}.\hfill \end{array}$$

(49)

${\beta }_t\in [0,1]$ is a hyperparameter which determines how quickly the noise gets blended. This hyperparameter can be fixed across all time steps or it can vary according to a predefined schedule. While single-step updates are possible for adjacent time steps by applying Equation 48, it is usually helpful to be able to move directly from $t=0$ to $t=\mathrm{T}$ in a single step as well. This is important especially in the reverse process when denoising. The expression required to do this can be derived by generalizing Equation 48 over multiple time steps and is expressed as,

$${\alpha }_t=\prod _{i=1}^{t}1-{\beta }_i,z_t={\sqrt{\alpha }}_{t}x+{\sqrt{1-\alpha }}_t{ϵ}_t.$$

(50)

Derivation of this expression takes advantage of the properties of sums of Gaussian random variables. On the other hand, during the reverse process, the transition process is modelled as,

$$z_{t-1}={\displaystyle \frac{1}{\sqrt{1-{\beta }_t}}}{z}_t-{\displaystyle \frac{{\beta }_t}{{\sqrt{1-\alpha }}_t{\sqrt{1-\beta }}_t}}{ϵ}_{\theta }$$

(51)

with an associated transition probability $p(z_{t-1}\mid z_t)$. ${ϵ}_{\theta }$ is the noise predicted by the neural network, which is parameterized by $\theta$. The goal is to learn a mapping that makes ${ϵ}_t$ and ${ϵ}_{\mathit{\theta }}$ as close as possible, which leads to the following loss function,

$$\mathbb{E}\left[{∥{ϵ}_t-{ϵ}_{\mathit{\theta }}∥}^2\right].$$

(52)

In other words, the model tries to predict the amount of noise that was added to the data sample during the forward (diffusion) process.

The discussion thus far has primarily focused on denoising diffusion probabilistic models (DDPMs), often considering the unconditional generation setting where the data distribution is modeled without auxiliary information. However, general DMs are also capable of conditional generation, where the data distribution is modeled given some conditioning variable such as class labels or other contextual attributes. DMs have been proposed for forecasting [177,178,179], classification [180,181], anomaly detection [182,183] , imputation [184,185], and generation [186,187]. A comprehensive review of the use of diffusion models in time series has been conducted by Lin et al. [188], covering key models, methodologies, and applications across forecasting, imputation, and generation tasks. More recently, Yang et al. [189] provided an extensive survey encompassing both time series and spatio-temporal applications of diffusion models, categorizing existing approaches based on task type and model structure.

5.6. Autoregressive Models

Given a sequence, $\mathbf{x}=(x_1,x_2,\dots ,x_T)$, the joint probability distribution can be factorized using the chain rule of probability,

$$p_{\mathit{\theta }}\left(\mathbf{x}\right)=\prod _{t=1}^{\mathrm{T}}{p}_{\theta }(x_t\mid x_{<t}),$$

(53)

where $x_{<t}=(x_1,x_2,\dots ,x_{t-1})$.

This is called an Autoregressive model (ARM) and it allows for the exact computation of the likelihood [190]. To capture complex dependencies within the data, each conditional distribution in equation 53 is parameterized by a neural network [191]. Tomczak [192] provides good cases for ARMs. Generally, the conditional distribution can be modeled as it is, but in practice, it is usually computationally expensive as the sequence length increases. An alternative approach employs the first-order Markov assumption to simplify the conditional distribution at time t.Thus, the conditional distribution at time step t depends only on the immediate previous time step $t-1$, i.e.

$$p(x_t\mid x_{<t})\approx p(x_t\mid x_{t-1}).$$

(54)

An ARM is trained by maximizing the likelihood of the observed data. This is also equivalent to minimizing the negative log-likelihood (NLL) given by

$${\mathcal{L}}_{\mathrm{NLL}}=-\sum _{t=1}^{\mathrm{T}}lnp(x_t\mid x_{<t}),$$

(55)

Once trained, data generation follows an ancestral sampling procedure, where values are sequentially sampled based on the learned conditionals.

5.7. Energy-Based Models

Energy-Based Models (EBMs) offer a general framework for modeling probability densities by associating scalar energy values to configurations of the input space.EBMs have their origins in statistical physics.The underlying assumption is that observed, high-probability data should correspond to low-energy states, while implausible or unlikely configurations should be mapped to regions of higher energy.

Formally, the model defines an unnormalized density over $\mathcal{X}$ via a Boltzmann distribution,

$$p_{\mathit{\theta }}\left(x\right)={\displaystyle \frac{exp(-E_{\mathit{\theta }}\left(x\right))}{Z_{\mathit{\theta }}}},\mathrm{where}{Z}_{\mathit{\theta }}={\int }_{\mathcal{X}}exp(-E_{\mathit{\theta }}\left(x\right))dx.$$

(56)

Here, $E_{\mathit{\theta }}\left(x\right)\in \mathbb{R}$ is a learned energy function, typically parameterized by a deep neural network, and $Z_{\theta }$ is the partition function that ensures the density integrates to one. However, this normalization term is generally intractable in high dimensions, rendering exact likelihood evaluation and sampling computationally prohibitive. The parameters of the energy function are found through maximum likelihood estimation, which results in the following objective function,

$$\mathcal{L}\left(\mathit{\theta }\right)={\mathbb{E}}_{x\sim p_{\mathrm{data}}}\left[E_{\mathit{\theta }}\left(x\right)\right]+ln{Z}_{\mathit{\theta }},$$

(57)

with the gradient given by,

$${\nabla }_{\mathit{\theta }}\mathcal{L}\left(\mathit{\theta }\right)={\mathbb{E}}_{x\sim p_{\mathrm{data}}}\left[{\nabla }_{\mathit{\theta }}{E}_{\mathit{\theta }}\left(x\right)\right]-{\mathbb{E}}_{y\sim p_{\mathit{\theta }}}\left[{\nabla }_{\mathit{\theta }}{E}_{\mathit{\theta }}\left(y\right)\right].$$

(58)

The first term is readily computable using observed data, while the second requires samples from the model distribution, which is itself defined implicitly by the energy function. Consequently, much of the complexity in training EBMs arises from the need to draw samples from a model whose normalization is unknown. Markov Chain Monte Carlo (MCMC) methods, such as Langevin dynamics, are commonly used for this purpose [190].

Popular examples of EBMs in literature include Boltzmann Machines [193] and Hopfield Networks [62]. Similar to the other generative models, EBMs have also been applied to time series analysis. Brakel et al. [194] outlined the training strategy for training EBMs for data imputation.Yan et al. [195] proposed ScoreGrad, which a is multivariate probabilistic time series forecasting framework based on continuous energy-based generative models and has been found to achieve state of the art results on a number of datasets.

5.8. Summary of Generative Models

Table 3 presents a comparative summary of the generative models discussed.

6. Uncertainty Quantification

SHM data can be limited or noisy. Noise can emanate from different sources such measurement errors and environmental variability. Limited data might be due to sensor malfunction, missing data, and intermittent data collection due to huge costs associated with the process. Considering these factors, uncertainty quantification is critical to SHM for effective decision making. Uncertainties are generally classified as aleatoric and epistemic [83,196]. Aleatoric uncertainty is associated with inherent randomness and cannot be reduced. According to Wikipedia, epistemic uncertainty is associated with things which one could in principle know but does not know [197]. Thus, once knowledge is available, epistemic uncertainty can be reduced. Epistemic uncertainty is further grouped into two major categories: homoscedastic, meaning it is constant across all observations, or heteroscedastic, meaning it varies with covariates. While homoscedastic uncertainty is relatively straightforward to handle, heteroscedasticity introduces additional complexity into modeling tasks. If these uncertainties are not properly accounted for, model predictions are usually unreliable can may potentially lead to suboptimal decision-making.The goal of this section is therefore to discuss how uncertainties are handled in deep learning. We do this from a Bayesian context.

6.1. Bayesian Inference

In many practical scenarios, prior knowledge plays a critical role in shaping our understanding of a problem. However, conventional deterministic modeling approaches typically disregard this prior information. A significant limitation of such methods is their reliance on large volumes of data to achieve good performance. In reality, data is often limited, incomplete, or corrupted by noise. In the late 1980s and early 1990s, Denker et al. [198], Tishby et al. [199], Denker and LeCun [200], Buntine and Weigend [201], MacKay [202], Neal [203] and Neal [204] proposed the bayesian framework as an alternative learning method. This approach provides a way to incorporate prior knowledge into a model as well estimate uncertainties associated with its outputs. The central idea that governs this philosophy is Bayes’ rule. Mathematically, Bayes’ rule is expressed as

$$p(\mathit{\theta }\mid \mathcal{D})={\displaystyle \frac{p(\mathcal{D}\mid \mathit{\theta })p\left(\mathit{\theta }\right)}{p\left(\mathcal{D}\right)}},$$

(59)

where $p(\mathit{\theta }\mid \mathcal{D})$ is the posterior, $p(\mathcal{D}\mid \mathit{\theta })$ is the likelihood, $p\left(\mathit{\theta }\right)$ is the prior, and the evidence term $p\left(\mathcal{D}\right)$ is given by,

$$p\left(\mathcal{D}\right)=\int p(\mathcal{D}\mid \mathit{\theta })p\left(\mathit{\theta }\right)d\mathit{\theta }.$$

To make predictions for a new data point in a supervised mode, the posterior predictive distribution is used,

$$p(y_{new}\mid x_{new},\mathcal{D})={\int }_{\Theta }p(y_{new}\mid x_{new},\mathit{\theta })p(\mathit{\theta }\mid \mathcal{D})d\mathit{\theta },$$

(60)

With regards to this expression, Murphy [205] writes, "the posterior is our internal belief state about the world and the way to test if our beliefs are justified is to use them to predict objectively observable quantities."

Despite the simplicity of Bayes’ rule, the posterior distribution is often analytically intractable in real-world applications. For deep neural neural networks, the parameter vector $\theta$ may contain millions or even billions of parameters. For this reason, Bayesian methods were historically not favored for training deep models. However, this has changed with the development of efficient approximation techniques. Common methods for parameter estimation are discussed in the following sections.

6.1.1. Analytical Methods (Conjugacy)

For certain combinations of the prior and likelihood, the posterior can be computed analytically. In this case, the posterior distribution is of the same functional form as the prior, and the likelihood-prior pair is known as a conjugate pair [83].Besides providing a closed form solution for the posterior distribution, conjugacy enables sequential learning where the posterior at t becomes the prior at $t+1$. Some of the conjugate pairs include Gaussian-Gamma, Gaussian-Inverse-Chi-Squared, Gaussian-Inverse-Gamma, and Gaussian-Inverse-Wishart [206].

6.1.2. Maximum Likelihood Estimation

Although not Bayesian, Maximum Likelihood Estimation (MLE) is widely used to estimate model parameters. The goal of MLE is to find the set of parameters $\mathit{\theta }\in \Theta$ that maximizes the likelihood of the data,

$${\widehat{\theta }}_{\mathrm{MLE}}=arg\underset{\mathit{\theta }}{max}\prod _{i=1}^{\mathrm{D}}p(y_i\mid \mathit{\theta },x_i).$$

(61)

Equation 61 considers that the data is identically and independently distributed (i.i.d.). For large datasets, the computation of the likelihood can result in numerical instability. Large values of of $p(y_i\mid \mathit{\theta },x_i)$ result in overflow, while small values result in underflow. It is important to note that $p(y_i\mid \mathit{\theta },x_i)$ is a likelihood function, not a probability, and therefore can take values greater than 1.

To improve numerical stability and make the expression amenable to optimization, the log-likelihood is considered instead. Applying the properties of logarithms, the log-likelihood function is given by,

$$lnp(Y\mid \mathit{\theta },X)=\sum _{i=1}^{\mathrm{D}}lnp(y_i\mid \mathit{\theta },x_i).$$

(62)

The above equation is a non-convex function and would otherwise benefit from readily available convex optimization algorithms by considering the negative log likelihood. This, in theory, is possible by taking the negative of the function. Thus, it is common to encounter negative log-likelihood in optimization problems

$$-lnp(Y\mid \mathit{\theta },X)=-\sum _{i=1}^{\mathrm{D}}lnp(y_i\mid \mathit{\theta },x_i).$$

(63)

Maximization of the log-likelihood is similar to the minimization of the negative log-likelihood, and thus to find the optimal parameters,

$${\mathit{\theta }}^{*}=arg\underset{\mathit{\theta }}{min}-lnp(Y\mid \mathit{\theta },X).$$

(64)

A key limitation of MLE is its tendency to overfit, which can be mitigated using priors, leading us to Maximum A Posteriori estimation.

6.1.3. Maximum A Posteriori (MAP)

Bayes’ rule can be rewritten, omitting the marginal likelihood (which is constant w.r.t. $\mathit{\theta }$) as,

$$p(\mathit{\theta }\mid \mathcal{D})\propto p(\mathcal{D}\mid \mathit{\theta })p\left(\mathit{\theta }\right).$$

(65)

With this expression , it is possible to maximize the unnormalized posterior directly, i.e.,

$${\theta }_{\mathrm{MAP}}^{*}=arg\underset{\mathit{\theta }}{max}p(\mathcal{D}\mid \mathit{\theta })p\left(\mathit{\theta }\right).$$

(66)

Taking the log for numerical stability,

$${\mathit{\theta }}_{\mathrm{MAP}}^{*}=arg\underset{\mathit{\theta }}{max}lnp(\mathcal{D}\mid \mathit{\theta })+lnp\left(\mathit{\theta }\right).$$

(67)

The second term serves as a regularization term, and this form is structurally similar to penalized optimization.

6.1.4. Laplace Approximation

The Laplace approximation is a technique used to approximate a posterior distribution with a multivariate Gaussian centered at the maximum a posteriori (MAP) estimate or, in some cases, the maximum likelihood estimate (MLE). It provides a way to account for uncertainty in the maximum likelihood or maximum a posteriori parameter estmates. Given a model with parameters $\mathit{\theta }\in {\mathbb{R}}^d$, and a log-posterior function $lnp(\mathit{\theta }\mid \mathcal{D})$, the Laplace approximation takes the following form,

$$p(\mathit{\theta }\mid \mathcal{D})\approx \mathcal{N}({\mathit{\theta }}^{*},\Sigma ),$$

(68)

where ${\mathit{\theta }}^{*}$ is the MAP estimate given by equation 67 and the covariance matrix $\Sigma$ is approximated by the inverse of the negative Hessian of the log-posterior evaluated at ${\mathit{\theta }}^{*}$,

$$\Sigma \approx {\left[-{\nabla }^{2}lnp(\mathit{\theta }\mid \mathcal{D}){|}_{\mathit{\theta }={\mathit{\theta }}^{*}}\right]}^{-1}.$$

(69)

This approach captures the local curvature of the posterior distribution around the mode, thereby providing a second-order approximation of the distribution [83]. However, computing the full Hessian matrix is often computationally expensive, especially in high-dimensional parameter spaces. As a result, various approximations and low-rank techniques have been proposed to reduce this burden.

6.1.5. Expectation Maximization

The procedures for MLE and MAP discussed so far assume that the dataset is complete and that all variables associated with the model are fully observed and known. However, in many practical applications, models often involve partial data or incorporate latent variables that are not directly observable. In such scenarios, the complete log-likelihood must account for these hidden variables, necessitating marginalization over them to obtain the observed data likelihood. Our derivation of the E-M algorithm is based on the work of Gao [207].

Formally, let $\mathcal{D}={\left\{x_i\right\}}_{i=1}^N$ represent the observed dataset, $Z={\left\{z_i\right\}}_{i=1}^N$ the corresponding latent variables, and $\mathit{\theta }$ the parameters of the model. The complete-data log-likelihood is given by

$$logp(\mathcal{D},Z\mid \mathit{\theta }).$$

(70)

However, since the latent variables Z are unobserved, the quantity of interest becomes the marginal log-likelihood

$$logp(\mathcal{D}\mid \mathit{\theta })=log\int p(\mathcal{D},Z\mid \mathit{\theta })dZ,$$

(71)

where the integral is taken over all possible configurations of the latent variables. In models containing a large number of latent variables, this integral typically becomes intractable, either due to high dimensionality or the complexity of the joint distribution $p(\mathcal{D},Z\mid \mathit{\theta })$.

To address this intractability, an arbitrary probability density function $q\left(Z\right)$ over the latent variables is introduced. By doing so, the marginal log-likelihood can be rewritten as

$$logp(\mathcal{D}\mid \mathit{\theta })=log\int q\left(Z\right){\displaystyle \frac{p(\mathcal{D},Z\mid \mathit{\theta })}{q\left(Z\right)}}dZ.$$

(72)

This re-expression enables the application of Jensen’s inequality, exploiting the concavity of the logarithm function to derive a lower bound. Specifically, Jensen’s inequality yields

$$log{\mathbb{E}}_q\left[{\displaystyle \frac{p(\mathcal{D},Z\mid \mathit{\theta })}{q\left(Z\right)}}\right]\ge {\mathbb{E}}_q\left[log{\displaystyle \frac{p(\mathcal{D},Z\mid \mathit{\theta })}{q\left(Z\right)}}\right],$$

(73)

which leads to

$$logp(\mathcal{D}\mid \mathit{\theta })\ge {\mathbb{E}}_{q\left(Z\right)}\left[logp(\mathcal{D},Z\mid \mathit{\theta })-logq\left(Z\right)\right].$$

(74)

Defining the functional,

$$\mathcal{L}(q,\mathit{\theta })={\mathbb{E}}_{q\left(Z\right)}\left[logp(\mathcal{D},Z\mid \mathit{\theta })\right]+H\left(q\right),$$

(75)

where $H\left(q\right)=-{\mathbb{E}}_{q\left(Z\right)}[logq\left(Z\right)]$ denotes the entropy of $q\left(Z\right)$, the above inequality can be compactly written as,

$$logp(\mathcal{D}\mid \mathit{\theta })\ge \mathcal{L}(q,\mathit{\theta }).$$

(76)

$\mathcal{L}(q,\mathit{\theta })$ serves as a surrogate objective function that can be maximized in place of the intractable marginal log-likelihood.

Expectation-Maximization (EM) algorithm is a classical procedure for maximizing the marginal likelihood by iteratively optimizing the lower bound $\mathcal{L}(q,\mathit{\theta })$. EM algorithm proceeds in two alternating steps, referred to as the Expectation step (E-step) and the Maximization step (M-step).

In the E-step, given the current estimate of the model parameters ${\mathit{\theta }}^{\left(t\right)}$, the function $q\left(Z\right)$ is set to the posterior distribution of the latent variables conditioned on the observed data and the current parameters,

$$q\left(Z\right)=p(Z\mid \mathcal{D},{\mathit{\theta }}^{\left(t\right)}).$$

(77)

Under this choice, the lower bound becomes tight, satisfying

$$\mathcal{L}(q,{\mathit{\theta }}^{\left(t\right)})=logp(\mathcal{D}\mid {\mathit{\theta }}^{\left(t\right)}),$$

(78)

thus removing the looseness introduced by Jensen’s inequality. The E-step therefore consists of computing the expected complete log-likelihood,

$$Q(\mathit{\theta },{\mathit{\theta }}^{\left(t\right)})={\mathbb{E}}_{Z\sim p(Z\mid \mathcal{D},{\mathit{\theta }}^{\left(t\right)})}\left[logp(\mathcal{D},Z\mid \mathit{\theta })\right].$$

(79)

Following the E-step, the M-step seeks to update the model parameters by maximizing this expected complete log-likelihood

$${\theta }^{(t+1)}=arg\underset{\mathit{\theta }}{max}Q(\mathit{\theta },{\mathit{\theta }}^{\left(t\right)}).$$

(80)

6.1.6. Monte Carlo Integration

Bayesian inference often requires evaluating integrals like the marginal likelihood, $\int p(\mathcal{D}\mid \mathit{\theta })p\left(\mathit{\theta }\right)d\mathit{\theta }$, which is equivalent to an expectation, ${\mathbb{E}}_{p\left(\mathit{\theta }\right)}\left[p(\mathcal{D}\mid \mathit{\theta })\right]$. This can be approximated by sampling,

$${\mathbb{E}}_{p\left(\mathit{\theta }\right)}\left[p(\mathcal{D}\mid \mathit{\theta })\right]\approx {\displaystyle \frac{1}{\mathrm{N}}}\sum _{i=1}^{\mathrm{N}}p(\mathcal{D}\mid {\theta }_i),$$

(81)

where $\mathrm{N}$ is the number of samples.Monte Carlo integration is useful in high-dimensional spaces where analytical solutions are infeasible. The accuracy improves with more samples.

6.1.7. Importance Sampling

The typical monte carlo assumes that all samples are equally likely, which may not hold for heavily-tailed distributions such as the Gaussian. In such cases, a more effective technique is needed to assign higher weights to high-probability points and lower weights to low-probability points in order to provide accurate estimates of integrals. This can be achieved by modifying the marginal likelihood integral as follows,

$$p\left(\mathcal{D}\right)={\int }_{\Theta }p(\mathcal{D}\mid \theta ){\displaystyle \frac{p\left(\theta \right)}{q\left(\theta \right)}}q\left(\theta \right)d\theta ,$$

(82)

with $q\left(\theta \right)$ being a simple distribution from which samples can be drawn efficiently. This transformation allows for reweighting samples according to the ratio ${\displaystyle \frac{p\left(\theta \right)}{q\left(\theta \right)}}$. The integral is now written as an expectation with respect to $q\left(\theta \right)$

$$p\left(\mathcal{D}\right)={\mathbb{E}}_{q\left(\theta \right)}\left[p(\mathcal{D}\mid \theta ){\displaystyle \frac{p\left(\theta \right)}{q\left(\theta \right)}}\right].$$

(83)

We approximate this expectation using Monte Carlo sampling. Drawing ${\theta }_i\sim q\left(\theta \right)$ for $i=1,\dots ,\mathrm{N}$, the marginal likelihood is estimated as

$$\widehat{p}\left(\mathcal{D}\right)={\displaystyle \frac{1}{N}}\sum _{i=1}^{\mathrm{N}}p(\mathcal{D}\mid {\theta }_i){\displaystyle \frac{p\left({\theta }_i\right)}{q\left({\theta }_i\right)}}.$$

(84)

6.1.8. Variational Inference

Variational Inference (VI) is a method for approximating complex probability distributions in cases where exact computation is intractable. The goal is to approximate a posterior distribution $p(\mathit{\theta }\mid \mathcal{D})$ with a simpler variational distribution $q\left(\mathit{\theta }\right)$.Such a variational distribution is considered close to the posterior distribution. Closeness is usually achieved by minimization of the Kullback-Leibler (KL) divergence between the variational distribution and the true posterior,

$$\mathrm{KL}(q\left(\mathit{\theta }\right)\Vert p(\mathit{\theta }\mid \mathcal{D}))=\int q\left(\mathit{\theta }\right)ln{\displaystyle \frac{q\left(\mathit{\theta }\right)}{p(\mathit{\theta }\mid \mathcal{D})}}d\mathit{\theta }.$$

(85)

Minimizing this divergence is equivalent to maximizing the Evidence Lower Bound (ELBO), given by,

$$\mathcal{L}\left(q\right)={\mathbb{E}}_{q\left(\mathit{\theta }\right)}\left[lnp(\mathcal{D},\mathit{\theta })-lnq\left(\mathit{\theta }\right)\right].$$

(86)

The ELBO essentially transforms a bayesian inference problem into an optimization problem and can thus be optimized using gradient-based methods, allowing efficient approximation of the posterior $p(\mathit{\theta }\mid \mathcal{D})$ by $q\left(\mathit{\theta }\right)$. While VI is typically faster than sampling-based methods, it is not guaranteed to converge to the true posterior.

6.1.9. Markov Chain Monte Carlo

Markov Chain Monte Carlo (MCMC) methods use a Markov chain to generate samples from an unknown posterior distribution. The chain satisfies the Markov property, where the next state depends only on the current state,

$$p\left({\theta }_{1:\mathrm{T}}\right)=p\left({\theta }_1\right)\prod _{t=2}^{\mathrm{T}}p({\theta }_t\mid {\theta }_{t-1}).$$

(87)

A large number of samples are generated using the chain, which is designed so that its stationary distribution is the target posterior. Perturbations are added iteratively (as in a random walk), and new states are accepted with a certain probability. The stationary distribution $p^{*}\left({\mathit{\theta }}^{\prime }\right)$ satisfies,

$$p^{*}\left({\mathit{\theta }}^{\prime }\right)=\int p({\mathit{\theta }}^{\prime }\mid \mathit{\theta })p^{*}\left(\mathit{\theta }\right)d\mathit{\theta }.$$

(88)

If the chain is ergodic, this stationary distribution is unique and reachable from any initial distribution.

Metropolis-Hastings Algorithm

The Metropolis-Hastings (MH) algorithm is an MCMC approach that generates samples from a proposal (or transition) distribution $q({\mathit{\theta }}^{\prime }\mid \mathit{\theta })$, with acceptance based on a probability criterion. The MH algorithm was first proposed in 1950, and used a symmetric proposal distribution in its initial form. By symmetric, we mean that the probability of moving to state i is equal to the probability of returning to the starting state. This symmetric form is expressed as

$$q({\mathit{\theta }}^{\prime }\mid \mathit{\theta })=q(\mathit{\theta }\mid {\mathit{\theta }}^{\prime }).$$

(89)

Hastings [208] improved upon the algorithm by incorporating asymmetric proposal distributions.

$$q({\mathit{\theta }}^{\prime }\mid \mathit{\theta })\ne q(\mathit{\theta }\mid {\mathit{\theta }}^{\prime }).$$

(90)

The MH algorithm uses the following criteria to move to new states: if the proposed sample is greater than the previous sample, the new sample is accepted. And if the proposed sample is less than the current sample, the new sample is accepted with an acceptance probability

$$\alpha =min\left(1,{\displaystyle \frac{p\left({\mathit{\theta }}^{\prime }\right)q(\mathit{\theta }\mid {\mathit{\theta }}^{\prime })}{p\left(\mathit{\theta }\right)q({\mathit{\theta }}^{\prime }\mid \mathit{\theta })}}\right).$$

(91)

This recursive process continues until the chain reaches a stationary distribution. One drawback of the MH algorithm is that random sampling can result in slow convergence and correlated samples. The subsequent method, Hamiltonian Monte Carlo (HMC), addresses this issue by incorporating gradient information to enhance sampling efficiency.

Hamiltonian Monte Carlo (HMC)

The Hamiltonian Monte Carlo (HMC) is a sampling-based parameter estimation approach that uses a systematic approach to sampling. To enable faster and efficient exploration of the target distribution, HMC leverages its gradient information. HMC employs Hamiltonian mechanics, incorporating both position q and momentum p. The core principle is to simulate the trajectory of a particle moving through a potential energy landscape defined by the target distribution.

In HMC, the potential energy $U\left(q\right)$ is related to the target distribution $p\left(q\right)$. The kinetic energy $K\left(p\right)$ and Hamiltonian $H(q,p)$ are expressed as,

$$\begin{array}{cc}\hfill H(q,p)& =U\left(q\right)+K\left(p\right),\hfill \\ \hfill K\left(p\right)& ={\displaystyle \frac{1}{2}}ln\left({\left(2\pi \right)}^D\left|M\right|\right)+{\displaystyle \frac{1}{2}}{p}^T{M}^{-1}p,\hfill \end{array}$$

(92)

where M is the mass matrix. The Hamiltonian governs the dynamics of the system via Hamilton’s equations,

$${\displaystyle \frac{dq}{dt}}={\displaystyle \frac{\partial H(q,p)}{\partial p}},{\displaystyle \frac{dp}{dt}}=-{\displaystyle \frac{\partial H(q,p)}{\partial q}}.$$

(93)

To approximate these dynamics, HMC uses the leapfrog integrator. The updates for momentum and position are,

$${\displaystyle \frac{dq}{dt}}={\displaystyle \frac{\partial H(q,p)}{\partial p}},{\displaystyle \frac{dp}{dt}}=-{\displaystyle \frac{\partial H(q,p)}{\partial q}}.$$

(94)

Finally, the acceptance probability for the proposed state $(q_{\mathrm{new}},p_{\mathrm{new}})$ is calculated as,

$$P\left(\mathrm{accept}{q}_{\mathrm{new}}\right)=min\left(1,{\displaystyle \frac{exp(-H(q_{\mathrm{new}},p_{\mathrm{new}}))}{exp(-H(q,p\left)\right)}}\right).$$

(95)

More extensive treatments of HMC and its variations can be found in Marwala et al. [209]

6.2. Bayesian Neural Networks

6.2.1. Overview

Within a Bayesian context, model parameters are treated as distributions rather than point estimates [210]. The objective of a Bayesian Neural Network(BNN) is to learn these distributions. To this end, different methods have been developed to quantify uncertainty in neural networks. The discussions in this section will focus on 4 methods and the rationale for this choice is simplicity and similarity of each approach to the classic backpropagation. The first approach is discussed mainly because it has been extensively applied to SHM and can easily be scaled.

6.2.2. Tractable Approximate Gaussian Inference

The Tractable Approximate Gaussian Inference (TAGI) [211] is a probabilistic method that allows for the analytical inference of model parameters without relying on backpropagation. TAGI achieves this through a two-step process: propagation of the model uncertainty using moment functions and local linearization of the activation function. Parameters are then updated using a modified Rauch–Tung–Striebel (RTS) algorithm [211,212]. TAGI makes some assumptions about network parameters, layers, and units within layers. Parameters and units within a layer are considered i.i.d. and normally distributed. TAGI also uses the inherent conditional independence between hidden layers, i.e. ${\mathbf{Z}}^{(j-1)}\perp {\mathbf{Z}}^{(j+1)}\mid {\mathbf{z}}^{\left(j\right)}$.

Consider a neural network with an observation model described by

$$y=Z^{\left(0\right)}+v,\mathrm{where}v\sim \mathcal{N}(0,{\sigma }_V^2).$$

(96)

Mathematically, the i-th hidden state of the $(j+1)$-th layer can be written as

$$Z_i^{(j+1)}=\sum _{k=1}^A{\mathbf{W}}_{\mathbf{i},\mathbf{k}}^{\left(\mathbf{j}\right)}{A}_k^{\left(j\right)}+{\mathbf{B}}_{\mathbf{i}}^{\left(\mathbf{j}\right)},$$

(97)

where $Z_i^{(j+1)}$ is the hidden state of the $(j+1)$-th layer, ${\mathbf{W}}_{\mathbf{i},\mathbf{k}}^{\left(\mathbf{j}\right)}$ is the weight from the k-th neuron in layer j to the i-th neuron in layer $j+1$, $A_k^{\left(j\right)}$ is the activation of the k-th neuron in layer j, and ${\mathbf{B}}_{\mathbf{i}}^{\left(\mathbf{j}\right)}$ is the bias term for the i-th neuron in layer $j+1$.

To propagate uncertainty, TAGI uses the Gaussian Multiplicative Approximation (GMA) given by the following equations

$$\begin{array}{cc}\hfill \mathbb{E}\left[X_1{X}_2\right]& ={\mu }_1{\mu }_2+cov(X_1,X_2),\hfill \end{array}$$

(98)

$$\begin{array}{cc}\hfill cov(X_3,X_1{X}_2)& =cov(X_1,X_3){\mu }_2+cov(X_2,X_3){\mu }_1,\hfill \\ \hfill cov(X_1{X}_2,X_3{X}_4)& =cov(X_1,X_3)cov(X_2,X_4)+cov(X_1,X_4)cov(X_2,X_3)\hfill \\ \hfill & +cov(X_1,X_3){\mu }_2{\mu }_4+cov(X_1,X_4){\mu }_2{\mu }_3\hfill \end{array}$$

(99)

$$\begin{array}{cc}\hfill & +cov(X_2,X_3){\mu }_1{\mu }_4+cov(X_2,X_4){\mu }_1{\mu }_3,\hfill \\ \hfill var\left(X_1{X}_2\right)& ={\sigma }_1^2{\sigma }_2^2+cov{(X_1,X_2)}^2+2cov(X_1,X_2){\mu }_1{\mu }_2\hfill \end{array}$$

(100)

$$\begin{array}{cc}\hfill & +{\sigma }_1^2{\mu }_2^2+{\sigma }_2^2{\mu }_1^2.\hfill \end{array}$$

(101)

While the sum of Gaussian variables remains Gaussian, the same is not generally true for their product. Equation 97 poses this challenge, as it involves the product of two random variables, i.e., ${\mathbf{W}}_{\mathbf{i},\mathbf{k}}^{\left(\mathbf{j}\right)}$ and $A_k^{\left(j\right)}$.

However, under the GMA and with a sufficiently large number of hidden units, the resulting distribution can be shown to converge to a Gaussian, ${\mathbf{Z}}_{\mathbf{i}}^{(j+1)}\sim \mathcal{N}({\mu }_{{\mathbf{Z}}_{\mathbf{i}}^{(j+1)}},{\Sigma }_{{\mathbf{Z}}_{\mathbf{i}}^{(j+1)}}\mathrm{I})$.

Thus, TAGI is theoretically capable of approximating the moments of the hidden states. It is important to note that, although the propagation of information from the input layer to the output layer utilizes nonlinear activation functions, TAGI employs a locally linearized version of this transformation for the hidden state inputs. This linearization facilitates operations involving Gaussian random variables since a linear transformation of a Gaussian random variable results in another Gaussian random variable. TAGI achieves this by locally linearizing the nonlinear activation function about ${\mu }_Z$ using a first-order Taylor approximation. It is crucial to emphasize that this local linearization does not imply the use of a linear activation function. Instead, for each input covariate $x_i$, the linearization is performed at different ${\mu }_Z$, thereby preserving the nonlinear relationship between the inputs $x_i$ and the outputs $y_i$.

Since TAGI assumes that the model parameters and hidden states are independent, it is computationally efficient, operating with a complexity of $\mathcal{O}\left({\mathrm{A}}^2\right)$, where $\mathrm{A}$ is the number of hidden units per layer, and scales linearly with $\mathrm{L}$ hidden layers . So far, numerous enhancements have been made since the original implementation of TAGI. For instance, Deka et al. [213] extended TAGI to account for heteroscedastic uncertainty (Figure 13).

6.2.3. Learned Observation Noise in TAGI

Instead of treating ${\sigma }_V^2$ as a fixed term, it is now learned from the data using the following equations,

Figure 13. A compact representation of the network in Figure 3 coupled with a directed acyclic graph (DAG) for the estimation of the variance parameter [100].

Figure 13. A compact representation of the network in Figure 3 coupled with a directed acyclic graph (DAG) for the estimation of the variance parameter [100].

$$\begin{array}{cc}\hfill {\mu }_{V^2\mid y}& ={\mu }_{V\mid y}^2+{\sigma }_{V\mid y}^2,\hfill \end{array}$$

(102)

$$\begin{array}{cc}\hfill {\sigma }_{V^2\mid y}^2& =2{\sigma }_{V\mid y}^4+4{\sigma }_{V\mid y}^2{\mu }_{V\mid y}^2,\hfill \end{array}$$

(103)

$$\begin{array}{cc}\hfill {\mu }_{\overline{V^2}\mid y}& ={\mu }_{\overline{V^2}}+k({\mu }_{V^2\mid y}-{\mu }_{V^2}),\hfill \end{array}$$

(104)

$$\begin{array}{cc}\hfill {\sigma }_{\overline{V^2}\mid y}^2& ={\sigma }_{V^2}^2+k^2({\sigma }_{V^2\mid y}^2-{\sigma }_{V^2}^2),\hfill \end{array}$$

(105)

$$\begin{array}{cc}\hfill k& ={\displaystyle \frac{{\sigma }_{\overline{V^2}}^2}{{\sigma }_{V^2}^2}}.\hfill \end{array}$$

(106)

where

$$\begin{array}{cc}\hfill {\mu }_{\overline{V^2}}& =exp\left({\mu }_{\overline{V^2}}+0.5{\sigma }_{\overline{V^2}}^2\right),\hfill \\ \hfill {\sigma }_{V^2}^2& =exp\left(2{\mu }_{\overline{V^2}}+{\sigma }_{\overline{V^2}}^2\right)·\left(exp\left({\sigma }_{\overline{V^2}}^2\right)-1\right),\hfill \\ \hfill cov(\overline{V^2},\widetilde{V^2})& ={\sigma }_{\overline{V^2}}^2·{\mu }_{\overline{V^2}},\hfill \end{array}$$

and

$${\mu }_{V^2}={\mu }_{\overline{V^2}},{\sigma }_{V^2}^2=3{\sigma }_{\overline{V^2}}^2+2{\mu }_{\overline{V^2}}^2.$$

Note:

$$\begin{array}{cc}\hfill f\left(v\right)& =\mathcal{N}(v;0,{\mu }_{\overline{V^2}}),\hfill \\ \hfill f(v^2\mid \overline{v^2})& =\mathcal{N}(v^2;\overline{v^2},2{\overline{v^2}}^2),\hfill \\ \hfill \overline{V^2}& \sim \mathcal{N}(\overline{v^2};{\mu }_{\overline{V^2}},{\sigma }_{\overline{V^2}}^2).\hfill \end{array}$$

6.2.4. Further TAGI Extensions

TAGI has also been successfully extended to address problems involving sequential data, particularly in time series modeling. Recently, Vuong et al. [214] developed the theory for modeling uncertainty in LSTMs and Gated Recurrent Units (GRUs) using TAGI as the foundational learning algorithm.

6.2.5. Monte Carlo Dropout

Dropout [215] is a method to mitigate overfitting in neural networks. During model training, dropout randomly freezing a subset of nodes (typically 50%) by setting their values to zero. The overall effect of this is that the network is less dependent on any individual weight. Once the model is trained, dropout is no longer applied.

Monte Carlo (MC) [216] dropout reinterprets dropout as a method for performing approximate Bayesian inference. A similar dropout procedure is applied during inference, and multiple forward passes are performed on the same input. This repeated stochastic sampling allows estimation of the predictive uncertainty by treating each forward pass as a sample from the posterior predictive distribution.

In principle, the final output of an MC dropout is considered an ensemble of outputs from different stochastic realizations of the network.

MC dropout is attractive due to its simplicity in implementation. However, it is computationally expensive during inference, as multiple forward passes must be performed. Additionally, the uncertainty estimates produced by MC dropout are often poorly calibrated.

6.2.6. Bayes by Backpropagation

As already highlighted, the goal of a BNN is to learn distributions of the network weights. The posterior distribution of these weights, $p(\theta \mid D)$, is generally intractable since neural networks contain a huge number of weights. This renders the marginal likelihood $p\left(D\right)$ intractable as well. As already discussed, the standard variational approach employs a more tractable distribution $q\left(\theta \right)$ to approximate $p(\theta \mid D)$; with the Gaussian distribution being a common choice. As is always the case in VI, the goal is to make $q\left(\theta \right)$ as close as possible to $p(\theta \mid D)$. This is achieved through minimization of the Kullback-Leibler divergence. We won’t discuss the full details of the derivation of the loss function used in Bayes by backprop as it can be easily derived using same principles used in Section 5.2 and Section 6.1.8. Bayes by backprop [217] uses a similar technique to VAE (reparamaterization trick) to express the gradient of an expectation as an expectation of the gradient. This allows us to compute gradients based on the reparameterization parameters and use standard backpropagation to update the network weights.

6.2.7. Probabilistic Backpropagation

Due to its reliance on sampling during both training and inference, the Bayes by backprop is not scalable. An alternative, scalable approach is the probabilistic backpropagation (PBP) [218], which was developed by researchers at Havard. PBP models each weight as a univariate Gaussian. These weights are learned via forward and backward passes through the network, just like in classic backpropagation. During the forward pass, input is propagated through the network. However, since the weights are random, activations at each layer are also random, though this time more complex and intractable. PBP approximates these intractable distributions with tractable Gaussians using moment matching. At the end of the forward pass, PBP uses the marginal probability of the target variable to measure performance and the gradient of this expression is computed with respect to the means and variances of the approximate distributions. The gradient is then propagated back to update the means and variances of the weights. Unlike classic backpropagation, PBP uses Assumed Density Filtering (ADF) for the parameter updates. ADF is a principled sequential way of minimizing the Kullback-Leibler divergence and is useful in PBP considering that we are building the posterior sequentially from more tractable distributions.

7. Applications

Deep learning for time series have been employed in various applications within SHM. We conducted a systematic literature search to identify studies on deep learning for time series analysis in SHM. The following search querry was used to retrieve articles from the scopus database:

(TITLE-ABS-KEY ("deep learning" OR "neural network" OR "variational autoencoder" OR "VAE" OR "recurrent neural network" OR "RNN" OR "long short-term memory" OR "LSTM" OR "gated recurrent unit" OR "GRU" OR "convolutional neural network" OR "CNN" OR "temporal convolutional network" OR "TCNN" OR "diffusion model" OR "normalizing flow" OR "physics-informed neural network" OR "PINN" OR "Bayesian neural network" OR "BNN" OR "multilayer perceptron" OR "MLP" OR "transformer") AND TITLE-ABS-KEY ("time series*" OR "temporal" OR "forecast*" OR "causal*") AND TITLE-ABS-KEY ("structural health monitoring" OR "SHM")) AND PUBYEAR > 2019 AND PUBYEAR < 2025 AND (LIMIT-TO(SRCTYPE, "j")) AND (LIMIT-TO(PUBSTAGE, "final")) AND (LIMIT-TO(SUBJAREA, "ENGI") OR LIMIT-TO(SUBJAREA, "MULT") OR LIMIT-TO(SUBJAREA, "ENVI") OR LIMIT-TO(SUBJAREA, "EART")) AND (LIMIT-TO(DOCTYPE, "ar")) AND (LIMIT-TO(LANGUAGE, "English"))

Figure 14 summarizes the entire process, from the initial database search, through exclusion of irrelevant studies, to the identification of articles included in the final review.The exclusion criteria were: (a) articles that do not address time series in SHM for civil infrastructure; (b) articles not written in English; and (c) articles published before 2020 or in 2025.

Based on the identified studies, the applications were categorized into five thematic areas: damage assessment, structural response prediction, structural load prediction, data reconstruction, and anomaly detection. Each of these areas is explored in greater detail in the following sections.

7.1. Damage Assessment

Damage assessment includes the detection, localization, and quantification of structural deterioration. A wide range of approaches have been proposed to address this challenge. A hierarchical CNN-GRU framework was developed to exploit both spatial and temporal information for structural damage detection in bridge datasets [219]. Recognizing the importance of feature extraction, another study proposed a channel-spatial-temporal attention network to refine sample-specific features across multiple dimensions [220]. Offshore jacket platforms have also been a subject of interest, where a CNN-BiLSTM-Attention model, optimized via particle swarm optimization, demonstrated more than 95% accuracy in damage detection [221].

In vibration-based damage diagnosis, several studies have focused on improving attention mechanisms and deep feature learning. A multi-head self-attention LSTM autoencoder was introduced for unsupervised diagnosis and quantification of damage in ambient vibration data [222]. The Transformer architecture was leveraged for post-disaster damage state classification using seismic data [223], and a custom 1DCNN-BiLSTM model with an Inception module was used to detect minor changes in RC beams [224].

Additional studies emphasized the localization of structural damage using various techniques. A time-varying damage index was proposed in conjunction with a 1D-CNN to enhance Lamb wave-based localization [225], while an autoencoder-based unsupervised model enabled accurate damage localization using raw acceleration data [226]. Localization has also been explored using acoustic emissions and arrival time distributions via LSTM and Bayesian methods [227], and through Transformer-based learning of multivariate vibration signals for submerged offshore wind turbine structures [228]. Damage classification using CNN and time-frequency analysis has been achieved through continuous wavelet transform representations of acceleration data [229], and Trident, a ConvLSTM3D model, achieved high-resolution damage identification in bridges [230]. Further work addressed unsupervised real-time detection using statistical modeling [231], sensor signal fusion via a 1D-CNN-LSTM hybrid [232], and the integration of Burg Autoregressive features with stacked autoencoders [233]. Multiclass damage detection and classification have also been explored. A 1D CNN method using windowed time series was validated for a full-scale bridge [234], while another study introduced a CNN model with synchrosqueezing transform to detect damage in the Z24 bridge [235]. Others investigated transformer-based models [236], metaheuristic optimization [237], and GAN-based signal augmentation [238]. Stochastic configuration networks [239] and hybrid CNN-RNN methods [240] have further enriched this task.

7.2. Structural Response Prediction

Predicting structural responses such as displacement, strain, stress, and vibration is fundamental in SHM. Deep time series have shown exceptional ability in capturing nonlinear temporal patterns, making them well-suited for this task. Seismic and environmental responses have been a primary focus. A Dung Beetle Optimization-enhanced BiLSTM model, integrated with discrete wavelet transforms, was used for seismic response prediction in slope reinforcement structures [241]. A dual-stage attention LSTM (DALSTM) encoder-decoder model was used for displacement forecasting in arch dams [242]. Another study used STL-extra-trees-LSTM modeling for a similar task [243].

Tunnel and bridge response prediction has also received extensive attention. LSTM-based approaches have been adopted for modeling stress and strain during shield tunneling [244], deflection predictions under vehicle and thermal loads [245,246] , and buffeting response forecasting in a bridge under aerodynamic effects [247]. An ensemble RNN with short-sequence LSTMs has been proposed for high-frequency health monitoring [248], while LSTM and CNN hybrid models have been introduced for rail track systems [249] and long-span bridge deck wind-induced lateral displacement responses [250]. Temperature-induced structural responses have been modeled in numerous studies. These include CNN-based prediction of long-term strain using weather data [251,252], a neural network based on hydrostatic-temperature-time (HTT) for predicting concrete dam displacement [253,254], and HTT enhancements via the DeepLift framework for dams in extreme climates [255]. BiLSTM models have also been used to map temperature fields to strain distributions [256]. Multi-factor response models have been implemented using attention-enhanced BiLSTM frameworks [257] and transformer-based architectures for long-term structural state forecasting and strain prediction [258,259]. Physics-informed models such as Phy-Seisformer [260] and simplified rheological neural networks [261] further emphasize integration with engineering principles. The use of CNNs for spatiotemporal deformation modeling of arch dams[262], and structural performance modeling of underwater tunnels via GC-GRU [263] demonstrate the breadth of this thematic area. Studies have also addressed low-cost strain prediction techniques [264], FEM-integrated RNN approaches for high-rise buildings [265], and global-partial BiLSTM models correlating cable tension and bridge deflection [266].

7.3. Structural Load Prediction

Load-related studies have focused on temperature prediction using ANN and LSTM models. For example, environmental data from the Shanghai Tower was used to predict curtain wall temperature fields using ANN models [10], while a similar deep learning model was used for the Jinping-I arch dam [267]. LSTM models have also been developed for bridge temperature forecasting using meteorological big data [268].

Electric arc furnace temperature, as a proxy for structural load, was forecasted using GRU models with multivariate inputs [269,270], and further enhanced by attention mechanisms in encoder-decoder frameworks [271].

Wind, another significant load, is often predicted using deep learning models integrated with SHM and meteorological data, including adaptive residual CNNs [272], LSTM with empirical mode decomposition [273], and other hybrid methods [274].

7.4. Data Reconstruction

Data loss is a persistent challenge in SHM, and many studies have focused on reconstructing missing sensor data through advanced deep learning models. BiLSTM and encoder-decoder models with attention mechanisms have been central to this effort. BiLSTM was used to reconstruct acceleration responses in a long span cable-stayed bridge [275] and impute missing dam sensor data using attention and transfer learning [276]. Hybrid models that combined CNN with BiGRU or GRU were developed for strain recovery [277], accelerometer data restoration [278], and recovery of multi-modal heterogeneous signals [279].

GAN-based architectures have also played a key role. A BiLSTM-GAN model was used for bidirectional data imputation [280], while another study applied GANs to restore missing strain data on a concrete bridge [281]. Data augmentation via GANs was also explored for AE signal enhancement [238] , while other studies reconstructed vibration responses for bridge damage detection using ED-DCNNs and attention modules [282].

Additional frameworks used AOA-TCN [283], EMD-LSTM [284], ESMD-PE-BiGRU [285], SSA-optimized BiLSTM [286], and Kriging-based interpolation [287]. External feedback was integrated in an RNN-based approach for response recovery [288]. Meanwhile, CycleGAN [289] and domain adaptation strategies were used for synthetic-to-real transformations, improving model generalizability.

7.5. Anomaly Detection

Anomaly detection serves as an early warning mechanism for structural faults, sensor failures, or environmental interferences. Various models have been proposed to distinguish between normal and anomalous behavior.

The Temporal Fusion Transformer was used to detect anomalies in historic buildings through vibrational analysis [290], while time-series attention models and encoder-decoder frameworks were used to flag abnormal seismic responses [291]. Another encoder-decoder framework integrated temperature-driven PCA for anomaly detection in long-span bridges [292].

Other techniques include semi-supervised anomaly detection using MixMatch technique [293], pattern-recognition neural networks for multi-type data detection anomalies [294], and multi-class imbalanced anomaly detection using 1D CNNs [295]. CNNs were also used to detect abnormal temperature patterns by transforming time-series data into Gramian Angular Field images [296].

Physics-informed models, such as Koopman neural ODEs, were developed for anomaly detection alongside structural forecasting [297], while sparse Bayesian ELMs were used for outlier detection in railway track systems [298], and simple ANNs for SHM anomaly detection were explored in the context of IFC-BIM data integration [299].

Table 4. Deep Learning Applications in Structural Health Monitoring.

Application	Scope (Specific Papers)	Deep Learning Architectures Used
Damage Assessment	Damage detection: [300,301,302,303,304,305,306,307,308,309,310,311,312,313]Damage localization: [314,315,316,317,318,319,320]Damage classification: [321,322,323,324]Damage progression prediction: [325]	CNN, GRU, LSTM, BiLSTM, Autoencoder, Transformer, CNN-RNN hybrids
Structural Response Prediction	Strain prediction: [326,327]Displacement/Deflection: [328,329,330,331,332,333,334,335]Seismic and vibration: [336,337,338,339,340]Thermal-induced: [341]Tunnel responses: [342,343,344]Mechanical/Stress: [326,334,340,345,346,347]Cable tension: [348]	LSTM, BiLSTM, CNN, GRU, Attention, Transformer-based models
Structural Load Prediction	Dynamic load: [349,350]	CNN, Bayesian optimization, Autoencoder, CNN-BiLSTM hybrids
Data Reconstruction	Wind: [287,351,352,353]Vibration: [354,355,356]Temperature: [357]Dam monitoring: [358]	CNN, BiLSTM, GAN, Autoencoder, CNN-GRU hybrids, VMD, EMD
Anomaly Detection	Sensor faults: [359,360,361,362,363]Outliers: [363,364,365]	CNN, LSTM, BiLSTM, FCN, Transformer, PCA
Sensor Placement	Optimized placement: [336,349]	Attention-based RNN, CNN-BiLSTM
Data Augmentation/Generation	Synthetic data: [366]	GAN, CycleGAN, CNN, BiLSTM
Other SHM Tasks	Traffic classification: [367]Data compression: [368]Structural state ID: [369]Defect diagnosis: [370]	CNN, Autoencoder, Transformer, CNN-RNN hybrids

Table 4. Deep Learning Applications in Structural Health Monitoring.

Application	Scope (Specific Papers)	Deep Learning Architectures Used
Damage Assessment	Damage detection: [300,301,302,303,304,305,306,307,308,309,310,311,312,313]Damage localization: [314,315,316,317,318,319,320]Damage classification: [321,322,323,324]Damage progression prediction: [325]	CNN, GRU, LSTM, BiLSTM, Autoencoder, Transformer, CNN-RNN hybrids
Structural Response Prediction	Strain prediction: [326,327]Displacement/Deflection: [328,329,330,331,332,333,334,335]Seismic and vibration: [336,337,338,339,340]Thermal-induced: [341]Tunnel responses: [342,343,344]Mechanical/Stress: [326,334,340,345,346,347]Cable tension: [348]	LSTM, BiLSTM, CNN, GRU, Attention, Transformer-based models
Structural Load Prediction	Dynamic load: [349,350]	CNN, Bayesian optimization, Autoencoder, CNN-BiLSTM hybrids
Data Reconstruction	Wind: [287,351,352,353]Vibration: [354,355,356]Temperature: [357]Dam monitoring: [358]	CNN, BiLSTM, GAN, Autoencoder, CNN-GRU hybrids, VMD, EMD
Anomaly Detection	Sensor faults: [359,360,361,362,363]Outliers: [363,364,365]	CNN, LSTM, BiLSTM, FCN, Transformer, PCA
Sensor Placement	Optimized placement: [336,349]	Attention-based RNN, CNN-BiLSTM
Data Augmentation/Generation	Synthetic data: [366]	GAN, CycleGAN, CNN, BiLSTM
Other SHM Tasks	Traffic classification: [367]Data compression: [368]Structural state ID: [369]Defect diagnosis: [370]	CNN, Autoencoder, Transformer, CNN-RNN hybrids

8. State of Deep Times Series in Reviewed Literature

We performed a preliminary visual inspection of the literature to identify prevalent modeling trends and challenges.These are presented in the subsequent sections.

8.1. Models

Advanced recurrent neural networks (LSTM and GRU) are among the most commonly used models for sequence-related processing in structural health monitoring (SHM). However, they are often combined with other architectures, such as convolutional neural networks (CNNs), or enhanced with attention mechanisms. To date, the Mamba architecture has not been employed in any of the studies reviewed in this paper.

The reviewed literature was also examined for the application of deep generative models. So far, autoencoders (AEs), generative adversarial networks (GANs), and autoregressive flows have been used. GANs have been applied in data recovery [280,281], as well as data augmentation[238,366].VAEs and AEs have been used in data compression [297,368] and structural response recovery[279]. However, it is evident that normalizing flows, diffusion models, and energy-based models are not reported at all. This observation aligns with Luleci and Catbas [33].

Furthermore, our analysis reveals that most SHM models are deterministic, with only a few studies such as those of Li et al. [298], Fernández et al. [371], Hlaing et al. [372], Pereira and Glisic [373] explicitly addressing uncertainty quantification.

8.2. Challenges

Several studies have highlighted challenges in deep time series modeling for SHM. Common issues include data related challenges, computational complexity,environmental variability, sensor placement and model transferability and generalization.

Data-related challenges prominently feature in SHM studies, especially regarding missing or incomplete data. Numerous researchers have directed efforts toward developing robust data reconstruction and imputation methods for strain data [277,286], acceleration response [275], wind [351,352] and structural temperature data [357], to mention a few. Additionally, limited availability of labeled data and imbalance between damaged and undamaged structural states pose considerable difficulty for supervised learning frameworks, necessitating innovative approaches in few-shot [303] and unsupervised learning [311] paradigms.

The inherent complexity and variability of structural behavior under diverse environmental and operational conditions further complicate modeling. Structural responses often intertwine with environmental factors such as temperature, wind, and humidity, making it challenging to accurately detect and localize damage. Several studies have thus focused on developing methods to distinguish between genuine structural anomalies and normal variations [292,296]. Other studies have made progress at predicting structural responses under nonlinear and extreme conditions [223,334]. Despite the progresses made so far, this area is still an active area of research and might benefit from uncertainty quantification.

Computational limitations also represent a critical constraint. Studies have highlighted the challenge of balancing model complexity with computational efficiency. Studies have highlighted the challenge of balancing model complexity with computational efficiency.Honarjoo et al. [236], Cao et al. [347] provided in-depth computational complexity analysis of models.

Sensor placement emerges as another practical yet challenging dimension. Optimal sensor placement is vital for maximizing monitoring efficiency while minimizing instrumentation costs and complexity. Li et al. [336] took into account the optimization of sensor locations for seismic displacement response of a building.

Issues surrounding model transferability and generalization pose a significant barrier to deploying machine learning models trained on limited experimental or numerical data to real-world structures. Studies have emphasized the necessity of robust domain adaptation and generalization strategies, highlighting methods such as domain adversarial training and physics-informed approaches [289], and out-of-distribution representation learning [321] to overcome these limitations. The capacity to generalize learned models across different structural types, environments, and operational conditions remains a critical research frontier.

9. Conclusions

The initial objective of this review was to provide a mathematically grounded discussion of deep learning models for time series. Reaching this far, we believe such a discussion has been provided. The paper further discussed the application of these models in structural health monitoring of civil structures. Five key thematic were identified, the current status of deep learning was presented, and critical challenges were pinpointed. What is key from the findings of this review is that many studies in SHM rarely quantify uncertainties in their models. The use of domain knowledge to enhance model performance is also limited. These issues might serve as the basis for future research in SHM.