2.1. DDMD and sparse regression

In general, sparse (symbolic) regression is a type of machine learning algorithm which searches a predefined feature space in order to output mathematical expressions which best describe the input dataset. The objective is to use sparse regression to infer the underlying partial differential equations (PDEs) from the data. The inference of a PDE is done by regression on a predefined line up of candidate PDE terms – the feature space. The method used for this work is suited to be applied to time-series data with spatio-temporal variation, in particular data sets containing the evolution of the system dynamics in the nonlinear regime. The feature space $\varTheta$ is defined as a matrix in which each column represents a scaled feature: a numerical estimate of a variable, its derivatives, linear combinations of variables and derivatives or polynomial combinations of terms up to a specified order, over a data sample over time. The algorithm discovers both the model structure and the model parameters. In practice prior knowledge and physical intuition about the systems dynamics informs the design of candidate terms – features – included in the feature space. However, in theory, prior knowledge about the system dynamics is not used. This makes the approach potentially suitable for both model validation as well as probing for new underlying physics.

The balance between the discovered model accuracy and complexity is achieved via Pareto analysis that prevents overfitting. In practice this leads to the optimization favouring models with less selected features (equation terms) over the ones with more features at the same accuracy. The cost function being minimized is of the following general form:

(2.1)\begin{equation} ||\partial_t \langle F \rangle - \langle \varTheta \rangle \xi||^2_2 + \lambda ||\xi||_0, \end{equation}

where the variable $F$ is a quantity whose dynamics is governed by a PDE that we seek to discover. The terms of this equation are the features in $\varTheta$ selected by multiplication with the sparsest vector $\xi$ which contains the model parameters (coefficients in front of the corresponding terms in the PDE). The vector elements of $\xi$ are computed using a thresholded least-squares algorithm (Brunton et al. Reference Brunton, Proctor and Kutz2016), while its level of sparsity is obtained by means of a 30-fold cross-validation procedure. Vector $\lambda$ contains $100$ regularization parameters used in sequentially thresholded least-squares optimization of $\xi$. The values of the regularization parameters are logarithmically increased from $10^{-6}$ up to $1$ to obtain increasingly sparser solutions of $\xi$. Each candidate feature in $\varTheta$ and the time derivatives vector $\partial _t F$ are approximated by central finite differencing and integrated over volumes on the sampled data. The angle brackets denote averaging over the sampled volumes. Therefore the governing PDE is discovered in its integral form, this methodology has been described in detail in Alves & Fiuza (Reference Alves and Fiuza2020) where it is applied to a hierarchy of basic plasma physics models.

2.2. The HW model

In our particular case the input data are generated through direct numerical simulation of resistive drift-wave turbulence according to two-dimensional (2-D) HW and mHW models (Hasegawa & Mima Reference Hasegawa and Mima1978; Wakatani & Hasegawa Reference Wakatani and Hasegawa1984; Hasegawa & Wakatani Reference Hasegawa and Wakatani1987) in the $(x,y)$ plane perpendicular to the magnetic field along the $z$ axis. These models are a good test bed for our model discovery method for two main reasons. Firstly, these are relevant models providing a description of drift-wave turbulence near the tokamak plasma edge. Secondly, they are simple consisting of two coupled PDEs with a relatively small number of terms and therefore allow straightforward interpretation and checking of obtained results. The HW couples the continuity and vorticity equations in reduced MHD (Wakatani & Hasegawa Reference Wakatani and Hasegawa1984; Hasegawa & Wakatani Reference Hasegawa and Wakatani1987; Biskamp & Zeiler Reference Biskamp and Zeiler1995)

(2.2)\begin{gather} \frac{\partial n}{\partial t}={-}[\phi,n]+\alpha (\phi - n)-\kappa \frac{\partial \phi}{\partial y}-D_n \nabla^{4}n, \end{gather}

(2.3)\begin{gather}\frac{\partial \omega}{\partial t}={-}[\phi,\omega]+\alpha (\phi - n)-D_{\omega}\nabla^{4}\omega, \end{gather}

where $n$ and $\omega$ represent the electron density and vorticity,Footnote ² respectively, $\phi$ stands for the electrostatic potential and the constants $\alpha$ and $\kappa$ are the electron adiabaticity and the background density gradient driver. The Poisson bracket for two functions of phase space and time, f and g, is defined in the usual manner $[f,g]=\hat {z} \times \boldsymbol {\nabla } f \boldsymbol {\cdot } \boldsymbol {\nabla } g$. Similarly, mHW couples the continuity and vorticity equations but captures the electron adiabatic response with the substitutions $n \rightarrow \tilde {n}$ and $\omega \rightarrow \tilde {\omega }$ where the new variables are defined by

(2.4)\begin{equation} \tilde{f} =f-\frac{1}{L_y}\int f\,{\rm d}y \end{equation}

and the equations being numerically solved become

(2.5)\begin{gather} \frac{\partial n}{\partial t}={-}[\phi,n]+\alpha (\tilde{\phi} - \tilde{n})-\kappa \frac{\partial \phi}{\partial y}-D_n \nabla^{4}n, \end{gather}

(2.6)\begin{gather}\frac{\partial \omega}{\partial t}={-}[\phi,\omega]+\alpha (\tilde{\phi} - \tilde{n})-D_{\omega}\nabla^{4}\omega. \end{gather}

Inclusion of the electron adiabatic response enables (for certain parameter values) the simulation of elongated asymmetric vortex modes, so-called, zonal flows.Footnote ³ The introduced substitutions as defined by (2.4) relaxes the coupling between the electric potential $\phi$ and the density $n$ for these modes as they are not subject to the coupling arising from the parallel momentum equation. Zonal flows are driven by nonlinear energy transfer from drift waves and represent a self-regulation agent for drift-wave transport and turbulence (Diamond et al. Reference Diamond, Itoh, Itoh and Hahm2005). Zonal flows play an important role in fusion plasmas because they suppress plasma turbulence and thus reduce the transport driven by it in turn improving plasma confinement.

It should be noted here that the standard drift-wave normalizations (Wakatani & Hasegawa Reference Wakatani and Hasegawa1984; Biskamp & Zeiler Reference Biskamp and Zeiler1995) were used: $\phi \longrightarrow (e\phi /T_e)L_n/\rho _s$,$n \longrightarrow (n/n_0)L_n/\rho _s$, $t \longrightarrow tc_s/L_n$, $\boldsymbol {\nabla }_{\perp } \longrightarrow \rho _s \boldsymbol {\nabla }_{\perp }$ and $\boldsymbol {\nabla }_{\parallel }\longrightarrow L_{\parallel }\boldsymbol {\nabla }_{\parallel }$, where $\rho _s$ is the ion Larmor radius, $n_0$ the background electron density, e the electron charge, $c_s$ the ion acoustic velocity, $\nu _{ei}$ the electron ion collision rate and the $L_n$ the scale length of the density background. The parallel scale length is defined by $L_{\parallel }=(L_n T_e/m_e c_s \nu _{ei})^{1/2}$.

2.3. Simulation set-up

The simulation outputs the changes in plasma density, vorticity and potential in the ($x,y$) plane. The size of the simulated domain was $16{\rm \pi} \times 16{\rm \pi}$. The temporal resolution was varied while keeping the total simulated time $t$ constant and sufficient to allow for nonlinear effects to develop and evolve. Once the input data sets are available what follows is setting up the data sampling routine. In our particular case the data were sampled with varying sampling density both spatially (distribution and number of sampling volumes) and temporally (time-series length). Uniform sampling was applied only in space while the time sampling was non-uniform/random within the specified time interval. This time sampling allows enables capturing both fast and long-time scale features. The sampling is a routine which may require adjustment tailored to the input data structure and the expected underlying physics. For some problems the uniform sampling might not be appropriate and one might choose to use threshold-based sampling or targeted sampling restricted to predefined regions of the phase space.

The features (the potential terms in the sought after PDE) need to be chosen and evaluated on the sampled data points. The choice of features to include in the $\varTheta$ matrix relies on the physical intuition about the problem at hand. Once the feature matrix is constructed, sequentially thresholded least-squares regression is applied to the sampled data in order to obtain the sparsity vector $\xi$ which reveals the retained features and the numerical values of their corresponding coefficients (for code structure diagram see figure 1). Plotting the fraction of variance unexplained (FVU) as a function of the number of retained features reveals the discovered PDE. The FVU is in this work defined as the ratio of the model's mean squared error and the variance of the partial time derivative of the dependent variable (density $n$ or vorticity $\omega$). Therefore, for the density equation we have the following definition of FVU:

(2.7)\begin{equation} {\rm FVU} = \frac{\text{MSE(model)}}{{\rm var}\left(\dfrac{\partial n}{\partial t}\right)}, \end{equation}

and similarly for the vorticity we have

(2.8)\begin{equation} {\rm FVU} = \frac{{\rm MSE(model)}}{{\rm var}\left(\dfrac{\partial \omega}{\partial t}\right)}. \end{equation}

A significant reduction in FVU is observed for the optimal number of features necessary for describing the dynamics underlying the input data.

Figure 1. Code structure: the arrows denote the flow of information and the blocks depict the sequence of steps performed by the sparse regression (SR) algorithm. The colour coding refers to parts of the code which might need to be adapted to a particular physics problem depending on the input data structure and content (red) and the parts which essentially remain the same regardless of the input data structure and content (blue).

2.4. Structural similarity index

Successful regression requires that the sampled data set contains sufficient change in the system dynamics in the nonlinear regime. To quantify this change in our sample we use the so-called, structural similarity index (SSI). The index was originally defined for the purpose of image quality assessment based on degradation of structural information (Wang et al. Reference Wang, Bovik, Sheikh and Simoncelli2004). It is calculated for two given images and has the value in the range $[1,-1]$ where the extremes of the range indicate that the two images are identical (${\rm SSI}=1$) or completely different (opposite) (${\rm SSI}=-1$). The SSI for two images is defined by the following expression:

(2.9)\begin{equation} {\rm SSI}(\boldsymbol{x,y})=[l(\boldsymbol{x},\boldsymbol{y})]^{\alpha}\boldsymbol{\cdot}[c(\boldsymbol{x},\boldsymbol{y})]^{\beta} \boldsymbol{\cdot}[s(\boldsymbol{x},\boldsymbol{y})]^{\gamma}, \end{equation}

where $x$ and $y$ are two images of size $N \times N$ which are being compared based on the luminance $l$, contrast $c$ and structure $s$ functions weighed by the exponents $\alpha$, $\beta$ and $\gamma$. The luminance function $l(\boldsymbol {x},\boldsymbol {y})$ is given by

(2.10)\begin{equation} l(\boldsymbol{x},\boldsymbol{y}) = \frac{2\mu_x \mu_y +C_1}{\mu_x^2+\mu_y^2+C_1}, \end{equation}

where $\mu$ denotes the average of all pixel values (mean luminance intensity) and the constant $C_1$ prevents instability for denominators approaching $0$. Similarly, the contrast function $c(\boldsymbol {x},\boldsymbol {y})$ is given by

(2.11)\begin{equation} c(\boldsymbol{x},\boldsymbol{y})=\frac{2\sigma_x \sigma_y +C_2}{\sigma_x^2+\sigma_y^2+C_2}, \end{equation}

where $\sigma$ denotes the standard deviation of all pixel values and $C_2$ is a constant. Finally, the structure function $s(\boldsymbol {x},\boldsymbol {y})$ is given by

(2.12)\begin{equation} s(\boldsymbol{x},\boldsymbol{y})=\frac{\sigma_{xy}+C_3}{\sigma_x \sigma_y+C_3} \end{equation}

where the correlation coefficient $\sigma _{xy}$ is equal to

(2.13)\begin{equation} \sigma_{xy}=\frac{1}{N-1}\sum_{i=1}^{N} (x_i-\mu_x)(y_i-\mu_y) \end{equation}

Equations (2.9)–(2.13) are given here for the sake of completeness, for more details on the derivation and the SSI metric we refer the reader to the original reference (Wang et al. Reference Wang, Bovik, Sheikh and Simoncelli2004). Our data is a series of snapshots of the configuration space which can be converted to a series of images (see figure 2) to which a sequence of SSI can be attributed.

Figure 2. The HW (a,b) and mHW (c,d) density and vorticity data.