# Upper bound on the Guessing probability using Machine Learning

Sarnava Datta<EMAIL_ADDRESS>Hermann Kampermann  Dagmar Bruß  Institut
für Theoretische Physik III
Heinrich-Heine-Universität Düsseldorf

###### Abstract

The estimation of the guessing probability has paramount importance in quantum
cryptographic processes. It can also be used as a witness for nonlocal
correlations. In most of the studied scenarios, estimating the guessing
probability amounts to solving a semi-definite programme, for which potent
algorithms exist. However, the size of those programs grows exponentially with
the system size, becoming infeasible even for small numbers of inputs and
outputs. We have implemented deep learning approaches for some relevant Bell
scenarios to confront this problem. Our results show the capabilities of
machine learning for estimating the guessing probability and for understanding
nonlocality.

††preprint: APS/123-QED

##  I Introduction

Whenever the statistics of a measurement on a composite quantum state
contradict the assumptions of local realism, thus violating a Bell-type
inequality, the correlations are referred to as nonlocal [1]. These nonlocal
correlations are used to certify private randomness in device-independent
quantum key distribution (DIQKD) [2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13] and
device-independent randomness generation (DIRNG) [14, 15, 16, 17, 18, 19, 20,
21, 22]. For quantifying randomness, estimating the guessing probability is
often an important task. The guessing probability is the probability with
which an adversary can guess an outcome of another party’s measurement. If the
guessing probability is less than 1, the adversary cannot predict the outcome
with certainty. This implies the presence of intrinsic randomness in the
system. However, bounding the guessing probability is not an easy task.
Typically it is not possible to explicitly compute the guessing probability,
but one can only provide an upper bound by solving a semi-definite
optimization problem. Usually, one bounds the guessing probability from a
given Bell inequality, and the corresponding quantum violation [7, 14]. Here,
one needs to use the hierarchical structure of the quantum correlations [23,
24] to solve the semi-definite optimization problem. The complexity of this
optimization problem is increasing and becoming computationally demanding with
the number of settings and outcomes.

In this paper, motivated by the outstanding recent progress in utilizing
machine learning in the field of quantum information [25, 26, 27, 28, 29, 30,
31, 32, 33, 34], we develop deep learning (DL) models that predict the
guessing probability along with the optimal Bell inequalities (used to upper
bound the guessing probability) from an observed probability distribution
using supervised machine learning. A crucial element of supervised machine
learning is to generate sample data input and output to train the model. Here,
we sample random quantum probability distributions and use them as the input
of the training data. With this data, using the two-step method of Ref. [35],
we estimate the upper bound of the guessing probability and the optimal Bell
inequality, and use it as the output of the training data. After sufficient
training, our DL approach can recognize the pattern and predict the guessing
probability and the optimal Bell inequality with high accuracy and low average
statistical error.

We organize this work as follows. We start in Sec. II by explaining the
generalized Bell set-up, types of correlations and Bell inequalities. We
introduce the guessing probability and show how to estimate it by solving a
semi-definite programme in Sec. III. We introduce our deep learning approach
in Sec. IV. We discuss how to sample quantum probability distributions from
the quantum correlation space, which are then used as input for supervised
learning. We build several deep learning models for predicting the guessing
probability and the Bell inequality for various Bell scenarios and measure
their efficiency to show the model’s utility.

##  II Generalized Bell Set-up

In this section, we introduce a generalized Bell set-up. In each measurement
round, two parties, Alice and Bob, share a quantum state $\rho_{AB}$ acting on
$\mathcal{H}_{A}\otimes\mathcal{H}_{B}$. In the presence of an eavesdropper
Eve, her side information $E$ is described via the purification of the joint
system $\rho_{ABE}$ acting on
$\mathcal{H}_{A}\otimes\mathcal{H}_{B}\otimes\mathcal{H}_{E}$ where
$\mathrm{Tr}_{E}\left(\rho_{ABE}\right)=\rho_{AB}$. Each party selects locally
an input (a measurement setting) which produces an output (a measurement
outcome). We refer to this scenario as a Bell scenario.

Figure 1: Schematic description of a Bell scenario consisting of two parties,
Alice and Bob. For further explanation, see the main text.

Alice performs measurements specified by her input $x\in X=\\{1,\cdots,m\\}$,
where each input has $k$ possible outcomes $a\in A=\\{1,\cdots,k\\}$.
Similarly, Bob performs measurements specified by his input $y\in
Y=\\{1,\cdots,m\\}$ and produces the outputs $b\in B=\\{1,\cdots,k\\}$. We
denote this scenario as $[m,k]$ Bell scenario, i.e. $m$ measurement settings
with $k$ outcomes each; see Fig. 1 for visualization. After many repetitions,
the conditional probability $P(ab|xy)$ can be estimated. The Bell scenario is
completely characterized by the set
$\textbf{P}:=\\{P(ab|xy)\\}\subset\mathds{R}^{m^{2}k^{2}}$ of all joint
conditional probabilities which we refer to as a behavior [36]. Thus, the
following constraints are imposed: positivity $P(ab|xy)\geq 0$ $\forall$
$a,b,x,y$ and the normalization $\sum_{a,b=1}^{k}P(ab\arrowvert xy)=1$ for all
$x$ and $y$. We say the behavior is no-signaling if the input-output
correlation obeys

$\displaystyle\sum_{b=1}^{k}P(ab|xy)$ $\displaystyle=P(a|x)\qquad\forall
a,x,y\quad\textup{and}$ (1) $\displaystyle\sum_{a=1}^{k}P(ab|xy)$
$\displaystyle=P(b|y)\qquad\forall b,x,y\,.$

The set of all correlations satisfying the no-signaling constraints forms a
convex polytope $\mathcal{NS}$. A behavior is said to be local if it can be
written as a convex mixture of deterministic strategies [37, 38]. The set of
all local correlations forms a convex polytope $\mathcal{P}$. There exist
inequalities of the form [36]

$\sum_{a,b,x,y}C_{abxy}P(ab|xy)\leq\mathcal{I}_{L}\,,$ (2)

which separate the set of all local correlations (in other words, the convex
polytope $\mathcal{P}$) from the nonlocal behaviors. These inequalities are
called Bell inequalities. A Bell inequality is specified by the coefficients
$C_{abxy}\in\mathds{R}$. We denote a Bell inequality as $B$, and
$\sum_{a,b,x,y}C_{abxy}P(ab|xy)$ as the Bell value $B[\textbf{P}]$ in this
paper. Here, $\mathcal{I}_{L}$ is the classical bound, which is the maximal
value over all local behaviors. Thus, a behavior with a classical origin, i.e.
$\\{P(ab|xy)\\}\in\mathcal{P}$, cannot violate this inequality.

The Born rule of quantum theory postulates that a behavior is quantum if there
exists a quantum state $\rho_{AB}$ acting on a joint Hilbert space
$\mathcal{H}_{A}\otimes\mathcal{H}_{B}$ of arbitrary dimension and measurement
operators (POVM elements) $\\{M_{a|x}\\}$ with $M_{a|x}\geq 0$ and
$\sum_{a}M_{a|x}=\mathds{1}$ $\forall x$, and $\\{M_{b|y}\\}$ with analogous
properties such that

$P(ab|xy)=\mathrm{Tr}(\rho_{AB}M_{a|x}\otimes M_{b|y})\,.$ (3)

The set of all quantum correlations forms a convex set $\mathcal{Q}$. If a
behavior $\\{P(ab|xy)\\}\in\mathcal{Q}\setminus\mathcal{P}$, it violates at
least one Bell inequality of the form in Eq. (2).

Figure 2: A pictorial representation for the set of correlations. All
classical probabilities form a convex polytope $\mathcal{P}$, which is
embedded in the set $\mathcal{Q}$ of quantum correlations, which in turn is a
subset of the no-signaling polytope $\mathcal{NS}$. $v_{1}$ and $v_{2}$ are
vertices of the local polytope. $B$ (blue dashed line) represents the Bell
inequality which separates the classical polytope from the quantum and no-
signaling set.

The sets $\mathcal{P}$, $\mathcal{Q}$ and $\mathcal{NS}$ obey the following
relation: $\mathcal{P}\subsetneq\mathcal{Q}\subsetneq\mathcal{NS}$; see Fig. 2
for a pictorial representation.

##  III Guessing probability

In an adversarial black box scenario framework, the adversary Eve tries to
guess some outcomes obtained by Alice and Bob. The probability that Eve can
correctly guess the outcome is called the guessing probability. Here, we
denote the guessing probability as $P_{g}(a|x,E)$, which is the guessing
probability of Eve about Alice’s outcome $a$ corresponding to her measurement
setting $x$. In Ref.[7], it is shown that $P_{g}(a|x,E)$ can be upper bounded
by a function $G_{x}$ of the observed Bell value $B[\textbf{P}]$ of a
particular Bell inequality $B$ by semi-definite programming, i.e.
$P_{g}(a|x,E)\leq G_{x}(B[\textbf{P}])$. One crucial element to bound the
guessing probability $P_{g}(a|x,E)$ is to choose a suitable Bell inequality.
We follow the two-step procedure of [35] where the Bell inequality is
constructed from the input-output probability distribution P that leads to the
maximum Bell violation for that particular measurement statistics.

This is achieved by solving the linear program:

$\displaystyle\underset{\textbf{h},c}{\text{max}}$
$\displaystyle\textbf{h}^{T}\textbf{P}-c\,,$ (4) subject to
$\displaystyle\textbf{{h}}^{T}\textbf{v}_{p}\leq c\quad\forall\quad
p\in\\{1,\cdots,k^{2m}\\}\,,$ $\displaystyle\textbf{h}^{T}\textbf{P}>c\,,$
$\displaystyle-1\leq h_{i}\leq 1\quad\forall\quad
i\in\\{1,\cdots,m^{2}k^{2}\\}\,.$

Here h is the hyperplane specifying the Bell inequality $B$, P denotes the
measurement data, $\textbf{v}_{p}$ corresponds to the $p$ vertices of the
classical polytope $\mathcal{P}$ and $c$ is the classical bound. Thus the Bell
inequality $B$ found by the optimization of Eq. (4) and specified by the
hyperplane vector h is given as:

$\sum_{a,b,x,y}h_{abxy}P(ab|xy)\leq c\,,$ (5)

where $a\in A$, $b\in B$, $x\in X$, $y\in Y$. We will use the Bell inequality
$B$ and corresponding Bell value
$B[\textbf{P}]=\sum_{a,b,x,y}h_{abxy}P(ab|xy)$ to upper bound the guessing
probability $P_{g}(a|x,E)$ by solving the following semidefinite program [7]:

$\displaystyle\underset{\rho_{AB},\\{A(a|x)\\},\\{B(b|y)\\}}{\text{max}}\quad
P_{g}(a|x,E)$ (6) $\displaystyle\text{subject
to:}\quad\mathrm{Tr}(\rho_{AB}\mathcal{G})=B[\textbf{P}]\,.$

In the optimization problem of Eq. (6), the guessing probability is bounded
using the NPA-hierarchy [23, 24] up to level 2. The optimization is performed
using standard tools YALMIP [39], CVX [40, 41] and QETLAB [42]. Note that,
$A(a|x)$ and $B(b|y)$ are the measurement operators of Alice and Bob,
respectively, and $\rho_{AB}$ is the state shared between them. $\mathcal{G}$
is the Bell operator defined as:

$\mathcal{G}=\sum_{a,b,x,y}h_{abxy}A(a|x)B(b|y)\,.$ (7)

Let us denote $P^{*}_{g}(a|x,E)$ as the upper bound of the guessing
probability, which is the solution to the optimization problem of Eq. (6).

##  IV Machine learning approach

Providing an upper bound for the guessing probability by solving a semi-
definite program is a computationally demanding task. It becomes arduous when
the Bell scenario raises its complexity, i.e. for an increased number of
measurement inputs and/or outputs in the Bell scenario.

Thus, in this paper, we approach solving the problem via machine learning (ML)
(see Ref. [43] for detailed discussions on the concepts of machine learning)
so that the trained model can estimate the guessing probability
$P^{*}_{g}(a|x,E)$ from the input-output probability distribution
$\\{P(ab|xy)\\}$. We are going to use the supervised learning technique. In a
supervised ML approach, the first step is generating the training points. We
use random bipartite quantum probability distributions as the supervised ML
model’s input (features), after generating them from facet Bell inequalities
using the weighted vertex sampling method [44]. Since the guessing probability
for local behaviors is always 1 (i.e. Eve can guess the right outcome with
probability 1), we do not need to train the machine to perform well on those.
Thus we only take samples from the nonlocal part of the no-signaling set, i.e.
$\mathcal{NS}\setminus\mathcal{P}$. To single out the input-output correlation
with a quantum realization, we reduce the samples using the NPA hierarchy to
approximate the quantum realizable probability distribution.

Explicitly, we generate samples from the quantum set $\mathcal{Q}$ as follows.
For the $[m,k]$ Bell scenario (i.e. $m$ measurements, $k$ outcomes each), the
classical polytope $\mathcal{P}$ is specified by $k^{2m}$ local vertices. The
classical polytope can also be described by its facets, which represent the
hyperplanes (or Bell inequalities) that separate any non-classical (quantum
and no-signaling) behavior from the classical ones. These facets are called
facet Bell inequalities or tight Bell inequalities [36]; see Fig. 2 for a
pictorial representation. For the $[2,2]$ scenario, eight facet Bell
inequalities exist, all equivalent to the CHSH inequality [45]. For the
$[3,2]$ Bell scenario, there are 648 facet Bell inequalities. These facet Bell
inequalities are found using the formulation of Ref. [46] 111Using Ref. [46],
one can calculate all the facets of a convex polytope given its vertices. The
transformation from the vertex representation to the facet representation of a
polytope is known as facet enumeration or convex hull problem, which uses
Gaussian and Fourier-Motzkin elimination. The list of facets consists of
positivity constraints and the facet Bell inequalities. Here, we only focus on
the facet Bell inequalities alone.. Note that all the 648 facet Bell
inequalities correspond to two classes of independent facet Bell inequalities,
i.e. the CHSH inequalities and the $I3322$ inequalities [47, 48]. We consider
all facet Bell inequalities for the [2,2] and [3,2] Bell scenario while
generating training points for the supervised machine learning problem. For
the $[4,2]$ Bell scenario, there are 174 independent facet Bell inequalities
[49]. Since there are many (>10000) equivalent facets [50], we will only
consider the independent ones. These facet Bell inequalities are spanned by
some of the local vertices of the classical polytope 222For the $[2,2]$ Bell
scenario, all the facet inequalities are spanned by eight local vertices. For
the $[3,2]$ Bell scenario, facet Bell inequalities, equivalent to the CHSH
inequality, are spanned by thirty-two vertices. Twenty vertices span the
inequalities equivalent to the $I3322$ inequality.. These vertices provide the
maximum classical bound of the corresponding facet Bell inequality. Consider
the case that $n$ local vertices span a facet Bell inequality, where we denote
the set of $n$ vertices as
$\left\\{P_{i}^{\mathcal{L}}(ab|xy)\right\\}_{i=1}^{n}$. We denote the PR-box
of the corresponding facet Bell inequality as $P^{\mathrm{PR}}(ab|xy)$, see
Fig. 2 for visualization. The PR-Box $P^{\text{PR}}(ab|xy)$ can be defined as
the probability distribution that provides the maximal no-signaling bound of
the corresponding facet Bell inequality [51, 52]. We take uniformly random
weighted mixtures of the $n+1$ vertices ($n$ vertices that span the facet Bell
inequality and the corresponding PR-box) with an $n$-fold weight on the PR-
box. Formally, the sample behavior from the set
$\mathcal{NS}\setminus\mathcal{P}$ can be generated as:

$\textbf{P}:=P(ab|xy)=\frac{nw_{0}P^{\mathrm{PR}}(ab|xy)+\sum_{i=1}^{n}w_{i}P_{i}^{\mathcal{L}}(ab|xy)}{nw_{0}+\sum_{i=1}^{n}w_{i}}$
(8)

where the $w_{i}\in[0,1]$ with $i=0,1,\cdots,n$, are uniformly drawn random
numbers. This process is done for all facet inequalities. From this set of
samples, we only select the ones with a $\mathcal{Q}_{2}$ realization (the
second level of NPA hierarchy [23, 24]). Here we work under the assumption
that $\mathcal{Q}_{2}$ provides a good approximation for the original quantum
set $\mathcal{Q}$.

We store the probability distribution {$P(ab|xy)$} and use it as the input
(features) of the supervised machine learning problem, i.e.

$\textbf{X}:=\\{P(ab\arrowvert xy)\\}^{a,b=1,\cdots,k}_{x,y=1,\cdots,m}\,.$
(9)

We calculate the guessing probability of each input P using the two-step
method (see Sec. III for details), and use it as the output (target), i.e.

$y=P^{*}_{g}(a|x,E)\,.$ (10)

Without loss of generality, we have always calculated the guessing probability
of Alice’s first measurement setting. We use a deep neural network to assess
the dataset and make predictions. We fed the input-output pair
$\\{\textbf{X},y\\}$ (see Eq. (9) and Eq. (10)) into an artificial neural
network (ANN) to learn the best possible fit. For an elaborate explanation of
an artificial neural network, see Ref. [43]. Following the standard approach,
we divide the dataset into two parts. The first part of the dataset is for
training and validation (80%), and the second is for testing (20%). We choose
a ’linear’ 333Here, linear means that there is no branching in the hidden
layers of the neural network architecture. ANN with several layers as our
model; see Fig. 3 for visualization.

Figure 3: Schematic description of a ’linear’ neural network. It consists of
an input layer, several hidden layers and an output layer without branching.
The hidden layers and the output layer are dense layers, meaning that the
neurons of the layer are connected to every neuron of its preceding layer.

The input layer has $m^{2}k^{2}$ neurons corresponding to the elements in
{$P(ab\arrowvert xy)$}. The output (last) layer has only one neuron since we
only have to predict one element: the guessing probability $P^{*}_{g}(a|x,E)$.
We perform 100 rounds of training using the optimizer ADAM [53], of which the
first 50 rounds have a fixed learning rate of 0.001. For the next 50 rounds,
we reduce the learning rate by 90% in every tenth round. We choose the
activation function ReLu (Rectified linear unit) 444Relu activation function:
$\phi(x)=max(0,x)$ in the input and the hidden layers while using the sigmoid
activation function555sigmoid activation function:
$\phi(x)=\frac{1}{1+e^{-x}}$ in the output layer. The ReLu activation function
introduces non-linearity and the sigmoid activation function keeps the output
within the range of [0,1]. The mean squared error (MSE)
666$\text{MSE}(y,\hat{y})=\frac{1}{N}\sum_{i=0}^{N}\left(y_{i}-\hat{y}_{i}\right)^{2}$,
where $y$ is the original output and $\hat{y}$ is the estimated output of the
model. is used as our loss function, which is minimized during the training
process.

Metrics Bell Scenario | [2,2] | [3,2] | [4,2]
---|---|---|---
MSE | 0.00007 | 0.007 | 0.009
MAE | 0.0003 | 0.01 | 0.025
Table 1: Performance measures for different Bell scenarios when estimating
$P^{*}_{g}(a|x,E)$ from the probability distribution $P(ab\arrowvert xy)$.
MAE: Mean Absolute Error (see Eq. (11)), MSE: Mean Squared Error (see Eq.
(12))

The trained model generates the predicted value of the guessing probability
$P^{\text{pred}}_{g}(a|x,E)$. To check the efficiency of our approach, we have
used the mean absolute error (MAE)

$\displaystyle\text{MAE}\left[P^{*}_{g}(a|x,E),P^{\text{pred}}_{g}(a|x,E)\right]$
(11)
$\displaystyle=\frac{1}{N_{\mathrm{test}}}\sum_{i=1}^{N_{\mathrm{test}}}\left|P^{*}_{g}(a|x,E)_{i}-P^{\text{pred}}_{g}(a|x,E)_{i}\right|\,,$

and the mean squared error (MSE)

$\displaystyle\text{MSE}\left[P^{*}_{g}(a|x,E),P^{\text{pred}}_{g}(a|x,E)\right]$
(12)
$\displaystyle=\frac{1}{N_{\mathrm{test}}}\sum_{i=1}^{N_{\mathrm{test}}}\left(P^{*}_{g}(a|x,E)_{i}-P^{\text{pred}}_{g}(a|x,E)_{i}\right)^{2}\,,$

as a performance measure. $N_{\mathrm{test}}$ is the number of data points in
the test set. We analyze the results for different bipartite Bell scenarios
and list the errors in Table 1. The average error is in the order of $10^{-4}$
to $10^{-2}$. Such high accuracy and small error without knowing the Bell
inequality are truly remarkable.

Time per sample Bell Scenario | [2,2] | [3,2] | [4,2]
---|---|---|---
Mosek | 95ms | 496ms | 1568 ms
Neural Network | 27$\mu$s | 35$\mu$s | 49$\mu$s
Table 2: Runtime per sample comparison for SDP solver Mosek and the neural
network method for estimating the guessing probability for different Bell
scenarios.

We also compare the runtime performance of the neural network model with the
frequently used SDP solver Mosek [54] (that can be used to upper bound the
guessing probability by solving the optimization problem of Eq. (6)) in Table
2. The Mosek task is generated and solved using the Ncpol2sdpa [55]. The
results are evaluated over 10000 unknown samples and performed on a personal
computer 777Specifications: Intel(R) Core(TM) i7-10510U Processor, 2.30GHz
Frequency, 16.0 GB RAM under comparable conditions. Once the neural network is
trained, we get a speed-up of $10^{3}-10^{5}$ for obtaining a prediction about
a new instance, compared to the runtime of the usual method for solving the
optimization problem; see Table 2. This follows from the fact that the number
of variables in the optimization process of Eq. (6) increases exponentially
with the number of measurement settings (or outcomes per measurement) in the
Bell scenario. Thus, it takes more computational time to perform the SDP using
a classical solver like Mosek. A trained neural network only calculates the
functional output using the optimized weights and biases. Only the neural
network size affects the computational time needed to complete the prediction
task.

However, the upper bound on the guessing probability calculated from a trained
machine learning model only provides an estimate of its real value. Thus, we
cannot use this estimate to bound the secret key rate. The predicted Bell
inequality on the other side that generates a non-zero Bell violation (for a
particular measurement statistics) can be used to bound the guessing
probability (see Sec. III for details) and the secret key rate. That’s why in
the next step, we use deep learning to predict the associated optimal Bell
inequality $B$, which is then used to upper bound the guessing probability
(see Sec. III for details). For this purpose, we again use the neural network
architecture where supervised learning is incorporated. We start by preparing
the dataset where our input features are

$\textbf{X}:=\\{P(ab\arrowvert xy)\\}^{a,b=1,\cdots,k}_{x,y=1,\cdots,m}\,.$
(13)

Note that, the input is identical to Eq. (9). The outputs are now the
coefficients of the optimal Bell inequality $B$ (specified by
$\\{h_{abxy}\\}^{a,b=1,\cdots,k}_{x,y=1,\cdots,m}$, see Eq. (4)) and the
guessing probability $P^{*}_{g}(a|x,E)$, i.e.

$\textbf{y}:=\left[\\{h_{abxy}\\}^{a,b=1,\cdots,k}_{x,y=1,\cdots,m},P^{*}_{g}(a|x,E)\right]\,.$
(14)

Here we use two types of neural network architectures. The first neural
network is a usual ’linear’ feed forward neural network (see [43] for details,
schematically represented in Fig. 3). For $[m,k]$ Bell scenario, the input
layer has $m^{2}k^{2}$ neurons (corresponds to the elements of
$\\{P(ab\arrowvert xy)\\}$). The input layer is followed by several hidden
layers. Unlike in the previous scenario, the output layer has $m^{2}k^{2}+1$
neurons in this case, where $m^{2}k^{2}$ neurons correspond to the
coefficients of the Bell inequality $h_{abxy}$, and one neuron corresponds to
the guessing probability $P^{*}_{g}(a|x,E)$. In this paper, we denote this
construction of the ’linear’ deep neural network as $\mathrm{NN}_{1}$.
Following the standard approach, we divide the dataset
$\\{\textbf{X},\textbf{y}\\}$ (see Eq. (13) and Eq. (14)) into two sets; the
first part of the dataset is for training and validation (80%), and the second
part is for testing (20%). Similar to the training of the previous network, we
perform 100 rounds (first 50 rounds with a 0.001 learning rate and then reduce
the learning rate by 90% in every tenth round) of training using the gradient
solver ADAM. Similar to the previous scenario, we use the activation function
ReLu in the input and the hidden layers. In the output layer, the linear
activation function 888Linear activation function: $\phi(x)=x$ is used for
$m^{2}k^{2}$ neurons that correspond to the optimal Bell inequality and the
sigmoid activation function is incorporated for the neuron that corresponds to
the guessing probability. As the cost function, we use the Mean Squared Error
(MSE) which is minimized during the training process.

In addition, we use another neural network architecture with two parallel sub-
models (by using branching) to interpret parts of the output that share the
same input. In this construction, the input layer has $m^{2}k^{2}$ neurons
corresponding to the elements of the probability distribution $\\{P(ab|xy)\\}$
of the $[m,k]$ Bell scenario. The input layer is followed by hidden layers
consisting of multiple neurons. Then we bifurcate one hidden layer to create
two branches. Several hidden layers then follow both branches; see Fig. 4 for
visualization. The first branch of the network is for predicting the
coefficients of the optimal Bell inequality
$\\{h_{abxy}\\}^{a,b=1,\cdots,k}_{x,y=1,\cdots,m}$ and thus has $m^{2}k^{2}$
neurons. The second branch of the network is for predicting the guessing
probability. Thus, the output layer will have only one neuron corresponding to
$P^{*}_{g}(a|x,E)$. In this paper, we refer to this neural network as
$\mathrm{NN}_{2}$ which is built using the Keras functional API [56].

Figure 4: Schematic diagram of a neural network where a hidden layer is
bifurcated into two different arms which predict different parts of the output
separately. In our scenario, input layer: $\\{P(ab|xy)\\}$ ($m^{2}k^{2}$
neurons), first output: $\\{h_{abxy}\\}$ ($m^{2}k^{2}$ neurons) and second
output: $P^{*}_{g}(a|x,E)$ (1 neuron).

In $\mathrm{NN}_{2}$, we use the ReLu activation function in the input and all
the hidden layers. The linear activation function is used in the output layer
of the first branch (which predicts the coefficients of the Bell inequality)
while the sigmoid activation function is used in the second branch (which
predicts the guessing probability). The other details of the training steps
are the same as for the $\mathrm{NN}_{1}$ neural network stated previously.
Both $\mathrm{NN}_{1}$ and $\mathrm{NN}_{2}$ predict the Bell inequality
$B^{\text{pred}}$ (specified by the predicted coefficients
$\\{h_{abxy}^{\text{pred}}\\}^{a,b=1,\cdots,k}_{x,y=1,\cdots,m}$) and the
guessing probability $P^{\text{pred}}_{g}(a|x,E)$.

Since the neural networks $\mathrm{NN}_{1}$ and $\mathrm{NN}_{2}$ predict two
separate entities (the optimal Bell inequality and the guessing probability),
we evaluate their performance separately. We use the mean absolute error (see
Eq. (11)) and mean squared error (see Eq. (12)) as our performance measure of
predicting the guessing probability. The errors for different bipartite Bell
scenarios are listed in Table 3.

ANN | Metrics | [2,2] | [3,2] | [4,2]
---|---|---|---|---
$\mathrm{NN}_{1}$ | MSE | 8.2$\times 10^{-6}$ | 0.002 | 0.009
MAE | 0.001 | 0.02 | 0.07
$\mathrm{NN}_{2}$ | MSE | 1.9$\times 10^{-7}$ | 0.001 | 0.002
MAE | 0.0003 | 0.013 | 0.027
Table 3: Statistical errors of the predicted guessing probability
$P^{\text{pred}}_{g}(a|x,E)$ with respect to the guessing probability
$P^{*}_{g}(a|x,E)$ for different Bell scenarios for $\mathrm{NN}_{1}$ and
$\mathrm{NN}_{2}$. Here the neural network is trained for predicting the
guessing probability $P^{*}_{g}(a|x,E)$ and the Bell inequality $B$ from the
probability distribution $P(ab|xy)$.

Note that, for estimating the guessing probability, $\mathrm{NN}_{2}$ yields
lower statistical errors than $\mathrm{NN}_{1}$. The reason lies in the
structure of the neural network architectures. Since we create a branch in the
neural network only to estimate the guessing probability, the
$\mathrm{NN}_{2}$ neural network assigns more nodes to only estimate the
guessing probability than $\mathrm{NN}_{1}$. In the case of predicting the
optimal Bell inequality $B$ (characterized by its coefficients
$\\{h_{abxy}\\}^{a,b=1,\cdots,k}_{x,y=1,\cdots,m}$), we use the performance
measure MSE, which reads:

$\displaystyle\text{MSE}\left[B,B^{\text{pred}}\right]$ (15)
$\displaystyle=\frac{1}{m^{2}k^{2}}\frac{1}{N_{\mathrm{test}}}\sum_{i=1}^{N_{\mathrm{test}}}\sum_{a,b=1}^{k}\sum_{x,y=1}^{m}\left(\left(h_{abxy}\right)_{i}-(h_{abxy}^{\text{pred}})_{i}\right)^{2}\,,$

and MAE, which reads:

$\displaystyle\text{MAE}\left[B,B^{\text{pred}}\right]$ (16)
$\displaystyle=\frac{1}{m^{2}k^{2}}\frac{1}{N_{\mathrm{test}}}\sum_{i=1}^{N_{\mathrm{test}}}\sum_{a,b=1}^{k}\sum_{x,y=1}^{m}\left|(h_{abxy})_{i}-(h_{abxy}^{\text{pred}})_{i}\right|\,.$

The errors are listed in Table 4.

ANN | Metrics | [2,2] | [3,2] | [4,2]
---|---|---|---|---
$\mathrm{NN}_{1}$ | MSE | 0.0004 | 0.0007 | 0.014
MAE | 0.001 | 0.002 | 0.067
$\mathrm{NN}_{2}$ | MSE | 0.0003 | 0.0005 | 0.015
MAE | 0.0009 | 0.002 | 0.069
Table 4: Statistical errors of the coefficients of the predicted Bell
inequality $B^{\text{pred}}$ (predicted by the trained deep learning models
$\mathrm{NN}_{1}$ and $\mathrm{NN}_{2}$),
$\\{h_{abxy}^{\text{pred}}\\}^{a,b=1,\cdots,k}_{x,y=1,\cdots,m}$ with respect
to the original coefficients
$\\{h_{abxy}\\}^{a,b=1,\cdots,k}_{x,y=1,\cdots,m}$ for different Bell
scenarios.

Another way to evaluate the quality of the predicted Bell inequality is to use
it for upper bounding the guessing probability problem (see Eq. (6)). First,
we estimate the probability of $P^{*}_{g}(a|x,E)<1$, where $P^{*}_{g}(a|x,E)$
is calculated from the predicted Bell inequality $B^{\text{pred}}$ and the
input-output probability distribution $\\{P(ab|xy)\\}$ of the test set. We
present the results in Table 5.

ANN | [2,2] | [3,2] | [4,2]
---|---|---|---
$\mathrm{NN}_{1}$ | 99.5% | 98.7% | 93.4%
$\mathrm{NN}_{2}$ | 99.6% | 99.4% | 94.6%
Table 5: Probability of $P_{g}(a|x,E)<1$ when using the Bell inequality
$B^{\text{pred}}$ (predicted by the trained deep learning models
$\mathrm{NN}_{1}$ and $\mathrm{NN}_{2}$).

We also look into the statistical errors between the original guessing
probability $P^{*}_{g}(a|x,E)$ from the test set and the guessing probability
calculated from the predicted Bell inequality $B^{\text{pred}}$. We use MAE
and MSE as the performance measures listed in Table 6.

ANN | Metrics | [2,2] | [3,2] | [4,2]
---|---|---|---|---
$\mathrm{NN}_{1}$ | MSE | 1.7 $\times 10^{-8}$ | 0.002 | 0.006
MAE | 6.3 $\times 10^{-5}$ | 0.014 | 0.038
$\mathrm{NN}_{2}$ | MSE | 1.5 $\times 10^{-8}$ | 0.002 | 0.004
MAE | 4.1 $\times 10^{-5}$ | 0.014 | 0.031
Table 6: Statistical errors between the guessing probability calculated from
the predicted Bell inequality and the original guessing probability from the
test set for various Bell scenarios and neural network constructions.

The high probability of generating $P^{*}_{g}(a|x,E)<1$ with the predicted
Bell inequalities (see Table 5) and the small statistical errors (see Table 6)
demonstrate the quality and accuracy of the predicted Bell inequality.

We again compare the computational runtime of predicting the optimal Bell
inequality using the standard linear optimization of Eq. (4) with the neural
network $\mathrm{NN}_{1}$ and $\mathrm{NN}_{2}$. The runtime for different
methods is shown in Table 7.

Time per sample Bell Scenario | [2,2] | [3,2] | [4,2]
---|---|---|---
LP | 95 ms | 159 ms | 379 ms
$\mathrm{NN}_{1}$ | 29 $\mu$s | 40 $\mu$s | 57 $\mu$s
$\mathrm{NN}_{2}$ | 31 $\mu$s | 41 $\mu$s | 55 $\mu$s
Table 7: Runtime per sample comparison for estimation of the optimal Bell
inequality using linear programming of Eq. (4) and the neural network method
for different Bell scenarios. LP stands for linear programming which is
performed using the Mosek solver.

Similar to the previous scenario, the runtimes are evaluated over 10000
unknown samples and performed on the same personal computer in the same
condition. The linear programming of Eq. (4) is performed with the Mosek
solver using PICOS [57] python interface. We notice a significant speed-up
when using the trained neural network models compared to the Mosek solver.
This again follows from the fact that the number of variables in the
optimization process of Eq. (4) increases with the number of measurement
settings (or outcomes per measurement) in the Bell scenario while the
computational time for the neural networks only depends on its size.

##  V Discussion & Conclusion

Estimating the guessing probability is a cornerstone for device-independent
quantum key distribution and device-independent randomness generation. This
paper introduces a novel method to estimate the guessing probability using
trained deep learning models to bypass the computationally complex and
cumbersome semi-definite optimization process. Computation with the trained
deep learning models is significantly faster than using a conventional solver.
With current technology, Bell test event rates are around 100 kHz, which
results in new data every 10$\mu$s [58]. This frequency is too high for
conventional SDP solvers on a single CPU. For those cases, our deep learning
approach improves the computation significantly. In principle, optimizing the
size of a deep neural network that can process each event as the experiment is
being conducted is possible.

The deep learning model only provides an estimation of the upper bound of the
guessing probability. But it will not provide a certification. Thus,
additionally, our DL model provides an estimation of the optimal Bell
inequality for which the Bell violation using the measurement statistics
certifies the nonlocality of input-output correlations and guarantees that the
guessing probability will be less than one. Our trained deep learning models,
which significantly speed up the prediction of the Bell inequality compared to
a conventional linear program solver, predict a Bell inequality that can
generate $P_{g}(a|x,E)<1$ with a very high probability. The mean average error
between the guessing probability calculated from the predicted Bell inequality
and the optimal Bell inequality (calculated using Eq. (4)) is in the order of
$10^{-5}-10^{-2}$ (mean squared error is in the order of $10^{-8}-10^{-3}$)
which shows the quality of this approach such that it can efficiently be used
in a DIQKD or DIRNG protocol.

We also demonstrate a method for sampling random quantum correlations
(correlations which have a realization of NPA hierarchy level of 2) using the
facet Bell inequalities, which is then used as input in the supervised machine
learning process. Note that, while generating probability distributions, we
consider all facet Bell inequalities for the $[2,2]$ and $[3,2]$ Bell
scenario. However, since there are more than 10000 facet Bell inequalities for
the $[4,2]$ Bell scenario, we only restrict ourselves to generating
probability distributions using the independent facet Bell inequalities.

To illustrate the benefits of our method, we have applied it to several
relevant Bell scenarios. Note that we design and train our neural networks to
minimize statistical errors. However, we do not claim that our choice of the
trained neural network is optimal for estimating the guessing probability and
the associated optimal Bell inequality from the measurement statistics. Other
constructions of neural networks will lead to different results.

We observed that the statistical errors in the estimation of the guessing
probability and the optimal Bell inequality increase with the complexity of
the Bell scenario (i.e., the increase in the number of measurements per
party). Since there are more inputs and outputs, our neural network
architecture might not be able to generalize the extensive system with a
limited number of hidden layers and nodes in each layer. To decrease the
errors, one can take two steps. First, one can generate a larger dataset to
train the model. Second, one can build a more extensive neural network
architecture (i.e., more hidden layers or nodes in every layer). However,
using a larger dataset for training or/and training a more extensive neural
network will result in significantly more computational time. There is also
the possibility of overfitting in an extensive network. A larger neural
network architecture will also take more time to predict new instances.
Therefore, one has to change the network architecture to optimize the speed
and precision of a specific scenario.

Note that while comparing the runtime for the Mosek optimization solver with
the trained neural network for the estimation of the guessing probability (see
Table 2), we implement the NPA hierarchy of level 2. The difference in
computational runtime between the methods will be much more pronounced with
increasing hierarchy.

Our research demonstrates the applicability of deep learning techniques for
Bell nonlocality and upper bounding the guessing probability. We believe that
this strategy will create several research lines. The logical next step is to
apply our approach to Bell scenarios with a higher number of measurement
settings and outcomes. It is also possible to expand our framework to a
multipartite scenario. Another direction worth exploring for future work is
investigating other neural network constructions. Beyond the advantage in
speed, one could use neural network architectures to search for new Bell
inequalities. Also, recall that our methodology does not account for
uncertainty or offers certification of the output. It remains for future work
to use techniques like probabilistic modeling [59] that can certify the
correctness of the model’s output.

##  VI Acknowledgement

The authors acknowledge support from the Federal Ministry of Education and
Research BMBF (Project Q.Link.X). We thank Lucas Tendick for helpful
discussions.

## References

  * [1] J. S. Bell, “On the Einstein Podolsky Rosen paradox,” Physics Physique Fizika, vol. 1, no. 3, p. 195, 1964.
  * [2] A. Acín, N. Brunner, N. Gisin, S. Massar, S. Pironio, and V. Scarani, “Device-independent security of quantum cryptography against collective attacks,” Physical Review Letters, vol. 98, no. 23, p. 230501, 2007.
  * [3] S. Pironio, A. Acin, N. Brunner, N. Gisin, S. Massar, and V. Scarani, “Device-independent quantum key distribution secure against collective attacks,” New Journal of Physics, vol. 11, no. 4, p. 045021, 2009.
  * [4] R. Arnon-Friedman, R. Renner, and T. Vidick, “Simple and tight device-independent security proofs,” SIAM Journal on Computing, vol. 48, no. 1, pp. 181–225, 2019.
  * [5] R. Arnon-Friedman, F. Dupuis, O. Fawzi, R. Renner, and T. Vidick, “Practical device-independent quantum cryptography via entropy accumulation,” Nature Communications, vol. 9, no. 1, pp. 1–11, 2018.
  * [6] J. Barrett, L. Hardy, and A. Kent, “No signaling and quantum key distribution,” Physical Review Letters, vol. 95, no. 1, p. 010503, 2005\.
  * [7] L. Masanes, S. Pironio, and A. Acín, “Secure device-independent quantum key distribution with causally independent measurement devices,” Nature Communications, vol. 2, p. 238, 2011.
  * [8] U. Vazirani and T. Vidick, “Fully device-independent quantum key distribution,” Physical Review Letters, vol. 113, no. 14, p. 140501, 2014\.
  * [9] L. Masanes, R. Renner, M. Christandl, A. Winter, and J. Barrett, “Full security of quantum key distribution from no-signaling constraints,” IEEE Transactions on Information Theory, vol. 60, no. 8, pp. 4973–4986, 2014\.
  * [10] A. Acin, S. Massar, and S. Pironio, “Efficient quantum key distribution secure against no-signalling eavesdroppers,” New Journal of Physics, vol. 8, no. 8, p. 126, 2006.
  * [11] G. Murta, S. B. van Dam, J. Ribeiro, R. Hanson, and S. Wehner, “Towards a realization of device-independent quantum key distribution,” Quantum Science and Technology, vol. 4, no. 3, p. 035011, 2019.
  * [12] T. Holz, H. Kampermann, and D. Bruß, “Genuine multipartite bell inequality for device-independent conference key agreement,” Physical Review Research, vol. 2, no. 2, p. 023251, 2020.
  * [13] E. Hänggi and R. Renner, “Device-independent quantum key distribution with commuting measurements,” arXiv preprint arXiv:1009.1833, 2010.
  * [14] S. Pironio, A. Acín, S. Massar, A. B. de La Giroday, D. N. Matsukevich, P. Maunz, S. Olmschenk, D. Hayes, L. Luo, T. A. Manning, et al., “Random numbers certified by bell’s theorem,” Nature, vol. 464, no. 7291, pp. 1021–1024, 2010.
  * [15] O. Nieto-Silleras, C. Bamps, J. Silman, and S. Pironio, “Device-independent randomness generation from several bell estimators,” New Journal of Physics, vol. 20, no. 2, p. 023049, 2018.
  * [16] S. Pironio and S. Massar, “Security of practical private randomness generation,” Physical Review A, vol. 87, no. 1, p. 012336, 2013.
  * [17] J.-D. Bancal, L. Sheridan, and V. Scarani, “More randomness from the same data,” New Journal of Physics, vol. 16, no. 3, p. 033011, 2014.
  * [18] O. Nieto-Silleras, S. Pironio, and J. Silman, “Using complete measurement statistics for optimal device-independent randomness evaluation,” New Journal of Physics, vol. 16, no. 1, p. 013035, 2014.
  * [19] F. Bischof, H. Kampermann, and D. Bruß, “Measurement-device-independent randomness generation with arbitrary quantum states,” Physical Review A, vol. 95, no. 6, p. 062305, 2017.
  * [20] A. Acín, S. Massar, and S. Pironio, “Randomness versus nonlocality and entanglement,” Physical Review Letters, vol. 108, no. 10, p. 100402, 2012\.
  * [21] A. Acín and L. Masanes, “Certified randomness in quantum physics,” Nature, vol. 540, no. 7632, pp. 213–219, 2016.
  * [22] P. Skrzypczyk and D. Cavalcanti, “Maximal randomness generation from steering inequality violations using qudits,” Physical Review Letters, vol. 120, no. 26, p. 260401, 2018.
  * [23] M. Navascués, S. Pironio, and A. Acín, “Bounding the set of quantum correlations,” Physical Review Letters, vol. 98, no. 1, p. 010401, 2007\.
  * [24] M. Navascués, S. Pironio, and A. Acín, “A convergent hierarchy of semidefinite programs characterizing the set of quantum correlations,” New Journal of Physics, vol. 10, no. 7, p. 073013, 2008.
  * [25] J. Carrasquilla and R. G. Melko, “Machine learning phases of matter,” Nature Physics, vol. 13, no. 5, pp. 431–434, 2017.
  * [26] P. Broecker, J. Carrasquilla, R. G. Melko, and S. Trebst, “Machine learning quantum phases of matter beyond the fermion sign problem,” Scientific Reports, vol. 7, no. 1, pp. 1–10, 2017.
  * [27] A. Canabarro, S. Brito, and R. Chaves, “Machine learning nonlocal correlations,” Physical Review Letters, vol. 122, no. 20, p. 200401, 2019\.
  * [28] G. Carleo and M. Troyer, “Solving the quantum many-body problem with artificial neural networks,” Science, vol. 355, no. 6325, pp. 602–606, 2017.
  * [29] D.-L. Deng, “Machine learning detection of bell nonlocality in quantum many-body systems,” Physical Review Letters, vol. 120, no. 24, p. 240402, 2018.
  * [30] X. Gao and L.-M. Duan, “Efficient representation of quantum many-body states with deep neural networks,” Nature Communications, vol. 8, no. 1, pp. 1–6, 2017.
  * [31] Y.-C. Ma and M.-H. Yung, “Transforming bell’s inequalities into state classifiers with machine learning,” NPJ Quantum Information, vol. 4, no. 1, pp. 1–10, 2018.
  * [32] P. Mehta, M. Bukov, C.-H. Wang, A. G. Day, C. Richardson, C. K. Fisher, and D. J. Schwab, “A high-bias, low-variance introduction to machine learning for physicists,” Physics Reports, 2019.
  * [33] G. Torlai, G. Mazzola, J. Carrasquilla, M. Troyer, R. Melko, and G. Carleo, “Neural-network quantum state tomography,” Nature Physics, vol. 14, no. 5, pp. 447–450, 2018.
  * [34] J. Venderley, V. Khemani, and E.-A. Kim, “Machine learning out-of-equilibrium phases of matter,” Physical Review Letters, vol. 120, no. 25, p. 257204, 2018.
  * [35] S. Datta, H. Kampermann, and D. Bruß, “Device-independent secret key rates via a postselected bell inequality,” Phys. Rev. A, vol. 105, p. 032451, Mar 2022.
  * [36] N. Brunner, D. Cavalcanti, S. Pironio, V. Scarani, and S. Wehner, “Bell nonlocality,” Reviews of Modern Physics, vol. 86, no. 2, p. 419, 2014.
  * [37] I. Pitowsky, “Resolution of the einstein-podolsky-rosen and Bell paradoxes,” Physical Review Letters, vol. 48, no. 19, p. 1299, 1982.
  * [38] I. Pitowsky, “Correlation polytopes: their geometry and complexity,” Mathematical Programming, vol. 50, no. 1-3, pp. 395–414, 1991.
  * [39] J. Löfberg, “Yalmip: A toolbox for modeling and optimization in matlab,” in In Proceedings of the CACSD Conference, (Taipei, Taiwan), 2004\.
  * [40] M. Grant and S. Boyd, “Graph implementations for nonsmooth convex programs,” in Recent Advances in Learning and Control (V. Blondel, S. Boyd, and H. Kimura, eds.), Lecture Notes in Control and Information Sciences, pp. 95–110, Springer-Verlag Limited, 2008.  http://stanford.edu/~boyd/graph_dcp.html.
  * [41] M. Grant and S. Boyd, “CVX: Matlab software for disciplined convex programming, version 2.1.” http://cvxr.com/cvx, Mar. 2014.
  * [42] N. Johnston, “QETLAB: A MATLAB toolbox for quantum entanglement, version 0.9.” http://qetlab.com, Jan. 2016.
  * [43] I. Goodfellow, Y. Bengio, and A. Courville, Deep Learning.  MIT Press, 2016.  http://www.deeplearningbook.org.
  * [44] T. Kriváchy, Y. Cai, J. Bowles, D. Cavalcanti, and N. Brunner, “High-speed batch processing of semidefinite programs with feedforward neural networks,” New Journal of Physics, vol. 23, no. 10, p. 103034, 2021.
  * [45] J. F. Clauser, M. A. Horne, A. Shimony, and R. A. Holt, “Proposed experiment to test local hidden-variable theories,” Physical Review Letters, vol. 23, no. 15, p. 880, 1969.
  * [46] K. Fukuda, “Cddlib reference manual,” Report version 093a, McGill University, Montréal, Quebec, Canada, 2003.
  * [47] D. Collins and N. Gisin, “A relevant two qubit bell inequality inequivalent to the chsh inequality,” Journal of Physics A: Mathematical and General, vol. 37, no. 5, p. 1775, 2004.
  * [48] K. F. Pál and T. Vértesi, “Maximal violation of a bipartite three-setting, two-outcome bell inequality using infinite-dimensional quantum systems,” Physical Review A, vol. 82, no. 2, p. 022116, 2010.
  * [49] E. Z. Cruzeiro and N. Gisin, “Complete list of tight bell inequalities for two parties with four binary settings,” Physical Review A, vol. 99, no. 2, p. 022104, 2019.
  * [50] T. Cope and R. Colbeck, “Bell inequalities from no-signaling distributions,” Phys. Rev. A, vol. 100, p. 022114, Aug 2019.
  * [51] J. Barrett, N. Linden, S. Massar, S. Pironio, S. Popescu, and D. Roberts, “Nonlocal correlations as an information-theoretic resource,” Physical review A, vol. 71, no. 2, p. 022101, 2005.
  * [52] S. Popescu and D. Rohrlich, “Quantum nonlocality as an axiom,” Foundations of Physics, vol. 24, no. 3, pp. 379–385, 1994.
  * [53] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” arXiv preprint arXiv:1412.6980, 2014.
  * [54] M. ApS, “Mosek optimization toolbox for matlab,” User’s Guide and Reference Manual, Version, vol. 4, 2019.
  * [55] P. Wittek, “Algorithm 950: Ncpol2sdpa—sparse semidefinite programming relaxations for polynomial optimization problems of noncommuting variables,” ACM Transactions on Mathematical Software (TOMS), vol. 41, no. 3, pp. 1–12, 2015.
  * [56] F. Chollet, “keras.” https://github.com/fchollet/keras, 2015.
  * [57] G. Sagnol and M. Stahlberg, “Picos: A python interface to conic optimization solvers,” Journal of Open Source Software, vol. 7, no. 70, p. 3915, 2022\.
  * [58] P. Bierhorst, E. Knill, S. Glancy, Y. Zhang, A. Mink, S. Jordan, A. Rommal, Y.-K. Liu, B. Christensen, S. W. Nam, et al., “Experimentally generated randomness certified by the impossibility of superluminal signals,” Nature, vol. 556, no. 7700, pp. 223–226, 2018.
  * [59] Z. Ghahramani, “Probabilistic machine learning and artificial intelligence,” Nature, vol. 521, no. 7553, pp. 452–459, 2015.