Exploring helical dynamos with machine learning

The forcing function is taken to be homogeneous and isotropic with explicit control over the fractional helicity determined by the $\sigma$ parameter. It is defined as:

where the wavevectors $\text{bi}k\left(t\right)$ and the phase $|\varphi \left(t\right)|<\pi$ are random at each time step. The wavevectors are chosen within a band to specify a range of forcing centered around $\text{bi}{k}_f$. The normalization factor is defined such that the forcing term matches the physical dimensions of the other terms in the Navier Stokes equation, $N=f_0{c}_s(|\text{bi}k|c_s/\delta t{)}^{1/2}$, where $f_0$ is the forcing amplitude, $c_s$ is the sound speed, $\delta t$ is the time step. The Fourier space forcing has the form [24]:

where $\sigma$ is a measure of the helicity of the forcing; for positive maximum helicity, $\sigma =1$, and $\text{bi}{f}_{\text{bi}k}^{nohel}=\text{bi}k\times \widehat{\text{bi}e}/\sqrt{k^2-(\text{bi}k\cdot \widehat{\text{bi}e}{)}^2}$.

The initial velocity and density are set to zero while the magnetic field is initialized through a vector potential with a small amplitude (${10}^{-3}$) Gaussian noise. We define the units [25]:

such that the magnetic field, $\text{bi}B$, is in the units of Alfv́en speed ($c_s\sqrt{{\mu }_0{\rho }_0}$). The domain size is $L^3=(2\pi {)}^3$, which leads to $k_1=2\pi /L=1$. Moreover, we use the following dimensionless numbers:

where $k_f/k_1=10$ is the forcing wavemode. We provide a summary of the simulations in Table 1.

Numerical simulations of large scale astrophysical flows with magnetic Reynolds number well in the excess of ${10}^6$ are prohibitive. Moreover, state-of-the-art simulations requires O(${10}^6$) CPU hours to compute for moderate Reynolds numbers. To make modeling such systems more tractable, it is therefore necessary to look for simpler surrogate models that can capture the most significant properties of the high fidelity simulations (DNS).

Large scale dynamo theory concerns itself with the evolution equation of the filtered induction equation,

where EMF is given as $E=\overline{\text{bi}v\times \text{bi}b}$, $\eta$ is the turbulent resistivity, and $\overline{Q}$ represents a filtered quantity. The filter could be ensemble, temporal or spatial depending on the problem. The first term on the right-hand side is typically ignored unless there is a strong non-constant mean velocity field such as in the case of galaxies (the mean differentially rotating flow for the galactic dynamo is considered in [26]).

In keeping with previous work that has used horizontal averaging as the filter in equation 8 [25, 7], our data has been averaged over the whole $xy$ plane:

This averaging scheme has the interesting property that $\nabla \times (\overline{\text{bi}V}\times \overline{\text{bi}B})=0$ since ${\overline{V}}_z={\overline{B}}_z=0$ (due to incompressibility and periodic boundary conditions). We only consider the fields reduced to 1D by furthermore averaging them in two different ways:

1. $z$ averaged: $⟨\overline{Q}⟩\left(t\right)=1/(z_2-z_1){\int }_{z_1}^{z_2}\overline{Q}dz$ where $z_1$, $z_2$ are chosen such that they do not cover the entire z-domain since that will correspond to volume averaged quantities that are zero.

2. $t$ averaged: $\left[\overline{Q}\right]\left(z\right)=1/(t_2-t_1){\int }_{t_1}^{t_2}\overline{Q}dt$ where $t_1$, $t_2$ are chosen to correspond to either the kinematic or saturation regime for different runs.

For horizontally averaged ($xy$-plane) fields with periodic boundary conditions, ${\overline{B}}_z=0$ due to incompressibility. This leaves two independent magnetic fields and corresponding EMF components. Analytical models that use $E=\alpha \overline{\text{bi}B}$ are termed ${\alpha }^2$ dynamos and have been extensively studied in the literature. In principle, a filter that retains $3$ independent components of the data while using $E=\alpha \overline{\text{bi}B}$ will be termed ${\alpha }^3$ dynamo model.

In Fig. 1, we show a three-dimensional phase-space visualization of $E_x$ vs. ${\overline{B}}_x$ vs. ${\overline{B}}_y$ that shows the relation between the EMF and the mean magnetic fields (situation is the same when considering $E_y$ too). The circle in ${\overline{B}}_x-{\overline{B}}_y$ plane represents the fact that the total magnetic energy ${\overline{B}}_x^2+{\overline{B}}_y^2\sim \text{constant}$. The slight tilt with respect to the EMF (i.e., towards $z$ direction) is a result of the linear correlation with the EMF. The increasing disorder with increasing $Rm$ is characteristic of highly turbulent systems.

Ordinary Least Squares (OLS) assumes that the data is linear, features/terms are independent, variance for each point is roughly constant (homoscedasticity), features/terms do not interact with one another. Moreover, causal inference based on linear regression can be misleading [27], a problem shared by all curve fitting methods. We should highlight the distinction between linearity of the regression as opposed to linearity of the features. For example, one can construct a non-linear basis (example: $x,x^2,x^3,...$ or $\sin x,\cos x,...$) and do linear regression on it. Linear regression on non-linear data can potentially represent non-linear data well assuming that the non-linear basis terms are independent and have low variance. Constructing the correct interacting/non-linear terms requires domain expertise. For this reason, it is important to also consider non-linear regression methods such as random forests that are robust to noise and are capable of modeling interactions and non-linearities.

In this subsection, we apply various methods on time averaged data that only has vertical dependence ($256$ grid points).

In Fig. 2, we show the linear (Pearson) correlation coefficients among the various fields. The correlation between the current densities and the magnetic fields is nearly $1$ indicating they are linearly dependent. Strong linear correlation does not rule out the physical importance of a variable - it only implies that from a curve fitting point of view, it has redundant information. Since we are using helical forcing leading to a helical state, $\overline{\text{bi}J}=\nabla \times \overline{\text{bi}B}\sim \overline{\text{bi}B}$. At high $Rm$, one might expect that the power in the velocity and magnetic field spectrum will be distributed across multiple modes implying weaker correlations. However, Fig. 2 (right) indicates the correlation between $\overline{\text{bi}B}$ and $\overline{\text{bi}J}$ is still perfect. This means that we can eliminate current density as an independent variable. Nonetheless, in keeping with previous work, for linear basis we do consider the current density as an independent variable but eliminate it for the non-linear (polynomial) basis where its inclusion will lead to several redundant terms.

We consider two sets of basis in this subsection:

1. Linear: ${\overline{B}}_x$, ${\overline{B}}_y$, ${\overline{J}}_x$, ${\overline{J}}_y$.

2. Polynomial of $O\left(3\right)$: ${\overline{B}}_x$, ${\overline{B}}_y$, ${\overline{B}}_x{\overline{B}}_y$, ${\overline{B}}^2$, ${\overline{B}}^2{\overline{B}}_x$, ${\overline{B}}^2{\overline{B}}_y$.

We note that for kinetically forced simulations, $v^2\sim \text{constant}$, which means that the total energy conservation equation yields, ${\overline{B}}_x^2+{\overline{B}}_y^2={\overline{B}}^2\sim \text{constant}$ in the saturated state. This makes ${\overline{B}}_x^2$ and ${\overline{B}}_y^2$ linearly dependent terms and thus we do not consider them separately. Another motivation to use ${\overline{B}}^2$ is because work on dynamo quenching often considers $\alpha \sim 1/(1+{\overline{B}}^2/B_{eq}^2)$, which when expanded for the $\text{bi}E$ leads to cubic terms of the form ${\overline{B}}^2{\overline{B}}_i$ ($i=x,y$).

We will be using the publicly available python library scikit-learn (see [38], [39] for introductions) for modeling the data.

Previous work on dynamo theory has often relied on linear regression methods [12, 40] with some effort to incorporate non-linear effects through magnetic helicity conservation ([9, 41]; see also non-linear test field method [42]). In this section, we use linear regression with a regularization term for both the linear and polynomial basis discussed previously.

In Fig. 3, we show the fits using random forests regression and LASSO. Because our data is small (only $256$ points), we use two-fold cross validation to determine various hyperparameters of LASSO and random forests. We found that bootstrapping leads to spurious results so we turned it off. A similar problem was found with random train/test splits: the spatial structure of the fields is like one full sinusoidal wave and randomly picking a fraction of the points leads to strange patterns. Moreover, while deep trees can lead to overfitting, in our case we found that a maximum depth of $8$ was optimal through 2-fold cross validation. The end results are still not as good as LASSO perhaps owing to the fact that the training set is only $204$ entries ($80\%$ of the data) - this is too low for an algorithm like random forest to really excel. That is the key lesson we have learnt while modeling this dataset throughout: “small” datasets mean simpler algorithms will outperform more sophisticated algorithms like random forests.

Feature Ranking

In Table 2, we compare the coefficients from OLS, LASSO, randomized LASSO (stability), Recursive Feature Elimination (RFE) from OLS and $RF{E}_{r}f$ from random forests, Mutual Information Criterion (MIC), Pearson correlation (Corr) and the average of all these predictions in the last column. The meaning of these numbers is different. For OLS and LASSO, the numbers corresponding to the fields ${\overline{B}}_x$ represent coefficients from linear regression and can be negative or positive. For random forests, the numbers represent the relative feature importance based on ‘mean decrease impurity’ that is only positive. These feature importance numbers are scaled to sum to unity. Stability selection is a general term that refers to checking robustness of model predictions by testing it on different random (bootstrapped) subsets of data with different hyperparameter values. In our particular application, we are using randomized LASSO that randomly varies the the LASSO coefficient [43] ⁵⁵5influenced by: https://blog.datadive.net/selecting-good-features-part-iv-stability-selection-rfe-and-everything-side-by-side/.

RFE is a model-agnostic (wrapper) feature selection method that starts with a number of features and removes the most important features recursively. This is an example of ‘backward’ feature selection [44, 28] and we used it with OLS and random forests in this case. Maximal Information Coefficient (MIC) is a measure of mutual information between two (random) variables and is capable of accounting for non-linear relationships that the Pearson correlation coefficient cannot (example: correlation between $x$ and $x^2$ will be nearly zero, but MIC will yield an answer close to one). We point out that we scaled the whole set of features $X_{t}rain$ using standard scaling before calculating coefficients/feature importances. Moreover, because LASSO, OLS, Corr returned negative as well as positive values and RFE returned values ranging from $1$ to $6$, we used minmax scaling to re-scale all values between $0-1$.

Ensembles of machine learning models can lead to improved accuracy [45]. Here we take the simplest possible approach: we take the unweighted average of various predictors to compute the overall score for each feature. While stability selection shows that $B_x$ is the dominant feature, the mean of all predictors gives $B^2{B}_x$ a slightly higher score. This could be because both $B_x$ and $B^2{B}_x$ have the same sinusoidal form.

In contrast to the previous section, we square and then average the fields vertically and study the temporal dependence. The magnetic fields and the EMF show exponential growth in time and thus it was important to do some pre-processing where we took the logarithm of squared averaged fields. Because of this particular pre-processing, we do not consider polynomial expansions like we did before and instead focus on either the linear $E\propto \overline{\text{bi}B}$ form, or the quenched form $E\propto 1/(1+{\overline{B}}^2/B_e{q}^2)$.

For both the vertical and temporal profiles, various models suggest that $E\propto \overline{\text{bi}B}$ is a good fit. But if $d\overline{\text{bi}B}/dt\propto \alpha \overline{\text{bi}B}$, this will lead to exponential growth without any saturation. Is there an inconsistency here? In the induction equation, the fields at time ‘$t$’ are related to the fields at earlier times. The machine learning models we use here are not capable of storing long term memory. The fits shown in this paper are comparing $E$ and $\overline{\text{bi}B}$ in some nearest neighbor sense and because both $E$ and $\overline{\text{bi}B}$ are growing exponentially in time, they appear to be linearly correlated. Capturing the dynamical information requires more sophisticated algorithms that we leave for future work.

We summarize the properties of the different machine learning algorithms considered in this work in Table 3. Linear regression (both OLS and LASSO) are by definition capable of dealing with linear data only. However, they can handle some non-linear properties using a non-linear basis (for example, polynomial regression). The biggest advantage of linear models is the interpretability: the sensitivity of the target variable (EMF) with respect to the features (${\overline{B}}_x,{\overline{B}}_y$) is as easy as determining the coefficients of these terms. Linear models are also good at dealing with data of all sizes but matrix inverses are expensive for large data. Linear models might also be able to model trends in data if the features they consider have the same trend as the target variable.

Random forests are based on ensembles of decision trees. While random forests are quite robust to outliers and missing data, the feature importances computed from random forests can be misleading [32]. Moreoever, random forests are not as interpretable as linear regression or single decision trees. A typical random forest fit can involve $10-100$ trees implying that interpretation is difficult. But computation of feature importances and partial dependence plots help in making sense of random forest results [31, 28]. Random forests typically bootstrap data and use a random subset of features for each decision tree to determine which one leads to the lowest error. These two sources of randomness imply that random forests are best suited for intermediate to big data. Indeed, in our work we find that bootstrap had to be turned off and the best fit models used less than $10$ decision trees.

Bayesian methods offer an interesting complementary toolkit to various machine learning techniques. They also enable us to properly take the noisy nature of DNSs into account. Here the big caveat is, that the model(s) needs to be known beforehand. This, however, further arguments in favour of our hybrid approach applied in this paper. We have used machine learning techniques to detect important features and then applied our physical knowledge to build valid models out from those. Biggest caveat here is that Bayesian model optimization becomes more and more computationally demanding with multidimensional (and often multimodal) data. This puts a limitation on the real-life usability of the method when dealing with increasingly complex data.

The data considered in this paper is either spatially or temporally averaged reducing to $O(100-1000)$ entries. Big data typically refers to the case where the observations/rows are $\ge {10}^6$. “Small” data faces problems of overfitting and is more likely to be sensitive to outliers. Because of these issues, simpler models like linear regression with regularization tend to do well. It is, therefore, no surprise that LASSO has outperformed random forests. For larger datasets with irregular non-linear patterns, random forests and gradient boosting tend to do much better [46].

We used the simplest model: isotropic with no non-local effects in either space and/or time. While earlier work has focused primarily on linear models of dynamos where $\text{bi}E\sim \overline{\text{bi}B}$ using linear regression and the test field method [7], non-linear models have also been considered [42]. We assumed isotropy, which could be a misleading assumption in the presence of shear (see [12, 47]). Furthermore, in our model, we did not incorporate non-local effects that were considered by [12, 13, 14]. Combining machine learning algorithms with anisotropic, local, non-linear models might lead to new insights, and we leave this for future work.

The ${\alpha }^2$ dynamo considered in this work is an idealized system but has the advantage that due to maximally helical forcing and periodic boundary conditions, magnetic helicity is both conserved and is gauge invariant (see [41]). This simplistic setup allows developing analytical theory to explain the origin and saturation of magnetic fields [9, 48]. Numerical results seem to be consistent with theory [25, 49]. In this work, we considered existing analytical formulations and used machine learning tools to test whether the EMF is linear or non-linear in the mean magnetic fields. Machine learning models can complement analytical theory in looking for reduced descriptions of high dimensional systems [18].