2.1.1 Spline interpolation

The simplest way of constructing a model for a system, be it a real-world system or some complex simulation, is to use a set of $n$ training points and to predict that every unknown position will have the same value as the nearest neighboring training point (see Figure 2(a)). A straightforward extension with slightly higher computational requirements is to connect the training points with straight lines, resulting in a piecewise linear function. Both methods, however, are not continuously differentiable, and for instance less suited for the integration of the model into an optimization process (see Section 4). A more advanced approach to interpolating the training points is to use splines (see Figure 2(b)), which require the piecewise-interpolated functions to agree up to a certain derivative order. For instance, cubic splines are continuously differentiable up to the second derivative. While higher-order spline interpolation works well in 1D cases, it becomes increasingly difficult in multi-dimensional settings. Furthermore, the interpolation approach does not allow for incorporating uncertainty or stochasticity, which is present in real-world measurements. Therefore, it will treat noise components as a part of the system, including, for instance, outliers.

Figure 2 Illustration of standard approaches to making predictive models in machine learning. The data were sampled from the function $y = x\left(1+\sin {x}^2\right)+\epsilon$ with random Gaussian noise, $\epsilon$ , for which $\left\langle {\epsilon}^2\right\rangle = 1$ . The data have been fitted by (a) nearest neighbor interpolation, (b) cubic spline interpolation, (c) linear regression of a third-order polynomial and (d) Gaussian process regression.

2.1.2 Regression

In some specific cases the shortcomings of interpolation approaches can be addressed by using regression models. For instance, simple systems can often be described using a polynomial model, where the coefficients of the polynomial are determined via a least-squares fit (see Figure 2(c)), that is, minimizing the sum of the squares of the difference between the predicted and observed data. The results are generally improved by including more terms in the polynomial but this can lead to overfitting – a situation where the model describes the noise in the data better than the actual trends, and consequently will not generalize well to unseen data. Regression is not restricted to polynomials, but may use all kinds of mathematical models, often motivated by a known underlying physics model. Crucially, any regression model requires the prior definition of a function to be fitted, thus posing constraints on the relationships that can be modeled. In practice this is one of the main problems with this approach, as complex systems can scarcely be described using simple analytical models. Before using these models, it is thus important to verify their validity, for example, using a measure for the goodness of fit such as the correlation coefficient ${R}^2$ , the ${\chi}^2$ test or the mean-squared error (MSE).

2.1.3 Probabilistic models

The field of probabilistic modeling relies on the assumption that the relation between the observed data and the underlying system is stochastic in nature. This means that the observed data are assumed to be drawn from a probability distribution for each set of input parameters to a generative model. Inversely, one can use statistical methods to infer the parameters that best explain the observed data. We will discuss such inference problems in more detail in the context of solving (ill-posed) inverse problems in Section 3.2.

Probabilistic models can generally be divided into two methodologies, frequentist and Bayesian. At the heart of the frequentist approach lies the concept of likelihood, $p\left(y|\theta \right)$ , which is the probability of some outcomes or observations $y$ given the model parameters $\theta$ (see, e.g., Ref. [Reference Anderson and Burnham86] for an introduction). A model is fitted by maximum likelihood estimation (MLE), which means finding the model parameters $\widehat{\theta}$ that maximize the likelihood, $p\left(y|\widehat{\theta}\right)$ . When observations $y = \left({y}_1,{y}_2,\dots \right)$ are independent, the probability of observing $y$ is the product of the probabilities of each observation, $p(y) = {\prod}_{i = 1}^np\left({y}_i\right)$ . As sums are generally easier to handle than products, this is often expressed in terms of the log-likelihood function that is given by the sum of the logarithms of each observation’s probability, that is

(1) $$\begin{align}\log p\left(y|\theta \right) = \sum \limits_{i = 1}^n\log p\left({y}_i|\theta \right).\end{align}$$

Probabilities have values between 0 and 1, so the log-likelihood is always negative and, the logarithm being a strictly monotonic function, minimizing the log-likelihood maximizes the likelihood:

(2) $$\begin{align}{\widehat{\theta}}_{\mathrm{MLE}} = \underset{\theta }{\arg \max}\left\{p\left(y|\theta \right)\right\} = \underset{\theta }{\arg \min}\left\{\log p\left(y|\theta \right)\right\}.\end{align}$$

The optimum can be found using a variety of optimization methods, for example, gradient descent, which are described in more detail in Section 4. A simple example is the use of MLE for parameter estimation in regression problems. In the case in which the error $\left( Ax-y\right)$ is normally distributed, this turns out to be equivalent to the least-squares method we discussed in the previous section.

The MLE is often seen as the simplest and most practical approach to probabilistic modeling. One advantage is that it does not require any a priori assumptions about the model parameters, rather only about the probability distributions of the data. However, this can also be a drawback if useful prior knowledge of the model parameters is available. In this case one would turn to Bayesian statistics. Here one assesses information about the probability that a hypothesis about the parameter values $\theta$ is correct given a set of data $x$ . In this context the probabilistic model is viewed as a collection of conditional probability densities $p\left(\theta |y\right)$ for each set of observed data $y$ , with the aim of finding the posterior distribution  $p\left(\theta |y\right)$ , that is, the probability of some parameters given the data observations. This can be done by applying Bayes’ rule, as follows:

(3) $$\begin{align}p\left(\theta |y\right) = \frac{p\left(y|\theta \right)p\left(\theta \right)}{p(y)},\end{align}$$

where $p\left(y|\theta \right)$ is known from above as the likelihood function and $p\left(\theta \right)$ denotes the prior distribution, that is, our a priori knowledge about the parameters before we observe any data $y$ . The denominator, $p(y) = \int p\left(y|{\theta}^{\prime}\right)\cdot p\left({\theta}^{\prime}\right)\mathrm{d}{\theta}^{\prime }$ , is called the evidence and ensures that both sides of Bayes’ rule are properly normalized by integrating over all possible input parameters $\theta$ . Once the posterior distribution is known, we can maximize it and get the maximum a posteriori (MAP) estimate:

(4) $$\begin{align}{\widehat{\theta}}_{\mathrm{MAP}} = \underset{\theta }{\arg \max}\left\{p\left(\theta |y\right)\right\}.\end{align}$$

As mentioned above, a particular strength of the Bayesian approach is that we can encode a priori information in the prior distribution. Taking the example of polynomial regression, we could, for instance, set a priori distribution $p\left(\theta \right)$ for the regression coefficients $\theta$ that favors small coefficients, thus penalizing high-order polynomials. Another advantage of using the Bayesian framework is that one can quantify the uncertainty in the model result. This is particularly simple to compute in the special case of Gaussian distributions and their generalization, Gaussian processes (GPs), which we are going to discuss in the next section.

2.1.4 Gaussian process regression

A popular version of Bayesian probabilistic modeling is GP regression^[ Reference MacKay^87, Reference Rasmussen^88]. This kind of modeling was pioneered in geostatistics in the context of mining exploration and is historically also referred to as kriging, after the South African engineer Danie G. Krige, who invented the method in the 1950s. Conceptually, one can think of it as loosely related to the spline interpolation method, as it is also locally ‘interpolates’ the training points. Compared to splines and conventional regression methods, kriging has a number of advantages. Being a regression method, kriging can deal with noisy training points, as seen in experimental data. At the same time, the use of GPs involves minimal assumptions and can in principle model any kind of function. Lastly, the ‘interpolation’ is done in a probabilistic way, that is, a probability distribution is assigned to the function values at unknown positions (see Figure 2(d)). This allows for quantifying the uncertainty of the prediction. These features make kriging an attractive method for the construction of surrogate models for complex systems for which only a limited number of the function evaluations is possible, for example, due to the long runtime of the system or the high costs of the function evaluations. GP regression forms the backbone of most implementations of Bayesian optimization (BO), which we will discuss in Section 4.5, including examples for potential use cases.

Mathematically speaking, a GP is an infinite collection of normal random variables, where each finite subset is jointly normally distributed. The mean vector and covariance matrix of the multivariate normal distribution are thereby generalized to the mean function $\mu (x)$ and the covariance function $\sigma \left(x,{x}^{\prime}\right)$ , respectively, where we use the short-hand notation $x = \left({x}_1,{x}_2,\dots \right)$ to denote the function inputs as a vector of real numbers.

A GP can be written in the following form:

(5) $$\begin{align}f(x)\sim \mathcal{GP}\left(\mu (x),\sigma \left(x,{x}^{\prime}\right)\right),\end{align}$$

denoting that the random function $f(x)$ follows a GP with mean function $\mu (x)$ and covariance function $\sigma \left(x,{x}^{\prime}\right)$ .

The mean function $\mu (x)$ is defined as the expectation of the GP, that is, $\mu (x) = \left\langle f(x)\right\rangle$ , whereas the covariance function is defined as $\sigma \left(x,{x}^{\prime}\right) = \left\langle f(x)-\mu (x)\left(f\left({x}^{\prime}\right)-\mu \left({x}^{\prime}\right)\right)\right\rangle$ . Note that in the special case of the constant mean function $\mu (x) = 0$ and a constant diagonal covariance function $\sigma \left(x,{x}^{\prime}\right) = {\sigma}^2\delta \left(x-{x}^{\prime}\right)$ , the GP simply reduces to a set of a normal random variable with zero mean and variance ${\sigma}^2$ commonly referred to as white noise. The covariance function is also referred to as the kernel. Its value at locations $x$ and ${x}^{\prime }$ is proportional to the correlation between the function values $f(x)$ and $f\left({x}^{\prime}\right)$ . A common choice for the covariance function is the squared exponential function, also known as the radial basis function (RBF) kernel:

(6) $$\begin{align}\sigma \left(x,{x}^{\prime}\right) = \exp \left(-\frac{{\left(x-{x}^{\prime}\right)}^2}{2{\mathrm{\ell}}^2}\right),\end{align}$$

where $\mathrm{\ell}$ is the length scale parameter, which controls the rate at which the correlation between $f(x)$ and $f\left({x}^{\prime}\right)$ decays. This kernel hyperparameter can usually be optimized by using the training data to minimize the log marginal likelihood.

It is possible to encode prior knowledge by choosing a specialized kernel that imposes certain restrictions on the model. A variety of such kernels exist. For instance, the periodic kernel (also known as exp-sine-square kernel) given by the following:

(7) $$\begin{align}\sigma \left(x,{x}^{\prime}\right) = \exp \left(-\frac{2\sin^2\left(\pi d\left(x,{x}^{\prime}\right)/\lambda \right)}{{\mathrm{\ell}}^2}\right),\end{align}$$

with the Euclidean distance $d\left(x,{x}^{\prime}\right)$ , length scale $\mathrm{\ell}$ and periodicity $\lambda$ as free hyperparameters, is particularly suitable to model systems that show an oscillatory behavior.

A useful property of kernels is the generation of new kernels through the addition or multiplication of existing kernels^[ Reference Genton and Mach^89]. This property provides another way to leverage prior information about the form of the function to increase the predictive accuracy of the model. For instance, measurement errors incorporated into the model by adding a Gaussian white noise kernel, as for instance done in Figure 2.

Figure 3 visualizes the covariance functions and their influence on the result of GP regression. This includes correlation matrices for the three covariance functions discussed above (white noise, RBF and periodic kernels). Once the mean and the covariance functions are fully defined, we can use training points for regression, that is, fit the GP and kernel parameters to the data and obtain the posterior distribution (see the right-hand panel of Figure 3).

Figure 3 Gaussian process regression: illustration of different covariance functions, prior distributions and (fitted) posterior distributions. Left: correlation matrices between two values $x$ and ${x}^{\prime }$ using different covariance functions (white noise, radial basis function and periodic). Center: samples of the prior distribution defined by the prior mean $\mu (x) = 0$ and the indicated covariance functions. Note that the sampled functions are depicted with increasing transparency for visual clarity. Right: posterior distribution given observation points sampled from $y = {\cos}^2x+\epsilon$ , where $\epsilon$ is random Gaussian noise with ${\sigma}_{\epsilon} = 0.1$ . Note how the variance between observations increases when no noise term is included in the kernel (top row). Within the observation window the fitted kernels show little difference, but outside of it the RBF kernel decays to the mean $\mu = 0$ dependent on the length scale $\mathrm{\ell}$ . This can be avoided if there exists prior knowledge about the data that can be encoded in the covariance function, in this case periodicity, as can be seen in the regression using a periodic kernel.

Figure 4 Sketch of a random forest, an architecture for regression or classification consisting of multiple decision trees, whose individual predictions are combined into an ensemble prediction, for example, via majority voting or averaging.

2.1.5 Decision trees and forests

A decision tree is a general-purpose machine learning method that learns a tree-like model of decision rules from the data to make predictions^[ Reference Breiman, Friedman, Olshen and Stone^90]. It works by splitting the dataset recursively into smaller groups, called nodes. Each node represents a decision point, and the tree branches out from the node according to the decisions that are made. The leaves of the tree represent the final prediction; see Figure 4 for an illustration. This can be either a categorical value in classification tasks (see Section 6) or a numerical value prediction in regression tasks. An advantage of this approach is that it can learn nonlinear relationships in data without having to specify them.

To generate a decision tree one starts at the root of the tree and determines a decision node such that it optimizes an underlying decision metric, such as the MSE in regression settings or entropy and information gain in a classification setting. At each decision point the dataset is split and subsequently the metric is re-evaluated for the resulting groups, generating the next layer of decision nodes. This process is repeated until the leaves are reached. The more decision layers are used, called the depth of the tree, the more complex relationships can be modeled.

Decision trees are easy to implement and can provide accurate predictions even for large datasets. However, with an increasing number of decision layers they may become computationally expensive and may overfit the data, the latter being in particular a problem with noisy data. One method to address overfitting is called pruning, where branches from the tree that do not improve the performance are removed. Another effective method is to use decision-tree-based ensemble algorithms instead of a single decision tree.

Figure 5 Example of gradient boosting with decision trees. Firstly, a decision tree ${g}_1$ is fitted to the data. In the next step, the residual difference between training data and the prediction of this tree is calculated and used to fit a second decision tree ${g}_2$ . This process is repeated $n$ times, with each new tree ${g}_n$ learning to correct only the remaining difference to the training data. Data in this example are sampled from the same function as in Figure 2 and each tree has a maximum depth of two decision layers.

One example of such an algorithm is the random forest, an ensemble algorithm that uses bootstrap aggregating or bagging to fit trees to random subsets of the data and the predictions of individually fitted decision trees are combined by majority vote or average to obtain a more accurate prediction. Another type of ensemble algorithms is boosting, where the trees are trained sequentially and each tries to correct its predecessor. A popular implementation is AdaBoost^[ Reference Freund, Schapire and Jpn^91], where the weights of the samples are changed according to the success of the predictions of the previous trees. Gradient boosting methods^[ Reference Friedman^92, Reference Friedman^93] also use the concept of sequentially adding predictors, but while AdaBoost adjusts weights according to the residuals of each prediction, gradient boosting methods fit new residual trees to the remaining differences at each step (see Figure 5 for an example).

Compared to other machine learning algorithms, a great advantage of a decision tree is its explicitness in data analysis, especially in nonlinear high-dimensional problems. By splitting the dataset into branches, the decision tree naturally reveals the importance of each variable regarding the decision metric. This is in contrast to ‘black-box’ models, such as those created by neural networks, which must be interpreted by post hoc analysis.

2.1.6 Neural networks

Neural networks offer a versatile framework for fitting of arbitrary functions by training of a network of connected nodes or neurons, which are loosely inspired by biological neurons. In the case of a fully connected neural network, historically called a multilayer perceptron, the inputs to one node are the outputs of all of the nodes from the preceding layer. The very first layer is simply all of the inputs for the function to be modeled, and the very last layer is the function outputs. One of the very attractive properties of neural networks is the capacity to have many outputs and inputs, each of which can be multi-dimensional. The weighting of each connection ${w}_i$ and the bias for each node $b$ are the parameters of the network, which must be trained to provide the best approximation of the function to be modeled. Each node gets activated depending on both weights and biases according to a pre-defined activation function. The nonlinear properties of activation functions are widely seen as the key properties to allow neural networks to model arbitrary functions and achieve general learning properties.

One of the simplest, but also most common, activation functions is the rectified linear unit (ReLU), which is defined as follows:

(8) $$\begin{align}\mathrm{ReLU}(x) = \left\{\!\!\begin{array}{ll}x,& \mathrm{if}\;x\ge 0,\\ {}0,& \mathrm{otherwise}.\end{array}\right.\end{align}$$

The function argument $x = \sum {x}_i{w}_i+b$ is the sum of the weights of incoming connections, multiplied with their values ${x}_i$ , and the bias of the node. The ReLU function is easy to compute, and it has the main advantage that it is linear for $x>0$ , which greatly simplifies the gradient calculation in training (see the next paragraph). Drawbacks are that it is not differentiable at $x = 0$ , it is unbounded, and since it yields a constant value below $0$ , it has the potential to produce ‘dead’ neurons. The last issue is solved in the so-called leaky ReLU, which uses a reduced gradient for $x<0$ , giving an output of $\alpha x$ where $0<\alpha <1$ . Other common activation functions include the logistic or sigmoid function, $\kern0.1em \mathrm{sig}\kern0.1em (x) = {(1+{e}^{-x})}^{-1}$ , and the hyperbolic tangent function, $\tanh (x) = ({e}^x-{e}^{-x})/({e}^x+{e}^{-x})$ .

The training is performed iteratively by passing corresponding input–output pairs to the network and comparing the true outputs to those given by the network to calculate the loss function. This is often chosen to be the MSE between the model output and the reference from the training data, but many other types of loss functions are also used (see Section 4.1.1) and the choice of loss function strongly affects the model’s training. The loss is then used to modify the weights and biases via an algorithm known as ‘back-propagation’^[ Reference Wythoff^99]. Here the gradient of the loss function is calculated with respect to the weights and biases via the chain rule. One can then optimize the parameters using, for instance, stochastic gradient descent (SGD) or Adam (adaptive moment estimation)^[ Reference Kingma and Ba^100].

A number of hyperparameters are used to control the training process. Typical ones include the following: the number of epochs – the number of times the model sees each data point; the batch size – the number of data points the model sees before updating its gradients; and the learning rate – a factor determining the magnitude of the update to the gradient. Each factor can have a critical effect on the training of the model.

Figure 6 Simplified sketch of some popular neural network architectures. The simplest possible neural network is the perceptron, which consists of an input, which is fed into the neuron that processes the input based on the weights, an individual bias and its activation function. Multiple such layers can be stacked within so-called hidden layers, resulting in the popular multilayer perceptron (or fully connected network). Besides the direct connection between subsequent layers, there are also special connections common in many modern neural network architectures. Examples are the recurrent connection (which feeds the output of the current layer back into the input of the current layer), the convolutional connection (which replaces the direct connection between two layers by the convolutional operation) and the residual connection (which adds the input to the output of the current layer; note that the above illustration is simplified and the layers should be equal in size).

In the training process, the weights and biases are optimized in order to best approximate the function for converting the inputs into the corresponding outputs. The resulting model will yield an approximation for the unknown function $f(x)$ . However, learning via this method only optimizes how well the model reproduces the training data and, in the worst case, the network essentially ‘memorizes’ the training data, and learns little about the desired function $f(x)$ . To quantify this, a subset of the data is kept separate and used solely for validating the model. Ideally, the loss should be similar on both training and validation data, and if the model has a significantly smaller loss on the training data than on the validation set it is said to overfit, which signifies that the model has learnt the random noise of the training set. The validation loss thus quantifies how well this model generalizes and, with it, how useful it is for predictions on unseen data.

Common methods to combat overfitting include early stopping, that is, terminating the training process once the validation loss stagnates; or the use of dropout layers to randomly switch off some fraction of the network connections for each batch of the training process, thereby preventing the network from learning random noise by reducing its capacity. A further approach is to incorporate variational layers into the network. In these layers, pairs of nodes are used that represent the mean $\mu$ and standard deviation $\sigma$ of a Gaussian distribution. During training, output values are sampled from this distribution and then passed on the subsequent layers, thereby requiring the network to react smoothly to these small random variations to achieve a small training loss. Both dropout layers and variational layers provide regularization that smooths the network response and prevents single nodes from dominating, resulting in improved interpolation performance and better validation loss^[ Reference Kristiadi, Hein and Hennig^101, Reference Wager, Wang and Liang^102] (see Section 3.3 for a brief general explanation of regularization).

Figure 7 Real-world example of a multilayer perceptron for beam parameter prediction. (a) The network layout^[ Reference Kirchen, Jalas, Messner, Winkler, Eichner, Hübner, Hülsenbusch, Jeppe, Parikh, Schnepp and Maier^29] consists of 15 input neurons, two hidden layers with 30 neurons and three output neurons (charge, mean energy and energy spread). The input is derived from parasitic laser diagnostics (laser pulse energy ${E}_{\mathrm{laser}}$ , central wavelength ${\lambda}_0$ and spectral bandwidth $\Delta \lambda$ , longitudinal focus position ${z}_{\mathrm{foc}}$ and Zernike coefficients of the wavefront). Neurons use a nonlinear ReLU activation and 20% of neurons drop out for regularization during training. The (normalized) predictions are compared to the training data to evaluate the accuracy of the model, in this case using the mean absolute ${\mathrm{\ell}}_1$ error as the loss function. In training, the gradient of the loss function is then propagated back through the network to adjust its weights and biases. (b) Measured and predicted median energy ( $\overline{E}$ ) and (c) measured and predicted energy spread ( $\Delta$ E), both for a series of 50 consecutive shots. Sub-figures (b) and (c) are adapted from Ref. [Reference Kirchen, Jalas, Messner, Winkler, Eichner, Hübner, Hülsenbusch, Jeppe, Parikh, Schnepp and Maier29].

Beyond the classical multilayer perceptron network, there exists a plethora of different neural network architectures, as partially illustrated in Figure 6, which are suited to different tasks. In particular, convolutions layers, which extract relevant features from one and two dimensions by learning suitable convolution matrices, are commonly used in the analysis of physical signals and images. Architecture selection and hyperparameter tuning are central challenges in the implementation of neural networks, and are often performed by an additional machine learning algorithm, for example, using BO (see Section 4.5).

The great strength of neural networks is their flexibility and relatively straightforward implementation, with many openly accessible software platforms available to choose from. However, trained networks are effectively ‘black-box’ functions, and do not, in their basic form, incorporate uncertainty quantification. As a result, the networks may make over-confident predictions about unseen data while giving no explanation for those predictions, leading to false conclusions. Various methods exist for incorporating uncertainty quantification into neural networks (see, for example, Ref. [Reference Gawlikowski, Tassi, Ali, Lee, Humt, Feng, Kruspe, Triebel, Jung, Roscher, Shahzad, Yang, Bamler and Zhu103]), such as by including variational layers (discussed above) and training an ensemble of networks on different training data subsets. There are several approaches to try and make neural network models explainable and a review of methods for network interpretation is given, for instance, by Montavon et al. ^[ Reference Montavon, Samek and Müller^104].

Another strength of neural networks is that the performance of a model on a new task can be improved by leveraging knowledge from a pre-trained model on a related task. This so-called transfer learning is typically used in scenarios where it is difficult or expensive to train a model from scratch on a new task, or when there is a limited amount of training data available. For example, transfer learning has been used to successfully train models for image classification and object detection tasks (see Section 6), using pre-trained models that have been trained on large image datasets such as ImageNet^[ Reference Huh, Agrawal and Efros^105]. In the context of laser–plasma physics, transfer learning could be used to improve the performance of a neural network by leveraging knowledge from a pre-trained model that has been trained on data from previous experiments or simulations.

2.1.7 Physics-informed machine learning models

The ultimate application of machine learning for modeling physics systems would arguably be to create an ‘artificial intelligence physicist’, as coined by Wu and Tegmark^[ Reference Wu and Tegmark^113]. One prominent idea at the backbone of how we build physical models is Occam’s razor, which states that given multiple hypotheses, the simplest one that is consistent with the data is to be preferred. In addition to this guiding principle, it is furthermore assumed that a physical model can be described using mathematical equations. A program should therefore be able to automatically create such equations, given experimental data. While a competitive artificial intelligence (AI) physicist is still years away,Footnote ² first steps have been undertaken in this direction. For instance, in 2009 Schmidt and Lipson^[ Reference Schmidt and Lipson^116] presented a genetic algorithm approach (Section 4.4) that independently searched and identified governing mathematical representations such as the Hamiltonian from real-life measurements of some mechanical systems. More recently, research has concentrated more on the concept of Occam’s razor and the underlying idea that the ‘simplest’ representation can be seen as the sparsest in some domain. A key contribution was SINDy (sparse identification of nonlinear dynamics), a framework for discovering sparse nonlinear dynamical systems from data^[ Reference Brunton, Proctor and Kutz^117]. As in CS (Section 3.4), the sparsity constraint was imposed using an ${\mathrm{\ell}}_1$ -norm regularization, which results in the identification of the fewest terms needed to describe the dynamics.

An important step towards combining physics and machine learning was undertaken in physics-informed neural networks (PINNs)^[ Reference Karniadakis, Kevrekidis, Lu, Perdikaris, Wang and Yang^118]. A PINN is essentially a neural network that uses physics equations, which are often described in form of ordinary differential equations (ODEs) or partial differential equations (PDEs), as a regularization to select a solution that is in agreement with physics. This is achieved by defining a custom loss function that includes a discrepancy term, the residuals of the underlying ODEs or PDEs, in addition to the usual data-based loss components. As such, solutions that obey the selected physics are enforced. In contrast to SINDy, there is no sparsity constraint imposed on the network weights, meaning that the network could still be quite complex. Early examples of PINNs were published by Raissi et al. ^[ Reference Raissi, Perdikaris and Karniadakis^119– Reference Raissi, Yazdani and Karniadakis^122] and Long et al. ^[ Reference Long, Lu, Ma and Dong^123] in 2017–2020. Since then, the architecture has been applied to a wide range of problems in the natural sciences, with a quasi-exponential growth in publications^[ Reference Cuomo, Di Cola, Giampaolo, Rozza, Raissi and Piccialli^124]. Applications include, for instance, (low-temperature) plasma physics, where PINNs have been successfully used to solve the Boltzmann equation^[ Reference Kawaguchi and Murakami^125], and quantum physics, where PINNs were used to solve the Schrödinger equation of a quantum harmonic oscillator^[ Reference Stiller, Bethke, Böhme, Pausch, Torge, Debus, Vorberger, Bussmann and Hoffmann^126]. The work by Stiller et al. ^[ Reference Stiller, Bethke, Böhme, Pausch, Torge, Debus, Vorberger, Bussmann and Hoffmann^126] uses a scalable neural solver that could possibly also be extended to solve, for example, the Vlasov–Maxwell system governing laser–plasma interaction.

2.2.1 Classical models

To model a time series one usually starts with a set of assumptions regarding its structure, that is, the interdependence between values at different times. A simple, approximate assumption is that the observed values in a discrete time series are linearly related. An important, wide-spread class of such models is the so-called autoregressive models, which assert that the next value ${y}_t$ in a time series is given as a linear function of the $p$ prior values ${y}_{t-1},{y}_{t-2},\dots, {y}_{t-p}$ . In addition, each value is assumed to be corrupted by additive white noise ${\varepsilon}_t\sim \mathcal{N}\left(0,{\sigma}^2\right)$ , representing, for example, measurement errors or inherent statistical fluctuations of the underlying process, which endows the models with stochasticity. In its most common form, the autoregressive model of order $p$ , denoted by $\mathrm{AR}(p)$ , is defined via the recurrence relation:

(9) $$\begin{align}{y}_t = \sum \limits_{i = 1}^p{\varphi}_i{y}_{t-i}+{\varepsilon}_t,\end{align}$$

where ${\varphi}_1,\dots, {\varphi}_p$ are the model’s parameters to be estimated. In this form, the problem of finding the model’s parameters is a linear regression problem that can be solved via the least-squares method:

(10) $$\begin{align}\left\{{\widehat{\varphi}}_1,\dots, {\widehat{\varphi}}_p\right\} = \mathop{\mathrm{argmin}}_{{\varphi_1,\dots, {\varphi}_p}}{\left({y}_t-\sum \limits_{i = 1}^p{\varphi}_i{y}_{t-i}\right)}^2.\end{align}$$

In another approach we might assume that the next value ${y}_t$ is instead given by a linear combination of (external) statistical fluctuations, sometimes called shocks, from the past $q$ points in time, again encoded in white noise terms ${\varepsilon}_t\sim \mathcal{N}\left(0,{\sigma}^2\right)$ . The corresponding model, called the moving average model and denoted by $\mathrm{MA}(q)$ , is defined as follows:

(11) $$\begin{align}{y}_t = \mu +\sum \limits_{i = 1}^q{\vartheta}_i{\varepsilon}_{t-i}+{\varepsilon}_t,\end{align}$$

where $\mu$ is the mean of the time series.

The parameters of the moving average process cannot be inferred by linear methods such as least squares, but rather have to be estimated by means of maximum likelihood methods.

Contrary to the $\mathrm{AR}(p)$ model, which is only stationaryFootnote ³ for certain parameters ${\varphi}_1,\dots, {\varphi}_p,$ the $\mathrm{MA}(q)$ model is always (strictly) stationary per definition. Put loosely, stationarity can be understood as the stochastic properties (mean, variance, …) of the time series being constant over time.

Another distinction between the autoregressive and the moving average model is how far into the future the effects of statistical fluctuations (shocks) are propagated. In the moving average model the white noise terms are only propagated $q$ steps into the future, while in the autoregressive model the effects are propagated infinitely far into the future. This is because the white noise terms are part of the prior values, which themselves are part of the future values in Equation (9).

If the time series in question cannot be explained by the $\mathrm{AR}(p)$ or the $\mathrm{MA}(q)$ model alone, both models can be combined into what is called an autoregressive-moving-average model, denoted by $\mathrm{ARMA}\left(p,q\right)$ :

(12) $$\begin{align}{y}_t = \mu +\sum \limits_{i = 1}^p{\varphi}_i{y}_{t-i}+\sum \limits_{j = 1}^q{\vartheta}_j{\varepsilon}_{t-j}+{\varepsilon}_t.\end{align}$$

Note that for the special cases of $p = 0$ and $q = 0$ , the $\mathrm{ARMA}\left(p,q\right)$ reduces to the $\mathrm{MA}(q)$ and $\mathrm{AR}(p)$ , respectively.

When fitting autoregressive models, special care should be taken that the time series to be modeled is stationary before fitting, otherwise spurious correlations are introduced. Spurious correlations are apparent correlations between two or more time series that are not causally related, thereby potentially leading to fallacious conclusions as warned of in the well-known adage ‘correlation does not imply causation’.

If the time series $y$ is non-stationary in general but stationary with respect to its mean, that is, the variations relative to the mean value are stationary, it might be possible to transform it into a stationary time series. For this we introduce a new series ${z}_t = {\nabla}^d{y}_t$ by differencing the original series $d$ times, where we defined the differencing operator $\nabla {y}_t\equiv {y}_t-{y}_{t-1}$ . Applying the differencing operator once ( $d = 1$ ), for example, removes a linear trend from ${y}_t$ , applying it twice removes a quadratic trend and so forth. If $z$ is stationary after differencing $y$   $d$ -times it can be readily modeled by Equation (12), which leads us to the autoregres sive-integrated-moving-average model $\mathrm{ARIMA}\left(p,d,q\right)$ :

(13) $$\begin{align}{z}_t = \sum \limits_{i = 1}^p{\varphi}_i{z}_{t-i}+\sum \limits_{j = 1}^q{\vartheta}_j{\varepsilon}_{t-j}+{\varepsilon}_t.\end{align}$$

Note that in contrast to the $\mathrm{ARMA}\left(p,q\right)$ model in Equation (12), the time series ${z}_t$ appearing here is a $d$ -times differenced version of the original series ${y}_t$ . Further extensions that will only be noted here for completeness are seasonal ARIMA models, called SARIMA models, that allow the modeling of time series that exhibit seasonality (periodicity), and exogeneous ARIMA models, called ARIMAX models, that allow the modeling of time series that are influenced by a separate, external time series. Both extensions can be combined to yield the so-called SARIMAX model, which is general enough to cover a large class of time series problems. However, we are still in any case limited to problems in which the time series is generated by a linear stochastic process. To cover more general nonlinear problems, we introduce SSMs.

2.2.2 State-space models

SSMs offer a very general framework to model time series data^[ Reference Commandeur and Koopman^127]. In this framework, presuming that the time series is based on some underlying system, it is assumed that there exists a certain true state of the system ${x}_t$ of which we observe a value ${y}_t$ subject to measurement noise. The true state ${x}_t$ , usually inaccessible and hidden to us, and the observed state ${y}_t$ are modeled by the state equation and the observation equation, respectively. A prominent example of SSMs in machine learning is hidden Markov models (HMMs)^[ Reference Zucchini and MacDonald^128], in which the hidden state ${x}_t$ is modeled as a Markov process. In general there are no restrictions on the functional form of the state and observation equations; however, the most common type of SSM is the linear Gaussian SSM, also referred to as the dynamic linear model, in which the state and observation equations are modeled as linear equations and the noise is assumed to be Gaussian. The prototypical example of such linear Gaussian SSMs is the Kalman filter ^[ Reference Kalman^129], ubiquitous in control theory, robotics and many other fields^[ Reference Auger, Hilairet, Guerrero, Monmasson, Orlowska-Kowalska and Katsura^130], including for instance adaptive optics^[ Reference Poyneer and Véran^131] and particle accelerator control^[ Reference Ushakov, Echevarria and Neumann^132]. The term ‘filter’ originates from the fact that it filters out noise to estimate the underlying state of a system. The state equations can be written as follows:

(14) $$\begin{align}{x}_t &= {A}_t{x}_{t-1}+{B}_t{u}_t+{C}_t{\varepsilon}_t,\nonumber\\ {}{y}_t &= {D}_t{x}_t+{E}_t{\eta}_t,\end{align}$$

where the system noise ${\varepsilon}_t$ and observation noise ${\eta}_t$ are sampled from a normal distribution ${\varepsilon}_t,{\eta}_t\sim \mathcal{N}\left(0,1\right)$ . Note that the separation of these two noise terms differs from the so-far discussed models. This permits the modeling of systems that are driven by noise, while also accounting for observation noise of possibly differing amplitude. Another difference is the addition of a second process ${u}_t$ , commonly called the control input in signal processing, which is a deterministic external process driving the system state. This permits the modeling of systems with known external inputs (e.g., a controller-driven motor) or with known deterministic drivers (e.g., system temperature). Note that we can recover the autoregressive models discussed earlier as a special case of the Kalman filter. For example, the $\mathrm{AR}(1)$ process is obtained by setting ${A}_t = \mathrm{const}. = {\varphi}_1$ , ${B}_t = 0$ , ${C}_t = \mathrm{const}. = {\sigma}^2$ , ${D}_t = 1$ and ${E}_t = 0$ . For a more detailed discussion on Kalman filters we refer the interested reader to the literature^[ Reference Welch and Bishop^133]. Other (nonlinear) ways to perform inference in SSMs exist, for instance using sequential Monte Carlo estimation^[ Reference Ristic, Arulampalam and Gordon^134].

2.2.3 Forecasting networks

While predictions based on autoregression or SSMs can suffice for many applications, the nonlinear nature of neural networks can be harnessed to model time series with complex, nonlinear dependencies^[ Reference Fawaz, Forestier, Weber, Idoumghar and Muller^135, Reference Lim and Zohren^136]. A particularly relevant architecture is the recurrent neural networks (RNNs). These are similar to the ‘traditional’ neural networks we discussed in Section 2.1.6, but they have a ‘memory’ that allows them to retain information from previous inputs. This is implemented in a somewhat similar way to the linear recurrence relations discussed in the previous sections, with a key difference being that the state ${x}_t$ of RNNs is updated according to an arbitrary nonlinear function. The current state ${x}_t$ is calculated from the previous states ${x}_{t-1},{x}_{t-2},\dots$ through the recurrence equation:

(15) $$\begin{align}{x}_t = f\left({w}_{\mathrm{rec}}{x}_{t-1}+{b}_{\mathrm{rec}}\right),\end{align}$$

where $f$ is a nonlinear activation function, ${w}_{\mathrm{rec}}$ is a matrix of recurrent weights and ${b}_{\mathrm{rec}}$ is a bias vector. This equation allows the update to each element of the current state vector, ${x}_t$ , to be dependent on the whole of the previous state vector, ${x}_{t-1}$ . The output ${y}_t$ of the network is calculated from the current state ${x}_t$ through the output equation:

(16) $$\begin{align}{y}_t = f\left({w}_{\mathrm{out}}{x}_t+{b}_{\mathrm{out}}\right),\end{align}$$

where ${w}_{\mathrm{out}}$ contains the output weights and ${b}_{\mathrm{out}}$ is another bias vector. Thus, the entire network state is updated from the previous state through a recurrent equation, and the current output is calculated from the current state through an output equation. The network state therefore contains all previous information about the time series. However, as the network state is updated by a simple multiplication of the weights ${w}_{\mathrm{rec}}$ with the previous state ${x}_{t-1}$ , the so-called vanishing gradient problem may occur^[ Reference Hochreiter^137, Reference Hochreiter^138]. This problem is solved in a seminal work by Hochreiter and Schmidhuber^[ Reference Hochreiter and Schmidhuber^139], who introduced a special type of RNN, the long short-term memory (LSTM) network. In contrast to the simple update equation of RNNs, LSTMs use a special type of memory cell that can learn long-term dependencies. This includes a so-called forget gate  ${f}_t$ to update the previous state ${x}_{t-1}$ to the current state ${x}_t$ :

(17) $$\begin{align}{x}_t = {f}_t\odot {x}_{t-1}+{i}_t\odot {h}_t,\end{align}$$

where $\odot$ denotes the element-wise Hadamard product, ${\left[A\odot B\right]}_{ij} = {A}_{ij}{B}_{ij}$ . In addition, the LSTM uses two more gates, the input gate ${i}_t$ and output gate ${o}_t$ , where the latter is used to calculate the hidden state ${h}_t$ from the previous hidden state via ${h}_t = {o}_t\odot {h}_{t-1}$ . The three gates ${f}_t$ , ${i}_t$ and ${o}_t$ are calculated from the current input ${x}_t$ , the previous hidden state ${h}_{t-1}$ and the previous gate states ${f}_{t-1}$ , ${i}_{t-1}$ and ${o}_{t-1}$ using individually set weights and biases. The gates determine how much of the previous state ${x}_{t-1}$ and the current hidden state ${h}_t$ are used to calculate the current state ${x}_t$ . The output of the LSTM network is then calculated from the current state ${x}_t$ through the same output equation, Equation (16), as used in the RNN. Despite it being introduced in the early 1990s, the LSTM architecture remains one of the most popular network architectures for predictive tasks to date.

A more recently introduced type of neural network that can learn to interpret and generate sequences of data is the transformer. Transformers are similar to RNNs, but instead of processing the time series in a sequential manner, they use a so-called attention mechanism^[ Reference Vaswani, Shazeer, Parmar, Uszkoreit, Jones, Gomez, Kaiser and Polosukhin^140] to capture dependencies between all elements of the time series simultaneously and, thus, focus on specific parts of an input sequence. Assuming again a time series ${x}_t,{x}_{t-1},\dots$ , a transformer maps each point ${x}_t$ to a representation ${h}_t$ using a linear layer, for example, ${h}_t = {wx}_t+b$ . The transformer network then calculates a new representation for each point ${x}_t$ using the attention mechanism:

(18) $$\begin{align}{\tilde{h}}_t = \sum {\alpha}_{ti}{h}_i, \end{align}$$

where the attention weights ${\alpha}_{ti}$ are calculated using the point representations ${h}_i$ and the previous representation ${\tilde{h}}_t$ . The attention weights ${\alpha}_{ti}$ are used to calculate the new representation ${\tilde{h}}_t$ of each point ${x}_t$ from all previous representations ${\tilde{h}}_i$ , $i<t$ , of the previous point ${x}_i$ , $i<t$ . The new representations are then used to calculate the output ${y}_t$ of the network as in Equation (16). Thus, the attention mechanism of a transformer network can be interpreted as a nonlinear function that updates the representation of each point ${x}_t$ from all previous representations of the previous points ${x}_i$ , $i<t$ . It can be applied multiple times and enables transformers to learn complex patterns in data, outperforming RNNs on a variety of tasks, such as machine translation and language modeling. It has also been shown that transformers are more efficient than RNNs, meaning they can be trained on larger datasets in less time. One recent example of a transformer developed for time series prediction is the Temporal Fusion Transformer^[ Reference Lim, Ark, Loeff and Pfister^141]. Beyond RNNs, LSTMs and transformers, there exist many other network architectures that can be used for forecasting, for example, fully convolutional networks (FCNs), as discussed by Wang et al. ^[ Reference Wang, Yan and Oates^142].

For longer-term forecasting, predictive networks are usually employed iteratively, meaning that a single-step forecast uses only real data, while a two-step forecast will use the historical data and the most recent prediction value, and so forth. In the context of laser–plasma physics and experiments, predictive neural networks could be used to model the time series data from diagnostic measurements in order to make predictions about the future performance, for example, for predictive steering of laser or particle beams.

4.1.1 Objective functions

Most optimization algorithms are based on maximizing or minimizing the value of the objective function, which has to be constructed in a way that it represents the actual, user-defined objective of the optimization problem. Typically, the objective function produces a single value, where higher (or lower) values represent a more optimized state.Footnote ⁸ The optimization problem is then a case of finding the parameters that maximize (or minimize) the objective function.

When the objective is to construct a model, the objective function encodes some measure of similarity between the model and what it is supposed to represent, for example, a measure of how well the model fits some training data. For example, deep learning algorithms minimize an objective function that encodes the difference between desired output values for a given input and actual results produced by the algorithm. A common, basic metric for this is the MSE cost function, in which the difference between predicted and actual values is squared. We already mentioned this metric at several points of this review, for example, in Section 3.1. The MSE objective function belongs to a class of distance metrics, the $\mathrm{\ell}$ -norms, which are defined by the following:

(24) $$\begin{align}{\mathrm{\ell}}_p\left(x,y\right) = {\left(\sum \parallel x-y{\parallel}^p\right)}^{1/p},\end{align}$$

where $x$ and $y$ are vectors, and we use the notation ${{\mathrm{\ell}}_2 = \mathrm{MSE}}$ . One can also use other $\mathrm{\ell}$ -norms, for example, ${\mathrm{\ell}}_1$ , which is more robust to outliers than ${\mathrm{\ell}}_2$ since it penalizes large deviations less severely. Note that the number of non-zero values is sometimes referred to as the ‘ ${\mathrm{\ell}}_0$ ’ norm (see Section 3.4), even though it does not follow from Equation (24).

Another popular similarity measure is the Kullback–Leibler (KL) divergence, sometimes called relative entropy, which is used to find the closest probability distribution for a given model. It is defined as follows:

(25) $$\begin{align}\mathrm{KL}\left(p|q\right) = \sum p(x)\log \frac{p(x)}{q(x)},\end{align}$$

where $p(x)$ is the probability of observing the value $x$ according to the model, and $q(x)$ is the actual probability of observing $x$ . The KL divergence is minimized when the model prediction $p(x)$ is as close to the actual distribution $q(x)$ as possible. It is a relative measure, that is, it is only defined for pairs of probability distributions. The KL divergence is closely related to the cross-entropy cost function:

(26) $$\begin{align}\mathrm{CE}\left(p,q\right) = -\sum p(x)\log q(x) = \mathrm{KL}\left(p|q\right)-H(p),\end{align}$$

where $p(x)$ and $q(x)$ are probability distributions, and $H(p) = -\sum p(x)\log p(x)$ is the Shannon entropy of the distribution $p(x)$ .

Other objective functions may rely on the maximization or minimization of certain parameters, for example, the beam energy or energy-bandwidth of particle beams produced by a laser–plasma accelerator. Both of these are examples of what is sometimes referred to as summary statistics, as they condense information from more complex distributions, in this case the electron energy distribution. While simple at the first glance, these objectives need to be properly defined and there are often different ways to do so^[ Reference Irshad, Karsch and Döpp^198]. In the example above, energy and bandwidth are examples for the central tendency and the statistical dispersion of the energy distribution, respectively. These can be measured using different metrics, such as the weighted arithmetic or truncated mean, the median, mode, percentiles, for the former; and full width at half maximum, median absolute deviation, standard deviation, maximum deviation, for the latter. Each of these seemingly similar measures emphasizes different features of the distribution they are calculated from, which can affect the outcome of optimization tasks. Sometimes one might also want to include higher-order momenta as objectives, such as the skewness, or use integrals, for example, the total beam charge.

4.1.2 Pareto optimization

In practice, optimization problems often constitute multiple, sometimes competing objectives ${g}_i$ . As the objective function should only yield a single scalar value, one has to condense these objectives in a process known as scalarization. Scalarization can, for instance, take the form of a weighted product $g = \prod {g}_i^{\alpha_i}$ or sum $g = \sum {\alpha}_i{g}_i$ of the individual objectives ${g}_i$ with the hyperparameters ${\alpha}_i$ describing the weight. Another common scalarization technique is $\epsilon$ -constraint scalarization, where one seeks to reformulate the problem of optimizing multiple objectives into a problem of single-objective optimization conditioned on constraints. In this method the goal is to optimize one of the ${g}_i$ given some bounds on the other objectives. All of these techniques introduce some explicit bias in the optimization, which may not necessarily represent the desired outcome. Because of this, the hyperparameters of the scalarization may have to be optimized themselves by running optimizations several times.

Figure 12 Pareto front. Illustration of how a multi-objective function $f(x) = y$ acts on a 2D input space $x = \left({x}_1,{x}_2\right)$ and transforms it to an objective space $y = \left({y}_1,{y}_2\right)$ on the right. The entirety of possible input positions is uniquely color-coded on the left and the resulting position in the objective space is shown in the same color on the right. The Pareto-optimal solutions form the Pareto front, indicated on the right, whereas the corresponding set of coordinates in the input space is called the Pareto set. Note that both the Pareto front and Pareto set may be continuously defined locally, but can also contain discontinuities when local maxima become involved. Adapted from Ref. [Reference Irshad, Karsch and Döpp199].

A more general approach is Pareto optimization, where the entire vector of individual objectives $g = \left({g}_1,\dots, {g}_N\right)$ is optimized. To do so, instead of optimizing individual objectives, it is based on the concept of dominance. A solution is said to dominate other solutions if it is both at least as good on all objectives and strictly better than the other solutions on at least one objective. Pareto optimization uses implicit scalarization by building a set of non-dominated solutions, called the Pareto front, and maximizing the diversity of solutions within this set. The latter can, for instance, be quantified by the hypervolume of the set or the spread of solutions along each individual objective. As it works on the solution for all objectives at once, Pareto optimization is commonly referred to as multi-objective optimization. An illustration of both the Pareto front and set is shown in Figure 12, where a multi-objective function $f$ ‘morphs’ the input space into the objective space. In this example $f$ is a modified version of the Branin–Currin function^[ Reference Dixon and Szego^200, Reference Currin, Mitchell, Morris, Ylvisaker and Amer^201], exhibiting a single, global maximum in ${y}_2$ but multiple local maxima in  ${y}_1$ . The individual 2D outputs ${y}_1 = {f}_1\left({x}_1,{x}_2\right)$ , with ${f}_1$ being the Branin function, and ${y}_2 = {f}_2\left({x}_1,{x}_2\right)$ , with ${f}_2$ being the Currin function, are depicted with red and blue colormaps at the bottom.

6.1.1 Support vector machines

Support vector machines (SVMs) are a popular set of machine learning methods used primarily in classification and recognition. For a simple binary classification task, a SVM draws a hyperplane that divides all data points of one category from all data points of the other category. While there could be many hyperplanes, an SVM looks for the optimal hyperplane that best separates the dataset by maximizing the margin between the hyperplane and the data points in both categories. The points that locate right on the margin are called ‘supporting vectors’. For a dataset with more than two categories, classification is performed by dividing the multi-classification problem into several binary classification problems and then finding a hyperplane for each. In practice, the data points in two categories can mix together so that they cannot be clearly divided by a linear hyperplane. For such nonlinear classification tasks, the kernel trick is used to compute the nonlinear features of the data points and map them to a high-dimensional space, so that a hyperplane can be found to separate the high-dimensional space. The hyperplane will then be mapped back to the original space.

The ideal application scenario for an SVM is for datasets with small samples but high dimensions. The choice of various kernel functions also adds to the flexibility of this method. However, an SVM would be very computationally expensive for large datasets. Besides, its accuracy can significantly decrease when analyzing datasets with large noise levels, as the hyperplanes cannot be defined precisely. Therefore, especially in the context of high-power laser experiments, one has to be cautious when applying SVMs if there are considerable stability issues.

6.1.2 Convolutional neural networks

CNNs are a type of neural network (cf. Section 2.1.6) that is particularly well-suited for image classification^[ Reference Rawat and Wang^253], but are also used in various other problems.

Such a network is composed of sequential convolutional layers, in each of which an $N\times N$ kernel (or ‘filter’ matrix) is convolved with the output of the previous layer. This operation is done by sliding the kernel over the input image, and each pixel in the output layer is calculated by the dot product of the kernel with a sub-section of the input image centered around the corresponding pixel. Resultingly, convolutional layers are capable of detecting local patterns in the $N\times N$ region of the kernel. The kernel is parameterized with weights, which are learned via back-propagation as in a regular neural network. Within a layer, there can also be multiple channels of the output, practically thought of as multiple kernels being passed over the image, allowing for different features to be detected. The early layers of a CNN detect simple structures such as edges, but by adding multiple layers with varying kernel sizes, the network can perform high-level abstraction in order to detect complicated patterns.

The convolutional layer makes the CNN very efficient for image classification. It allows the network to learn translation-invariant features – a feature learned at a certain position of an image can be recognized at any position on the same image.

In order to detect patterns that are non-local in the image, pooling is often applied. There exist many schemes of pooling, but the general concept is to take a set of pixels in the input and to apply some operation that turns them into 1 pixel. Examples include max-pooling (taking the maximum of the set of pixels) and average pooling. This operation decreases the dimensions of the image, and therefore allows a subsequent convolutional layer with the same $N\times N$ kernel to detect features that were much further apart in the original image. Typical CNNs will use multiple pooling layers to decrease the dimensions until a kernel can nearly span the whole image to detect any non-local patterns. The output is then flattened and several fully connected layers can be used to manipulate the data for the relevant (i.e., regression or classification) task.

While the use of deeper CNNs with an increasing number of layers tends to improve performance^[ Reference Szegedy, Liu, Jia, Sermanet, Reed, Anguelov, Erhan, Vanhoucke and Rabinovich^254], architectures can suffer from unstable gradients in training via back-propagation^[ Reference Balduzzi, Frean, Leary, Lewis, Ma and McWilliams^255], showing that some deeper architectures are not as easy to optimize. One particularly influential solution to this problem was the introduction of the residual shortcut, where the input to a block is added to the output of the block. In back-propagation, this enables the ‘skipping’ of layers, effectively simplifying the network and leading to better convergence. This was first proposed and implemented in ResNet ^[ Reference He, Zhang, Ren and Sun^256], which has since become a standard for deep CNN architectures^[ Reference Szegedy, Ioffe, Vanhoucke and Alemi^257], with any number of variations. It should be mentioned though that a number of competitive networks without residual shortcuts exist, for example, AlexNet ^[ Reference Krizhevsky, Sutskever and Hinton^258].