aFidelity Center of Applied Technology,
FMR LLC, Boston, MA 02210, USA
b IonQ Inc, 4505 Campus Dr, College Park, MD 20740, USA

# Generative Quantum Learning of Joint Probability Distribution Functions

Elton Yechao Zhua, Sonika Johrib, Dave Baconb, Mert Esencana, Jungsang Kimb,
Mark Muira, Nikhil Murgaia, Jason Nguyenb, Neal Pisentib, Adam Schouelaa,
Ksenia Sosnovab, Ken Wrightb

###### Abstract

Modeling joint probability distributions is an important task in a wide
variety of fields. One popular technique for this employs a family of
multivariate distributions with uniform marginals called copulas. While the
theory of modeling joint distributions via copulas is well understood, it gets
practically challenging to accurately model real data with many variables. In
this work, we design quantum machine learning algorithms to model copulas. We
show that any copula can be naturally mapped to a multipartite maximally
entangled state. A variational ansatz we christen as a ‘qopula’ creates
arbitrary correlations between variables while maintaining the copula
structure starting from a set of Bell pairs for two variables, or GHZ states
for multiple variables. As an application, we train a Quantum Generative
Adversarial Network (QGAN) and a Quantum Circuit Born Machine (QCBM) using
this variational ansatz to generate samples from joint distributions of two
variables for historical data from the stock market. We demonstrate our
generative learning algorithms on trapped ion quantum computers from IonQ for
up to 8 qubits and show that our results outperform those obtained through
equivalent classical generative learning. Further, we present theoretical
arguments for exponential advantage in our model’s expressivity over classical
models based on communication and computational complexity arguments.

##  1 Introduction

Understanding the statistical relationship between several random variables is
critical to all data-based analysis and decision-making. A few examples of its
diverse applications include risk management [1], portfolio optimization [2],
reliability analysis [3], recommender systems[4], climate research [5], and
medical imaging [6]. Traditionally, single-parameter quantities such as the
Pearson correlation or Spearman’s correlation have been used to model
dependence between variables. However, such measures are good for monotonic
dependence, which is frequently too simplistic for real data. Data, such as
that from the financial markets, engineering reliability studies,
earth/atmospheric sciences tends to exhibit tail dependence, which means they
do not appear to have much correlation but exhibit dependence in extreme
deviations, as in the case of a black swan event [7, 8, 9].

Due to the reasons above, the relationship between random variables is now
commonly modeled using a dependence function between uniformly distributed
variables, called a ‘copula’. Sklar’s theorem[10], which will be explained in
the text, states that any multivariate joint distribution can be written in
terms of univariate marginal distributions and a copula that describes the
dependence structure between the variables. This provides the theoretical
foundations for the use of copulas. Since then, copulas have found many
applications in quantitative finance [11], engineering [12] and medicine [13].

Some of the commonly used copulas include elliptical and Archimedean. However,
empirical copulas (copulas from real data) tend to be a mixture of copulas and
are commonly modeled using parametric methods like maximum likelihood
estimation [14]. As a rule, the more complex the copula and the more
completely it describes the data, the more computationally challenging it
becomes to extend it to higher dimensional data.

More recently, generative models have been proposed for statistical modeling
that learn to generate data with the same statistics as a given training
dataset, effectively learning its distribution. The model can be used to
output new samples that could plausibly have belonged to the original dataset
[15, 16]. One of the most successful of these is the Generative Adversarial
Network (GAN), in which two neural networks compete with one another in a
minimax game. In Ref. [15], the authors showed theoretically that a GAN will
learn the data distribution if given enough capacity, data and training time.
Since then, GANs have found applications in areas like data augmentation [17],
fashion design [18], and super resolution [19], etc. However, there is still
debate over whether GANs will successfully learn a given distribution [20]. In
particular, problems like vanishing or unstable gradient, mode collapse and
non-convergence, sensitivity to hyperparameters, and either the generator or
discriminator overpowering the other, can make them challenging to train [21].

In this context, the question arises whether quantum computers can provide any
advantage over classical models for generating samples that reproduce the
inter-dependence between multiple variables of a given dataset. Quantum
wavefunctions naturally represent probability distributions which can be hard
to sample from by classical algorithms. Quantum computers can also be used to
generate correlations which cannot be efficiently reproduced by classical
means [22, 23]. Therefore, intuitively, it seems that learning dependence
between multiple random variables is something at which a quantum computer may
excel. Further, an application that requires sampling from a probability
distribution is likely one of the near term applications of quantum computers
as evidenced by this already being the basis of the quantum supremacy
demonstration [24].

Quantum generative models such as the Quantum Circuit Born Machine (QCBM) [25,
26] and Quantum Generative Adversarial Network (QGAN) [27, 28] have been
proposed as quantum algorithms for the learning of both classical and quantum
data. Experimental demonstrations include training a QCBM to produce sample
from the bars and stripes dataset [29] and to approximately prepare many-body
states [26]. QGANs have been used in state preparation [30], learning of qubit
channels [31], and image generation [32]. In Ref. [33], the output from a
quantum circuit is used as the prior to a classical GAN. Since quantum
algorithms are known to achieve speedup over classical algorithms in a variety
of tasks [34, 35], it has been hypothesized that quantum generative models can
do better than classical generative models in terms of expressivity, stability
during training, and shorter training time [28].

In the existing literature, generative quantum learning has either been
restricted to learning quantum data, univariate distributions, or multivariate
distributions without separating the correlation structure from the marginal
distributions. Here, we fundamentally extend this prior work using a natural
connection between entanglement and the copula function. We prove that a
variational ansatz based on this insight leads to an exponential advantage
over some classical methods of generating multivariate correlations. Our tests
on quantum simulators and hardware show advantages for the stability of the
training process as well which is relevant for practical applications.

This work is especially relevant as quantum hardware has also been advancing
rapidly implying that useful applications of quantum algorithmic advances may
not be far off. A variety of noisy intermediate scale quantum (NISQ) computers
are now available over the cloud, of which we use the trapped ion quantum
processing units (QPU) from IonQ. We also use a latest generation QPU at IonQ,
not yet available over the cloud, with substantially increased ion lifetime
and performance. This meaningfully increases the overall fidelity of
operations and decreases the amount of time needed for calibration, and makes
it possible to run a large number of quantum circuits using the same ions over
the course of a few hours.

The outline of our paper is as follows. We begin by introducing the basic
concepts of GAN, QGAN and QCBM necessary for this work in Section 2. In
Section 3, we introduce copulas and show how they can be learned using
maximally entangled quantum states. Section 4 gives a brief overview of the
quantum hardware used for our experiments. In Sections 5 and 6, we devise
training routines for our quantum generative models that are implementable on
a NISQ device. We then demonstrate via numerical simulation and executions on
the IonQ QPUs that our generative learning algorithms can learn from empirical
data. Section 7 compares the performance of our algorithms with a classical
GAN of similar number of parameters as well as classical parametric methods,
revealing advantages for the quantum methods. We argue via communication and
computational complexity that our quantum learning algorithms have a
computational advantage over classical algorithms for learning joint
distributions in Section 8. Section 9 provides an overview of future work. We
present a discussion and conclusion in Section 10.

##  2 Classical and Quantum Generative Adversarial Network, Quantum Circuit
Born Machine

Generative modeling is an unsupervised learning task that involves the
learning of patterns in the training data such that the model can generate new
examples as if they are drawn from the original data. A generative adversarial
network is a framework to train a generative model by transforming it to a
supervised learning task, in which two neural networks are trained
alternatively in a minimax game.

GAN employs two sub-models, a generator $G$ and a discriminator $D$.
$G_{\theta_{g}}:Z\to X$ represents a mapping from a latent space $Z$ to the
data space $X$, with $\theta_{g}$ denoting the parameters of $G$.
$D_{\theta_{d}}:X\to[0,1]$ represents a mapping from the data space $X$ to a
scalar, with $\theta_{d}$ denoting the parameters of $D$. For each $x\in X$,
$D(x)$ can be interpreted as the probability that $x$ belongs to the original
training data. $D$ is trained to maximize the probability of differentiating
true data from samples generated by $G$. Simultaneously, $G$ is trained to
’fool’ the discriminator by generating realistic data samples. In other words,
the two models $D$ and $G$ are trained in a two-player minimax game with value
function

$\min_{G}\max_{D}V(D,G)=\mathbb{E}[\log D(x)]+\mathbb{E}[\log(1-(D(G(z))))].$
(1)

After training, the generator is retained to generate new data samples for
downstream use.

In practice, the minimax optimization is realized with two loss functions, one
for the generator and one for the discriminator. For a batch of $m$ samples
from the real data $\\{x^{(l)}\\}$ and $m$ noise vectors $\\{z^{(l)}\\}$, the
loss function of the generator is

$L_{G}(\theta_{g},\theta_{d})=-\frac{1}{m}\sum_{l=1}^{l=m}[\log
D(G(z^{(l)}))]$ (2)

and the loss function of the discriminator is

$L_{D}(\theta_{g},\theta_{d})=-\frac{1}{2m}\sum_{l=1}^{l=m}[\log
D(x^{(l)})+\log(1-D(G(z^{(l)})))].$ (3)

A typical GAN training schedule is show in Algorithm 1.

Result: Trained Generator $G$

Initialize Network;

for  _$i\leftarrow 1$ to $iterations$_ do

Sample minibatch of $m$ noise vectors $\\{z^{(l)}\\}$ ;

Sample minibatch of $m$ data examples $\\{x^{(l)}\\}$ ;

Calculate Discriminator Loss $L_{D}$ ;

Update the discriminator by descending its stochastic gradient
$\nabla_{\theta_{d}}L_{D}(\theta_{g},\theta_{d})$ ;

Calculate Generator Loss $L_{G}$ ;

Update the generator by descending its stochastic gradient
$\nabla_{\theta_{g}}L_{G}(\theta_{g},\theta_{d})$ ;

end for

Algorithm 1 GAN Training Loop

The stochastic gradient descent is usually performed using adaptive algorithms
like Adam[36] with a suitable learning rate. Some researchers also update the
discriminator/generator multiple times in a training iteration or use
different learning rates, depending on the problem.

Intuitively, the generator is trying to get better at generating data samples
and thus fool the discriminator. The discriminator is trying not to be fooled.
The discriminator provides positive feedback to the generator, such that the
generator continues to improve. However, since it is a minimax optimization,
it is very common for GAN training to suffer from non-convergence (GAN fails
to learn well from the training data) and mode collapse (GAN generates data
with limited variety). While many techniques have been proposed to improve the
stability of GAN training, these are still common problems researchers face.

Recently, quantum generative adversarial network (QGAN) [27, 28] has been
proposed as a generalization of classical GAN that can run on a quantum
computer. QGAN can be used to learn either classical or quantum data. For the
learning of quantum data, each training sample $x$ becomes a quantum state
$\ket{x}\in\mathcal{H}_{X}$.
$G_{\theta_{g}}:\mathcal{H}_{Z}\to\mathcal{H}_{X}$ represents a quantum
circuit from the latent Hilbert space $\mathcal{H}_{Z}$ to the data space
$\mathcal{H}_{X}$, and $D_{\theta_{d}}:\mathcal{H}_{X}\to[0,1]$ can be either
a quantum circuit or a classical neural network. Note that if $D$ is a
classical neural network, initial measurement is needed to convert $\ket{x}$
to some classical data.

For the learning of classical data, one way of realizing the quantum generator
$G_{\theta_{g}}:Z\to X$ is to encode a noise vector $z$ as some quantum state
$\ket{z}$, apply a quantum circuit $U_{G}$, and then measure the expectation
of some observable $O_{G}$, i.e.,
$G(z)=\bra{z}U_{G}^{\dagger}O_{G}U_{G}\ket{z}$. Another way, as noted in Ref.
[30], is to represent $G_{\theta_{g}}$ as a quantum circuit $U_{G}$ followed
by some POVM (positive operator-valued measure) $M_{x}$ in $X$ space. In this
case, the dependence of $G$ on a latent space $Z$ is optional as POVM
measurement is inherently probabilistic and data samples are drawn from the
probability distribution

$\displaystyle q(x)=\text{tr}(M_{x}U_{G}\ket{0}\bra{0}U_{G}^{\dagger}).$ (4)

A Quantum Circuit Born Machine (QCBM) is a quantum generative model similar to
the generator part of QGAN, except the training objective is different. In the
QCBM case, if we suppose that the target distribution is $p$, a QCBM seeks to
minimize the distance between distributions $p$ and $q$ directly, through
either Kullback–Leibler (KL) divergence [26] or maximum mean discrepancy [37].
Here, we use the KL divergence as a cost function, defined as

$\displaystyle d_{\text{KL}}=\sum_{x}p(x)\log\frac{p(x)}{q(x)},$ (5)

##  3 Copula and Quantum Entanglement

Given random variables $(\mathcal{X}_{1},\ldots,\mathcal{X}_{d})$, probability
integral transform states that the marginal cumulative distribution functions
defined as $F_{i}(x)=\text{Pr}[\mathcal{X}_{i}\leq x]$ satisfies
$\text{Pr}[F_{i}(\mathcal{X}_{i})\leq u]=u$. Therefore the random variables

$(\mathcal{U}_{1},\ldots,\mathcal{U}_{d})=(F_{1}(\mathcal{X}_{1}),\ldots,F_{d}(\mathcal{X}_{d}))$
(6)

have marginals which are uniformly distributed on $[0,1]$.

A copula $C$ of $(\mathcal{X}_{1},\ldots,\mathcal{X}_{d})$ is defined as the
joint cumulative distribution function of
$(\mathcal{U}_{1},\ldots,\mathcal{U}_{d})$:

$C(u_{1},\ldots,u_{d})=\text{Pr}[\mathcal{U}_{1}\leq
u_{1},\ldots,\mathcal{U}_{d}\leq u_{d}].$ (7)

Reversing the above steps gives a way to generate samples from multivariate
distributions. Given a procedure to generate samples
$(\mathcal{U}_{1},\ldots,\mathcal{U}_{d})$ from the copula function $C$,

$(\mathcal{X}_{1},\ldots,\mathcal{X}_{d})=(F_{1}^{-1}(\mathcal{U}_{1}),\ldots,F_{d}^{-1}(\mathcal{U}_{d}))$
(8)

would be samples from the original multivariate distribution.

Sklar’s theorem states that every multivariate cumulative distribution
function can be expressed in terms of its marginals and a copula, i.e.,

$\displaystyle\text{Pr}[\mathcal{X}_{1}\leq x_{1},\ldots,\mathcal{X}_{d}\leq
x_{d}]=C(F_{1}(x_{1}),\ldots,F_{d}(x_{d})).$ (9)

Copulas are used to simulate correlated variables because they remove the
structure of the marginal distributions, and capture just the point-wise
correlation between variables. Simple copula formulas are extensively used in
finance, engineering and medicine. However these simple formulae often fall
short in accurately capturing the relationship between variables. Instead more
sophisticated or empirical copulae can be used but these become hard to model
and sample from as the number of variables increases. Very recently,
generative models have been successfully used to model copulas [38, 39].

For our study, we test our generative models using the following framework

  * •

Use Eq. (6) to transform training data into copula space

  * •

Fit known copula model/GAN/QGAN/QCBM in copula space

  * •

Sample from known copula model/GAN/QGAN/QCBM in copula space

  * •

Use Eq. (8) to transform synthetic data from copula space back to original
space

Next, we show that every copula with density can be represented by a maximally
entangled state. Suppose $C:[0,1]^{d}\to[0,1]$ is a $d$-dimensional copula
with density $c:[0,1]^{d}\to[0,1]$. $c$ can be obtained as

$\displaystyle
c(u_{1},\ldots,u_{d})=\frac{\partial^{d}C(u_{1},\ldots,u_{d})}{\partial
u_{1}\cdots\partial u_{d}}.$ (10)

This can be represented as a quantum state

$\displaystyle\ket{c}=\sum_{u_{1},\ldots,u_{d}\in[0,1]^{d}}\sqrt{c(u_{1},\ldots,u_{d})}\ket{u_{1}}\cdots\ket{u_{d}}.$
(11)

Since $\mathcal{U}_{1},\ldots,\mathcal{U}_{d}$ have uniform marginals,

$\displaystyle\int\dots\int c(u_{1},\ldots,u_{d})du_{1}\cdots
du_{i-1}du_{i+1}\cdots du_{d}=1~{}~{}~{}\forall i.$ (12)

This means the reduced density matrix $\rho_{i}$ for partition $i$ is
completely mixed and $\ket{c}$ is a maximally entangled state with respect to
each of the $d$ partitions.

Here we construct a quantum circuit that can prepare such maximally entangled
states up to relative phases, which we christen as a ‘qopula’ circuit. The
quantum circuit for two random variables is shown in Fig. 1 (a). The first
part of the circuit consists of the creation of $N_{q}/2$ Bell pairs which are
distributed between two registers A and B. This results in the creation of a
state of the form

$\displaystyle\frac{1}{2^{N_{q}/4}}\sum_{i=0}^{2^{N_{q}/2}-1}\ket{i_{A}}\ket{i_{B}}.$
(13)

In this state, the reduced density matrix of the states in each register is
completely mixed ($\rho_{A}=I/|A|$). In the next step, unitaries act locally
on registers A and B. Since $U^{\dagger}IU=I$, the action of the unitaries
does not change the reduced density matrix, but it does change the
correlations between the registers A and B. Note that the unitaries on each
register can be set by parameters that are independent of the other.

If each register is measured in the computational basis with output
$x_{1},\ldots,x_{n},y_{1},\ldots,y_{n}$, the measured bitstrings can represent
discretized values of each random variable in the domain $[0,1]$ through the
following linear transformation

$\displaystyle x=\frac{1}{2^{1+n}}+\sum_{i=0}^{n}\frac{1}{2^{n}}x_{n}.$ (14)

Therefore this ansatz can model the correlations between uniform random
variables and capture the copula dependence structure.

The number of qubits in each register determines the discretization of the
probability distribution. The more the qubits, the more finely the point-wise
correlation can be learnt. If the number of qubits is restricted, for
practical application, one can use the available qubits to learn the
significant bits and then append random bits in the position of the less
significant bits.

(

(

Figure 1: (a) The ‘qopula’ circuit for 2 random variables. The top (bottom)
half of the circuit consists of qubits that belong to the register which
provides samples of the first (second) random variable. (b) The ansatz for $U$
corresponding to 3 qubits. All gates are parametrized by angles which are
optimized during the learning process. The structure can be repeated for
multiple layers each with different parameters. Here the 2-qubit RXX gates
represent $\exp(i\theta X_{i}X_{j})$ acting on qubits $i$ and $j$, which can
be executed as $\text{CNOT}(i,j)\exp(i\theta X_{i})\text{CNOT}(i,j)$ using
standard gates available on gate-model quantum computers. The RZ and RX gates
represent rotations around the Z and X axes.

a) @C=1em @R=1.7em —0⟩ & H 4 3U_A
—0⟩ H 4 U_A
⋮
—0⟩ H 4 U_A
—0⟩ 3U_B
—0⟩ U_B
⋮
—0⟩ U_B
b) @C=1em @R=1.7em & RZ RX RZ RXX1 RXX2 …
RZ RX RZ RXX RXX1 …
RZ RX RZ RXX RXX …

The circuit can be extended to more than 2 variables by preparing GHZ states
instead of Bell pairs as the starting point on which the operators $U$ act.
The exact structure of $U$ is not restricted and a subject of future research.
Here we use a structure that consists of layers alternating between a ‘driver’
layer that consists of parametrized single qubit rotations RZ and RX and an
‘entangler’ layer that consists of $RXX$ gates as shown in Fig. 1 (b) [40]. We
also note that the uniform marginal condition can be satisfied by many states
other than the simple Bell or GHZ states [41].

The dataset that we use is the daily return of AAPL and MSFT between
2010-2018. Fig. 2 (a) represents the hypothetical growth of capital if $10,000
is invested in both stocks (assuming dividends are re-invested). The daily
returns are calculated as the percentage change of close prices between
adjacent trading days. Suppose on Day 1 a stock’s close price is $P_{1}$ and
its close price is $P_{2}$ on Day 2, then its daily return on Day 2 is
$r=(P_{2}-P_{1})/P_{1}$. The daily return data is plotted both in
$\mathcal{X}$ space (real data space, Fig. 2 (b)) and $\mathcal{U}$ space
(copula space, Fig. 2(c)). In model training, no time series model is assumed
and hence the daily returns are assumed to be independent events.

Figure 2: (a) Hypothetical Growth of $10,000 invested in AAPL and MSFT between
2010-2018 (b) Scatter Plot of Daily Returns (c) Scatter Plot of Data after
Probability Integral Transform Figure 3: Breakdown of fractional time spent on
various system tasks for currently available computers on the cloud (left) and
next generation IonQ system (right). The next generation system has a marked
improvement in the time spent loading ions ( 40x reduction) and in the amount
of time spent running algorithms ($>2$X improvement).

##  4 Quantum Hardware

The experiments were run on IonQ trapped ion quantum processing units (QPUs),
which utilize stable ground states of trapped ytterbium atomic ions as qubits.
Two ground states of these atoms with hyperfine splitting are used as the two
qubit states, and they are initialized, manipulated and detected with
adequately tailored laser lights applied to them. The quantum circuit is
executed by applying quantum logic gates driven by Raman transitions using up
to 32 individually addressed beams of 355 nm light. These Raman transitions
are programmed to provide a universal gate set by applying individual
rotations of a given qubit’s internal states to realize a single-qubit gate,
as well as two-qubit gates by inducing Molmer-Sorenson type interaction to
generate entanglement between a pair of qubits in the system. IonQ’s QPUs
feature all-to-all connectivity, where a two-qubit gate can be applied
directly to any pair of qubits in the system regardless of their physical
location [42]. IonQ’s QPUs available on the cloud are self-calibrating, which
means that the system monitors the quality of the quantum logic gates
constantly and makes sure they are fully calibrated before the computational
tasks are executed. The hardware used in this work utilizes IonQ QPUs that are
available over the cloud, as well as a next generation QPU with significant
hardware improvements in vacuum quality, optical stability, and crosstalk in
ion addressing. These major improvements led to about an order-of-magnitude
improvements in gate errors, and reduction in the fraction of time spent in
ensuring that the system is ready to execute computational tasks. Figure 3
shows the improvement made on the fraction of time spent running computational
tasks between the current and next generation systems.

##  5 Quantum Learning through QGAN

In our QGAN implementation, the quantum generator comprises of 6 qubits, and
each unitary uses 1 layer of ansatz as described in Fig. 1. Therefore the
quantum generator has 24 trainable parameters. A measurement in the
computational basis is performed at the end of the circuit and Eq. (14) is
applied to get the 2-dimensional sample in copula space. To deal with the
issue of discretization, the measurement output is extended from
$(x_{0},x_{1},x_{2},y_{1},y_{2},y_{3})$ to $(x_{0},x_{1},\cdots,x_{22},$
$y_{0},y_{1},\cdots,y_{22})$ with the additional 40 bits uniformly randomly
generated, and then the linear transformation is applied.

The discriminator is a feed-forward neural network with input dimension 2, a
hidden layer of dimension 32 and an output layer of dimension 1. The input and
hidden layer comprises of a linear layer and Leaky Rectified Linear Unit
(ReLU), whereas the output layer is a linear layer followed by Sigmoid
function.

To compare the expressivity of QGAN vs GAN, we design the generator of
classical GAN to have the same number of trainable parameters. This can be
achieved by a 2-layer feed-forward neural network with input dimension 6 and
output dimension 2. The input layer consists of a linear layer followed by
batch norm and ReLU. The output layer consists of a linear layer followed by a
sigmoid function. A sigmoid function is needed to transform samples back to
the unit square.

Neural-network models are typically trained using gradient-based optimizers.
In the case where the cost function has a quantum circuit component, the
gradient with respect to circuit parameters can be calculated using the
parameter-shift rule [43]. However, this procedure requires two executions of
the original circuit (each with shifted parameters), and thus the number of
circuit executions scales two times as the number of circuit parameters. For
example, in our case, it would take 48 circuit executions to compute the
gradient with respect to the 24 parameters in the quantum generator. This is
executed sequentially in most software for quantum computers. In addition to
slower evaluation of the cost function itself, since commercially available
quantum computers (including IonQ’s) are typically cloud-based, sequential
evaluation will also be slow due to network latency. For this reason, we
avoided the computation of quantum gradients and used the simultaneous
perturbation stochastic approximation (SPSA) algorithm [44, 45]. SPSA
algorithm updates the parameter in an iterative process, with the $k$-th step

$\displaystyle\theta_{k+1}=\theta_{k}-a_{k}\hat{g}_{k}(\theta_{k}),$ (15)

where the $i$-th component of the gradient estimator $\hat{g}_{k}(\theta_{k})$
is

$\displaystyle\left(\hat{g}_{k}(\theta_{k})\right)_{i}=\frac{C(\theta_{k}+c_{k}\Delta_{k})-C(\theta_{k}-c_{k}\Delta_{k})}{2c_{k}(\Delta_{k})_{i}}.$
(16)

Here $\Delta_{k}$ is a random perturbation vector with each element drawn
uniformly at random from $\\{-1,1\\}$. For initial learning rate $a$ and step
size $c$, $a_{n}$ and $c_{n}$ are updated as $a_{k}=a/k$ and
$c_{k}=c/k^{\gamma}$, where $\gamma$ is chosen from $[1/6,1/2]$. Through
numerical simulation, we observe that running SPSA algorithm with 3-5
iterations in place of a gradient descent step for the generator, yields good
performance in our QGAN training schedule. Since each iteration step only
requires 2 circuit executions, an optimization step using SPSA would only
require 6-10 circuit executions in total.

Therefore, our QGAN training schedule will be as shown in Algorithm 2.

Result: Trained Generator $G$

Initialize Network and angles of the Quantum Circuit ansatz;

for  _$i\leftarrow 1$ to $iterations$_ do

Run the quantum generator with $m$ shots to generate $m$ measurements ;

Sample minibatch of $m$ data examples $\\{x^{(l)}\\}$ ;

Calculate Discriminator Loss $L_{D}$ ;

Update the discriminator by descending its stochastic gradient
$\nabla_{\theta_{d}}L_{D}(\theta_{g},\theta_{d})$ ;

Calculate Generator Loss $L_{G}$ ;

Update the Quantum Generator by using SPSA algorithm ;

end for

Algorithm 2 QGAN Training Loop

We first test our model in simulation before running it on IonQ’s cloud QPUs.
All simulations and experiments of QGAN are trained with 1000 iterations and
use the parameter $a=0.008$, $c=0.01$, number of iterations $n_{iter}=5$,
$\gamma=0.101$ for SPSA, learning rate 0.0015 for the discriminator, batch
size $m=2048$. To reduce the experiment time, we used the maximum learning
rates possible without breaking down the learning. To better understand the
effect of hardware noise, random initialization with the same seed is used for
comparing simulation and experiment results. However, other random
initializations were also used in simulation and results are summarized in Sec
7.

(a) Loss Plot from Simulator (b) Loss Plot from Experiment

Figure 4: Plot of the loss functions from QGAN training for $N_{q}=6$ qubits.

Results from the simulator, Fig. 4 (a) indicate that in our training process,
the losses from the generator and discriminator first diverge before
converging around iteration 250. This implies that, the discriminator learns
faster than the generator in the beginning, but then the generator experiences
an Aha! moment and quickly catches up.

As quantum computers at the moment still have significant noise, it is
essential to investigate how noise would affect our training algorithm. We ran
the same training algorithm on a mixed state simulator with the addition of up
to $4\%$ depolarizing noise following each two-qubit gate, which is an
appropriate model for noise on IonQ’s cloud QPUs [46]. We observe that the
losses of the generator and discriminator diverge a bit further. This is to be
expected as noise makes the generator a weaker learner and therefore it needs
more time to learn. However, in most of the cases, our QGAN converges
eventually (anywhere between iterations 300-500), indicating robustness
against noise and suitability to run on NISQ machines. We observe that in our
experiment run on IonQ’s cloud QPUs, the losses converged at around iteration
400, agreeing with results from noisy simulations.

QGAN-generated data is visualized in both $\mathcal{X}$ space (real data
space, Fig. 5(a)) and $\mathcal{U}$ space (copula space, Fig. 5(b)). We
observe that the data in copula space (Fig. 5 (b)) are concentrated around a
few clusters and are less evenly distributed than the original data (Fig. 2
(c)). This is because we used only 3 qubits to represent each dimension and
therefore the output is very discrete. One would only see 64 unique data
points had we not added the random bits after measurements. The random bits
have the effect of creating small fluctuations around the original discrete
value, hence the clusters.

Figure 5: (a) Scatter Plot of QGAN-generated Daily Returns after Reverse
Probability Integral Transform (b) Scatter Plot of QGAN-generated Copula Data.
$N_{q}=6$ qubits were used.

##  6 Quantum Learning through QCBM

For the QCBM, we use the KL divergence as a cost function. As the target, we
bin the data in copula space into $2^{N_{q}/2}$ bins, where $N_{q}$ is the
number of available qubits, generating a discrete probability distribution
$p$. This is compared to the probability distribution $q$ over the
computational basis states measured at the output of the circuit. The KL
divergence cost function is

$\displaystyle d_{\text{KL}}=\sum_{i,j}q_{i,j}\log\frac{q_{i,j}}{p_{i,j}},$
(17)

where $i$ and $j$ are for the first and second variable respectively.
Numerically, we implement the clipped version of the KL divergence where
$p_{i,j}$ are set to $10^{-6}$ if below that value. The QCBM training schedule
is shown in Algorithm 3.

Result: Trained Quantum Circuit

Initialize angles of the Quantum Circuit ansatz;

for  _$i\leftarrow 1$ to $iterations$_ do

Run the quantum circuit with $m$ shots to generate $m$ measurements ;

Calculate the KL divergence $d_{KL}$ ;

Update the circuit angles by using SPSA algorithm ;

end for

Algorithm 3 QCBM Training Loop

Figure 6: QCBM Training (a) Average of the 10 minimum values of the cost
function reached during training as a function of the number of qubits for
optimization parameters $a=0.5$, $c=0.5$. The error bars are smaller than the
marker sizes. (b) Convergence of the cost function as a function of the number
of iterations for $a=0.1$, $c=0.1$. The panels indicate results from different
numbers of qubits $N_{q}$.

For converting the measured bits into samples, Eq. (14) is used along with the
procedure of appending the random bits as described in the QGAN section.

The results of training are shown in Fig. 6. The results are shown for a
series of system sizes consisting of 2, 4, and 6 qubits. A one-layer ansatz
was used in each case. Each system size was initialized with the optimal
angles from the previous smaller system size which can be thought of as a form
of transfer learning. The training was done on the simulator, IonQ’s cloud
QPUs, and its latest generation QPU which is not yet accessible through the
cloud.

The optimization algorithm used for training is the SPSA similar to the QGAN
training. 200 iterations were used for all experiments. Fig. 6 (a) shows
results for up to 8 qubits on the simulator and cloud quantum computer and up
to 6 qubits on the next generation quantum computer. The optimization
parameters used were $a=0.5$ and $c=0.5$, for which the training finds the
minimum within 20 iterations. As can be seen in this figure, which shows the
average over the 20 smallest values of the measured cost function, for 2 and 4
qubits, all three backends converge to a similar value. For 6 qubits, somewhat
surprisingly, we observe that the noisy hardware has a smaller value of the
final cost function. This observation of the noise being favorable for
converging to a low value of the cost function at the beginning of the
optimization is further confirmed in Fig. 6 (b) which shows the convergence
for $a=0.1$ and $c=0.1$, where it is slow and hence can be seen clearly. We
note that we do not expect that the noise will continue to remain beneficial
at larger system sizes however.

While vanishing gradient (barren plateau) has been observed in several
variational quantum algorithms [47, 48], we note that we see no evidence of it
even as the system size increases. Instead, the training landscape becomes
noisier even for the simulator as shown by the broadening of the cost function
values as the number of qubits increases. This is to be expected since only a
finite number of shots are used whereas the number of possible outputs
increases exponentially with the number of qubits. We note that this would
also be true for a classical generator. More discussion of the QCBM training
can be found in Appendix 12.

We further note that for the cases we tested, the QCBM technique requires
fewer function calls to train the model than the QGAN approach. Which
technique may end up performing better will depend on the details of the
application. Finally, we also note that an improved convergence is not
observed for the increased quantum gate fidelity of the next generation QPU
for either approach, implying that the noise plays some positive role in the
statistical machine learning models. The major advantage of the next
generation QPU in this application is its increased stability that leads to
reduced runtime of the entire optimization task by about a factor of 2.

##  7 Evaluation of Model Performance

We use the 2-dimensional 2-sample Kolmogorov-Smirnov (KS) test[49] to evaluate
our generative models together with classical GAN and a parametric model. The
null hypothesis of the 2-sample KS test is that the two samples are drawn from
the same distribution, and the alternative hypothesis is that they are drawn
from different distributions. It is typically performed on 1-dimensional
distributions but here we are performing on 2-dimensional distributions. We
use a significance level of 0.05 here. If the p-value produced by the test is
more than the significance level specified, then we accept the null
hypothesis. In the case where multiple models produce samples which are
accepted by the KS test, we use the associated KS statistics to compare among
the models. The KS statistic quantifies the distance between the empirical
distributions of two samples. For 1-dimensional data, it can be defined as

$\displaystyle D_{KS}(P,Q)=\sup_{x}|P(X\leq x)-Q(X\leq x)|,$ (18)

where $P$ and $Q$ are the empirical cumulative distributions. Please refer to
Ref. [49] for the 2-dimensional generalization.

For the classical models and QGAN, 2048 samples are collected from our
generative models, and for the QCBM, 2000 samples are collected from 4
circuits with minimal cost function that were run with 500 shots each. These
are compared with a bootstrapped version of the training data. The
architecture of classical GAN is described earlier. The parametric model we
used involves fitting the copula data with a Gaussian copula.

Model |  $D_{KS}$ (the smaller the better) | p-value (threshold 0.05)
---|---|---
Parametric model | 0.0449 | 0.117
Classical GAN | 0.0363 - 0.0508 | 0.0530 - 0.309
QGAN simulation | 0.0320 - 0.0396 | 0.226 - 0.473
QGAN experiment, QPU cloud | 0.0352 | 0.3570
QCBM simulation | 0.0425 - 0.0520 | 0.0511 - 0.1717
QCBM experiment, QPU cloud | 0.0373 - 0.0515 | 0.0548 - 0.3030
QCBM experiment, QPU Next Gen | 0.0330 - 0.0510 | 0.0578 - 0.4465
Table 1: KS statistics and p-value of KS test across multiple models. The
quantum models use $N_{q}=6$ qubits. Figure 7: KS metric averaged over the 10
best converged values of the cost function as a function of system size for
the QCBM.

For classical GAN and QGAN on simulator, we were able to run the model
multiple times with different initializations and thus there is a range of
results. For QGAN, we initialized the parameters uniformly at random from
$[0,2\pi]$ and note that our model converges regardless of initialization. For
classical GAN, we used the default initialization method from PyTorch [50]
with different seeds. We note that about $40\%$ of our model instances are
accepted by the 2D KS test with threshold 0.05. The rest failed to learn and
were rejected by the test. In Table 1, the metrics of those accepted instances
were included. For the QCBM, we first collected data from 5 sets of 4 jobs
each that corresponded to the minimum cost function values found during the
optimization. For each set, we calculated the performance by affixing less
significant random bits to the ones from the generator. We then also averaged
over 20 such sets of random bits.

As one can see, the parametric model does a good job in modeling the copula
but generative models are able to do better. The result from QGAN is
consistent and outperforms classical GAN with similar number of parameters.
The QCBM performs slightly worse than the QGAN but comparable to the classical
GAN. Fig. 7 shows the results for the $D_{KS}$ as a function of number of
qubits for the QCBM training. Upto 6 qubits, the quantum computers perform
better for QCBM training than the simulator as is to be expected from their
better performance in copula space as well. We note that $D_{KS}$ decreases
upto 6 qubits as expected but increases after that indicating the need for
more shots and iterations for the optimization to converge.

We note that it’s entirely possible for classical GAN to achieve better
results with a deeper neural network as generator. However, the same can be
said about our quantum generative models as they are only learning the
significant bits at the moment. With better hardware, more qubits, more layers
and better design for the ansatz $U$, our quantum generative models will also
do better. Nevertheless, our results provide compelling evidence that our
quantum generator has more expressivity than classical generative models.

We also note that we are able to train QGAN/QCBM at a much faster learning
rate and therefore conclude the training with much fewer iterations than
classical GAN. In classical GAN, the learning rate used is 0.0001 and model
training concludes after 20000 iterations. Attempts to increase the learning
rate failed due to non-convergence in model training. In QGAN, model training
concludes after 1000 iterations. In QCBM for 6 qubits, the training converges
to a good value for as little as 20 iterations.

Currently the gate speeds of commercially available QPUs are typically slower
than those of classical computers. Additionally, due to the current access
mode of IonQ QPUs over the cloud that includes queued access priority, network
latency and calibration time on the machine, our QGAN experiment took 2 weeks
to complete. At this time, this is much slower as compared to other models.
The training time for our classical GAN is 6 minutes on a CPU-only machine and
that of QGAN simulation is 4 minutes. For the QCBM on the next generation
quantum hardware, the training for 6 qubits with 200 iterations took
approximately 5 hours whereas it took 9 hours on the cloud system. We expect
many of the network latencies and calibration procedures to improve on QPUs in
the future, as well as dedicated access to be available. Further, as shown in
the next section, there will be certain joint probability distributions for
which the run time for the classical models and the simulations are expected
to grow exponentially with the number of qubits, while the QPU time will only
grow polynomially, ultimately making the classical and quantum technologies
competitive for learning applications.

##  8 Quantum Computational Advantage

We next present arguments for the presence of computational advantage in our
quantum machine learning algorithms over classical machine learning algorithms
such as neural networks.

One source of advantage during learning stems from Bell’s Theorem. Consider a
Bell pair of two qubits. One qubit is given to an entity (say Alice) and the
other to a different entity (say Bob). At this stage, they each also have
access to random variables which share arbitrary correlations. Next, suppose
they are given measurements $x$ and $y$ which are not shared, and based on
these they produce outcomes $a$ and $b$ respectively. Then quantum mechanics
implies that the joint conditional probability distribution, $P(a,b|x,y)$,
cannot be reproduced by classical means without communication between Alice
and Bob. Further, when $n$ Bell pairs are shared between Alice and Bob, some
quantum correlations quantified by $P(a,b|x,y)$, where $x$ and $y$ are
selected from a set of $\\{0,1\\}^{2^{n}}$ measurements, can only be
reproduced by classical means if $\mathcal{O}(2^{n})$ number of bits are
exchanged between Alice and Bob [51].

Let us analyze the implication of this for the quantum learning algorithms
described in the paper vs. classical learning schemes. The communication flow
in the quantum learning scheme is shown in the left part of Fig. 8. Once the
Bell pairs are created, the classical computer executing the optimization
scheme can separately send $x$ and $y$ to Alice and Bob respectively. The
measurements $a$ and $b$ are then collectively processed by the classical
optimizer to produce the input to the quantum computer for the next iteration.

Now consider how one would simulate the quantum learning algorithm
classically. Let us assume that, as is usually the case, the classical
algorithm is deterministic in nature, i.e., for every unique input, there is
one unique output. In general, the input to the classical generator is a
vector of random bits, which has to contain at least $n$ bits for each random
variable to produce as many unique outputs as produced in the quantum case.
Distribute these bits equally between Alice and Bob. Any classical operation
can be written as a combination of operations that only acts on the bits of
Alice or Bob separately, and operations that act on both sets of bits. The
latter can be further modeled as operations on the bits belonging to Alice,
followed by communication to Bob and then operations on just on the bits
belonging to Bob, and vice versa (right part of Fig. 8). Then, without loss of
generality, we can restrict Alice’s operation $C_{A}$ to be of the form
$C_{A}(x,d_{AB})$, i.e. parametrized by $x$ from the optimizer, and the
sequence of bits $d_{AB}$ communicated from Bob to Alice. Similarly, for Bob
we have $C_{B}(y,d_{BA})$. Note that $d_{AB}$ can contain information about
$y$ and $d_{BA}$ about $x$, so we have not restricted the form of the
classical generator in any fashion. Then, in order to simulate $P(a,b|x,y)$
for arbitrary $x$ and $y$ sent by the classical optimizer, it follows that the
amount of communication between Alice and Bob, $d_{BA}+d_{AB}$, will scale as
$\mathcal{O}(2^{n})$.

Figure 8: Comparison of quantum (left) vs. classical (right) learning. $U_{A}$
and $U_{B}$ are unitaries that act on entangled Bell pairs distributed between
Alice and Bob. $C_{A}$ and $C_{B}$ are operations on random classical bit
strings belonging to Alice and Bob.
$\vdots$$\vdots$$\chi_{1}$$\chi_{2}$$\chi_{3}$$\chi_{k}$$o_{1}$$o_{m}$InputlayerOutputlayer
Figure 9: Feedforward Neural Network with dimension $[k,m]$

Next, let us see the concrete implications of this argument by considering a
2-layer feedforward neural network architecture as shown in Fig 9. We assume
Alice and Bob share the neurons at each layer and their job is to jointly
compute $o=W\chi+b$. Writing it in blocks as

$\begin{bmatrix}o_{A}\\\
o_{B}\end{bmatrix}=\begin{bmatrix}W_{A\mathcal{A}}&W_{A\mathcal{B}}\\\
W_{B\mathcal{A}}&W_{B\mathcal{B}}\end{bmatrix}\begin{bmatrix}\chi_{\mathcal{A}}\\\
\chi_{\mathcal{B}}\end{bmatrix}+\begin{bmatrix}b_{A}\\\ b_{B}\end{bmatrix},$

where $W_{A\mathcal{A}},W_{B\mathcal{A}},\chi_{\mathcal{A}},b_{A}$ belongs to
Alice and $W_{A\mathcal{B}},W_{B\mathcal{B}},\chi_{\mathcal{B}},b_{B}$ belongs
to Bob, Alice locally stores $W_{A\mathcal{A}}\chi_{\mathcal{A}}$ and sends
$W_{B\mathcal{A}}\chi_{\mathcal{A}}$ to Bob. Similarly, Bob locally stores
$W_{B\mathcal{B}}\chi_{\mathcal{B}}$ and sends
$W_{A\mathcal{B}}\chi_{\mathcal{B}}$ to Alice. After the communication, Alice
receives $W_{A\mathcal{B}}\chi_{\mathcal{B}}$ and computes
$o_{A}=W_{A\mathcal{A}}\chi_{\mathcal{A}}+W_{A\mathcal{B}}\chi_{\mathcal{B}}+b_{A}$,
and Bob does so similarly. In this process, Alice sends $\mathcal{O}(|B|)$
bits to Bob and Bob sends $\mathcal{O}(|A|)$ bits to Alice. Since $|A|+|B|=m$,
$O(m)$ bits of communication occur.

Thus, if the generator is represented by a deep feedforward neural network
where each layer has at most $m$ neurons, then the number of layers has to
scale exponentially asymptotically as $\mathcal{O}(2^{n}/m)$ to reproduce
certain probability distributions generated by experiments on the Bell pairs.
This ultimately implies that a neural network will have to have exponentially
scaling depth in the precision of the output.

In the argument above, $x$ and $y$ belong to $\\{0,1\\}^{2^{n}}$. For the
quantum generator, they can be more compactly specified as angles to the
ansatzes specified by $U_{A}$ and $U_{B}$ respectively, which can even consist
of a constant number of layers. This is related to the sampling complexity
argument presented next.

The second source of advantage stems from the quantum supremacy argument of
“instantaneous Quantum Polynomial” (IQP) circuits [52]. The output
distribution of IQP circuits is of the form

$\displaystyle q^{IQP}(z)=|\bra{z}H^{\otimes n}U_{D}H^{\otimes n}\ket{0^{n}}|$
(19)

where $U_{D}=\exp(iD)$ and

$\displaystyle D=\sum_{k,l}J_{k,l}Z_{k}Z_{l}+\sum_{k}M_{k}Z_{k}.$ (20)

In Ref [52], it’s argued that this type of distribution cannot be efficiently
sampled classically, unless the polynomial hierarchy [53] collapses to the
third level.

We note that the output distribution of IQP circuits is a special case of the
conditional output distribution from our circuit. The output distribution of
our circuit can be written as

$\displaystyle q(z_{A},z_{B})=\left|\bra{z_{A},z_{B}}U_{A}\otimes
U_{B}\left(\frac{\ket{0_{A}0_{B}}+\ket{1_{A}1_{B}}}{\sqrt{2}}\right)^{\otimes
N_{q}/2}\right|.$ (21)

If we set all the parameters in $U_{B}$ to 0 such that $U_{B}$ is trivial,
then

$\displaystyle q\left(z_{A}|z_{B}=0\cdots
0\right)=|\bra{z_{A}}U_{A}\ket{0^{\otimes N_{q}/2}}|$ (22)

If we set $U_{A}$ to use 1 layer of the ansatz and set angles of the RZ
rotations to 0, then we can see that Eq. (22) is of the same form as Eq. (19),
as $HZH=X$. Hence if the output distribution of our circuit can be efficiently
sampled classically, it would also imply the collapse of the polynomial
hierarchy (see Appendix 13).

We note that while it is true that some states and measurements exhibit the
separations outlined above, these are arguments for the worst case and it may
not apply to the ones relevant for a particular learning problem. While
determining the generality of these exponential communication and
computational separations is a topic of ongoing research, the reasoning above
clearly shows that the quantum computer expands the space of joint probability
distributions that can be efficiently explored for learning.

##  9 Future Work

GANs are prone to mode collapse and non-convergence, and as such many
improvements have been made to mitigate these issues. While our pedagogical
case did not suffer from them, they may occur with more complex models and
problems. Therefore it would be interesting to see if our framework can
benefit from the GAN enhancements discussed in this section. For instance, one
may use a critic ($C$) in place of a discriminator. Instead of using Binary
Cross Entropy (BCE) as loss function, one may instead use Wasserstein Loss
(W-Loss) [54]

$\min_{G}\max_{C}V(G,C)=\mathbb{E}[C(x)]-\mathbb{E}[C(G(z)],$ (23)

which tends to make model training more stable. W-Loss solves the vanishing
gradient problem of the discriminator, in which the discriminator doesn’t
provide enough information for the generator to make progress in the BCE case.

WGAN can be implemented straightforwardly in our case as we can replace the
sigmoid function with a linear activation function such that the discriminator
becomes the critic and replace BCE with W-Loss as the loss function. We note
that this is different from Ref. [55] as we are learning classical data and
the discriminator is classical, whereas in Ref. [55] both the discriminator
and the generator are quantum.

Another way to enhance our model is the analogue of the addition of noise
vectors from the latent space as input to the classical GAN. Recall that in
classical GAN, once a model is trained, this input noise vector completely
determines the model output (see Section 2). A significant improvement with
the addition of input noise vectors is that it allows the controlled
generation of data, such as changing specific features of the output by
tweaking the noise vector. In our example, this would mean that one can
generate samples with desired characteristics such as where all absolute
returns are larger than a preset value, as in the case of a black swan event.
Controlled generation allows for the emergence of desired features by updating
the noise vector according to pre-trained feature classifiers. This would also
allow more complex GAN architectures such as styleGAN [56].

A proposed circuit for such a generator from Section 2 allows one to calculate
the expectation value of an observable with respect to the noise vector in the
latent space:

$\displaystyle G(z)=\bra{z}U_{G}^{\dagger}O_{G}U_{G}\ket{z}.$ (24)

$\ket{z}$ may be realized in the following way, where the $Z$ gates are
parametrized rotations defined by the noise vector $z$.

@C=1em @R=1.7em —0⟩ & H Z
—0⟩ H Z
⋮
—0⟩ H Z

It is important to note that this model might be more costly to implement on
quantum hardware. A circuit execution with 1000 shots corresponds to 1000 data
samples in our current model architecture, whereas if using expectation values
this would only correspond to one data sample. However, in the long run, such
a model might outperform our current model due to the mitigation of problems
mentioned above.

Moreover, it is not clear if such a model will be advantageous against
classical cases in terms of complexity. While we have argued that producing
the distribution in our case is difficult for a classical generator, it might
not be the case for producing expectation values as this procedure is
deterministic (up to small statistical noise from the measurement). For
example, it has been argued that even the lowest depth version of the Quantum
Approximate Optimization Algorithm (QAOA) can produce distributions which are
hard for classical computers [57]. However, the computational task of QAOA is
to optimize the expectation value
$\bra{\boldsymbol{\gamma},\boldsymbol{\beta}}C\ket{\boldsymbol{\gamma},\boldsymbol{\beta}}$
of some objective function $C$ and it’s unclear if the low depth version of
QAOA can outperform classical algorithms [58].

Another area of potential exploration is the development of higher-level
ansatzes. Currently, the number of trainable parameters in our ansatz scales
quadratically with the number of qubits. This is undesirable as we use more
qubits to represent each dimension of data. It may be beneficial to develop
ansatz analogue to building blocks in classical neural networks (e.g., linear,
batch norm), where operations are applied to the variables instead of the bits
representing those variables.

Lastly, as discussed in Section 5, instead of utilizing SPSA, the generator
can be updated using gradient-based optimizers with the parameter-shift rule.
This could allow for more efficient optimization through the cost function.
Implementing the parameter-shift rule in parallel would reduce the runtime for
quantum gradient estimation. However, this is only possible if the quantum
hardware has enough qubits and low crosstalk noise. A networked quantum
computer architecture would also make this possible.

##  10 Discussion & Conclusion

In our view, the theory of quantum machine learning (QML) has gone through
three development phases. The first phase consists of using quantum computers
to speed up subroutines of classical learning algorithms [59, 60]. These
algorithms most often used the quantum linear systems algorithm [61] as the
basis of quantum advantage. The resource requirements for this class of
algorithms are very high and it often requires loading classical data in a
superposition. The second stage, which started with the advent of cloud-based
machine learning consists of heuristic techniques utilizing variational
quantum circuits which are applied to discriminative or generative problems
[62, 63, 29]. While these have been experimentally demonstrated on very small
problem sizes, these are heuristic in nature without theoretical support for
their advantage over classical methods. Since they cannot be tested at scale,
often their efficacy cannot even be numerically predicted.

The third phase of QML development is now under way. It consists of quantum
algorithms with reasonable resource requirements. Moreover, the structure of
circuits can often be shown to have advantage over the corresponding classical
algorithms. For example, in Ref. [64], quantum advantage is proven for certain
problems where the objective is achieving a specified worst-case prediction
error. Ref. [65] uses the hardness of the discrete log problem to show that a
quantum classifier can achieve high accuracy and is robust to additive error
where a classical classifier can do no better than random guessing. Perhaps
the simplest example of advantage is the proof of separation in expressive
power between commonly-used Bayesian networks and their minimal quantum
extension shown in Ref. [66] based on results from quantum foundations. Our
work falls in this third phase. We anticipate that the future of practical QML
will follow this path.

For the particular algorithm presented here, we can project when it will
become competitive with classical techniques as follows. Typical structured
datasets have dimension of order 10-100. With each dimension represented by
5-10 qubits, we expect our model to have production value if running on
quantum hardware of 50 qubits or more.

In conclusion, we have demonstrated quantum generative learning for
multivariate data using a QGAN and a QCBM. Our learning models utilize a
variational ansatz that we have designed which has provable advantage in
capturing correlations between the data over classical methods. We note that
our training of the quantum generative models does not show the barren plateau
problem and is resilient to noise. The outcome of the QGAN training reaches a
better performance and is more stable than training of the equivalent
classical GAN. We anticipate that our observation of a fundamental connection
between modelling correlated distributions and quantum entanglement will
become the standard technique for modeling multivariate data with quantum
computers, and will be used as a basis for a wide variety of learning
applications.

##  11 Acknowledgement

This work is a collaboration between Fidelity Center of Applied Technology,
Fidelity Investments and IonQ Inc.

##  12 QCBM Training

Here we discuss details of the QCBM training. Fig. 10 shows the effect of
changing the number of shots $N_{shots}$ on the simulator. For $N_{q}=2$, the
cost function converges faster for fewer shots while this is reversed for
$N_{q}=6$. This further supports our observation in the text that for the QCBM
at small system sizes, a noisy estimate of the cost function can play a
favorable role in the optimization. For $N_{q}=6$, one can also observe that
the spread of the values towards the end of the optimization decreases due to
the increased accuracy available from larger number of shots.

Figure 10: Effect of increasing shots on the convergence of the QCBM for
different system sizes

Fig. 11 shows the marginal distributions obtained at the minimum value of the
cost function achieved during optimization. Theoretically, the marginals
should be uniform within statistical error. This is seen to be true even on
the quantum hardware, revealing one reason for the robustness of the algorithm
to noise.

Figure 11: Marginal distributions in the copula space on the three different
backends used. The red circles correspond to the first variable and the blue
triangles correspond to the second. The width of the green shaded bands is the
theoretical statistical uncertainty which is $1/\sqrt{N_{shots}}$ and they are
centered at $1/2^{N_{q}/2}$, where $N_{q}=$ 2, 4 and 6 qubits from top to
bottom in each subfigure. Here $N_{shots}=500$.

##  13 Argument on Computational Complexity

Here we argue in detail the advantage via computational complexity by
introducing post-selection. A quantum circuit with post-selection has output
registers $O$ and post-selection registers $P$. Instead of sampling
measurement results $x$ directly from measurement on the output registers, we
consider only those runs of the process for which a measurement on the post-
selection registers $P$ yield $0\dots 0$. In this construction, we require the
circuit to have the property $\text{Pr}[P=0\dots 0]\neq 0$ so that the
conditional disstribution $\text{Pr}[O=x|P=0\dots 0]$ is well-defined.

Since post-selecting the “qopula” circuit with $B=0\dots 0$ allows us to
generate distributions identical to that of IQP circuits, we see that IQP
$\subseteq$ post-Qopula (the complexity class of “qopula” circuits with post-
selection) and therefore post-IQP $\subseteq$ post-Qopula $\subseteq$ post-BQP
(Bounded error Quantum Polynomial time). Since post-IQP = post-BQP [52], we
get post-Qopula = post-BQP. If the output distribution of “qopula” circuits
can be sampled efficiently classically, we would have Qopula $\subseteq$ BPP
(Bounded error Probabilistic Polynomial time) and therefore post-Qopula
$\subseteq$ post-BPP. This would imply post-BQP $\subseteq$ post-BPP, which
would result in the collapse of the polynomial hierarchy to the third level.

Since the collapse of the polynomial hierarchy is widely believed to be
unlikely, this means the output distribution of “qopula” circuits cannot be
sampled efficiently classically.

## References

  * [1] G. James, D. Witten, T. Hastie & R. Tibshirani, “An Introduction to Statistical Learning: with Applications in R”, Springer (2013)
  * [2] H. Markowitz, “Portfolio Selection”, The Journal of Finance 7, 77 (1952)
  * [3] S. Zacks, “Introduction to Reliability Analysis: Probability Models and Statistical Methods”, Springer-Verlag New York (1992)
  * [4] C. C. Aggarwal, “Recommender Systems: The Textbook”, 1st edition, Springer Publishing Company, Incorporated (2016)
  * [5] H. v. Storch & F. W. Zwiers, “Statistical Analysis in Climate Research”, Cambridge University Press (1999)
  * [6] C. Palmer, “Statistics of Medical Imaging”, Health Physics 104, (2013)
  * [7] R. Aloui, M. S. B. Aissa & D. K. Nguyen, “Global Financial Crisis, Extreme Interdependences, and Contagion Effects: The Role of Economic Structure?”, Journal of Banking and Finance 35, 130 (2011)
  * [8] D. Kenourgios, A. Christopoulos & D. I. Dimitriou, “Asset Markets Contagion During the Global Financial Crisis”, Mutual Funds , (2013)
  * [9] F. Serinaldi, A. Bárdossy & C. G. Kilsby, “Upper tail dependence in rainfall extremes: would we know it if we saw it?”, Stochastic Environmental Research and Risk Assessment 29, 1211 (2015)
  * [10] M. Sklar, “Fonctions de repartition a n dimensions et leurs marges”, Publications de l’Institut Statistique de l’Université de Paris (1959)
  * [11] R. W. van den Goorbergh, C. Genest & B. J. Werker, “Bivariate option pricing using dynamic copula models”, Insurance: Mathematics and Economics 37, 101 (2005), Papers presented at the DeMoSTAFI Conference, Québec, 20-22 May 2004
  * [12] R. T. Kilgore & D. B. Thompson, “Estimating Joint Flow Probabilities at Stream Confluences by Using Copulas”, Transportation Research Record 2262, 200 (2011)
  * [13] J. Lapuyade-Lahorgue, J.-H. Xue & S. Ruan, “Segmenting Multi-Source Images Using Hidden Markov Fields With Copula-Based Multivariate Statistical Distributions”, IEEE Transactions on Image Processing 26, 3187 (2017)
  * [14] H. Joe, “Multivariate Models and Multivariate Dependence Concepts”, Chapman and Hall/CRC (1997), New York
  * [15] I. Goodfellow, J. Pouget-Abadie, M. Mirza, B. Xu, D. Warde-Farley, S. Ozair, A. Courville & Y. Bengio, “Generative Adversarial Nets”, in “Advances in Neural Information Processing Systems”
  * [16] D. P. Kingma & M. Welling, “Auto-Encoding Variational Bayes”, in “2nd International Conference on Learning Representations, ICLR 2014, Banff, AB, Canada, April 14-16, 2014, Conference Track Proceedings”, ed: Y. Bengio & Y. LeCun
  * [17] A. Antoniou, A. Storkey & H. Edwards, “Data Augmentation Generative Adversarial Networks”, arXiv:1711.04340
  * [18] J. Wu, C. Zhang, T. Xue, B. Freeman & J. Tenenbaum, “Learning a Probabilistic Latent Space of Object Shapes via 3D Generative-Adversarial Modeling”, in “Advances in Neural Information Processing Systems 29: Annual Conference on Neural Information Processing Systems 2016, December 5-10, 2016, Barcelona, Spain”, 82–90
  * [19] C. Ledig, L. Theis, F. Huszár, J. Caballero, A. Cunningham, A. Acosta, A. Aitken, A. Tejani, J. Totz, Z. Wang & W. Shi, “Photo-Realistic Single Image Super-Resolution Using a Generative Adversarial Network”, in “2017 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)”, 105-114
  * [20] S. Arora, A. Risteski & Y. Zhang, “Do GANs learn the distribution? Some Theory and Empirics”, in “International Conference on Learning Representations”
  * [21] “Common problems — generative adversarial networks — google developers”, Google , , https://developers.google.com/machine-learning/gan/problems
  * [22] S. Aaronson & A. Arkhipov, “The Computational Complexity of Linear Optics”, in “Proceedings of the Forty-Third Annual ACM Symposium on Theory of Computing”, 333–342, Association for Computing Machinery (2011), New York, NY, USA
  * [23] S. Boixo, S. V. Isakov, V. N. Smelyanskiy, R. Babbush, N. Ding, Z. Jiang, M. J. Bremner, J. M. Martinis & H. Neven, “Characterizing quantum supremacy in near-term devices”, Nature Physics 14, 595 (2018)
  * [24] F. Arute, K. Arya, R. Babbush, D. Bacon, J. C. Bardin, R. Barends, R. Biswas, S. Boixo, F. G. S. L. Brandao, D. A. Buell, B. Burkett, Y. Chen, Z. Chen, B. Chiaro, R. Collins, W. Courtney, A. Dunsworth, E. Farhi, B. Foxen, A. Fowler, C. Gidney, M. Giustina, R. Graff, K. Guerin, S. Habegger, M. P. Harrigan, M. J. Hartmann, A. Ho, M. Hoffmann, T. Huang, T. S. Humble, S. V. Isakov, E. Jeffrey, Z. Jiang, D. Kafri, K. Kechedzhi, J. Kelly, P. V. Klimov, S. Knysh, A. Korotkov, F. Kostritsa, D. Landhuis, M. Lindmark, E. Lucero, D. Lyakh, S. Mandrà, J. R. McClean, M. McEwen, A. Megrant, X. Mi, K. Michielsen, M. Mohseni, J. Mutus, O. Naaman, M. Neeley, C. Neill, M. Y. Niu, E. Ostby, A. Petukhov, J. C. Platt, C. Quintana, E. G. Rieffel, P. Roushan, N. C. Rubin, D. Sank, K. J. Satzinger, V. Smelyanskiy, K. J. Sung, M. D. Trevithick, A. Vainsencher, B. Villalonga, T. White, Z. J. Yao, P. Yeh, A. Zalcman, H. Neven & J. M. Martinis, “Quantum supremacy using a programmable superconducting processor”, Nature 574, 505 (2019)
  * [25] S. Cheng, J. Chen & L. Wang, “Information Perspective to Probabilistic Modeling: Boltzmann Machines versus Born Machines”, Entropy 20, (2018)
  * [26] M. Benedetti, D. Garcia-Pintos, O. Perdomo, V. Leyton-Ortega, Y. Nam & A. Perdomo-Ortiz, “A generative modeling approach for benchmarking and training shallow quantum circuits”, npj Quantum Information 5, 45 (2019)
  * [27] P.-L. Dallaire-Demers & N. Killoran, “Quantum generative adversarial networks”, Phys. Rev. A 98, 012324 (2018)
  * [28] S. Lloyd & C. Weedbrook, “Quantum Generative Adversarial Learning”, Phys. Rev. Lett. 121, 040502 (2018)
  * [29] D. Zhu, N. M. Linke, M. Benedetti, K. A. Landsman, N. H. Nguyen, C. H. Alderete, A. Perdomo-Ortiz, N. Korda, A. Garfoot, C. Brecque, L. Egan, O. Perdomo & C. Monroe, “Training of quantum circuits on a hybrid quantum computer”, Science Advances 5, (2019)
  * [30] C. Zoufal, A. Lucchi & S. Woerner, “Quantum Generative Adversarial Networks for learning and loading random distributions”, npj Quantum Information 5, 103 (2019)
  * [31] L. Hu, S.-H. Wu, W. Cai, Y. Ma, X. Mu, Y. Xu, H. Wang, Y. Song, D.-L. Deng, C.-L. Zou & L. Sun, “Quantum generative adversarial learning in a superconducting quantum circuit”, Science Advances 5, (2019)
  * [32] H.-L. Huang, Y. Du, M. Gong, Y. Zhao, Y. Wu, C. Wang, S. Li, F. Liang, J. Lin, Y. Xu, R. Yang, T. Liu, M.-H. Hsieh, H. Deng, H. Rong, C.-Z. Peng, C.-Y. Lu, Y.-A. Chen, D. Tao, X. Zhu & J.-W. Pan, “Experimental Quantum Generative Adversarial Networks for Image Generation”, arXiv:2010.06201
  * [33] M. S. Rudolph, N. B. Toussaint, A. Katabarwa, S. Johri, B. Peropadre & A. Perdomo-Ortiz, “Generation of High-Resolution Handwritten Digits with an Ion-Trap Quantum Computer”, arXiv:2012.03924
  * [34] L. K. Grover, “A Fast Quantum Mechanical Algorithm for Database Search”, in “ANNUAL ACM SYMPOSIUM ON THEORY OF COMPUTING”, 212–219, ACM (1996)
  * [35] P. W. Shor, “Polynomial-Time Algorithms for Prime Factorization and Discrete Logarithms on a Quantum Computer”, SIAM Journal on Computing 26, 1484 (1997)
  * [36] D. P. Kingma & J. Ba, “Adam: A Method for Stochastic Optimization”, in “ICLR (Poster)”
  * [37] J.-G. Liu & L. Wang, “Differentiable learning of quantum circuit Born machines”, Phys. Rev. A 98, 062324 (2018)
  * [38] N. Patki, R. Wedge & K. Veeramachaneni, “The Synthetic Data Vault”, in “2016 IEEE International Conference on Data Science and Advanced Analytics (DSAA)”, 399-410
  * [39] N. Tagasovska, D. Ackerer & T. Vatter, “Copulas as High-Dimensional Generative Models: Vine Copula Autoencoders”, in “Advances in Neural Information Processing Systems”, ed: H. Wallach, H. Larochelle, A. Beygelzimer, F. d' Alché-Buc, E. Fox & R. Garnett, Curran Associates, Inc. (2019)
  * [40] D. Maslov, “Basic circuit compilation techniques for an ion-trap quantum machine”, New Journal of Physics 19, 023035 (2017)
  * [41] M. Enríquez, I. Wintrowicz & K. Życzkowski, “Maximally Entangled Multipartite States: A Brief Survey”, Journal of Physics: Conference Series 698, 012003 (2016), https://doi.org/10.1088/1742-6596/698/1/012003
  * [42] K. Wright, K. M. Beck, S. Debnath, J. M. Amini, Y. Nam, N. Grzesiak, J.-S. Chen, N. C. Pisenti, M. Chmielewski, C. Collins, K. M. Hudek, J. Mizrahi, J. D. Wong-Campos, S. Allen, J. Apisdorf, P. Solomon, M. Williams, A. M. Ducore, A. Blinov, S. M. Kreikemeier, V. Chaplin, M. Keesan, C. Monroe & J. Kim, “Benchmarking an 11-qubit quantum computer”, Nature Communications 10, (2019)
  * [43] M. Schuld, V. Bergholm, C. Gogolin, J. Izaac & N. Killoran, “Evaluating analytic gradients on quantum hardware”, Phys. Rev. A 99, 032331 (2019)
  * [44] J. C. Spall, “A Stochastic Approximation Technique for Generating Maximum Likelihood Parameter Estimates”, in “1987 American Control Conference”, 1161-1167
  * [45] J. Spall, “Multivariate stochastic approximation using a simultaneous perturbation gradient approximation”, IEEE Transactions on Automatic Control 37, 332 (1992)
  * [46] S. Johri, S. Debnath, A. Mocherla, A. Singh, A. Prakash, J. Kim & I. Kerenidis, “Nearest Centroid Classification on a Trapped Ion Quantum Computer”, arXiv:2012.04145
  * [47] J. R. McClean, S. Boixo, V. N. Smelyanskiy, R. Babbush & H. Neven, “Barren plateaus in quantum neural network training landscapes”, Nature Communications 9, 4812 (2018)
  * [48] M. Cerezo, A. Sone, T. Volkoff, L. Cincio & P. J. Coles, “Cost function dependent barren plateaus in shallow parametrized quantum circuits”, Nature Communications 12, 1791 (2021)
  * [49] J. A. Peacock, “Two-dimensional goodness-of-fit testing in astronomy”, Monthly Notices of the Royal Astronomical Society 202, 615 (1983)
  * [50] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, A. Desmaison, A. Kopf, E. Yang, Z. DeVito, M. Raison, A. Tejani, S. Chilamkurthy, B. Steiner, L. Fang, J. Bai & S. Chintala, “PyTorch: An Imperative Style, High-Performance Deep Learning Library”, in “Advances in Neural Information Processing Systems 32”, ed: H. Wallach, H. Larochelle, A. Beygelzimer, F. d' Alché-Buc, E. Fox & R. Garnett, Curran Associates, Inc. (2019), 8024–8035
  * [51] G. Brassard, R. Cleve & A. Tapp, “Cost of Exactly Simulating Quantum Entanglement with Classical Communication”, Phys. Rev. Lett. 83, 1874 (1999)
  * [52] M. J. Bremner, R. Jozsa & D. J. Shepherd, “Classical simulation of commuting quantum computations implies collapse of the polynomial hierarchy”, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 467, 459 (2011)
  * [53] S. Arora & B. Barak, “Computational Complexity: A Modern Approach”, 1st edition, Cambridge University Press (2009), USA
  * [54] M. Arjovsky, S. Chintala & L. Bottou, “Wasserstein Generative Adversarial Networks”, in “Proceedings of the 34th International Conference on Machine Learning”, ed: D. Precup & Y. W. Teh, 214–223, PMLR (2017)
  * [55] S. Chakrabarti, H. Yiming, T. Li, S. Feizi & X. Wu, “Quantum Wasserstein Generative Adversarial Networks”, in “Advances in Neural Information Processing Systems”
  * [56] T. Karras, S. Laine & T. Aila, “A Style-Based Generator Architecture for Generative Adversarial Networks”, in “Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)”
  * [57] E. Farhi & A. W. Harrow, “Quantum Supremacy through the Quantum Approximate Optimization Algorithm”, arXiv:1602.07674
  * [58] M. Hastings, “Classical and quantum bounded depth approximation algorithms”, 2019 Quantum Information & Computation 19, 1116 (2018)
  * [59] N. Wiebe, D. Braun & S. Lloyd, “Quantum Algorithm for Data Fitting”, Phys. Rev. Lett. 109, 050505 (2012)
  * [60] P. Rebentrost, M. Mohseni & S. Lloyd, “Quantum Support Vector Machine for Big Data Classification”, Phys. Rev. Lett. 113, 130503 (2014)
  * [61] A. W. Harrow, A. Hassidim & S. Lloyd, “Quantum Algorithm for Linear Systems of Equations”, Phys. Rev. Lett. 103, 150502 (2009)
  * [62] I. Cong, S. Choi & M. D. Lukin, “Quantum convolutional neural networks”, Nature Physics 15, 1273 (2019)
  * [63] E. Farhi & H. Neven, “Classification with Quantum Neural Networks on Near Term Processors”, arXiv:1802.06002 [quant-ph]
  * [64] H.-Y. Huang, R. Kueng & J. Preskill, “Information-Theoretic Bounds on Quantum Advantage in Machine Learning”, Phys. Rev. Lett. 126, 190505 (2021)
  * [65] Y. Liu, S. Arunachalam & K. Temme, “A rigorous and robust quantum speed-up in supervised machine learning”, Nature Physics , (2021)
  * [66] X. Gao, E. R. Anschuetz, S.-T. Wang, J. I. Cirac & M. D. Lukin, “Enhancing Generative Models via Quantum Correlations”, arXiv:2101.08354 [quant-ph]