# From Causal Pairs to Causal Graphs

Rezaur Rashid  Department of Computer Science
UNC Charlotte
Charlotte, USA
<EMAIL_ADDRESS>Jawad Chowdhury  Department of Computer Science
UNC Charlotte
Charlotte, USA
<EMAIL_ADDRESS>Gabriel Terejanu  Department of Computer Science
UNC Charlotte
Charlotte, USA
<EMAIL_ADDRESS>
###### Abstract

Causal structure learning from observational data remains a non-trivial task
due to various factors such as finite sampling, unobserved confounding
factors, and measurement errors. Constraint-based and score-based methods tend
to suffer from high computational complexity due to the combinatorial nature
of estimating the directed acyclic graph (DAG). Motivated by the ‘Cause-Effect
Pair’ NIPS 2013 Workshop on Causality Challenge, in this paper, we take a
different approach and generate a probability distribution over all possible
graphs informed by the cause-effect pair features proposed in response to the
workshop challenge. The goal of the paper is to propose new methods based on
this probabilistic information and compare their performance with traditional
and state-of-the-art approaches. Our experiments, on both synthetic and real
datasets, show that our proposed methods not only have statistically similar
or better performances than some traditional approaches but also are
computationally faster.

###### Index Terms:

causal pairs, causal graphs, DAG, causal structure learning

##  I Introduction

Machine learning methods, deep learning in particular, have achieved
unparalleled predictive performance in the past two decades. Nevertheless,
these correlation-based models exhibit significant limitations when applied to
out-of-distribution data and prescriptive analytics, which is grounded in
causal inference. Learning the underlying causal structure is an important
task to both constrain a model by reducing spurious correlations [1] and
perform What-If analysis [2]. Learning causal relationships from observational
data, also known as causal discovery, remains an active and challenging
research topic [2, 3, 4].

Several causal discovery methods have been proposed in the literature.
Constraint-based approaches learn the causal skeleton using conditional
independence test using the joint probability distribution of the data and
identify edge directions up to their Markov equivalence class [3, 5, 6, 7, 8].

Score-based approaches learn the causal graph $\mathcal{G}$ by optimizing a
score function generally computed with respect to observational data [9, 10,
11]. Unfortunately, these methods suffer from super-exponential computational
complexity in the number of nodes. Tsamardinos et al. [12] propose a hybrid
method where they use a constraint-based approach to reduce the search space
in score-based methods. However, this method relies on local heuristics and
lacks a standard way of choosing score functions and search strategies [13].

A promising direction, NOTEARS [14], formulates a smooth characterization of
acyclicity that can be incorporated into a continuous optimization and solved
using well-known numerical methods. NOTEARS was later extended to parametric
nonlinear models and nonparametric models [15]. GOLEM[16] also adopts a
continuous optimization framework, however, it makes use of a linear DAG
learning model and doesn’t capture non-linear relationships. A different
approach looks at identifying cause-effect pairs using the statistical
techniques from observational data [17, 18]. Singh et al.[19] use deep
convolutional neural network (CNN) models to determine the directions of
pairwise causal edges from observational data. Hassanzadeh et al.[20]
formulate the pairwise causal discovery techniques as binary causal problems
where they try to answer if there exist any causal relations between two
variables in the context of Natural Language Processing (NLP). Nevertheless,
they have not studied how the predicted edge directions can be used to provide
a solution to causal graph identification.

Motivated by applications in biological networks, Medvedovsky et al.[21]
propose an approximation algorithm to orient a graph by maximizing the number
of pairs that admit a directed path from known pairs of sources and targets.
Nevertheless, given the peculiarities of the application based on prior
knowledge constrains, this approach falls short from identifying the entire
causal graph. Therefore, discovering the causal graph, including pairwise
causal relations from observational data remains a challenging task due to
various factors such as finite sampling and measurement errors. In general,
identifying cause-effect relationships requires controlled experimentation
which is expensive and/or technically and ethically impossible to perform
[22].

We propose a probabilistic approach to discover causal structures using the
cause-effect pairs features proposed in response to the ‘Cause-Effect Pair’ at
the NIPS 2013 Workshop on Causality challenge. The paper introduces the
following novel contributions: (1) generate a probability distribution over
all the edges of a digraph using various statistical and information-theoretic
features that describe the relationships between any two variables in the
dataset; (2) generate the most likely probability distribution of directed
acyclic graphs (DAG) using the maximum spanning DAG; (3) generate an
approximate solution to the causal graph problem by estimating the digraph and
DAG using maximum likelihood estimate with the probability distributions in
(1) and (2) respectively; (4) overall, our proposed methods are comparable
with traditional ones (PC, GES), while benefiting from polynomial time
complexity as compared to super-exponential time complexity; and (5) finally,
by comparing with state-of-the-art methods such as NOTEARS-MLP, we show that
future improvements are possible by further leveraging global graph
information.

Section II introduces the problem formulation and details our causal discovery
approach based on causal-effect pairs. The empirical evaluation of our methods
is described in Section III. Lastly, in Section IV we present our findings in
brief and opportunities for future improvements.

##  II Methodology

Given $n$ i.i.d. observations in the data matrix
$\mathbf{X}=[\mathbf{x_{1}}\ldots\mathbf{x_{d}}]\in\mathbb{R}^{n\times d}$,
the goal of causal discovery is to estimate the underlying causal relations
encoded by the directed acyclic graph (DAG), $\mathcal{G}_{\text{DAG}}=(V,E)$.
$V$ comprises of nodes corresponding to the observed random variables $X_{i}$
for $i=1\ldots d$ and the edges in $E$ correspond to the causal relations
encoded by $\mathcal{G}_{\text{DAG}}$. Namely, the existence of the edge $i\to
j$ corresponds to a direct causal relationship between $X_{i}$ (the cause) and
$X_{j}$ (the effect).

Our approach is to leverage the work on cause-effect pairs which uses machine
learning to predict the probability distribution $p(e_{ij}|f)$ of causal
relation between two variables $X_{i}$ and $X_{j}$ given the observational
dataset $[\mathbf{x_{i}},\mathbf{x_{j}}]\in\mathbb{R}^{n\times 2}$.

$\displaystyle p(e_{ij}|f)$ $\displaystyle=$ $\displaystyle
f([\mathbf{x_{i}},\mathbf{x_{j}}]),~{}\text{for}~{}i<j$ (1)

Here, we assume that $f(\cdot)$ is a pre-trained machine learning model and
$e_{ij}\in[-1,0,1]$.

$e_{ij}=\begin{cases}-1:&\text{$j\to i$, causal relation exists from $X_{j}$
to $X_{i}$}\\\ \;\;\>0:&\text{$i\not\to j$ and $j\not\to i$, no direct causal
relation}\\\ &\text{between $X_{i}$ and $X_{j}$}\\\ \;\;\>1:&\text{$i\to j$
causal relation exists from $X_{i}$ to $X_{j}$}\end{cases}$

After calculating the probability distributions of causal relations between
all the pairs in the dataset, a naive approach to construct the probability
distribution of a digraph $\mathcal{G}$ is to assume that the causal-pairs are
independent. In Section II-B, we show that our proposed approach of enforcing
DAGness does correlate these causal-pairs and provides us with additional
global information to constrain the graph probability distribution and allows
us to move beyond the initial edge independence.

$p(\mathcal{G}|f)=\prod_{i<j}p(e_{ij}|f)$ (2)

Given this rich probabilistic information on all the causal relationships in
the dataset, one may choose to generate the maximum likelihood digraph.

$\mathcal{G}_{\text{ML}}=\arg\max_{\mathcal{G}}p(\mathcal{G}|f)$ (3)

Note that the samples from the probability distribution, Eq. 2, and the
maximum likelihood estimate, Eq. 3, are digraphs with no guarantees that they
are acyclic. In the following, we propose to generate the most likely
probability distribution of directed acyclic graphs (DAG) using the maximum
spanning DAG approach [23] as well as estimate a representative DAG using the
maximum likelihood estimate.

###  II-A Developing causal-pair models

The model $f(\cdot)$ in Eq. 1 can be trained using synthetic datasets or real
datasets with known causal relations. Given a set of labeled datasets
$\\{([\mathbf{x_{i}},\mathbf{x_{j}}],y_{ij})_{k}\\}$, in this paper we take
the approach of engineering features from this dataset using various
statistical and information-theoretic measures such as: minimum or maximum
value of a variable; number of unique samples of a variable; entropy, mutual
information, uniform divergence; slope-based information geometric causal
inference (IGCI), Hilbert Schmidt independence criterion (HSIC); Pearson R
coefficient; Spearman’s rank coefficient; moments and mixed moments such as
skewness and kurtosis. Therefore, the machine learning model $f(\cdot)$ can be
trained on the new engineered dataset with features previously introduced. We
note that our proposed methodology is agnostic to the features deployed and it
works with any causal-pairs model. Also, the computational runtime of
calculating any of these features such as HSIC, Pearson R coefficient has no
effect asymptotically on the computational complexity - it just increases the
constant of the polynomial runtime. We also note that additional improvements
might be brought by developing deep neural network architectures capable to
extract informative representations for predicting the target
$y_{ij}\in[-1,0,1]$ directly from the sample dataset, but it is left as future
work.

###  II-B Enforcing DAGness

In this section, we propose to derive a probability distribution that
guarantees that the sample graphs are DAGs. This takes the form of the
probability distribution in Eq. 4 which unlike Eq. 2 contains the DAGness
condition.

$p(\mathcal{G}|f,\text{DAG})=\sum_{\pi}p(\mathcal{G}|f,\text{DAG},\pi)p(\pi|f)$
(4)

Due to computational intractability, we have chosen to build this conditional
distribution not using the Bayes rule and utilizing Eq. 2 as prior, but rather
as the law of total probability where we integrate out the topological
ordering $\pi$ of the vertices. We note that for a topological ordering where
node $i$ comes before node $j$, it implies that a directed edge can only
happen from $i$ to $j$. We also note that the possibility of no causal
relation between $i$ and $j$ is not excluded in this context. Both causal,
$i\rightarrow j$, and noncausal, $i\not\to j$ and $j\not\to i$, are possible.

To generate a representative DAG one can use the maximum likelihood estimate.
This however is intractable, and we do assume that there is a topological
sorting of vertices that also covers our maximum likelihood DAG.

$\displaystyle\mathcal{G}_{\text{DAG}}$ $\displaystyle=$
$\displaystyle\arg\max_{\mathcal{G}}p(\mathcal{G}|f,\text{DAG})$ (5)
$\displaystyle\approx$
$\displaystyle\arg\max_{\mathcal{G}}p(\mathcal{G}|f,\text{DAG},\pi_{\text{ML}})$
(6)

We propose to approximate the topological ordering, $\pi_{\text{ML}}$, by the
topological sorting of the Maximum Spanning DAG (MSDAG) [23] of the induced
weighted graph by the probability of causal relations.

$\displaystyle\pi_{\text{ML}}$ $\displaystyle=$
$\displaystyle\arg\max_{\pi}p(\pi|f)$ (7) $\displaystyle\approx$
toposort(MSDAG($\mathcal{G}_{A}$)) (8)

We build the following weighted adjacency matrix $A\in\mathcal{R}^{d\times
d}$, which contains the probability of all directed edges as weights.

$\displaystyle A[i,j]=p(e_{ij}=1|f)$ (9) $\displaystyle A[j,i]=p(e_{ij}=-1|f)$

Let $\mathcal{G}_{A}$ be a weighted graph induced by the adjacency matrix $A$.
The goal is to find the topological sorting of the MSDAG of $\mathcal{G}_{A}$.
The motivation is to accommodate as many directed edges with large
probabilities as possible. We use the approach introduced by [24] to
approximate the MSDAG by first constructing the maximum spanning tree and
greedily adding edges in the descending order of the weights as long as no
cycles are formed. Note that a topological sorting derived this way still
accommodates the possibility of no edges to account for their probability in
the maximum likelihood DAG, Eq. 6.

$p(\mathcal{G}|f,\text{DAG},\pi_{\text{ML}})=\prod_{\pi^{-1}_{\text{ML}}[i]<\pi^{-1}_{\text{ML}}[j]}p(e_{i\rightarrow
j}|f)$ (10)

Given the maximum likelihood of topological ordering, one can easily calculate
the probability distribution in Eq. 6 by constraining the direction of the
edge based on the node ordering, see Eq. 10. In this context, we are left with
only two possibilities when node $i$ appears before $j$ in $\pi_{\text{ML}}$.
Either there is an edge from $i$ to $j$ or there is no edge between them.

$e_{i\rightarrow j}=\begin{cases}1:&\text{$i\to j$, causal relation exists
from $X_{i}$ to $X_{j}$}\\\ 0:&\text{$i\not\to j$ and $j\not\to i$, no direct
causal relation}\\\ &\text{between $X_{i}$ and $X_{j}$}\end{cases}$

By constraining the direction of edges we need to re-normalizing the edge
probabilities as follows.

$p(e_{i\rightarrow j}=1|f)=\frac{p(e_{ij}=1|f)}{p(e_{ij}=1|f)+p(e_{ij}=0|f)}$
(11)

Enforcing DAGness in Eq. 10 and consequently Eq. 6 provides us with additional
global information to constrain the graph probability distribution and allows
us to move beyond the initial edge independence in Eq. 2 and Eq. 3
respectively, which was derived from just pair-wise (local) information.

##  III Experiments

We have used the following labels for our approaches: PG given by Eq. 2, MLG
given by Eq. 3, PDAG given by Eq. 10, and MLDAG given by Eq. 6. We show the
empirical results of our approaches applied on both synthetic and real-world
datasets and we compare our performance with two traditional approaches: the
PC algorithm [3] and the GES algorithm [9] and with a state-of-the-art
approach, NOTEARS-MLP [15]. Both PC and GES have public implementations and we
used the code111PC: https://github.com/keiichishima/pcalg,222GES:
https://github.com/juangamella/ges from their git repositories. In addition,
for the NOTEARS-MLP[15] approach, we followed the implementation stated in the
paper and git repository333NOTEARS-MLP: https://github.com/xunzheng/notears.
We used the default settings and default hyper-parameters for all these three
implementations. However, Reisach et. al. [25] in their paper highlight that
continuous score-based approaches i.e. NOTEARS-MLP[15] in particular suffer
highly from data scaling which was addressed from a theoretical perspective.
Therefore, we standardize the features for both the synthetic and real-world
datasets by removing the mean and scaling to unit variance.

###  III-A Data Sets

#### Cause-Effect Pair Train Data

To train our model, we have used the Cause-effect pairs dataset444Cause-Effect
Pair Dataset: https://www.kaggle.com/c/cause-effect-pairs/data from the NIPS
2013 Workshop on Causality. The train data contains 4050 samples of pairs
(attribute A and attribute B) of real variables with known causal
relationships and every node within a pair has the same number of data
samples. In addition, these known ground truths are derived from expert
domains such as chemistry, ecology, engineering, medicine, physics, sociology,
etc., and these pairs are intermixed with controls such as pairs of
independent variables and pairs of dependent variables but not causally
related.

#### Synthetic Test Data

To evaluate the performance of our methods on causal graph estimation, we have
generated synthetic data for testing. We have considered 16 types of different
graph combinations having similar criteria: number of nodes, $d=[10,20]$,
number of edges, $e=[1d,4d]$, number of data samples per node, $n=[200,1000]$,
and graph models from Erdos-Renyi(ER) and Scale-Free (SF). We have generated
non-linear data samples for the graph nodes similar to data generation
utilities available in the NOTEARS-MLP implementation. In addition, for each
of these 16 graph types, we have generated 10 random graph structures with
ground truths to test our methods. The outputs are then summarized over these
10 graph structures to report our results for all 16 graph combinations.

###  III-B Metrics

We consider three performance metrics to evaluate the causal graphs: True
Positive Rate (TPR), False Positive Rate (FPR) and Structural Hamming Distance
(SHD). A lower SHD and FPR indicate a better performance whereas a higher TPR
is better. However, since both PC and GES may generate outputs with undirected
edges, we treated an undirected edge as a true edge with a probability of
$0.5$ and a false edge with the same probability. SHD, TPR, and FPR were
implemented by their definition for PC, GES, NOTEARS-MLP, and our two maximum
likelihood estimates (MLG given by Eq. 3 and MLDAG given by Eq. 6).

As for our probabilistic approaches (PG given by Eq. 2 and PDAG given by Eq.
10), given a true graph $\mathcal{G}_{true}$ and an adjacency matrix
$\mathcal{A}$ with edge probabilities, we calculate SHD using Eq. 12, TPR
using Eq. 14 and FPR using Eq. 16.

$\displaystyle SHD=$ $\displaystyle\sum_{\begin{subarray}{c}{i<j}\\\ {(i,j)\in
E(\mathcal{G}_{true})}\\\ {(j,i)\not\in
E(\mathcal{G}_{true})}\end{subarray}}(1-\mathcal{A}[i,j])+\sum_{\begin{subarray}{c}{i<j}\\\
{(i,j)\not\in E(\mathcal{G}_{true})}\\\ {(j,i)\in
E(\mathcal{G}_{true})}\end{subarray}}(1-\mathcal{A}[j,i])$ (12)
$\displaystyle+\sum_{\begin{subarray}{c}{i<j}\\\ {(i,j)\not\in
E(\mathcal{G}_{true})}\\\ {(j,i)\not\in
E(\mathcal{G}_{true})}\end{subarray}}(\mathcal{A}[i,j]+\mathcal{A}[j,i])$
$\displaystyle TP=\sum_{\begin{subarray}{c}{i<j}\\\ {(i,j)\in
E(\mathcal{G}_{true})}\\\ {(j,i)\not\in
E(\mathcal{G}_{true})}\end{subarray}}\mathcal{A}[i,j]+\sum_{\begin{subarray}{c}{i<j}\\\
{(i,j)\not\in E(\mathcal{G}_{true})}\\\ {(j,i)\in
E(\mathcal{G}_{true})}\end{subarray}}\mathcal{A}[j,i]$ (13) $\small
TPR=\frac{TP}{max(|E(\mathcal{G}_{true})|,1)}$ (14)

$\displaystyle FP=$ $\displaystyle\sum_{\begin{subarray}{c}{i<j}\\\ {(i,j)\in
E(\mathcal{G}_{true})}\\\ {(j,i)\not\in
E(\mathcal{G}_{true})}\end{subarray}}\mathcal{A}[j,i]+\sum_{\begin{subarray}{c}{i<j}\\\
{(i,j)\not\in E(\mathcal{G}_{true})}\\\ {(j,i)\in
E(\mathcal{G}_{true})}\end{subarray}}\mathcal{A}[i,j]$ (15)
$\displaystyle+\sum_{\begin{subarray}{c}{i<j}\\\ {(i,j)\not\in
E(\mathcal{G}_{true})}\\\ {(j,i)\not\in
E(\mathcal{G}_{true})}\end{subarray}}(\mathcal{A}[i,j]+\mathcal{A}[j,i])$
$\small FPR=\frac{FP}{max((M-|E(\mathcal{G}_{true})|),1)}$ (16)

Here, $M=\frac{d(d-1)}{2}$ is the number of possible edges of the graph
$\mathcal{G}_{true}$ and $d$ is the total number of nodes in the graph.

Note that since we are calculating these metrics over 160 different graph
structures of various sizes in our test data, we report a normalized SHD over
the number of graph nodes (SHD/d) as our SHD measure. TPR and FPR are
normalized by definition.

TABLE I: Edge probability model trained on cause-effect pairs data provided at
the NIPS 2013 Workshop on Causality. The means and standard errors of the
performance metrics are based on the 80 Erdos-Renyi (ER) graph structures in
the test data.
Metrics |  |  PG (Eq. 2)
---
|  MLG (Eq. 3)
---
|  PDAG (Eq. 10)
---
|  MLDAG (Eq. 6)
---
PC | GES | NOTEARS-MLP
SHD/$d$ |  2.38$\pm$0.14 |  2.32$\pm$0.17 |  2.30$\pm$0.15 |  2.18$\pm$0.16 |  2.40$\pm$0.21 |  1.78$\pm$0.13 |  1.33$\pm$0.10
TPR |  0.39$\pm$0.02 |  0.15$\pm$0.02 |  0.38$\pm$0.02 |  0.28$\pm$0.02 |  0.17$\pm$0.02 |  0.48$\pm$0.02 |  0.58$\pm$0.02
FPR |  0.72$\pm$0.10 |  0.07$\pm$0.01 |  0.61$\pm$0.09 |  0.29$\pm$0.05 |  0.22$\pm$0.04 |  0.87$\pm$0.15 |  0.32$\pm$0.06
TABLE II: Edge probability model trained on cause-effect pairs data provided
at the NIPS 2013 Workshop on Causality. The means and standard errors of the
performance metrics are based on the 80 Scale-Free (SF) graph structures in
the test data.
Metrics |  |  PG (Eq. 2)
---
|  MLG (Eq. 3)
---
|  PDAG (Eq. 10)
---
|  MLDAG (Eq. 6)
---
PC | GES | NOTEARS-MLP
SHD/$d$ |  2.02$\pm$0.12 |  1.97$\pm$0.13 |  1.96$\pm$0.12 |  1.88$\pm$0.13 |  1.93$\pm$0.15 |  1.43$\pm$0.11 |  1.36$\pm$0.11
TPR |  0.31$\pm$0.01 |  0.12$\pm$0.01 |  0.30$\pm$0.01 |  0.20$\pm$0.01 |  0.17$\pm$0.02 |  0.51$\pm$0.03 |  0.47$\pm$0.02
FPR |  0.26$\pm$0.02 |  0.03$\pm$0.01 |  0.21$\pm$0.02 |  0.09$\pm$0.01 |  0.08$\pm$0.01 |  0.26$\pm$0.04 |  0.12$\pm$0.02

###  III-C Simulation

For our implementation, we first extract features from the data-pairs using
the feature extraction method of Team-Jarfo[26], the second winner from the
NIPS 2013 Workshop on Causality challenge. We have also used the causal pairs
model of the third winning team from the above-mentioned competition and we
have found that it doesn’t have better performance than Team-Jarfo[26] as
expected. We couldn’t run a performance analysis with the model proposed by
the No.1 team due to the lack of availability of their implemented code. We
have extracted 130 features from the pairs for all training and testing
datasets using the code555Team-Jarfo Code: https://github.com/jarfo/cause-
effect of Team-Jarfo[26]. This feature set contains some standard statistical
features as well as new measures based on variable measures of the conditional
distribution. We train a multi-classifier model based on LightGBM [27] which
is a Gradient Boosting Decision Tree (GBDT) algorithm developed by Microsoft.
LightGBM speeds up the training process by using a histogram-based algorithm
[28, 29] and combines weak learners into strong ones using an iterative
approach [30] to optimize parallel learning. We create the LightGBM classifier
from the Python library PyCaret [31] and select the hyper-parameters by tuning
the model using the PyCaret ‘tune_model()’ function optimized over AUC.

We train our classifier model on cause-effect pairs data using the 130
extracted features and make predictions on the synthetic testing data. We
predict a probability distribution over the three classes of edge directions
(backward edge, no edge, forward edge) for all edges in the testing dataset.
Finally, using the methods described in Section II we calculate the metrics of
our causal graph estimation methods from the predicted probability
distribution and compare the results with the benchmark approaches.

Table I and Table II show the empirical results of our methods applied to 80
different Erdos-Renyi graph structures and 80 different Scale-Free graphs,
respectively. From these two tables, we see that estimating causal graphs is
more challenging for Erdos-Renyi graph structures than Scale-Free graphs for
all the methods. We also note that the performance of NOTEARS-MLP is superior
to all other methods in terms of SHD and TPR in particular and it does have
implications for further developing our proposed methods as detailed later.

Regardless of the graph structure, we make the following observations. (1) PG
performs better than PC in terms of SHD and TPR and is better than GES in
terms of FPR. (2) MLG does improve the FPR but at the cost of degrading the
TPR. (3) We do see a significant improvement by enforcing DAGness. Namely,
PDAG does improve FPR over PG, but not sufficient to be statistically
better/similar to PC. (4) However, this does happen when we take the maximum
likelihood of the conditional probability. Namely, MLDAG performs at the same
FPR level as PC.

From these results, it becomes clear that using global information by
constraining the graphs to be DAGs does improve the performance compared to
using the pair-wise probabilities naively. However, when we compare it with
NOTEARS-MLP, it also becomes clear that this is not sufficient and that the
next iteration of methods needs to develop features that intrinsically exploit
global information.

Nevertheless, one of the distinguishable advantages our methods have over PC
and GES is that they not only perform statistically better/similar but also
they have significantly low computational complexity. While both PC and GES
have exponential time complexity due to their combinatorial approach, our
methods run in polynomial time $\mathcal{O}(d^{2})$ in the number of nodes $d$
as they exploit local node pairs information.

###  III-D Real-World Data

For real data, we consider the dataset published by [32], which is based on
the expression level of proteins. This protein network is largely used in the
scientific community as a real application due to the consensus ground truths
of the graph structure. It has $11$ different protein cells which are
considered as graph nodes $d$. We considered $e=17$ edges as ground truths
representing the protein signaling network. We aggregated the $9$ different
data files, resulting in a sample size $n=7466$ in our experiment.

TABLE III: Comparison of our probabilistic methods with GES and NOTEARS-MLP
that were applied on protein network dataset using cause-effect pairs as
training data.
Metrics |  |  PG (Eq. 2)
---
|  MLG (Eq. 3)
---
|  PDAG (Eq. 10)
---
|  MLDAG (Eq. 6)
---
GES | NOTEARS-MLP
Predicted Edges | 36.14 | 9.82 | 33.16 | 18.48 | 34 | 42.23
Correct Edges | 6.7 | 3.04 | 7.42 | 4.91 | 5.5 | 5.83
Reversed Edges | 7.77 | 4.26 | 6.62 | 5.41 | 9.5 | 7.18
TABLE IV: Comparison of our probabilistic methods with GES and NOTEARS-MLP
(results reported from the original manuscript[15]). These methods were
applied on non-standardized protein network dataset.
Metrics |  |  PG (Eq. 2)
---
|  MLG (Eq. 3)
---
|  PDAG (Eq. 10)
---
|  MLDAG (Eq. 6)
---
GES | NOTEARS-MLP
Predicted Edges | 38.01 | 10.41 | 34.81 | 20.60 | 34 | 13
Correct Edges | 6.21 | 1.52 | 6.47 | 4.71 | 5.5 | 7
Reversed Edges | 8.26 | 4.04 | 7.49 | 6.32 | 9.5 | 3

Unlike synthetic data sets, we have considered three different metrics for
this protein dataset: the total number of predicted edges, the number of
correct edge predictions, and the number of reversed edge predictions. Since
this protein network uses a consensus over the number of true edges ($17$
known edges as ground truth), we do not know the actual true graph for the
entire network. Therefore, a metric such as SHD becomes meaningless.
Furthermore, since GES algorithms generate graphs that contain bidirectional
edges, similarly to the synthetic results, we considered the bidirectional
edge as increasing the number of corrected edges with $0.5$ and the number of
reversed edges with $0.5$ as well.

In Table III, we show the performance evaluation of our methods applied to the
protein network dataset. We note that both PG and PDAG perform better than GES
and NOTEARS-MLP in terms of predicting the number of correct edges. In
addition, we find that taking the maximum likelihood estimates, MLG and MLDAG,
significantly reduce the number of false edge predictions but their
performances in predicting correct edges also degrade. The degradation is less
severe between PDAG and MLDAG than PG and MLG.

PG, PDAG, as well as GES have the best performance in terms of the total
number of predicted edges than NOTEARS-MLP, but MLDAG has the best total
number of predicted edges which is close to the number of true edges in the
dataset. We also note that MLG and MLDAG have better performance than NOTEARS-
MLP in terms of fewer reversed edges. Finally, PGDAG shows an overall better
performance than PG, further demonstrating the impact of enforcing DAGness.

Sensitivity to data scaling. Table IV shows the results of our methods as well
as GES and NOTEARS-MLP [15] applied on protein dataset before scaling them. We
observe that our methods: PG, PDAG, MLG, MLDAG have almost similar results in
both Table III (after data scaling) and Table IV (before data scaling). GES is
not affected by the data scaling. However, as expected, we see that NOTEARS-
MLP has very different results in these two tables where it suffers highly
from data-scaling, which is consistent with the sensitivity to scaling results
in Ref.[25]. It is to mention that in Table IV, NOTEARS-MLP results are
reported from the original manuscript [15], whereas in Table III, we use the
same implementation3 of NOTEARS-MLP [15] with default parameters and applied
on the standardized protein dataset by removing the mean and scaling to unit
variance.

##  IV Conclusions

In this paper, we have introduced a novel approach to causal discovery by
leveraging the probabilistic information of pairwise causal edges. We have
proposed to go beyond the naive approach to generate graph probabilities from
causal pair probabilities by enforcing the graph to be acyclic and
approximating its solution using the maximum spanning directed acyclic graph
approach. Enforcing acyclicity clearly improves the performance on both
synthetic and real datasets compared with the naive approach.

We have shown that our methods have statistically better and/or similar
performances than some traditional methods. More importantly, this performance
comes with just polynomial run-time as compared with the exponential run-time
of traditional methods that are combinatorial in nature. These promising
results prompt us to further look into improving the causal pair feature
generation to intrinsically capture global information.

## Acknowledgement

Research was sponsored by the Army Research Office and was accomplished under
Grant Number W911NF-22-1-0035. The views and conclusions contained in this
document are those of the authors and should not be interpreted as
representing the official policies, either expressed or implied, of the Army
Research Office or the U.S. Government. The U.S. Government is authorized to
reproduce and distribute reprints for Government purposes notwithstanding any
copyright notation herein.

## References

  * [1] B. M. Lake, T. D. Ullman, J. B. Tenenbaum, and S. J. Gershman, “Building machines that learn and think like people,” Behavioral and brain sciences, vol. 40, 2017.
  * [2] J. Pearl, Causality. Cambridge university press, 2009.
  * [3] P. Spirtes, C. N. Glymour, R. Scheines, and D. Heckerman, Causation, prediction, and search. MIT press, 2000\.
  * [4] J. Peters, D. Janzing, and B. Scholkopf, Elements of causal inference: foundations and learning algorithms. The MIT Press, 2017.
  * [5] J. Pearl and T. S. Verma, “A theory of inferred causation,” in Studies in Logic and the Foundations of Mathematics. Elsevier, 1995, vol. 134, pp. 789–811.
  * [6] D. Colombo and M. H. Maathuis, “Order-independent constraint-based causal structure learning,” J. Mach. Learn. Res., 2014, 15(1), pp.3741-3782.
  * [7] X. Sun, D. Janzing, B. Scholkopf, and K. Fukumizu, “A kernel-based causal learning algorithm,” in Proceedings of the 24th international conference on Machine learning, 2007, pp. 855–862.
  * [8] K. Zhang, J. Peters, D. Janzing, and B. Scholkopf, “Kernel-based conditional independence test and application in causal discovery,” arXiv preprint arXiv:1202.3775, 2012.
  * [9] D. M. Chickering, “Optimal structure identification with greedy search,” Journal of machine learning research, vol. 3, no. Nov, pp. 507–554, 2002.
  * [10] J. Ramsey, M. Glymour, R. Sanchez-Romero, and C. Glymour, “A million variables and more: the fast greedy equivalence search algorithm for learning high-dimensional graphical causal models, with an application to functional magnetic resonance images,” International journal of data science and analytics, vol. 3, no. 2, pp. 121–129, 2017.
  * [11] P. Nandy, A. Hauser, and M. H. Maathuis, “High-dimensional consistency in score-based and hybrid structure learning,” The Annals of Statistics, vol. 46, no. 6A, pp. 3151–3183, 2018.
  * [12] I. Tsamardinos, L. E. Brown, and C. F. Aliferis, “The max-min hill-climbing bayesian network structure learning algorithm,” Machine learning, vol. 65, no. 1, pp. 31–78, 2006.
  * [13] S. Zhu, I. Ng, and Z. Chen, “Causal discovery with reinforcement learning,” arXiv preprint arXiv:1906.04477, 2019\.
  * [14] X. Zheng, B. Aragam, P. K. Ravikumar, and E. P. Xing, “Dags with no tears: Continuous optimization for structure learning,” Advances in Neural Information Processing Systems, vol. 31, 2018.
  * [15] X. Zheng, C. Dan, B. Aragam, P. Ravikumar, and E. Xing, “Learning sparse nonparametric dags,” in International Conference on Artificial Intelligence and Statistics. PMLR, 2020, pp. 3414–3425.
  * [16] I. Ng, A. Ghassami, and K. Zhang. On the role ofsparsity and DAG constraints for learning linear DAGs.InNeurIPS,2020
  * [17] I. Guyon, O. Goudet, and D. Kalainathan, “Evaluation methods of cause-effect pairs,” in Cause Effect Pairs in Machine Learning. Springer, 2019, pp. 27–99.
  * [18] J. A. Fonollosa, “Conditional distribution variability measures for causality detection,” in Cause Effect Pairs in Machine Learning. Springer, 2019, pp. 339–347.
  * [19] K. Singh, G. Gupta, L. Vig, G. Shroff, and P. Agarwal, “Deep convolutional neural networks for pairwise causality,” arXiv preprint arXiv:1701.00597, 2017.
  * [20] O. Hassanzadeh, D. Bhattacharjya, M. Feblowitz, K. Srinivas, M. Perrone, S. Sohrabi, and M. Katz, “Answering binary causal questions through large-scale text mining: An evaluation using cause-effect pairs from human experts.” in IJCAI, 2019, pp. 5003–5009.
  * [21] A. Medvedovsky, V. Bafna, U. Zwick, and R. Sharan, “An algorithm for orienting graphs based on cause-effect pairs and its applications to orienting protein networks,” in International Workshop on Algorithms in Bioinformatics. Springer, 2008, pp. 222–232.
  * [22] O. Stegle, D. Janzing, K. Zhang, J. M. Mooij, and B. Scholkopf, “Probabilistic latent variable models for distinguishing between cause and effect,” Advances in neural information processing systems, vol. 23, pp. 1687–1695, 2010.
  * [23] N. Schluter, “On maximum spanning DAG algorithms for semantic DAG parsing,” in Proceedings of the ACL 2014 Workshop on Semantic Parsing. Baltimore, MD: Association for Computational Linguistics, Jun. 2014, pp. 61–65. [Online]. Available: https://aclanthology.org/W14-2412
  * [24] R. McDonald and F. Pereira, “Online learning of approximate dependency parsing algorithms,” in 11th Conference of the European Chapter of the Association for Computational Linguistics. Trento, Italy: Association for Computational Linguistics, Apr. 2006, pp. 81–88. [Online]. Available: https://aclanthology.org/E06-1011
  * [25] A. Reisach, C. Seiler, and S. Weichwald, ”Beware of the simulated dag! causal discovery benchmarks may be easy to game.” Advances in Neural Information Processing Systems, 34, 2021, 27772-27784, in press.
  * [26] J. A. R. Fonollosa, ”Conditional distribution variability measures for causality detection.” Cause-Effect Pairs in Machine Learning. Springer, Cham, 2019. 339-347.
  * [27] G. Ke, Q. Meng, T. Finley, T. Wang, W. Chen, W. Ma, Q. Ye, and T. Y. Liu, “Lightgbm: A highly efficient gradient boosting decision tree,” Advances in neural information processing systems, vol. 30, pp. 3146–3154, 2017\.
  * [28] S. Ranka and V. Singh, “Clouds: A decision tree classifier for large datasets,” in Proceedings of the 4th knowledge discovery and data mining conference, vol. 2, no. 8, 1998\.
  * [29] R. Jin and G. Agrawal, “Communication and memory efficient parallel decision tree construction,” in Proceedings of the 2003 SIAM international conference on data mining. SIAM, 2003, pp. 119–129.
  * [30] Y. Hua, “An efficient traffic classification scheme using embedded feature selection and lightgbm,” in 2020 Information Communication Technologies Conference (ICTC). IEEE, 2020, pp. 125–130.
  * [31] M. Ali, PyCaret: An open source, low-code machine learning library in Python, July 2020, pyCaret version 2.3. [Online]. Available: https://www.pycaret.org
  * [32] K. Sachs, O. Perez, D. Pe’er, D. A. Lauffenburger, and G. P. Nolan, “Causal protein-signaling networks derived from multiparameter single-cell data,” Science, vol. 308, no. 5721, pp. 523–529, 2005.


# Fate of bubble clusters rising in a quiescent liquid

Tian Ma1,∗<EMAIL_ADDRESS>Hendrik Hessenkemper1,∗<EMAIL_ADDRESS>Dirk Lucas1  Andrew D. Bragg2<EMAIL_ADDRESS>1Helmholtz-Zentrum
Dresden – Rossendorf, Institute of Fluid Dynamics, 01328 Dresden, Germany
2Department of Civil and Environmental Engineering, Duke University, Durham,
NC 27708, USA

###### Abstract

We use experiments to study the evolution of bubble clusters in a swarm of
freely rising, deformable bubbles. A new machine learning-aided algorithm
allows us to identify and track bubbles in clusters and measure the cluster
lifetimes. The results indicate that contamination in the carrier liquid can
enhance the formation of bubble clusters and prolong the cluster lifetimes.
The mean bubble rise velocities conditioned on the bubble cluster size are
also explored, and we find a positive correlation between the cluster size and
the rise speed of the bubbles in the cluster, with clustered bubbles rising up
to $20\%$ faster than unclustered bubbles.

**footnotetext: These authors contributed equally to this work

##  1 Introduction

Bubbles rising though a liquid is a frequently occurring situation in nature
(e.g. bubble plumes rising from the bottom of a lake), daily life (e.g. bubble
chains rising in carbonated drinks) and technology (e.g. waste water
treatment). For the case of bubble swarms rising freely in a quiescent fluid,
the bubbles tend to distribute inhomogeneously, spontaneously forming
clusters. When the bubbles are large enough, the rising bubbles generate
turbulence in the liquid, a phenomena referred to as bubble-induced turbulence
(BIT). This BIT in turn influences the clusters and understanding this is
important both for its own sake, and also due to its impact on many aspects of
bubble motion, including their rise velocities, collisions, dispersion and the
intensity of the BIT generated (Takagi et al., 2008; Tagawa et al., 2013;
Lohse, 2018; Mathai et al., 2020).

Early numerical investigations (e.g. Smereka, 1993) on bubble clustering in
the presence of BIT assumed spherical bubbles and considered potential flow,
and found that the bubbles form large horizontal clusters as they rise.
However, subsequent experiments and more realistic numerical simulations did
not clearly observe such clusters, but 2-bubble clusters with a wide range of
orientations were identified as the most commonly occurring clusters (Zenit et
al., 2001; Bunner & Tryggvason, 2003; Esmaeeli & Tryggvason, 2005). The
prevalence of 2-bubble clusters has motivated the community to explore the
dynamics of bubble pairs with varying separations and orientations, usually
for the case of clean bubbles (Hallez & Legendre, 2011). The two extreme cases
are bubbles aligned side-by-side (Legendre et al., 2003) and in-line (Zhang et
al., 2021). The results indicate that the side-by-side configuration is more
stable than the in-line case, and this is because for the in-line
configuration a slight transverse movement of the trailing bubble relative to
the leading bubble makes it ‘feel’ a shear flow, that can drive the trailing
bubble out of the leading bubbles wake. Nevertheless, stable in-line bubble
chains are often observed in carbonated drinks. This contradiction was
recently explained by Atasi et al. (2023) as being due to the combined effects
of bubble deformation and contamination in such liquids that can result in a
reversal of the lift force and stable chain. In addition to exploring the
stability of nearby bubble pairs, the impact of neighbouring bubbles on their
rise velocity has been investigated and compared to the case of isolated
bubbles. Hallez & Legendre (2011) showed that the side-by-side configuration
maximizes the drag force acting on a pair of bubbles, while the in-line bubble
configuration minimizes the drag due to wake entrainment for the trailing
bubble.

These studies on bubble pair dynamics have provided much insight, however,
there are many open questions concerning the behaviour of bubble swarms where
two or more bubbles may be clustered together, whose motion may also be
affected by the wakes of other bubbles and bubble clusters in the flow.
Indeed, while the rise velocity of bubble pairs in a quiescent liquid is well
understood, its behaviour in the context of bubble swarms is debated. For
example, the experiments of Stewart (1995) and Brücker (1999) for large
deformable bubbles in a swarm found that the mean rise velocity was
considerably larger than that for a single bubble. However, this contradicts
other experimental (Ishii & Zuber, 1979) and numerical (Roghair et al., 2011)
studies for large bubbles which argue that the mean bubble rise velocity
decreases monotonically as the gas void fraction is increased.

Several fundamental questions remain mostly unexplored: what is the
probability to form clusters involving $N_{b}$ number of bubbles? What is the
lifetime of these clusters? How does the rise velocity of bubbles in a cluster
depend on $N_{b}$? How do the answers to these questions depend on
contaminants in the liquid? In this paper we explore these questions
experimentally by tracking thousands of deformable bubbles in a vertical
column, using a recently developed machine-learning algorithm to detect and
follow the evolution of bubble clusters, and explore how the bubble rise
velocities depend on $N_{b}$. We also consider the effect of surfactants to
provide a more complete picture for real systems where contaminants may cause
behavior that differs substantially from that of an idealized clean system.

##  2 Experimental method

###  2.1 Experimental set-up

The experimental apparatus is identical to that in Ma et al. (2022), and we
therefore refer the reader to that paper for additional details; here, we
summarize. The experiments were conducted in a rectangular bubble column
(depth $50\,\mathrm{mm}$ and width $112.5\,\mathrm{mm}$), with a water fill
height of 1,000 mm. Air bubbles are injected through 11 spargers which are
homogeneously distributed at the bottom of the column.

We use tap water in the present work as the base liquid and consider two
different bubble sizes by using spargers with different inner diameters. For
each bubble size, we manipulate the gas flow rate and ensure that all cases
are not in the heterogeneous regime of dispersed bubbly flows. In total, we
have six mono-dispersed cases (see supplementary movies 1-6) labelled as
SmTapLess, SmTap, SmPen+, LaTapLess, LaTap, and LaPen+ in table 1, including
some basic characteristic dimensionless numbers for the bubbles. Here, “Sm/La”
stand for smaller/larger bubbles, “Pen+” stands for corresponding cases added
by 1,000 ppm 1-Pentanol, and “Less” stands for lower gas void fraction than
**Tap/**Pen+ cases for smaller/larger bubbles, respectively. It should be
noted that the three cases with larger bubble sizes have higher gas void
fractions than the three cases with smaller bubbles. This is because in our
setup it is not possible to have the same flow rate for two different spargers
while also maintaining a homogeneous gas distribution for mono-dispersed
bubbles. Furthermore, the bubble size is slightly reduced when adding
1-Pentanol for both types of sparger. This is due to the influence of the
surfactants that reduce the surface tension and hence affect the bubble
formation at the rigid orifice.

To identify and track bubble clusters, we use planar shadow images obtained by
recording the flow with a high-speed camera and illuminating the setup with a
LED. The measurement resolution in time and space are 250 fps and 59.9 µm/Px,
respectively, with a field of view (FOV) of $90\,\mathrm{mm}\times
76\,\mathrm{mm}$. For each case, we record 1,000 sequences with each having
70-75 frames – approximately the time a bubble passing through the complete
image height.

Parameter | SmTapLess | SmTap | SmPen+ | LaTapLess | LaTap | LaPen+
---|---|---|---|---|---|---
$\alpha$ | 0.51% | 0.79% | 0.71% | 1.2% | 1.98% | 1.91%
$d_{b}\;(\mathrm{mm})$ | 3 | 3.1 | 2.7 | 4 | 4.3 | 3.8
$\chi$ | 1.9 | 1.9 | 1.2 | 1.9 | 2.0 | 1.3
$L_{b}\;(\mathrm{mm})$ | 14.1 | 12.7 | 11.4 | 14.1 | 13.0 | 11.2
$Ga$ | 512 | 538 | 437 | 788 | 879 | 730
$Eo$ | 1.29 | 1.38 | 1.05 | 2.30 | 2.66 | 2.08
$Re_{b}$ | 755 | 739 | 493 | 912 | 1022 | 782
$C_{D}$ | 0.61 | 0.70 | 1.04 | 0.98 | 0.97 | 1.14
Table 1: Selected characteristics of the six bubble swarm cases. Here,
$\alpha$ is the averaged gas void fraction, $d_{b}$ the equivalent bubble
diameter, $\chi$ the aspect ratio, $L_{b}$ the inter-bubble distance,
$\dot{V}$ the gas flow rate at injection,
$Ga\equiv\sqrt{\left|\pi_{\rho}-1\right|gd_{b}^{3}}/\nu$ the Galileo number,
$Eo\equiv\Delta\rho gd_{b}^{2}/\sigma$, the Eötvös number. The bubble Reynolds
number $Re_{b}$ and drag coefficient $C_{D}$ are based on $d_{b}$ and the
bubble to fluid relative velocity.

###  2.2 Bubble identification and tracking

In our study only one camera is used, however, as will be shown in § 2.3, we
are nevertheless able to perform quasi-3D tracking of bubbles. Independent on
the number of cameras used, the task to track bubbles in image sequences can
be done in a detect-to-track or in a track-to-detect fashion. While the former
links previously detected bubbles in each frame to form suitable tracks, the
latter uses extrapolations of already established tracks to detect bubbles in
follow-up images. We use the former detect-to-track strategy, which allows to
incorporate detections among multiple frames to establish tracks, with however
relying more strongly on an accurate detector that finds bubbles in individual
frames.

Even for low gas volume fractions, detecting bubbles in individual images is a
challenging task since bubbles can overlap in the images. Fully overlapping
bubbles cannot be detected, but partially overlapping bubbles can be dealt
with and deep-learning-based strategies for this have recently shown very
promising results (e.g. Kim & Park, 2021). In our previous work (Hessenkemper
et al., 2022), we developed such an approach that used a trained convolutional
neural network (CNN) to segment overlapping bubbles. Furthermore, the contour
of each detected bubble is reconstructed using 64 radial vectors pointing from
the segmentation centre to the boundary (figure 1a), and the radial vectors of
partly occluded bubbles are corrected using an additional multi-layer
perceptron (MLP).

Figure 1: Steps of the tracking algorithm: (a) detections represented with 64
radial vectors, (b) graph construction, (c) feature encoding with node encoder
MLP (dark blue rectangle) and edge encoder MLP (orange rectangle) together
with edge update MLP (green rectangle), (d) time-aware node update MLPs (light
blue rectangle), (e) predicted active edges. Cluster search region: view from
(f) the top, (g) the front and (h) the side. Note that (f,h) are only
schematic representations, showing possible bubble arrangements.

The subsequent tracking of multiple detected bubbles in close proximity poses
further challenges, as the tracker not only has to be robust against
inaccuracies of the detector, i.e. missing or false detections, but also has
to be able to track bubbles that are fully occluded even for multiple times
steps, while at the same time having numerous possible associations in the
near vicinity. To solve these issues, a graph-based tracking formalism is
used. Specifically, we follow the framework of Brasó & Leal-Taixé (2020),
utilizing multiple MLPs to predict valid connections of detections on graph
structured data. The four main aspects of this tracking framework are
described as follows. Details on network architectures, the created training
dataset as well as validation tests are provided in the Supplementary
Materials.

Graph construction: To track the bubbles, each sequence is modelled as a
graph, with detections (bubbles) being the nodes of the graph and possible
connections in time being the edges of the graph, i.e. a pair of detections
forward or backward in time that are possibly from the same bubble (figure
1b). The task is then to classify the edges in active and non-actives edges,
which at the end form a set of valid tracks that fulfil the so-called ‘flow
conservation constraints’– each node having an active edge to at most one node
forward in time and at most one node backward in time.

Feature encoding: For both the nodes and the edges of the graph, features are
encoded with two separate MLPs (figure 1c). The node embeddings represent the
appearance features of the detections, which are usually encoded with a CNN
(Brasó & Leal-Taixé 2020). However, monochromatic bubble images show few
distinct features, with the size and the shape of the bubble image being the
most relevant ones. Thus, we have chosen the 64 radial vectors from the bubble
detector as input for the node feature encoder, providing not only a more
accurate description of the features bubble size/shape, but also a better 2D
bubble contour in the case of overlapping bubbles due to the additional
correction of the radial locations. The edge embeddings represent tracking-
related features. For each detected pair of frames, the time increment,
relative coordinate and size are fed into the edge feature encoding MLP to
generate edge embeddings.

Message Passing Network: The core of the tracking algorithm is the Message
Passing Network (MPN) whose main purpose is to update node and edge embeddings
w.r.t. their surrounding nodes and edges in the graph, and this is done
iteratively using message passing steps. First, the edge embeddings are
updated by combining their embeddings with the embeddings of the adjacent pair
of nodes and feeding them into an edge-update MLP (figure 1c). Then, a time-
aware node update step is applied by aggregating edge embeddings of adjacent
edges, which already contain information of connected nodes due to the
previous edge update step. The time-awareness is introduced by at first
aggregating and updating separately incoming edges, i.e. connections backward
in time, and outgoing edges, i.e. connections forward in time, with individual
MLPs and then concatenating the outcome to finally update the node embeddings
with a node update MLP (figure 1d). For each iteration, information of nodes
one step further in time is passed through the network to the node/edge to be
updated. Thus, the number of iterations defines the time increment of the
information of other nodes that are supplied to the current node.

Edge classification and post-processing: After updating all node and edge
embeddings with the MPN, the edges are classified with a classifier MLP
(figure 1e). The predictions are post-processed and remaining violations of
the flow conservation constraints are treated with an exact rounding scheme
(Brasó & Leal-Taixé 2020). Lastly, missing links in the trajectories are
interpolated using bilinear interpolation and each trajectory is smoothed with
a uniform filter.

###  2.3 Identification and tracking of bubble clusters

The detection of bubble clusters at each time step follows a distance criteria
between neighbouring bubbles whose centres in the 2D image domain are below a
predefined threshold $2d_{b}$ from each individual case (figure 1g). This
value is mainly based on the work of Legendre et al. (2003) that a
considerable drag enhancement is observed for a bubble pair rising side-by-
side within this distance. Tests for different thresholds ($2d_{b}\pm
0.5d_{b}$) were conducted and the trends of the results in § 3 were found to
be insensitive to the choice of this parameter. Furthermore, since we attempt
to detect the bubble cluster in a quasi-3D manner, we keep the in-focus region
in the depth direction to also be $2d_{b}$ (figure 1f,h). To estimate this
depth distance to the centre plane we use the gray value gradient of the
detected bubbles and consider only sharp bubbles in the shallow Depth of Field
(DoF) region (see supplementary materials for more detail). In summary, we
utilize a cylindric search volume, $\pi(2d_{b})^{2}\times 2d_{b}$, for the
cluster identification, radially in the 2D image domain and linear in the
depth direction. For all the cases, the mean inter-bubble distance $L_{b}$
based on the global void fraction (table 1) is much larger than the search
radius $2d_{b}$, indicating that the bubble clusters to be discussed are
dynamically significant.

(a) (b)

Figure 2: Method of tracking clusters: (a) illustrates an example of two
clusters merging into one cluster; (b) illustrates an example of one cluster
separating into two clusters.

The cluster tracking strategy is inspired by the work of Liu et al. (2020) for
characterising the temporal evolution of inertial particle clusters in
turbulence. Considering two clusters identified in two consecutive time steps
($\Delta t=1/250s$), we take both to be successive realizations of the same
cluster if the number of bubbles they share is above a given threshold. The
shared bubbles across clusters in successive time steps are termed
connections. We consider forward-in-time and backward-in-time connections, and
apply thresholds on the fraction of connected bubbles over the total number of
bubbles in each cluster. We illustrate in figure 2(a) an example: Cluster A
(identified in time step 1) shares all its bubbles with cluster C (identified
in time step 2), while C shares 2/3 of its bubbles with A. Therefore, the
fractions of forward and backward connections between A and C are 1 and 2/3,
respectively. On the other hand, B shares 1/3 of its bubbles with C, and C
shares 1/3 of its bubbles with B. Thus, the forward and backward connections
between B and C are 1/3 and 1/3, respectively. Following Liu et al. (2020),
two clusters in consecutive time steps are identified as the same cluster when
the fractions of their backward and forward connections are both $\geq 1/2$.
In the example of figure 2(a), A and C are recognized as belonging to the same
cluster. The cluster lifetime is defined as the time elapsed between birth
(the first instance a cluster is identified) and death (the last time it is
recognized). Here, we explicitly include the lower threshold of $1/2$, since
many 2-bubble clusters appear and require an additional criterion for
tracking. In figure 2(b) we give an example where cluster A at time step 2
splits into B and C. To decide whether B or C should be regarded as the
continuation of cluster A for the purposes of tracking, we consider whether
cluster B or C persists longer into the future. In this example, while cluster
C survives until time n, B does not. Therefore, we regard C, D, and A as
belonging to the same cluster, while cluster B is considered to be a newborn
cluster at time step 2. This approach eliminates ambiguities since it ensures
that a cluster at any instant can only be associated with at most one cluster
either in the past or future.

##  3 Results

###  3.1 Probability to be clustered

Figure 3: (a) Percentage of bubbles in cluster for 3 cases with larger
bubbles. (b-g) Probability of number of bubbles within a cluster for the
different cases.

We first consider the percentage of bubbles in the flow that are clustered,
and the results in figure 3(a) show that this percentage increases in the
order of LaTapLess, LaTap to LaPen\+ for the larger bubbles. While the
increase from LaTapLess to LaTap is quite understandable due to the increase
of the gas void fraction, the result from LaPen\+ (whose $\alpha$ is sightly
less than that of LaTap) shows the surfactant promotes the formation of
clusters. We obtained similar results (not shown) for the three cases with
smaller bubbles.

In figure 3b-g we consider the probability to find a given number of bubbles
within a cluster for all cases, respectively. The results show that the
probability decreases with increasing $N_{b}$, and consistent with previous
studies, $N_{b}=2$ is the most common cluster size for all 6 cases (Zenit et
al., 2001; Bunner & Tryggvason, 2003). However, the results also show that
$N_{b}=3,4$ clusters occur with non-negligible probability, and there are even
rare events with $N_{b}=8$ clusters. The results also show that adding
contaminants decreases the probability to form $N_{b}=2$ clusters, and
increases the probability to form larger clusters, although the dependence is
not too strong. For a fixed contaminant level, increasing $\alpha$ has the
same effect.

###  3.2 Lifetime

(a) (b)

(c) (d)

Figure 4: Mean lifetime of 2-, 3-, 4- and 5-bubble clusters (a,b) and PDF of
the bubble cluster lifetime (c,d): the smaller bubble cases (a,c) and the
larger bubble cases (b,d).

We now turn to consider the mean lifetime of the clusters, $\langle
t_{life}\rangle$, as a function of $N_{b}$ (only the results for $N_{b}\leq 5$
are shown as the statistics for $N_{b}>5$ are not converged). Figure 4(a,b)
shows $\langle t_{life}\rangle$ normalized by a characteristic timescale of
BIT, $t_{BIT}\equiv d_{b}/U_{b}$, where $U_{b}$ is the mean vertical slip
velocity between the bubble and liquid at the column center. The values of
$\langle t_{life}\rangle/t_{BIT}$ are order unity, suggesting that $t_{BIT}$
is indeed a dynamically relevant timescale for the cluster lifetimes. The
results also reveal a systematic dependence on $N_{b}$ and the liquid
contamination. First, $\langle t_{life}\rangle$ decreases monotonically with
increasing $N_{b}$, such that larger bubble clusters are not only rarer (see §
3.1), but also more unstable. While this may not seem surprising, it is in
fact the opposite to what has been observed for inertial particles where the
cluster size and its lifetime are positively correlated (Liu et al., 2020).
The difference could be simply due to the fact that the most common values of
$N_{b}$ for our clusters are much smaller than those for the inertial
particles in Liu et al. (2020), and as a result relatively small changes in
the bubble configurations can result in the formation or destruction of a
given cluster. The other significant difference is that our bubbles
hydrodynamically interact, unlike the numerically simulated inertial particles
in Liu et al. (2020) where a one-way coupling assumption is used. Second,
increasing $\alpha$ not only leads to the formation of larger clusters, but
also slightly longer mean lifetimes for the clusters, although the lifetimes
for $N_{b}=2$ are the least sensitive to $\alpha$. Third, the mean lifetimes
of the bubble clusters notably increase with increasing contamination levels.
In a recent paper we showed that increased contamination leads to a reduction
of $Re_{b}$ and an increase in BIT (Ma et al., 2023). The reduction in
$Re_{b}$ causes the bubble trajectories to be less chaotic, and this may
explain why the cluster lifetimes increase with increasing contamination.

Figure 4(c,d) show the probability density functions (PDFs) for the cluster
lifetimes, which have been computed using clusters of all sizes. The general
dependence on the flow variables is similar to that observed for the mean
cluster lifetime, with the PDF tails becoming increasingly heavy in the order
TapLess, Tap and Pen\+ for both the small and large bubbles. The majority of
the bubble clusters survive for $t_{life}/t_{BIT}=O(1)$, however, there are
extreme cases where clusters survive for up to $t_{life}/t_{BIT}\approx 15$
for the Pen\+ cases. The central regions of the PDFs are well described by
stretched exponential functions with parameters that vary between the cases.

###  3.3 Mean $N_{b}$-bubble cluster rise velocity

(a) (b)

Figure 5: Mean bubble rise velocity as a function of the number of bubbles
$N_{b}$ in the cluster: (a) smaller bubbles and (b) larger bubbles. $N_{b}=1$
denotes unclustered bubbles.

We finally consider the role that clustering plays in the mean bubble rise
velocity. Figure 5 shows the mean rise velocity of bubbles in clusters $U_{b}$
conditioned on $N_{b}$, and the results for unclustered bubbles $N_{b}=1$ are
also shown for reference. Consistent with our previous results based on
averaging over all bubbles (Ma et al., 2023), the results show that for almost
all $N_{b}$, increasing the liquid contamination leads to a reduction in
$U_{b}$, due to the modification of the bubble boundary conditions. For the
larger bubbles we also observe a clear increase in $U_{b}$ with increasing
$N_{b}$, with an increase of up to $20\%$ when going from unclustered bubbles
($N_{b}=1$) to bubble pairs ($N_{b}=2$), while the increase is more moderate
when $N_{b}$ is increased beyond 2. The enhancement of $U_{b}$ when going from
$N_{b}=1$ to $N_{b}=2$ is also observed in the SmPen\+ case, with only slight
enhancements when $N_{b}$ is increased beyond 2. However, for the SmTapLess
and SmTap cases, $U_{b}$ varies only weakly with $N_{b}$, even in going from
$N_{b}=1$ to $N_{b}=2$.

What is the physical explanation for why the clustered bubbles rise faster
than unclustered bubbles? We begin by considering the case of bubble pairs
$N_{b}=2$ and plot in figure 6 the mean inclination angle $\theta$ of bubble
pair centreline with respect to the vertical direction (see sketch in the
figure). (It should be noted that since our measurements are only quasi-3D,
$\theta=0$ does not necessarily mean that the bubbles are in-line because they
may nevertheless be separated in the depth direction by up to a distance
$2d_{b}$.) For all cases an almost uniform distribution of $\theta$ is
observed, i.e. there is no preferential alignment for bubble pair. This is
consistent with visual inspection of the experimental images (see
supplementary movies) which show that the bubble pair orientations are not
persistent, but instead the bubbles continually trade places in a ‘leapfrog’
fashion. This observation was also found in many 3D experiments (Stewart,
1995; Riboux et al., 2010) and DNS for bubble swarm (Bunner & Tryggvason,
2003; Esmaeeli & Tryggvason, 2005). It is, however, strikingly different from
the behaviour observed for isolated bubble pairs where a stable configuration
is observed for two clean spherical bubbles where they rise side-by-side
(Hallez & Legendre, 2011), while for contaminated systems their stable
configuration is to be in-line due to the lift reversal experienced by the
trailing bubble (Atasi et al., 2023). One possible reason for this discrepancy
is that in our experiments the bubbles have oscillating and/or chaotic rising
paths, for which the probability that two bubbles will rise in a stable
arrangement is very low. By contrast, in Hallez & Legendre (2011) the bubbles
are fixed at various positions, and in Atasi et al. (2023) $Re_{b}$ is small
enough such that the bubbles have straight rising paths. Another reason is
that in our experiments the bubble pairs are not isolated and can experience
fluctuations and turbulence due to the wakes of other bubbles in the flow, and
this will readily suppress any preferential orientation that might have
occurred were the bubble pairs isolated.

Although the bubble pair orientation is almost random, the impact of their
interaction on $U_{b}$ will however depend on $\theta$, especially in the
present bubble regime (deformable bubbles with $Re_{b}\sim O$(100-1000)). For
example, in the side-by-side configuration the two bubbles are outside of each
others wake and the modification to the drag force on each bubble is minimal
(Kong et al., 2019). On the other hand, for the in-line configuration the
trailing bubble is sheltered by the leading bubble and the reduced pressure
behind the leading bubble causes the trailing bubble to be sucked towards it,
increasing the rise velocity of the trailing bubble while the leading bubble
is almost unaffected (Zhang et al., 2021). These effects mean that only the
rise velocity of trailing bubbles will be significantly affected by the
clustering, and hence when averaged over all orientations, the increased
vertical velocity of the trailing bubbles leads to an overall increase in
$U_{b}$. This explains the increased mean rise velocity for $N_{b}=2$ compared
to the $N_{b}=1$ results in figure 5. The increase is, however, minimal for
the cases SmTapLess and SmTap. This is most likely due to the bubble wakes
being weaker for these cases, a result of which is that the bubble interaction
and the associated effect on $U_{b}$ is also weaker.

The results in figure 5 show that $U_{b}$ further increases as $N_{b}$ is
increased beyond 2. The can be understood in terms of the enhanced opportunity
for bubbles to be sheltered by other bubbles as $N_{b}$ increases. However,
the increase is not as strong as when going from $N_{b}=1$ to $N_{b}=2$
because the greatest effect of sheltering will occur when two bubbles are in-
line; if the in-line bubbles are part of a cluster, the additional bubbles in
the cluster must be displaced in the horizontal direction due to the finite
size of the bubbles, and they will therefore be less effective in sheltering
the trailing bubble. It is interesting to note that for experiments on heavy
particles settling in a quiescent fluid, similar behaviour was also found,
with clustered particles falling faster (Huisman et al., 2016). In that case,
the enhanced settling velocity was also attributed to a sheltering effect,
i.e. reduced drag on particles that are falling in the wake of other
particles. However, in that context, the particle clusters were found to
exhibit strong alignment with the vertical direction, unlike our bubble
clusters whose orientations are almost random (at least for the $N_{b}=2$
case).

Figure 6: Orientation of bubble pair for the different cases (e.g. in-line
bubble pair for $\theta=0^{\circ}$ and side-by-side for $\theta=90^{\circ}$).

##  4 Conclusions

We conducted experiments on the temporal evolution of bubble clusters with the
aid of a new bubble tracking method for crowded swarms. Our results show that
2-bubble clusters are the most common, however, 3 and 4-bubble clusters also
often occur. The clusters persist on average for a time of order
$d_{b}/U_{b}$, although rare clusters persisting for an order of magnitude
longer are also observed. Furthermore, surfactants are observed to enhance the
cluster sizes and their lifetimes. A positive correlation between cluster size
and bubble rise speed is observed, with clustered bubbles rising up to $20\%$
faster than unclustered bubbles. Finally, while our cluster tracking method is
only quasi-3D, a fully 3D method for dense, deformable bubbles can be
developed by combining our bubble identification method with the recent
tracking algorithm of Tan et al. (2023) that currently applies to spherical
bubbles.

## References

  * Atasi et al. (2023) Atasi, O., Ravisankar, M., Legendre, D. & Zenit, R. 2023 The presence of surfactants controls the stability of bubble chains in carbonated drinks. Phys. Rev. Fluids 8, 053601\.
  * Brasó & Leal-Taixé (2020) Brasó, G. & Leal-Taixé, L. 2020 Learning a neural solver for multiple object tracking. In Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, pp. 6247–6257.
  * Brücker (1999) Brücker, C. 1999 Structure and dynamics of the wake of bubbles and its relevance for bubble interaction. Phys. Fluids 11 (7), 1781–1796.
  * Bunner & Tryggvason (2003) Bunner, B. & Tryggvason, G. 2003 Effect of bubble deformation on the stability and properties of bubbly flows. J. Fluid Mech. 495, 77–118.
  * Esmaeeli & Tryggvason (2005) Esmaeeli, A. & Tryggvason, G. 2005 A direct numerical simulation study of the buoyant rise of bubbles at o(100) Reynolds number. Phys. Fluids 17 (9), 093303.
  * Hallez & Legendre (2011) Hallez, Y. & Legendre, D. 2011 Interaction between two spherical bubbles rising in a viscous liquid. J. Fluid Mech. 673, 406–431.
  * Hessenkemper et al. (2022) Hessenkemper, H., Starke, S., Atassi, Y., Ziegenhein, T. & Lucas, D. 2022 Bubble identification from images with machine learning methods. Int. J. Multiphase Flow 155, 104169\.
  * Huisman et al. (2016) Huisman, S. G., Barois, T., Bourgoin, M., Chouippe, A., Doychev, T., Huck, P., Morales, C. E. B., Uhlmann, M. & Volk, R. 2016 Columnar structure formation of a dilute suspension of settling spherical particles in a quiescent fluid. Phys. Rev. Fluids 1 (7), 074204.
  * Ishii & Zuber (1979) Ishii, M. & Zuber, N. 1979 Drag coefficient and relative velocity in bubbly, droplet or particulate flows. AIChE J. 25 (5), 843–855.
  * Kim & Park (2021) Kim, Y. & Park, H. 2021 Deep learning-based automated and universal bubble detection and mask extraction in complex two-phase flows. Sci. Rep. 11 (1), 1–11.
  * Kong et al. (2019) Kong, G., Mirsandi, H., Buist, K.A., Peters, E.A.J.F., Baltussen, M.W. & Kuipers, J.A.M. 2019 Hydrodynamic interaction of bubbles rising side-by-side in viscous liquids. Exp. Fluids 60, 1–15.
  * Legendre et al. (2003) Legendre, D., Magnaudet, J. & Mougin, G. 2003 Hydrodynamic interactions between two spherical bubbles rising side by side in a viscous liquid. J. Fluid Mech. 497, 133–166.
  * Liu et al. (2020) Liu, Y., Shen, L., Zamansky, R. & Coletti, F. 2020 Life and death of inertial particle clusters in turbulence. J. Fluid Mech. 902, R1.
  * Lohse (2018) Lohse, D. 2018 Bubble puzzles: From fundamentals to applications. Phys. Rev. Fluids 3 (11), 110504.
  * Ma et al. (2022) Ma, T., Hessenkemper, H., Lucas, D. & Bragg, A. D. 2022 An experimental study on the multiscale properties of turbulence in bubble-laden flows. J. Fluid Mech. 936, A42.
  * Ma et al. (2023) Ma, T., Hessenkemper, H., Lucas, D. & Bragg, A. D. 2023 Effects of surfactants on bubble-induced turbulence. J. Fluid Mech. (under revision), arXiv:2212.03620.
  * Mathai et al. (2020) Mathai, V., Lohse, D. & Sun, C. 2020 Bubbly and buoyant particle–laden turbulent flows. Annu. Rev. Condens. Matter Phys. 11 11, 529–559.
  * Riboux et al. (2010) Riboux, G., Risso, F. & Legendre, D. 2010 Experimental characterization of the agitation generated by bubbles rising at high Reynolds number. J. Fluid Mech. 643, 509–539.
  * Roghair et al. (2011) Roghair, I, Lau, YM, Deen, NG, Slagter, HM, Baltussen, MW, Annaland, M Van Sint & Kuipers, JAM 2011 On the drag force of bubbles in bubble swarms at intermediate and high reynolds numbers. Chem. Eng. Sci. 66 (14), 3204–3211.
  * Smereka (1993) Smereka, P. 1993 On the motion of bubbles in a periodic box. J. Fluid Mech. 254, 79–112.
  * Stewart (1995) Stewart, C.W. 1995 Bubble interaction in low-viscosity liquids. Int. J. Multiphase Flow 21 (6), 1037–1046.
  * Tagawa et al. (2013) Tagawa, Y., Roghair, I., Prakash, V. N., van Sint Annaland, M., Kuipers, H., Sun, C. & Lohse, D. 2013 The clustering morphology of freely rising deformable bubbles. J. Fluid Mech. 721, R2.
  * Takagi et al. (2008) Takagi, S., Ogasawara, T. & Matsumoto, Y. 2008 The effects of surfactant on the multiscale structure of bubbly flows. Phil. Trans. R. Soc. A 366 (1873), 2117–2129.
  * Tan et al. (2023) Tan, S., Zhong, S. & Ni, R. 2023 3D lagrangian tracking of polydispersed bubbles at high image densities. Exp. Fluids 64 (4), 85.
  * Zenit et al. (2001) Zenit, R., Koch, D. L. & Sangani, A. S. 2001 Measurements of the average properties of a suspension of bubbles rising in a vertical channel. J. Fluid Mech. 429, 307–342.
  * Zhang et al. (2021) Zhang, J., Ni, M.-J. & Magnaudet, J. 2021 Three-dimensional dynamics of a pair of deformable bubbles rising initially in line. part 1. moderately inertial regimes. J. Fluid Mech. 920, A16.


# Real-time capable first principle based modelling of tokamak turbulent
transport

J. Citrin1,2, S. Breton2, F. Felici3, F. Imbeaux2, T. Aniel2, J.F. Artaud2, B.
Baiocchi4, C. Bourdelle2, Y. Camenen5, J. Garcia2 1FOM Institute DIFFER –
Dutch Institute for Fundamental Energy Research, P.O. box 6336, 5600 HH
Eindhoven, The Netherlands  2CEA, IRFM, F-13108 Saint Paul Lez Durance, France
3Eindhoven University of Technology, Department of Mechanical Engineering,
Control Systems Technology Group, P.O. Box 513, 5600MB, Eindhoven, The
Netherlands  4Istituto di Fisica del Plasma “P. Caldirola”, Associazione
Euratom-ENEA-CNR, Milano, Italy  5Aix-Marseille Universit´e, CNRS, PIIM UMR
7345, 13397 Marseille, France<EMAIL_ADDRESS>
###### Abstract

A real-time capable core turbulence tokamak transport model is developed. This
model is constructed from the regularized nonlinear regression of quasilinear
gyrokinetic transport code output. The regression is performed with a
multilayer perceptron neural network. The transport code input for the neural
network training set consists of five dimensions, and is limited to adiabatic
electrons. The neural network model successfully reproduces transport fluxes
predicted by the original quasilinear model, while gaining five orders of
magnitude in computation time. The model is implemented in a real-time capable
tokamak simulator, and simulates a 300 s ITER discharge in 10 s. This proof-
of-principle for regression based transport models anticipates a significant
widening of input space dimensionality and physics realism for future training
sets. This aims to provide unprecedented computational speed coupled with
first-principle based physics for real-time control and integrated modelling
applications.

Introduction.– Particle, momentum, and heat transport in the tokamak core is
dominated by turbulence driven by plasma microinstabilities [1, 2]. An
accurate predictive model for turbulent transport fluxes is thus vital for the
interpretation and optimization of present-day experiments, and extrapolation
to and control of future machines.

Direct numerical simulation with massively parallel nonlinear gyrokinetic
codes has provided tremendous insight to the underlying transport physics and
success in reproducing experimental fluxes in many regimes. However, the
computational cost – typically $10^{5}$ CPU hours for a local flux calculation
at a single radial point – precludes the routine use of such codes for
integrated tokamak transport simulations which demand $\sim 10^{3}$ flux
computations per 1 s of plasma evolution on JET scale devices.

Reduced turbulent transport models have been constructed to increase
tractability. They are based on the quasilinear approximation, which is proven
to be largely valid in the core of tokamak plasmas [3, 4, 5]. These rely on
nonlinear simulations for validating their ansatzes and normalizing factors.
These models have proven successful in reproducing experimental profiles in
many cases. Examples are TGLF [6] and QuaLiKiz [7]. A $\sim$6 orders of
magnitude speedup is gained in quasilinear calculations compared to nonlinear
simulations. However, while extremely useful, the tractability of such models
is still marginal for convenient large-scale scenario development over
discharge timescales. For example, 1 s of JET plasma evolution can take up to
10 hours with 10 processors, depending on the integrated modelling platform
used. This speed is also insufficient for applications such as trajectory
optimization, and simulations for developing real-time controllers.
Furthermore, any increase in physics fidelity in the models often results in a
trade off with further decrease in tractability.

This Letter illustrates an approach to overcome these challenges. The central
point is to relegate the expensive flux calculations to a stage precedent to
its use in a transport simulation. Instead, analytical formulae are to be used
in the simulation, based on a neural network (NN) nonlinear regression of
quasilinear fluxes previously compiled in a database. The advantage is
twofold: 1) the numerical resolution of analytical formulae is orders of
magnitude faster than original flux calculation; 2) the computation time
required for compiling the database is independent from the computation time
spent during the tokamak simulation itself, hence the training set for NN
regression can include results from more complete codes than used in
contemporary integrated transport modelling.

Neural networks have found multiple applications in tokamak research,
including: nonlinear regression for energy confinement scaling [8];
neoclassical transport [9]; rapid determination of equilibria [10], electron
temperature profiles [11], and charge exchange spectra [12]; classification of
disruption [13, 14, 15] and L-H transition onsets [16]. Most related to this
work is a regression of DIII-D heat fluxes from experimental power balance
databases [17].

Quasilinear transport model and training set.– The QuaLiKiz quasilinear
gyrokinetic transport model [7, 5, 18, 19, 20] was employed in this work.
QuaLiKiz solves a linear gyrokinetic dispersion relation for calculating
wavenumber spectra of instability growth rates and frequencies. Then,
integrating over the spectra, the transport fluxes are calculated via
quasilinear flux integrals and nonlinear saturation rules. The bulk of the
computational time is spent in the first stage, the dispersion relation
solver. QuaLiKiz has been coupled to the CRONOS [21] integrated modelling
suite, and has successfully reproduced temperature and density profiles of JET
and Tore-Supra discharges [22, 23]. Following recent upgrades [24], the
computational time for the QuaLiKiz eigenvalue solver at a single wavenumber
in QuaLiKiz is on the order of $\sim$1 s.

A database of QuaLiKiz solutions was constructed, in the ion temperature
gradient (ITG) instability regime. This instability is often the primary
driver of tokamak microturbulence. The code was run with adiabatic electrons
for simplicity, which also decreases the computational time to $\sim$300 ms.
The database covers four input parameters known to have significant impact on
ITG transport fluxes in this regime: the driving normalized logarithmic ion
temperature gradient $R/L_{Ti}$, the ion to electron temperature ratio
$T_{i}/T_{e}$, the safety-factor $q$, and the magnetic shear
$\hat{s}\equiv\frac{r}{q}\frac{dq}{dr}$. In addition, the input normalized
wavenumber $k_{\theta}\rho_{s}$was scanned, constricted to above ion-Larmor-
radius scales, where $\rho_{s}\equiv\sqrt{T_{e}m_{i}}/(Z_{i}q_{e}B)$. The
following parameters were maintained fixed: the normalized logarithmic density
gradient $R/L_{n}=3$, normalized radial location $r/a=0.5$. No Shafranov shift
was assumed in the geometry. The database consists of a dense grid of points
summarized in table 1, from which the training sets for the neural network
were sifted. The QuaLiKiz outputs we investigate are: growth rates and
frequencies, which correspond to 5D input space; ion heat flux, which
corresponds to 4D input space due to integration over wavenumbers. The
database includes cases corresponding to unstable modes, and cases where no
instabilities were found by QuaLiKiz, and the outputs are set to 0.

A regression of the ion heat flux has immediate application for transport
modelling. However, a regression of the more primitive linear output has its
own specific applications. For example, since the dispersion relation solver
is the slowest part of the code, a fast reproduction of growth rates,
frequencies and eigenfunctions would allow rapid tests of various saturation
rule formulations throughout parameter space. These saturation rules typically
evolve following continuous comparisons with nonlinear simulations and
experiments. In this sense, a database consisting of the complete outputs of
linear codes does not become obsolete, while a quasilinear flux database can.

Table 1: Summary of input parameters for the QuaLiKiz adiabatic electron ITG
database employed in this work

Parameter | Min value | Max value | No. of points
---|---|---|---
$R/L_{Ti}$ | 2 | 12 | 30
$T_{i}/T_{e}$ | 0.3 | 3 | 20
$q$ | 1 | 5 | 20
$\hat{s}$ | 0.1 | 3 | 20
$k_{\theta}\rho_{s}$ | 0.05 | 0.8 | 16
Total no. of points | 3 840 000

Neural networks.– The goal is to find analytical formulae which robustly
reproduce the various QuaLiKiz outputs. To this end, a multilayer perceptron
neural network is used, which is a nonlinear function with tunable variables
(weights and biases), with the property of universal approximation [25, 26].
For an overview, with an emphasis on applications for fusion, see Ref.[27].
Linear combinations of the inputs and biases are propagated through a series
of nonlinear transfer function vectors (named ‘hidden layers’), until
eventually linearly combined to an output layer. With two hidden layers and a
single output value (as used in this work), this is represented as:

$y=b_{3}+\sum_{i}^{N}{w^{2}_{i}}g\left(b^{2}_{i}+\sum_{j}^{M}w^{1}_{ij}g\left(b^{1}_{j}+\sum_{k}^{I}w^{in}_{jk}x_{k}\right)\right)$
(1)

Where $y$ is the output ‘neuron’ containing the output value (i.e. growth
rate, frequency, or ion heat flux), $x_{k}$ the vector of input values,
$b^{n}$ the bias vectors, $w^{in}$ the M$\times$I weight matrix connecting the
input vector to the 1st hidden layer, $w^{1}$ the N$\times$M weight matrix
connecting the two hidden layers, and $w^{2}$ the weight vector connecting the
2nd hidden layer to the output neuron. $g$ is the nonlinear transfer function,
defined as a sigmoid in this work:

$g(x)=\frac{2}{1+e^{-2x}}-1$ (2)

Following a series of optimization tests, two hidden layers, as shown in
equation 1, were employed here. The hidden layer sizes M and N were set to 40.
The input layer size, I, is 4 for ion heat fluxes, and 5 for growth rates and
frequencies.

The key stage is the determination of the optimized values of the weights and
biases. This is done by minimizing a cost function consisting of the average
squared error between the network output and known target output. This set of
target output, known as the ‘training set’, is a subset of the QuaLiKiz output
values from the database. The BFGS algorithm [28], implementing a quasi-Newton
method, was used for the weight and bias optimization. All NN training in this
work was carried out with the MATLAB neural network toolbox [29]. Following
training, the network output then emulates the original model within the
database input parameter envelope. This is validated by comparison to
validation sets sifted from the database, which are different from the
training set.

To avoid overfitting the data, regularization techniques were used in the
regression. This corresponds to adding a penalty term in the cost function
related to the sum of squares of the network weights and biases, leading to
smoother output. The use of regularization ensures that the NN response is
smooth (e.g. without strong oscillations) in sparse regions of training set
parameter space or when extrapolating beyond the training set envelope.

The analytic form of the nonlinear regression function allows for the
calculation of analytical gradients of the outputs with respect to the inputs.
This is vital for the efficient solution of fast implicit schemes in real-time
capable core transport simulators such as RAPTOR [30]. The regularization also
ensures smooth gradients throughout parameter space, important for the
stability of such implicit schemes.

Regression results.– We focus on the ion heat flux NN regression, due to its
direct relevance for transport modelling applications. Successful regressions
of the growth rates and frequencies were also obtained but for brevity not
discussed here.

To capture the instability thresholds with high fidelity, the regression was
only carried out for a training set corresponding to unstable modes. The NN
output for the stable regions in the validation set was then negative, since
the regularized network tends to smoothly extrapolate the trends observed
towards the training set envelope. For the final heat flux output, these
negative values were then set to zero to represent stability. This scheme
avoids having the regularized regression network attempt to directly fit the
discontinuous gradients at the instability thresholds, which would be
performed poorly due to to the regularization constraint. This is an important
point since tokamak transport often tends to be maintained near the critical
temperature gradient thresholds, especially in high temperature regimes.

The network was trained with a training set of 35000 points chosen randomly
from the set of unstable modes in the database. A comparison between the
regression NN and QuaLiKiz outputs for a validation set of 10000 unstable
cases (different from the training set) is shown in figure 1. The regression
network has an RMS error of 0.77 in gyroBohm units
($\chi_{iGB}=\frac{T_{i}^{3/2}m_{i}^{1/2}}{(Z_{i}q_{e}B)^{2}r}$) when compared
to the validation set. This RMS error is similar for the training set itself
and is primarily due to the regularization constraint. The impact of this
error on the simulated profiles is minor. This is due to stiffness, defined
here as the local gradient of $\chi_{i}$ with respect to the driving
$R/L_{Ti}$. To quantify this, a comparison was made between the $R/L_{Ti}$
values predicted by the NN and QuaLiKiz to balance a representative
$\chi_{iGB}=1$. This was done for all values of $q$, $\hat{s}$, and
$T_{i}/T_{e}$ in the database. The RMS error was
${\Delta}\left(R/L_{Ti}\right)=0.29$, which corresponded to an average
relative error of only 4.2%.

Figure 1: Comparison between normalized ion heat fluxes obtained directly from
QuaLiKiz (x-axis) and those from its NN regression (y-axis)

The typical quality of the fits can be seen in figure 2, displaying scans of
the 4 separate input parameters while the others remained fixed. Negative
outputs of the NN network are set to zero. Note the resulting excellent fit of
the instability thresholds. In addition, extrapolating the NN scans beyond the
range of the training set maintains the trend observed in the data, due to the
regularization. This is very encouraging with regard to extension of this
approach to more sparse datasets in higher dimensions. However, we do not
intend to routinely use NN models in poorly represented regions of parameter
space, as the quality of extrapolation cannot be determined a priori. Rather,
the training sets should be continuously expanded to cover such encountered
sparse or empty regions, and the NN then periodically retrained. Nevertheless,
the smoothness of the regularized NN response when extrapolating ensures its
robustness and stability during practical use as a transport model, including
during phases when such sparse regions are encountered.”

Figure 2: Comparison of NN parameter scans (blue solid lines) vs the original
QuaLiKiz ion heat flux calculations (red dots). The scans are in $R/L_{Ti}$
(top left panel), $T_{i}/T_{e}$ (top right panel), $q$ (bottom left panel) and
$\hat{s}$ (bottom right panel)

Each NN output is calculated on a sub 10 ${\mu}s$ timescale in MATLAB on a
Intel(R) Xeon(R) E5450 CPU @ 3.00GHz. This is a 5 order of magnitude speedup
in comparison to the original QuaLiKiz calculations.

Application in transport codes.–A transport model based on the trained neural
network was constructed, and implemented both in the CRONOS [21] and RAPTOR
[30] integrated modelling codes.

In CRONOS, the validity of the NN transport model was assured by a successful
comparison with a JET baseline H-mode shot 73342 with ion and electron heat
transport previously simulated [23] with the full QuaLiKiz model. For brevity
we will not focus further on this case. Rather, we focus on the real-time
simulation capabilities offered by coupling the NN model to RAPTOR.

Presently, RAPTOR only models electron heat transport. The NN model output was
thus modified to roughly approximate ITG regime electron heat transport with
kinetic electrons. This was done by assuming that heat fluxes in ITG kinetic
electron cases are higher by factor 3 compared with adiabatic electron cases,
and furthermore assuming an ion to electron heat flux ratio of
$q_{i}/q_{e}=3$. These approximations are based on typical nonlinear and
quasilinear observations in the ITG regime [31, 5].

In figure 3, we compare a RAPTOR simulation of an ITER hybrid scenario, using
the QuaLiKiz NN model for electron heat transport, with a simulation of the
same case originally carried out [32] using CRONOS and the GLF23 [33]
transport model. Using GLF23 allows to compare over ITER-scale discharge times
of $>$100s, which is less tractable using the original QuaLiKiz model. For
heat transport in a pure ITG regime, GLF23 and QuaLiKiz predictions are
expected to be similar, as illustrated in specific single-time-slice
comparisons [23].

The RAPTOR simulation uses all the same actuator (source) inputs and density
evolution as the CRONOS simulation. Ion temperatures were held fixed at
$T_{i}/T_{e}\sim 0.8$ in L-mode and $T_{i}/T_{e}\sim 0.9$ in H-mode. The NN
model was operational within a normalized toroidal flux coordinate ($\rho$)
range of 0.25 to 0.95. For $\rho>0.95$, $\chi_{e}$ was feedback controlled to
maintain a prescribed edge pedestal temperature of 4 keV. For $\rho<0.25$, a
constant $\chi_{e}$ was assumed to maintain a reasonable level of transport,
since GLF23 and QuaLiKiz both predicted stability within that region. A RAPTOR
simulation of an entire 300 s ITER discharge took 10 s on a single CPU,
corresponding to 30x faster than real-time. This combination of simulation
speed and first-principle modelling is unprecedented. With CRONOS/GLF23, the
same simulation took 24 hours.

Figure 3: Comparison between the $T_{e}$ predictions for an ITER hybrid
discharge carried out with CRONOS/GLF23 [32] (red curve) and a RAPTOR
simulation using the QuaLiKiz NN transport model (blue curve). A typical
H-mode profile (left panel) and time dependence at mid-radius (right panel)
are shown. The LH transition was set at 100 s.

Conclusions and outlook.– A neural network fit to a restricted subspace of
quasilinear gyrokinetic transport model calculations, relevant in the ITG
regime, was carried out and applied as a transport model for integrated
modelling. While the quasilinear model, QuaLiKiz, is 6 orders of magnitude
faster than nonlinear simulations, the NN regression leads to a further 5
order of magnitude speedup. This model is thus real-time capable while still
being based on first-principles, which is unprecedented. This model has been
coupled to the CRONOS integrating modelling suite, and validated against a
full QuaLiKiz simulation in the ITG regime. The model is also coupled to the
real-time capable RAPTOR tokamak simulation code, and can model a 300s ITER
discharge within 10s, with good agreement with previous modelling using CRONOS
and the GLF23 transport model.

This opens up many new possibilities for real-time controller design and
validation, scenario preparation and optimization, and real-time discharge
supervision. Such models can be used to design controllers for the plasma
profiles using model-based controller design methods (e.g. [34] or [35]). The
transport model can be used in closed-loop simulations to validate the
designed controllers. Recent work on plasma ramp-up trajectory optimization
[36] was carried out with an ad-hoc transport model, and can now be improved
using this first-principle-based transport model. Also, this transport model
can be used in real-time simulations to verify the measured plasma evolution
and warn a supervisory control system of any unexpected deviations during the
discharge [37]. Specifically for ITER, the faster-than-real-time opens up the
possibility of (on-line) real-time optimization of the discharge evolution in
response to such unexpected events.

While applications in the ITG transport regime are already feasible with this
model, there remains much scope for expanding the number of input dimensions
in the databases used for the fits, as well as employing slower yet more
complete linear gyrokinetic codes for populating the database. Neural network
topology complexity favourably scales linearly with input dimensionality.
However, we estimate that uniform density population of the input dimensions,
as carried out in this work, is feasible up to N$\sim$10\. This is due to
constraints on the NN training and quasilinear database calculation times. For
higher dimensionalities, a training set which captures the natural
correlations of the input parameters is then vital. This can be done by basing
the training sets on experimental parameters and reasonable extrapolations
thereof. This is a feasible goal, as also evidenced in Ref. [17], and this
work is ongoing.

Acknowledgments.– This work is part of the research programme ‘Fellowships for
Young Energy Scientists’ (YES!) of the Foundation for Fundamental Research on
Matter (FOM), which is financially supported by the Netherlands Organisation
for Scientific Research (NWO). The authors greatly thank O. Meneghini, S.
Smith, U. von Toussaint and P. Xanthopolous for inspiring discussions.

## References

## References

  * [1] W. Horton 1999 Rev. Mod. Phys. 71 735\.
  * [2] E.J. Doyle et al. Progress in the ITER Physics Basis Chapter 2: Plasma confinement and transport 2007 Nucl. Fusion 47 S18.
  * [3] Z. Lin et al. 2007 Phys. Rev. Lett. 99 265003\.
  * [4] T. Dannert and Frank Jenko 2005 Phys. Plasmas 12 072309\.
  * [5] A. Casati et al. 2009 Nucl. Fusion 49 085012\.
  * [6] G.M Staebler, J.E Kinsey and R.E Waltz 2007 Phys. Plasmas 14 055909\.
  * [7] C. Bourdelle et al. 2007 Phys. Plasmas 14 112501\.
  * [8] L. Allen and C.M. Bishop 1992 Plasma Phys. Control. Fusion 34 7 1291\.
  * [9] A. Wakasa et al. 2007 Jpn. J. Appl. Phys 46 3A 1157.
  * [10] J.B. Lister and H. Schnurrenberger 1991 Nucl. Fusion 31 129\.
  * [11] D.J. Clayton et al. 2013 Plasma Phys. Control. Fusion 55 095015\.
  * [12] J. Svensson, M. von Hellermann and R.W.T. König 1999 Plasma Phys. Control. Fusion 41 2 315.
  * [13] D. Wroblewski, G.L. Jahns and J.A. Leuer 1997 Nucl. Fusion 37 725\.
  * [14] B. Cannas, A. Fanni, E. Marongiu and P. Sonato 2004 Nucl. Fusion 44 68\.
  * [15] J. Vega et al. 2014 Nucl. Fusion 54 123001\.
  * [16] P Gaudio, A Murari, M Gelfusa, I Lupelli and J Vega 2014 Plasma Phys. Control. Fusion 56 114002\.
  * [17] O. Meneghini, C. J. Luna, S. P. Smith, and L. L. Lao 2014 Phys. Plasmas 21 060702\.
  * [18] J. Citrin et al. 2012 Phys. Plasmas 19 062305\.
  * [19] P. Cottier et al. 2014 Plasma Phys. Control. Fusion 56 015011\.
  * [20] C. Bourdelle. Turbulent Transport in Tokamak Plasmas: bridging theory and experiment. Physics. Aix Marseille Université, 2015. ¡tel-01113299¿.
  * [21] J.F. Artaud et al. 2010 Nucl. Fusion 50 043001\.
  * [22] A. Casati, Ph.D. dessertation, Universite de Provence (Aix-Marseille I), 2009.
  * [23] B. Baiocchi et al. 2015 Plasma Phys. Control. Fusion 57 035003\.
  * [24] J. Citrin et al. ITPA Transport & Confinement Topical Group Meeting, ITER IO, October 2014.
  * [25] Christopher M Bishop et al.  Neural networks for pattern recognition.  Clarendon press Oxford, 1995.
  * [26] Simon Haykin.  Neural Networks: A Comprehensive Foundation.  Prentice Hall, 1998.
  * [27] C. Bishop Rev. Sci. Instrum. 1994 65 1803\.
  * [28] John E Dennis Jr and Robert B Schnabel.  Numerical methods for unconstrained optimization and nonlinear equations, volume 16.  Siam, 1996.
  * [29] MATLAB and Neural Network Toolbox Release 2012b, The MathWorks, Inc., Natick, Massachusetts, United States.
  * [30] F. Felici and O. Sauter 2012 Plasma Phys. Control. Fusion 54 025002\.
  * [31] J.E Kinsey, R.E. Waltz and J. Candy 2005 Phys. Plasmas 12 062302\.
  * [32] J. Citrin et al., 2010 Nucl. Fusion 50 115007\.
  * [33] R.E. Waltz, G.M. Staebler, W. Dorland and G.W. Hammett 1997 Phys. Plasmas 4 2482\.
  * [34] E. Maljaars et al. 2015 Nucl. Fusion 55 023001\.
  * [35] D. Moreau et al. 2013 Nucl. Fusion 53 063020\.
  * [36] J. van Dongen, F. Felici, G.M.D Hogeweij, P. Geelen and E. Maljaars 2014 Plasma Phys. Control. Fusion 56 125008\.
  * [37] F. Felici et al. Proc. 24th Int. Conf. on Fusion Energy (San Diego, CA, 2012) EX/P3-12 www-web.iaea.org/napc/physics/FEC/FEC2012/papers/169_EXP312.pdf.