#  Parallel Quantum Computing Simulations via Quantum Accelerator Platform
Virtualization

Daniel Claudino  Dmitry I. Lyakh  Alexander J. McCaskey

###### Abstract

Quantum circuit execution is the central task in quantum computation. Due to
inherent quantum-mechanical constraints, quantum computing workflows often
involve a considerable number of independent measurements over a large set of
slightly different quantum circuits. Here we discuss a simple model for
parallelizing simulation of such quantum circuit executions that is based on
introducing a large array of virtual quantum processing units, mapped to
classical HPC nodes, as a parallel quantum computing platform. Implemented
within the XACC framework, the model can readily take advantage of its
backend-agnostic features, enabling parallel quantum circuit execution over
any target backend supported by XACC. We illustrate the performance of this
approach by demonstrating strong scaling in two pertinent domain science
problems, namely in computing the gradients for the multi-contracted
variational quantum eigensolver and in data-driven quantum circuit learning,
where we vary the number of qubits and the number of circuit layers. The
latter (classical) simulation leverages the cuQuantum SDK library to run
efficiently on GPU-accelerated HPC platforms.

###### keywords:

Quantum computing , Quantum software , Distributed computing

††journal: Future Generation Computer Systemst1t1footnotetext: This manuscript
has been authored by UT-Battelle, LLC, under Contract No. DE-AC0500OR22725
with the U.S. Department of Energy. The United States Government retains and
the publisher, by accepting the article for publication, acknowledges that the
United States Government retains a non-exclusive, paid-up, irrevocable, world-
wide license to publish or reproduce the published form of this manuscript, or
allow others to do so, for the United States Government purposes. The
Department of Energy will provide public access to these results of federally
sponsored research in accordance with the DOE Public Access Plan.

[inst1]organization=Quantum Information Science Section, Oak Ridge National
Laboratory,addressline=1 Bethel Valley Road, city=Oak Ridge, postcode=37831,
state=TN, country=USA

[inst2]organization=NVIDIA Corporation,city=Santa Clara, state=CA, country=USA

##  1 Introduction

Quantum computing is largely an umbrella term that refers to several
remarkably distinct architectural, logical, and computational approaches.
However, this picture is dominated by the digital, gate-based model of quantum
computation, and while there are theoretical guarantees that other competing
modes are equally valid, such as adiabatic quantum computing, [1] those remain
a much more niche technology. One trademark of the digital quantum computers
is that, despite their striking difference from conventional computers, they
retain the ability to enable some quantum operations to behave as analogs of
conventional logic gates (for example, the classical NOT gate and the quantum
X gate).

At the same time, quantum computers exhibit many unique characteristics that
are inherently quantum, and an important example of which is quantum
measurement. The state of a closed quantum system can be fully defined by its
wave function, a mathematical object that has information about the behavior
of the quantum system at every point in space. A critical complication arises
from the fact that one is denied direct access to such an object, and is
limited to drawing indirect conclusions based on quantum measurements. In
turn, these measurements yield classical outcomes according to the
corresponding observable operator through which the quantum state is measured,
but which ends up affecting the original quantum state. This is a process
known as the collapse of the wave function, which practically means that the
quantum state, initially in a superposition of the eigenstates of the
observable operator, upon measurement collapses into one of those eigenstates,
revealing the probabilistic and destructive nature of the quantum measurement.
Thus, in a general case where a quantum state is not prepared as an eigenstate
of the quantum operator to be observed/measured, the measurement leads to a
distribution of probabilities over its eigenstates. The prepared quantum state
then needs to be resolved in such an eigenbasis, in which case many
measurements need to take place in order to render a reliable statistics, each
requiring a quantum state to be prepared “from scratch”. It is important to
note the output information from this process is classical in nature. To
illustrate this point, in the common case of digital quantum computers that
operate on qubits which can store an arbitrary quantum superposition of two
basis states, the outcome of measurements on each qubit is still classical,
either 0 or 1. Thus, quantum algorithms ultimately still produce classical
information.

There are various ways the information in a quantum state can be encoded and
manipulated, resulting in a multitude of simulation techniques whose inner
workings follow those principles. The state can be naturally represented as a
vector in an exponentially large Hilbert space, with the density matrix
formalism requiring an even larger space, such that both can in principle be
exact. Alternatively, there are much more computationally amenable
representations of quantum states such as tensor network and stabilizer states
which are applicable to certain classes of cases. Recognizing that each
available mode of simulation has strengths and weaknesses means that one would
ideally have a flexible framework that would easily grant access to this
myriad of simulators while abstracting their intricacies. This is the
foundational design argument in the XACC quantum-classical computing
framework, [2, 3] which puts forth a modular infrastructure for hardware-
agnostic compilation and simulation of quantum circuits.

A common pitfall shared by all the different representations above is that
general simulation of quantum states by classical means demands an exponential
amount of resources, e.g., time and/or memory. Fortunately, this problem can
be partially alleviated by turning to high-performance computing (HPC), which
has become an invaluable tool in extending the scale and scope of how
classical computers can aid in quantum computing research. An example of this
is the validation of quantum protocols in which such simulations provide a
setting where the deleterious effects of noise can be carefully modeled, or
altogether removed from the computation. If, on the one hand, the nature of
quantum measurements implies repeating seemingly the same tasks many times
over, on the other hand, one can exploit the fact that all instances of these
tasks are independent from one another, alluding to carrying out those
measurements in an embarrassingly parallel fashion, which has become
commonplace in HPC. This realization is one of the main motivations behind our
proposed parallel quantum computing workflow.

In terms of programming and execution models, quantum computing has displayed
a similar trajectory to that of accelerated computing based on the graphical
processing units (GPU). That is, both approaches imply a heterogeneous
computing system of which some parts serve as accelerators called to handle a
limited subset of computational tasks that play to their strengths. While this
paradigm has already proven extremely successful, heterogeneity increases
complexity in various levels of the stack. This can be offset by adopting
suitable abstractions that model the interaction between heterogeneous
components, one example of which is compute resource virtualization [4].
Similar to the powerful quantum-classical programming framework CUDA-Q from
NVIDIA [5], in our approach each compute node (CPU or GPU) is modeled after a
virtual quantum processing unit (QPU). In other words, QPUs are simulated as
quantum virtual machines that are executed by classical compute resources,
i.e., CPUs and GPUs. In this picture, we leverage the standard message passing
interface (MPI) orchestrate classical communication between virtual QPUs, thus
being able to compose the entire quantum-classical parallel workflows.

One of the main technological limitations in quantum computing refers to the
narrow time window during which assertions regarding the state prepared by a
quantum computer can be made reliably, which is referred to as coherence time
and which prevents quantum circuits from implementing a large number of gates.
In various applications, this can be circumvented at the expense of an
increasing number of circuit repetitions. This trade-off can be exploited in a
variety of ways, and is at the heart of many quantum algorithms and
subroutines, which is one of the strongest arguments in favor of a parallel or
distributed quantum computing workflow. One of the most promising use cases of
quantum computing, and which was the motivating application behind Feynman’s
argument,[6] is the simulation of physical systems, with the leading algorithm
being the variational quantum eigensolver.[7] This algorithm alleviates the
circuit depth requirement of more precise alternatives (i.e. quantum phase
estimation) by performing classical optimization of quantum circuit
parameters, where the search for optimal parameters is informed by the
classical outcomes of a much larger number of measurements. Similarly, circuit
complexity can be lowered both in terms of the depth and width (number of
qubits) by increasing the number of measurements.[8, 9] One prominent example
is that of circuit cutting techniques,[10, 11, 12] which can lower the circuit
depth requirements by cutting circuits into smaller, simpler ones, which
incurs an exponential measurement overhead.

Against this backdrop, the design rationale behind XACC not only embodies the
versatility required to ensure performance across these different modalities,
granted in part by its robust plugin infrastructure, but also attains a high
degree of interoperability due to its layered compilation approach. Higher-
level quantum programming languages, and even domain-specific languages, can
be parsed into the XACC intermediate representation (IR). Developers then
target backend-specific languages with compiler plugins that map the XACC IR
onto the instruction set of the desired backend. Within this paradigm, users
can “mix and match” languages and backends, while programming in their
preferred languages and are freed from needing to possess explicit knowledge
regarding compilation and corresponding code optimizations. Looking forward,
this means that XACC is already equipped with the main features required to
manage hybrid, highly interoperable workflows, where computations carried out
by quantum accelerators are embedded in a largely classical computational
pipeline.

While the most immediate application of the proposed execution model would be
the  parallel simulation of quantum circuits with classical numerical
simulators, as discussed above,  over quantum processors would also be
attractive and fully supported by our approach, akin to the ensemble quantum
computing model.[13] With this in mind, our main contributions are:

  * 1.

We virtualize quantum compute resources, which we call virtual QPUs, by
abstracting away specifics of desired numerical simulators and actual quantum
hardware;

  * 2.

Provide a flexible and robust communication framework based on the message
passing interface (MPI) standard;

  * 3.

Demonstrate our approach’s scalability via important domain science problems
deploying simulations using both CPUs and GPUs, in two different HPC systems,
and with distinct numerical simulator backends.

The rest of the paper is organized as follows. We first review previous
related work in Section 2. In Section 3 we present the most relevant design
elements of our virtualization scheme and some basic aspects of our
implementation. In Section 4 we discuss the performance and scalability of our
approach in light of two relevant problems in the domain of quantum simulation
and quantum machine learning. Finally, in Section 5 we summarize our main
findings and consider potential future ramifications.

##  2 Related Work

The execution of an exponential number of calculations in parallel is at the
heart of the expected speedups in various instances of quantum computing. This
was realized early in the development of quantum computing as natural
execution model to be employed in many quantum algorithms, the quintessential
example being the Shor algorithm for prime number factorization. Later on,
parallel execution in quantum computing gained prominence primarily as a means
to circumvent the exponential demand in memory/time for classical simulation
of quantum systems. The latter modality is where our platform falls into and
the one for which we will discuss similar instances in the literature.

The first example of a parallel simulator in the literature comes from the
work of Obenland and Despain, which was implemented in C and enabled
parallelization via MPI. A distinct feature of this simulator is its
implementation of quantum gates in terms of experimental gate parameters in
order to mimic laser pulses in trapped ions.[14] Another example of a
simulator implemented in C is found in Ref.15, however with parallelization
based on a multi-threaded shared memory model. Additionally, the
implementation of a massively parallel simulator, initially implemented in
Fortran 90 and parallelized via MPI, has undergone several updates and offers
both shared and distributed memory models.[16, 17]

With the emergence of graphical processing units, their intrinsic
architectural paradigm can be naturally leveraged for parallel execution over
their many threads, and here we emphasize the initial works of Gutierrez et
al.[18] and Amariutei and Caraiman,[19] demonstrating parallelization within a
single GPU.

Recent works present simulators designed to be deployed in the HPC
environments. Two scalable, CPU-only simulators[20, 21] have demonstrated
performance for quantum chemistry applications [22, 23]. The NWQ suite of
simulators provides CPU and GPU implementations of the state vector and
density matrix representations of the quantum state.[24, 25, 26, 27] It is
also worth mentioning works in specialized simulators, with the ones based on
tensor networks standing out, such as ExaTN,[28] QTensor,[29] and ITensor.[30]

Perhaps the closest examples to our virtualization platform is CUDA-Q[5] and
QuEST,[31] later renamed to the Intel Quantum Simulator, and currently named
the Intel Quantum SDK. They offer both shared (OpenMP) and distributed (MPI)
memory parallelization and can target both CPU and GPU backends. However, the
Intel Quantum SDK differs from XACC in that the set of backends made available
by the Intel Quantum SDK is currently restricted to classical numerical
simulators, while XACC and CUDA-Q encompass those as well as actual quantum
hardware.

##  3 Design and Implementation

###  3.1 Design

In XACC, quantum hardware is viewed as accelerators embedded in a classical
computational workflow. This model draws inspiration from HPC nodes and how
they leverage GPUs, which is considered a staple of heterogeneous computation.
Within this picture, two important abstractions are in order, namely those
facilitated by the Accelerator and AcceleratorBuffer interfaces. The former
enables access to a target backend to carry out quantum circuit executions.
The AcceleratorBuffer serves two main purposes: it models the underlying
register of quantum bits (qubits) and provides data structures to persist
information resulting from circuit execution, such as bit-string counts, and
hardware specifics metadata like error rates. From a more practical
standpoint, XACC defers circuit execution to concrete implementations of the
Accelerator interface, which can be numerical simulators (IBM’s Aer, Google’s
qsim[32], Quantum++[33], TNQVM,[34, 35] etc.) or quantum computing vendor
hardware (IBM, Rigetti, Quantinuum, etc.).

One of the most distinctive features of XACC is it intermediate representation
(IR), which enables richly composable quantum workflows. A variety of
programming languages can be parsed into the XACC IR as well as circuit
construction via a builder class. At that point, the XACC puts forth a
Compiler interface whose main task is to map the instructions in the XACC IR
onto the native gate set of the desired backend. Once compiled into the set of
instructions and the format expected by the backend of choice, the circuits
are offloaded for execution. The primary result of circuit execution is the
distribution of observed bitstrings over the possible $2^{N}$ outcomes, which
may in turn be further processed in order to provide answers according to the
computational task at hand. The main components of this workflow, which also
serves as a basis for understanding the parallel pipeline, are illustrated in
Figure 1.

Figure 1: Summary of the XACC circuit execution workflow.

The salient features of the design behind our virtualization scheme are
summarized in the diagram in Figure 2, which follow the usual decorator
pattern. While the Accelerator interface prescribes other methods, we suppress
them here because they are not relevant in the current context. Here we
highlight that, for a collection of virtual QPUs over which we intend to
distribute/parallelize the computational workload, each virtual QPU is
assigned ownership of and has access to an independent instance of the chosen
Accelerator, with the two concrete implementations used in this work being
illustrated in Figure 2, namely qsim and Aer. Also of importance is that all
the necessary communication is facilitated by qpuComm attribute in each
virtual QPU, which will be discussed in more detail later on.

Figure 2: The class diagram representing the design of the HPCVirtDecorator.

Before we delve into the details of our implementation, we consider
appropriate to discuss some of the relevant features of MPI in light of the
programming and execution model we propose here. In this context, we judge
essential to avail ourselves of this standard to enable communication and
synchronization among virtual QPUs, which has a proven track-record of
maturity and robustness in HPC environments. One aspect of MPI that is crucial
to our model is the concept of communicators, which represent groups of
processes for MPI provide extensive capabilities that enables efficient and
streamlined communication. The communication functionality relevant to the
current work is encapsulated by ProcessGroup (Figure 2), which is a wrapper to
native MPI calls (in C) with the intention of simplifying message passing and
communication management. In our framework, each virtual QPU is modeled after
an MPI communicator, allowing processes associated with the same virtual QPU
to exchange data and coordinate execution. Another important concept is MPI
collective communication operations, of which we highlight MPI_Bcast() and
MPI_Allgatherv(), which are used for broadcasting data from one process to all
others or gathering data from all processes, respectively. These operations
are employed to consolidate execution results obtained from individual virtual
QPUs that become accessible to all processes. Because several MPI routines are
non-blocking in the most general sense of the term, synchronization is ensured
by MPI_Barrier(). By leveraging MPI’s communication primitives, our framework
achieves efficient and scalable parallel execution of quantum circuits across
heterogeneous computing resources.

In order to facilitate the discussion, we establish that data / operations
done at the virtual QPU level are referred to as local, while those that take
place across/are shared by all virtual QPUs are global. Once the Accelerator
offloads the quantum circuits to the quantum backend of choice (Figure 1),
barring trivial cases where the number of virtual QPUs is either set to 1 or
is less than the number of available MPI processes, all available MPI
processes are assigned to a certain virtual QPUs and identified by a color.
Furthermore, each process associated with a given color has access to the
corresponding qpuComm object (Figure 2). We assigned one process under each
qpuComm, conventionally the one with rank 0, to be responsible for managing
all the necessary communication. In a similar way as the MPI processes are
grouped according to colors, each virtual QPU is assigned a subset of the
entire collection of quantum circuits to be executed. A local
AcceleratorBuffer is instantiated to collect and persist the outcomes of
execution called by the corresponding virtual QPU and offloaded to the quantum
backend of choice. At this point, each virtual QPU is equipped with all the
required ingredients to carry out the required local circuit executions, whose
outcomes are generally probability distributions, but can also be expectation
values (of type double), in the case of numerical simulators. These results
are broadcast in order to be consolidated under an AcceleratorBuffer that is
shared among all virtual QPUs, thus enabling all MPI processes to have access
to all the execution data. This is summarized in Figure 3.

Figure 3: Main data structures and communication pattern across virtual QPUs.
The white arrows represent a call to the Accelerator instance in each virtual
QPU to offload quantum circuits to the backend, and the two rows of virtual
QPUs represent different states of the same virtual QPUs before (top) and
after (bottom) communication.

###  3.2 Implementation in the XACC framework

In this section we discuss some pertinent details of the virtualization
platform in the XACC framework, which is facilitated by turning to the example
in Figure 4 that illustrates its main features. We bring attention to the fact
that XACC is entirely developed in modern C++.

⬇

1using namespace xacc;

2// get reference to Accelerator

3auto accelerator = getAccelerator("qsim");

4// allocate a register of nBits qubits

5buffer = xacc::qalloc(nBits);

6// execute std::vector<CompositeInstruction*>

7// in serial manner

8accelerator->execute(buffer, circuits);

9// now turn to parallel execution

10// options is a HeterogeneousMap data structure

11// that maps n-virtual-qpus to # of virtual QPUs

12auto options = {{"n-virtual-qpus", 256}};

13// decorate Accelerator

14accelerator = getAcceleratorDecorator(

15 "hpc-virtualization",

16 accelerator, options);

17// decorated accelerator executes circuits

18accelerator->execute(buffer, circuits);

Figure 4: Code snippet showcasing how to execute a collection of circuits in
XACC using the qsim[32] simulator and how to enable parallel execution via the
HPCVirtDecorator.

In order to tie back to the previous section, the Accelerator interface
enables access to a target backend via the execute() method, which takes as
arguments a quantum circuit (or an std::vector thereof) and a
AcceleratorBuffer∗, whose function is examined in Section 3.1. In XACC, a
quantum circuit is embodied by an instance of CompositeInstruction, which
provides the abstraction for a cohesive set Instructions (e.g. quantum gates).
A collection of CompositeInstructions can be executed via
Accelerator::execute(), as shown in Figure 4.

Normally, the call to Accelerator::execute() would serialize the execution by
iterating over the contents of the
std::vector<CompositeInstruction∗> (line 8 in Figure 4). In order to enable
parallel quantum circuit execution, the original Accelerator::execute()
functionality is augmented based on the decorator pattern (HPCVirtDecorator).
We note that XACC makes extensive use of decorators as a very convenient means
to abstract complex pre- and post-processing protocols of quantum circuit
execution data. The usage of the proposed HPCVirtDecorator is illustrated in
Figure 4, lines 12-18. The options object (line 12) is a HeterogenousMap,
which is a flexible data structure that can map strings to any value type (via
std::any<T>), and through which one passes various parameters to customize
execution. Among these parameters, which the implementation elects to be
required or optional, we illustrate a map with just n-virtual-qpus by which we
pass to HPCVirtDecorator the number of virtual QPUs (256 in our example in
Figure 4, line 12). This parameter’s default value is 1, in which case the
trivial, serial circuit execution ensues.

In passing, we point that the getAcceleratorDecorator wrapper takes the name
of the decorator, and the Accelerator instance to be decorated, and a
HeterogeneousMap with parameters, e.g., number of virtual QPUs. Parallel
circuit execution follows by retrieving HPCVirtDecorator from the XACC service
registry via the string identifier hpc-virtualization and calling its
execute() method. Moreover, the xacc namespace provides the getAccelerator and
the
getAcceleratorDecorator wrappers (lines 8 and 14, respectively) that obviates
the use of templated functions explicitly.

As mentioned above, all the communication across virtual QPUs is handled by
the message passing interface (MPI) standard. With the XACC framework being
entirely developed in C++, we are faced with choice of either linking the
HPCVirtDecorator against the Boost implementation of MPI or to carry out all
the MPI calls in C. We chose the latter, as MPICH or OpenMPI are both adherent
to the MPI specification and updates tend to be more timely. That being said,
that modeling of each virtual QPU after an MPI communicator is facilitated by
a set of helper functions that wrap the relevant C calls under ProcessGroup
and thus enable utilization of C++ smart pointers (qpuComm in Figure 2).

At this point, the results from execution in each virtual QPU need to
consolidated, which in principle could be achieved via broadcast
(MPI_Bcast()). However, the results of execution of quantum circuits on
digital quantum computers are akin to a hash table of bit-strings and the
number of times they occur (counts), which in C++ can be simply embodied by
the std::map data structure, with no analog in C. This is a drawback of opting
to C MPI calls instead of relying on the Boost implementation. Each local
buffer (AcceleratorBuffer), when handed over to Accelerator::execute() stores
execution information of each unique circuit in a child buffer. For our
purposes, these children buffer persists three basic pieces of execution
information that we are concerned with: 1) a unique std::string identifier,
whose name defaults to the memory address the child buffer points to and is
often overridden by the string representation of the observable being
measured; 2) a hash table of bit-strings to counts, that is,
std::map<std::string, int>; 3) (Optional) the expectation value associated
with the measurement, which is relevant in certain applications and is of type
double.

We chose to circumvent the limitations arising from turning to a C
implementation of MPI by expressing all the information in the children
buffers in terms of MPI native data types. In other words, this means using
MPI_INT, MPI_DOUBLE, and MPI_CHAR in lieu of int, double, and std::string,
respectively. Apart from the communication involved in sending and receiving
execution data, this also entails explicit knowledge of the amount of data
being communicated. Once the necessary communication is carried out across the
virtual QPUs, each MPI communicator has access to all execution data and then
are in position to rebuild the global buffer (buffer in line 18) as though the
quantum circuits were executed serially. This is illustrated by the arrows and
the second row in Figure 3.

Finally, we comment on an implementation detail that may be relevant depending
on the intended application. By the end of HPCVirtDecorator::execute() call,
we impose synchronization across all virtual QPUs and the global buffer is now
accessible by all processes. In this case, the global buffer differs from that
of trivial execution by storing the corresponding MPI rank of the process.
This is important because in several contexts the post-processing that follows
from circuit execution can itself be data and compute intensive, thus by
enforcing post-processing to take place only from the data associated with a
certain rank, e.g., 0, we can avoid redundant computations and data transfer.

##  4 Performance evaluation

###  4.1 Multi-contracted variational quantum eigensolver

Many applications of quantum computing revolve around estimating the
expectation value of some quantum mechanical observable. The prime example of
this is the estimation of the energy eigenspectrum of a physical system, which
revolves around measurements of the quantum circuit via the Hamiltonian
operator $\hat{H}$, with emphasis in the ground state energy $E_{0}$. In the
context of quantum computing, the most natural representation of such
operators is in terms of Pauli spin operators. For complex quantum systems,
the Hamiltonian can be comprised of a large number of tensor products of such
Pauli spin operators, e.g., $\mathcal{O}(N^{4})$ for the Hamiltonian of a
collection of interacting electrons in a molecule. Each of these terms in
$\hat{H}$ requires a quantum circuit to be measured, an undertaking that can
greatly benefit from a parallel execution model.

A proposed means to circumvent current quantum hardware limitations (e.g.,
circuit depth) in ground state energy estimation is the ubiquitous variational
quantum eigensolver (VQE) algorithm.[7] Briefly, one choose a circuit
implemented by the unitary $\hat{U}(\vec{\theta})$ where $\vec{\theta}$ is a
series of parameters, generally single-qubit rotation gate angles, that are
varied until a cost function, which here is the expectation value of
$\hat{H}$, is minimized. This is summarized in Equation 1

$E_{0}\leq\langle\hat{H}\rangle(\vec{\theta})=\min_{\vec{\theta}}\langle
0|\hat{U}^{\dagger}(\vec{\theta})\hat{H}\hat{U}(\vec{\theta})|0\rangle.$ (1)

The multi-contracted variant of the VQE algorithm (MC-VQE)[36] specializes in
the absorption spectrum of a collection of weakly-interacting “monomers”,
mapping their ground and first excited states to the $|0\rangle$ and
$|1\rangle$ states of a qubit. The salient features of this model can be
successfully captured by the Hamiltonian $\hat{H}$ in Equation 2

$\displaystyle\hat{H}$
$\displaystyle=\mathcal{E}+\sum_{\text{A}}\mathcal{X}_{\text{A}}\hat{X}_{\text{A}}+\mathcal{Z}_{\text{A}}\hat{Z}_{\text{A}}$
$\displaystyle+\sum_{\text{A,
B}}\mathcal{XX}_{\text{AB}}\hat{X}_{\text{A}}\hat{X}_{\text{B}}+\mathcal{XZ}_{\text{AB}}\hat{X}_{\text{A}}\hat{Z}_{\text{B}}$
$\displaystyle+\mathcal{ZX}_{\text{AB}}\hat{Z}_{\text{A}}\hat{X}_{\text{B}}+\mathcal{ZZ}_{\text{AB}}\hat{Z}_{\text{A}}\hat{Z}_{\text{B}},$
(2)

where $\hat{X}$ and $\hat{Z}$ are the Pauli X and Z spin operators, and the
coefficients, e.g., $\mathcal{X}_{\text{A}}$, are real numbers.

The target state builds upon a classical configuration interaction singles
approximation ($\psi_{\text{CIS}}$).[37] At this point, the qubit register
holds a W state generalized to the desired number of qubits.[38] Entanglement
necessary to go beyond this approximation is introduced via entanglers
$\hat{U}_{\text{AB}}(\vec{\theta})$ that are comprised of 6 single-qubit
$R_{y}$ rotation gates and act between qubits associated with adjacent
monomers via two CNOT gates. The circuit and entangler gate layouts are shown
in Figure 5.

Figure 5: Layouts of the circuit and chosen entangler used in the MC-VQE
algorithm.

Each monomer has its own $\psi_{\text{CIS}}$ state preparation step, but they
all share the same entanglers $\hat{U}_{\text{AB}}(\vec{\theta})$. The central
task in the MC-VQE algorithm is to minimize the energy average over all
monomers, by minimizing the cost function in Equation 1 with the Hamiltonian
in 2. Because $\psi_{\text{CIS}}$ is often a good first approximation and the
final energy is often downhill from the energy of that initial state, it is
desirable to leverage gradient-based optimization. This can readily benefit
from application of the parameter-shift rule[39], which in this case involves
$\mathcal{O}(2N_{H}N_{\theta})$ circuit executions per parameter update, where
$N_{H}$ is the number of terms in the Hamiltonian (Equation 2), and
$N_{\theta}$ is the size of $\vec{\theta}$. Important metrics related to the
parameter update in this algorithm for a range of qubits are displayed in
Table 1. It is important to note that here we take advantage of the fact that
$R_{y}(\theta)R_{y}(\phi)=R_{y}(\theta+\phi)$ to merge adjacent $R_{y}$ gates,
thus reducing the number of variational parameters to be optimized.

Table 1: Number of monomers/qubits ($N_{q}$), number of terms in the Hamiltonian $N_{H}$, number of variational parameters $N_{\theta}$, and number of circuits in a single parameter update. $N_{q}$ | $N_{H}$ | $N_{\theta}$ | Number of circuit executions
---|---|---|---
16 | 92 | 76 | 13984
18 | 104 | 86 | 17888
20 | 116 | 96 | 22272
22 | 128 | 106 | 27136

Our scheme is implemented in C++ as is the XACC framework. Python users can
also take advantage of this functionality, along with much of what XACC has to
offer, via bindings facilitated by Pybind11[40]. We investigate the strong
scaling related to the computation of MC-VQE gradients using the parameter-
shift rule where qsim[32] serves as the XACC virtual backend (numerical,
noiseless statevector simulator). In this study, each virtual QPU comprises of
a compute node in the Oak Ridge Computing Facility (OLCF) Andes system, each
with 2 AMD EPYC 7302 16 Core Processor 3.0 GHz, totalling 32 cores per node.
Scaling plots based on the parameters in Table 1 are presented in Figure 6
with runtimes from both native C++ executables and also relying on the
corresponding Python bindings.

Figure 6: Strong scaling plots of the MC-VQE gradients for 16-22 qubits for a
single parameter update.

As the number of qubits increase, together with number of circuits, we
observed a trend that consistently better approximates the expected linear
behaviour. There is a small but noticeable disadvantage when relying on the
Python interface in comparison with exclusive C++, even though the overall
qualitative pattern emerges in either case, pointing to the introduction of
some overhead in the Python case. This overhead is becomes negligible as the
number of qubits increase, pointing to the fact the runtime in this case is
largely dominated by the execution of more demanding circuits.

In light of the discussion in Section 3, the vector containing the circuits in
Table 1 (third column) is split into as many segments as the desired number of
virtual QPUs, here determined by the $x$-axis in Figure 6 and which are in
turn offloaded to its instance of the qsim simulator. For a given choice of
the circuit parameters $\bar{\theta}$ in Equation 1 and Figure 5, and casting
Equation 2 as a sum of generic terms $\hat{P}_{i}$ weighed by scalars $c_{i}$,
the final outcome of quantum circuit execution in this case is to estimate the
expectation values $\langle
0|\hat{U}^{\dagger}(\bar{\theta})\hat{P}_{i}\hat{U}(\bar{\theta})|0\rangle$,
which are scalars. While in reality this requires each expectation value to be
resolved in the eigenspectrum of the corresponding $\hat{P}_{i}$, a numerical
simulator can directly evaluate such a scalar, which we do here. Thus, each
local buffer stores a collection of children buffers, each with the string
representation of the $\hat{P}_{i}$ it measured and the resulting expectation
value.

The local information now needs to be consolidated. We assigned the rank-0 MPI
process in each qpuComm communicator to handle the necessary communication.
The local information, i.e., expectation values and buffer names, are gathered
(via MPI_Allgather() and MPI_Allgatherv()) into vectors of the respective
types, which in turn each qpuComm broadcasts to its subordinate process (via
MPI_Bcast()). At this point all MPI processes have access to all the execution
data, which are used to populate the parent buffer that was originally passed
to the execute(). The copy of the parent buffer managed by each process also
stores its MPI rank (int) as metadata, and synchronization across all
communicators is enforced prior to the end of the execute() call.

###  4.2 Data-driven circuit learning

The use case presented in Subsection 4.1 is primarily concerned with the
computation of expectation values, which in the case of execution on a quantum
circuit simulator can be evaluated directly. More generally, the output of
execution of quantum circuits on digital quantum hardware is the probability
distribution over all possible $2^{N}$ bit-strings, which can give insight
into the prepared quantum state and can, in turn, be used to estimate
expectation values.

Various applications of quantum algorithms revolve around how well one can
match a target probability distribution. One of such approaches is the data-
driven circuit learning algorithm (DDCL).[41] The goal is to learn a set of
rotation angles $\vec{\theta}$ that parameterize a quantum circuit such that
the output distribution $P_{\theta}$ minimizes the Jensen–Shannon (JS)
divergence[42, 43] for a given target distribution $P$ as shown in Equation 3

$\text{JS}(P|P_{\theta})=\frac{1}{2}\sum_{b}\left(P(b)\log\frac{P(b)}{M}+P_{\theta}(b)\log\frac{P_{\theta}(b)}{M}\right),$
(3)

with $M=(P+P_{\theta})/2$ and $P(b)$ being the probability ascribed to the
bit-string $b$ in the distribution $P$.

Since $P_{\theta}$ can be estimated by measuring all qubits in the $Z$-basis,
it involves the measurement of a single Pauli string, and the total cost per
gradient estimate scales as $\mathcal{O}(2N_{\theta})$. The circuit starts
with each pair of neighboring qubits being entangled as a Bell pair. This is
followed by layers of gates where the angles $\vec{\theta}$ appear in the
layout portrayed in Figure 7, partially inspired by Ref. 43.

Figure 7: Layout of the layer used in the DDCL algorithm. Table 2: Number of qubits ($N_{q}$) and number of variational parameters $N_{\theta}$. The number of circuits implementations in a single parameter update is $2N_{\theta}$. $N_{q}$ | $N_{\theta}$ | Number of circuit executions
---|---|---
20 | 1200 | 2400
22 | 1320 | 2640
24 | 1440 | 2880
26 | 1560 | 3120

We again monitor the time associated with a single parameter update, whose
runtime is largely dominated by the execution of the $2N_{\theta}$ quantum
circuits required by the parameter-shift rule. In this case, our virtual QPUs
are powered by a locally modified version of the IBM Aer simulator which
leverages the efficient cuStateVec library available in the cuQuantum SDK from
NVIDIA [44]. Simulations were performed on the OLCF Summit HPC system, where
each virtual QPU maps onto a single NVIDIA Tesla V100 GPU out of each compute
node, with the results reported in Figure 8.

Figure 8: Strong scaling plots of the DDCL gradients for a single parameter
update with 20-26 qubits for the circuit in Figure 5 with 10 layers.

In Figure 8 we plot the average runtimes over 10 independent gradient
evaluations, along with error bars given by the standard deviation ($\sigma$).
It so happens that the standard deviation is negligible in the scale of the
plot, and not visible in Figure 8. In Table 3 we report those values, along
with the equation of the line from the linear fits from the corresponding data
points.

Table 3: Number of qubits ($N_{q}$), maximum standard deviation (max $\sigma$) and the equation corresponding to the linear fits of log10(runtime) vs. log2 (number of GPUs) shown in Figure 8. $N_{q}$ |  max $\sigma$ (s) | Linear fit
---|---|---
20 | 1.49 | $y=-0.132x+2.622$
22 | 0.97 | $y=-0.137x+2.726$
24 | 0.77 | $y=-0.159x+2.958$
26 | 2.96 | $y=-0.204x+3.398$

In the case of GPU execution of the DDCL algorithm, whose performance is
illustrated by Figure 8, we notice that the strong scaling deviates from the
linear trend as the number of GPUs increase. This can be attributed to the
fact that with a larger number of GPUs, each GPU execute fewer circuits, and
the time overhead associated with data transfer between CPU and GPU eventually
becomes a considerable part of the overall runtime.

The discussion on the messages exchanged in the previous example (end of
Section 4.1) largely holds in the case of the present use case, with one
crucial difference. Because the DDCL algorithm is concerned with matching a
target distribution, the relevant outcome of circuit execution in this case
are the bitstring counts, not expectation values. Since XACC stores those as a
hash table (std::map<std::string, int>), and such a data structure finds
analog in C, we start by storing the observed bitstrings and the corresponding
counts in appropriate local vectors. The communication pattern involving those
follow closely that regarding the other vectors considered in Section 4.1.

##  5 Conclusions

We put forward a parallel quantum computing software tool that is based on
virtualization of compute resources by abstracting details of the execution
and backend away from the end user. It takes advantage of the foundational
architectural feature that grants XACC a high degree of backend agnosticism to
be well positioned to stay relevant in the rapidly evolving field of quantum
computing. The independent nature of quantum circuit executions can be
explored by a parallel workflow, which we accomplish by extending the
capabilities found in the Accelerator interface following a decorator design
pattern and availing ourselves of the acclaimed MPI standard. We highlight the
performance of our scheme in the OLCF Andes and Summit leadership computers,
with both CPU-only and GPU-only runs, in the context the resource-intensive
gradient evaluation tasks in some representative quantum algorithms. The
strong scalings drawn from these simulations provide evidence that our
proposal constitute a sound and promising strategy in the pursue of
significant speedups in quantum computing workflows that revolved around
execution of a large number of circuits.

One of the main reasons our approach stands out is due to the modularity
afforded by XACC. An example of this is that by virtualizing compute
resources, the approach presented here is not constrained to classical
simulation of quantum circuits, but is equally applicable to workflows where
the accelerators are actual quantum hardware, and in principle even a
collection of radically different hardware technologies. This is line with
ensemble quantum computing,[13] with a previous demonstration presented in
Ref. 45.

Finally, in the present work we only consider a parallel workflow of execution
of quantum programs, where each virtual QPU independently executes an entire
program, which corresponds to a quantum circuit. A major development in
scalability would be achieved by an analog of MPI that would enable quantum
communication, making possible large scale quantum computation using an array
of several smaller quantum computers [46]. This scheme would allow for
distributed quantum computing by having this collection of QPUs execute a
single program, in which each individual QPU is assigned a subset of the
instructions, e.g. gates, that compose the quantum circuit in question.[47,
48] There is currently a specification for such a communication protocol, but
no available implementation, and interested readers are referred to extensive
reviews on the topic.[49, 50] As we anticipate this new paradigm in light of
what XACC is capable of, minor changes may required in order to accommodate
the communication primitives proposed by QMPI, such as ensuring that all
operations are reversible, which assumes uncomputation, as well as
modifications and/or additions in terms of abstraction layers that will enable
the virtualization platform to carry out distributed quantum computing. An
immediate example is how to determine which QPU implements a specific set of
gates. We are confident that XACC’s flexible design will carry over and be
crucial in moving toward this mode of quantum computing and is well positioned
to embrace those advancements.

##  6 Acknowledgments

DC and DL acknowledge funding by the US Department of Energy award ERKCG13
provided by the Office of Basic Energy Sciences. AJM acknowledges funding by
the US Department of Energy Office of Science Advanced Scientific Computing
Research (ASCR), Accelerated Research in Quantum Computing (ARQC). This
research used resources of the Oak Ridge Leadership Computing Facility, which
is a DOE Office of Science User Facility supported under Contract DE-
AC05-00OR22725.

Declaration of generative AI and AI-assisted technologies in the writing
process

During the preparation of this work the author(s) used ChatGPT in order to
improve readability. After using this tool/service, the author(s) reviewed and
edited the content as needed and take(s) full responsibility for the content
of the publication.

## References

  * [1] D. Aharonov, W. van Dam, J. Kempe, Z. Landau, S. Lloyd, O. Regev, Adiabatic quantum computation is equivalent to standard quantum computation, SIAM Review 50 (4) (2008) 755–787.  doi:10.1137/080734479.
URL http://dx.doi.org/10.1137/080734479

  * [2] A. McCaskey, E. Dumitrescu, D. Liakh, M. Chen, W. Feng, T. Humble, A language and hardware independent approach to quantum–classical computing, SoftwareX 7 (2018) 245 – 254.  doi:https://doi.org/10.1016/j.softx.2018.07.007.
  * [3] A. McCaskey, D. Lyakh, E. Dumitrescu, S. Powers, T. Humble, XACC: a system-level software infrastructure for heterogeneous quantum–classical computing, Quantum Science and Technology 5 (2) (2020) 024002.  doi:10.1088/2058-9565/ab6bf6.
  * [4] V. Sunderam, G. Geist, Heterogeneous parallel and distributed computing, Parallel Computing 25 (13) (1999) 1699–1721.  doi:https://doi.org/10.1016/S0167-8191(99)00088-5.
  * [5] NVIDIA, CUDA Quantum: The platform for hybrid quantum-classical computing.
URL https://github.com/NVIDIA/cuda-quantum

  * [6] R. P. Feynman, Simulating physics with computers, Int. J. Theor. Phys. 21 (6-7) (1982) 467–488.
  * [7] A. Peruzzo, J. McClean, P. Shadbolt, M.-H. Yung, X.-Q. Zhou, P. J. Love, A. Aspuru-Guzik, J. L. O’Brien, A variational eigenvalue solver on a photonic quantum processor, Nat. Commun. 5 (2014) 4213.
  * [8] S. Bravyi, G. Smith, J. A. Smolin, Trading classical and quantum computational resources, Phys. Rev. X 6 (2016) 021043.  doi:10.1103/PhysRevX.6.021043.
URL https://link.aps.org/doi/10.1103/PhysRevX.6.021043

  * [9] M. Steudtner, S. Wehner, Fermion-to-qubit mappings with varying resource requirements for quantum simulation, New J. Phys. 20 (6) (2018) 063010.
  * [10] T. Peng, A. W. Harrow, M. Ozols, X. Wu, Simulating large quantum circuits on a small quantum computer, Phys. Rev. Lett. 125 (2020) 150504.  doi:10.1103/PhysRevLett.125.150504.
URL https://link.aps.org/doi/10.1103/PhysRevLett.125.150504

  * [11] C. Tüysüz, G. Clemente, A. Crippa, T. Hartung, S. Kühn, K. Jansen, Classical splitting of parametrized quantum circuits, Quantum Mach. Intell. 5 (2) (Dec. 2023).
  * [12] A. W. Harrow, A. Lowe, Optimal quantum circuit cuts with application to clustered hamiltonian simulation (2024).  arXiv:2403.01018.
  * [13] G. L. Long, L. Xiao, Parallel quantum computing in a single ensemble quantum computer, Phys. Rev. A 69 (2004) 052303.  doi:10.1103/PhysRevA.69.052303.
URL https://link.aps.org/doi/10.1103/PhysRevA.69.052303

  * [14] K. M. Obenland, A. M. Despain, A parallel quantum computer simulator (1998).  arXiv:quant-ph/9804039.
  * [15] J. Niwa, K. Matsumoto, H. Imai, General-purpose parallel simulator for quantum computing, Phys. Rev. A 66 (2002) 062317.  doi:10.1103/PhysRevA.66.062317.
URL https://link.aps.org/doi/10.1103/PhysRevA.66.062317

  * [16] K. De Raedt, K. Michielsen, H. De Raedt, B. Trieu, G. Arnold, M. Richter, T. Lippert, H. Watanabe, N. Ito, Massively parallel quantum computer simulator, Comput. Phys. Commun. 176 (2) (2007) 121–136.
  * [17] H. De Raedt, F. Jin, D. Willsch, M. Willsch, N. Yoshioka, N. Ito, S. Yuan, K. Michielsen, Massively parallel quantum computer simulator, eleven years later, Comput. Phys. Commun. 237 (2019) 47–61.
  * [18] E. Gutierrez, S. Romero, M. A. Trenas, E. L. Zapata, Parallel quantum computer simulation on the CUDA architecture, in: Computational Science – ICCS 2008, Lecture notes in computer science, Springer Berlin Heidelberg, Berlin, Heidelberg, 2008, pp. 700–709.
  * [19] A. Amariutei, S. Caraiman, Parallel quantum computer simulation on the gpu, in: 15th International Conference on System Theory, Control and Computing, 2011, pp. 1–6.
  * [20] R. Li, B. Wu, M. Ying, X. Sun, G. Yang, Quantum supremacy circuit simulation on sunway TaihuLight, IEEE Trans. Parallel Distrib. Syst. 31 (4) (2020) 805–816.
  * [21] Z. Wang, Z. Chen, S. Wang, W. Li, Y. Gu, G. Guo, Z. Wei, A quantum circuit simulator and its applications on sunway TaihuLight supercomputer, Sci. Rep. 11 (1) (2021) 355.
  * [22] H. Shang, L. Shen, Y. Fan, Z. Xu, C. Guo, J. Liu, W. Zhou, H. Ma, R. Lin, Y. Yang, F. Li, Z. Wang, Y. Zhang, Z. Li, Large-scale simulation of quantum computational chemistry on a new sunway supercomputer, in: Proceedings of the International Conference on High Performance Computing, Networking, Storage and Analysis, SC ’22, IEEE Press, 2022.
  * [23] M. Morita, Y. Tomita, J. Koyama, K. Kimura, Simulator demonstration of large scale variational quantum algorithm on hpc cluster (2024).  arXiv:2402.11878.
  * [24] A. Li, O. Subasi, X. Yang, S. Krishnamoorthy, Density matrix quantum circuit simulation via the bsp machine on modern gpu clusters, in: Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis, 2020.
  * [25] A. Li, B. Fang, C. Granade, G. Prawiroatmodjo, B. Hein, M. Rotteler, S. Krishnamoorthy, Sv-sim: Scalable pgas-based state vector simulation of quantum circuits, in: Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis, 2021.
  * [26] M. Wang, F. Hua, C. Liu, N. Bauman, K. Kowalski, D. Claudino, T. Humble, P. Nair, A. Li, Enabling scalable vqe simulation on leading hpc systems, in: Proceedings of the SC ’23 Workshops of The International Conference on High Performance Computing, Network, Storage, and Analysis, SC-W ’23, Association for Computing Machinery, New York, NY, USA, 2023, p. 1460–1467.  doi:10.1145/3624062.3624221.
URL https://doi.org/10.1145/3624062.3624221

  * [27] I.-S. Suh, A. Li, Simulating quantum systems with nwq-sim on hpc (2024).  arXiv:2401.06861.
  * [28] D. I. Lyakh, T. Nguyen, D. Claudino, E. Dumitrescu, A. J. McCaskey, ExaTN: Scalable GPU-accelerated high-performance processing of general tensor networks at exascale, Front. Appl. Math. Stat. 8 (Jul. 2022).
  * [29] D. Lykov, A. Chen, H. Chen, K. Keipert, Z. Zhang, T. Gibbs, Y. Alexeev, Performance evaluation and acceleration of the qtensor quantum circuit simulator on gpus, in: 2021 IEEE/ACM Second International Workshop on Quantum Computing Software (QCS), IEEE, 2021.  doi:10.1109/qcs54837.2021.00007.
URL http://dx.doi.org/10.1109/QCS54837.2021.00007

  * [30] M. Fishman, S. R. White, E. M. Stoudenmire, The ITensor Software Library for Tensor Network Calculations, SciPost Phys. Codebases (2022) 4doi:10.21468/SciPostPhysCodeb.4.
URL https://scipost.org/10.21468/SciPostPhysCodeb.4

  * [31] T. Jones, A. Brown, I. Bush, S. C. Benjamin, Quest and high performance simulation of quantum computers, Scientific reports 9 (1) (2019) 10736.
  * [32] Q. A. team, collaborators, qsim (Sep. 2020).  doi:10.5281/zenodo.4023103.
URL https://doi.org/10.5281/zenodo.4023103

  * [33] V. Gheorghiu, Quantum++: A modern c++ quantum computing library, PLOS ONE 13 (12) (2018) e0208073.  doi:10.1371/journal.pone.0208073.
URL http://dx.doi.org/10.1371/journal.pone.0208073

  * [34] A. McCaskey, E. Dumitrescu, M. Chen, D. Lyakh, T. Humble, Validating quantum-classical programming models with tensor network simulations, PLOS ONE 13 (12) (2018) 1–19.  doi:10.1371/journal.pone.0206704.
URL https://doi.org/10.1371/journal.pone.0206704

  * [35] T. Nguyen, D. Lyakh, E. Dumitrescu, D. Clark, J. Larkin, A. McCaskey, Tensor network quantum virtual machine for simulating quantum circuits at exascale, ACM Transactions on Quantum Computing 4 (1) (oct 2022).  doi:10.1145/3547334.
URL https://doi.org/10.1145/3547334

  * [36] R. M. Parrish, E. G. Hohenstein, P. L. McMahon, T. J. Martínez, Quantum computation of electronic transitions using a variational quantum eigensolver, Phys. Rev. Lett. 122 (2019) 230401.  doi:10.1103/PhysRevLett.122.230401.
URL https://link.aps.org/doi/10.1103/PhysRevLett.122.230401

  * [37] J. B. Foresman, M. Head-Gordon, J. A. Pople, M. J. Frisch, Toward a systematic molecular orbital theory for excited states, The Journal of Physical Chemistry 96 (1) (1992) 135–149.  doi:10.1021/j100180a030.
URL http://dx.doi.org/10.1021/j100180a030

  * [38] W. Dür, G. Vidal, J. I. Cirac, Three qubits can be entangled in two inequivalent ways, Phys. Rev. A 62 (2000) 062314.  doi:10.1103/PhysRevA.62.062314.
URL https://link.aps.org/doi/10.1103/PhysRevA.62.062314

  * [39] M. Schuld, V. Bergholm, C. Gogolin, J. Izaac, N. Killoran, Evaluating analytic gradients on quantum hardware, Phys. Rev. A 99 (2019) 032331.  doi:10.1103/PhysRevA.99.032331.
URL https://link.aps.org/doi/10.1103/PhysRevA.99.032331

  * [40] W. Jakob, J. Rhinelander, D. Moldovan, pybind11 – Seamless operability between C++11 and Python, https://github.com/pybind/pybind11 (2017).
  * [41] 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 (1) (May 2019).  doi:10.1038/s41534-019-0157-8.
URL http://dx.doi.org/10.1038/s41534-019-0157-8

  * [42] V. Leyton-Ortega, A. Perdomo-Ortiz, O. Perdomo, Robust implementation of generative modeling with parametrized quantum circuits (2019).  doi:10.48550/ARXIV.1901.08047.
URL https://arxiv.org/abs/1901.08047

  * [43] A. Delgado, K. E. Hamilton, Unsupervised quantum circuit learning in high energy physics, Phys. Rev. D 106 (2022) 096006.  doi:10.1103/PhysRevD.106.096006.
URL https://link.aps.org/doi/10.1103/PhysRevD.106.096006

  * [44] H. Bayraktar, A. Charara, D. Clark, S. Cohen, T. Costa, Y.-L. L. Fang, Y. Gao, J. Guan, J. Gunnels, A. Haidar, A. Hehn, M. Hohnerbach, M. Jones, T. Lubowe, D. Lyakh, S. Morino, P. Springer, S. Stanwyck, I. Terentyev, S. Varadhan, J. Wong, T. Yamaguchi, cuQuantum SDK: A high-performance library for accelerating quantum science, in: 2023 IEEE International Conference on Quantum Computing and Engineering (QCE), Vol. 01, 2023, pp. 1050–1061.  doi:10.1109/QCE57702.2023.00119.
  * [45] S. Stein, N. Wiebe, Y. Ding, P. Bo, K. Kowalski, N. Baker, J. Ang, A. Li, EQC: Ensembled Quantum Computing for Variational Quantum Algorithms, in: Proceedings of the 49th Annual International Symposium on Computer Architecture, ISCA ’22, Association for Computing Machinery, New York, NY, USA, 2022, p. 59–71.  doi:10.1145/3470496.3527434.
URL https://doi.org/10.1145/3470496.3527434

  * [46] T. Häner, D. S. Steiger, T. Hoefler, M. Troyer, Distributed quantum computing with qmpi, in: Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis, SC ’21, Association for Computing Machinery, New York, NY, USA, 2021.  doi:10.1145/3458817.3476172.
URL https://doi.org/10.1145/3458817.3476172

  * [47] R. Van Meter, K. Nemoto, W. J. Munro, K. M. Itoh, Distributed arithmetic on a quantum multicomputer, SIGARCH Comput. Archit. News 34 (2) (2006) 354–365.  doi:10.1145/1150019.1136517.
URL https://doi.org/10.1145/1150019.1136517

  * [48] R. Beals, S. Brierley, O. Gray, A. W. Harrow, S. Kutin, N. Linden, D. Shepherd, M. Stather, Efficient distributed quantum computing, Proc. Math. Phys. Eng. Sci. 469 (2153) (2013) 20120686.
  * [49] M. Caleffi, M. Amoretti, D. Ferrari, D. Cuomo, J. Illiano, A. Manzalini, A. S. Cacciapuoti, Distributed quantum computing: a survey (2022).  arXiv:2212.10609.
  * [50] D. Barral, F. J. Cardama, G. Díaz, D. Faílde, I. F. Llovo, M. M. Juane, J. Vázquez-Pérez, J. Villasuso, C. Piñeiro, N. Costas, J. C. Pichel, T. F. Pena, A. Gómez, Review of distributed quantum computing. from single qpu to high performance quantum computing (2024).  arXiv:2404.01265.


capbtabboxtable[][]

# MAQA: A Quantum Framework for Supervised Learning

Antonio Macaluso
German Research Center for Artificial Intelligence
(DFKI)
Saarbruecken, Germany
<EMAIL_ADDRESS>
&Matthias Klusch
German Research Center for Artificial Intelligence
(DFKI)
Saarbruecken, Germany
<EMAIL_ADDRESS>
Stefano Lodi
University of Bologna
Bologna, Italy
<EMAIL_ADDRESS>
&Claudio Sartori
University of Bologna
Bologna, Italy
<EMAIL_ADDRESS>
corresponding author

###### Abstract

Quantum Machine Learning has the potential to improve traditional machine
learning methods and overcome some of the main limitations imposed by the
classical computing paradigm. However, the practical advantages of using
quantum resources to solve pattern recognition tasks are still to be
demonstrated.

This work proposes a universal, efficient framework that can reproduce the
output of a plethora of classical supervised machine learning algorithms
exploiting quantum computation’s advantages. The proposed framework is named
Multiple Aggregator Quantum Algorithm (MAQA) due to its capability to combine
multiple and diverse functions to solve typical supervised learning problems.
In its general formulation, MAQA can be potentially adopted as the quantum
counterpart of all those models falling into the scheme of aggregation of
multiple functions, such as ensemble algorithms and neural networks. From a
computational point of view, the proposed framework allows generating an
exponentially large number of different transformations of the input at the
cost of increasing the depth of the corresponding quantum circuit linearly.
Thus, MAQA produces a model with substantial descriptive power to broaden the
horizon of possible applications of quantum machine learning with a
computational advantage over classical methods. As a second meaningful
addition, we discuss the adoption of the proposed framework as hybrid quantum-
classical and fault-tolerant quantum algorithm.

_K_ eywords Quantum Machine Learning $\cdot$ Quantum Computing $\cdot$ Machine
Learning

##  1 Introduction

Quantum computers are machines that leverage the properties of quantum
mechanics to store and process information. Although a potential quantum
advantage has already been shown in different fields, such as quantum
chemistry [1], multi-agent systems [2, 3], it is still unclear whether quantum
computing can be efficiently used for machine learning (ML) tasks.

The intersection between ML and quantum computing (QC) is known as Quantum
Machine Learning (QML). There are two ways in which ML and QC can be combined:
one approach is to run the learning process predominantly in a quantum
computer so that the expensive subroutines can be executed efficiently. For
this purpose, a rich collection of quantum algorithms for basic linear algebra
subroutines have been proposed in literature [4, 5, 6]. Some popular examples
of this approach are QSVM [7] and QSplines [8], which obtain an exponential
speed-up with respect to their classical counterparts. However, the protocols
within this category usually assume the availability of a fault-tolerant
quantum computer.

Alternatively, variational quantum algorithms can be considered machine
learning models that can be trained using hybrid quantum-classical
optimization. In this case, a quantum algorithm is used to make a call to a
function that allows estimating the target variable of interest given the
input data and a set of rotation parameters [9, 10]. This approach requires a
parametrised quantum circuit and a classical optimisation procedure to find
the optimal set of parameters for a sequence of quantum gates. Although these
techniques represent the most promising attempt to leverage near-term quantum
technology, it is still unclear whether they can outperform classical
algorithms.

Despite the remarkable success of ML in numerous real-world applications, the
ever-increasing size of datasets and the high computational requirements of
modern algorithms indicate that the current computational tools will no longer
be sufficient in the future. In this work, we propose a novel and efficient
quantum framework to reproduce a plethora of machine learning models using
quantum computational advantages. The framework is called Multiple Aggregator
Quantum Algorithm (MAQA) due to its capability to combine multiple and diverse
functions to solve typical supervised learning tasks. Thanks to superposition,
entanglement and interference, the MAQA framework can compute the weighted
average of an exponentially large number of functions while increasing the
depth of the correspondent quantum circuit linearly. This allows for building
quantum models with incredible descriptive power that might be a credible
alternative to classical methods in the future.

##  2 Preliminaries

The objective of a supervised model is to find a useful approximation to the
function $f(x;\theta)$ that underlies the predictive relationship between the
input $x$ and output $y$ for a fixed set of parameters $\theta$. Assuming for
simplicity an additive error, the model of interest can be expressed as
follows:

$\displaystyle y=f(x;\theta)+\epsilon,$ (1)

where $\epsilon$ is a random variable whose conditioned probability
distribution given $x$ is centred in zero. Although Eq. (1) provides a general
mathematical formulation for supervised learning, several methods do not
estimate a single function but explicitly calculate multiple and diverse
functions. These functions belong to the same family but differ in either a
set of parameters or the training data. In all these cases, the final output
results from the weighted average of the estimated functions:

$\displaystyle y=f(x;\theta)$
$\displaystyle=\sum_{h=1}^{H}\beta_{h}g(x;\theta_{h}),$ (2)

where $f(x;\theta)$ is the final output and $g(x;\cdot)$ describes the
function component.

The calculation of $g(x;\cdot)$ corresponds to a specific transformation of
data $x$ based on $\theta_{h}$, whose contribution to the final output is
weighted by $\beta_{h}$. The estimation of a collection of functions
components allows producing an extremely flexible model, which is able to
approximate the behaviour of complex patterns. Different choices for $\beta$,
$g(x;\cdot)$ and $\theta_{h}$ determine different supervised models commonly
adopted in real-world applications.

For instance, a single-layer neural network (or Single Layer Perceptron - SLP)
assumes as function component $g(x;\cdot)$ the activation function
$\sigma_{\text{hidden}}$ applied to the linear combinations $L(x;\theta_{h})$
of the input vector $x$. In fact, an SLP with $H$ hidden neurons is a two-
stage model that takes as input training data $x$ and $H$ sets of linear
coefficients and estimates the target variable as follows:

$\displaystyle
f_{\text{SLP}}(x)=\sigma_{\text{output}}\left[\sum_{h=1}^{H}\beta_{h}\sigma_{\text{hidden}}\left(L(x;\Theta_{h})\right)\right],$
(3)

where $\sigma_{\text{output}}$ is the identity function when the task is the
function approximation. 111When considering a neural network with multiple
hidden layers, the only difference in Eq. (3) is that the function component
$g(x;\cdot)$ is, in turn, a neural network.

Another classical supervised learning approach that falls into the schema of
function aggregation is ensemble learning. In practice, ensemble methods
reduce to computing several predictions $g_{1}(x),g_{2}(x),\dots,g_{H}(x)$
using $H$ different training sets, which are then averaged to obtain a single
model:

$\displaystyle f_{\text{ens}}(x)=\frac{1}{H}\sum_{h=1}^{H}\beta_{h}g_{h}(x).$
(4)

In this case, the component functions $g(x;\cdot)$ are weak
classification/regression models and the choice of the weights depends on the
type of the ensemble in use (boosting, bagging, randomisation).

Other models that fit into the idea of multiple aggregations are Generalised
Additive Models [11], Support Vector Machines and Decision trees [12].

#### Contribution.

In this work, we propose a novel efficient quantum framework to reproduce the
idea of machine learning models as functions aggregators. The proposed
architecture, named Multiple Aggregator Quantum Algorithm (MAQA), can
potentially reproduce some of the most important classical supervised learning
algorithms introducing relevant computational advantages. In particular, MAQA
propagates an input state to multiple quantum trajectories in superposition,
and each trajectory describes a specific function $g(x;\cdot)$ that represents
the component function of the final model. The entanglement between the two
quantum registers involved (data and control) allows for efficient averaging
of those transformations, and the final result can be accessed by measuring
only a subset of qubits. The proposed approach has two main advantages: from a
classical perspective, it introduces an exponential scaling in the number of
aggregated functions while linearly increasing the time complexity of the
correspondent quantum algorithm. From a quantum perspective, the framework
opens the possibility of implementing a plethora of models not yet proposed in
the literature. Eventually, we discuss the adoption of MAQA to generalise some
existing QML algorithms, considering both fault-tolerant settings and hybrid
quantum-classical algorithms.

##  3 Multiple Aggregator Quantum Algorithm (MAQA)

In this section, we describe the MAQA framework that is able to reproduce the
classical model expressed in Eq. (2). The algorithm leverages the three main
properties of quantum computing (superposition, entanglement and interference)
to encode in a quantum state the sum of different input transformations
accessible by measuring a single quantum register. The proposed algorithm can
potentially reproduce all those models that refer to the idea of functions
aggregation and provide attractive computational advantages with respect to
the classical counterparts.

The quantum algorithm adopts two quantum registers: data and control. The
$data$ register encodes the model’s input data, and the $control$ register is
used to generate multiple trajectories in superposition, where each trajectory
represents a different transformation of data.

Starting from a $n$-qubit $data$ register and a $d$-qubit $control$ register
the Multiple Aggregator Quantum Algorithm (MAQA) involves four main steps:
state preparation, multiple trajectories in superposition, transformation via
interference and measurement.

#### (Step 1) State Preparation

State preparation consists of encoding the input in the $data$ register and
the initialisation of the $control$ register whose amplitudes depend on a set
of parameters $\beta=\\{\beta^{*}_{i}\\}_{i=1,\dots,2^{d}}$:

$\displaystyle\ket{\Phi_{0}}$ $\displaystyle=(S_{\beta}\otimes
S_{x})\ket{0}_{\text{control}}^{\otimes d}\otimes\ket{0}^{\otimes
n}_{\text{data}}$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{2^{d}}}\sum_{h=1}^{2^{d}}\beta^{*}_{h}\ket{h}\otimes\ket{x}.$
(5)

We refer to $S_{x}$ as a quantum routine to encode data into a quantum state,
and $S_{\beta}$ as a routine that transforms a $d$-qubit register from an all-
zero state to a quantum state which depends on a set of parameters $\beta$.
Importantly, the computational cost of this step is not considered classically
since any classical algorithm assumes the input $x$ to be directly accessible.

#### (Step 2) Multiple Trajectories in Superposition

The second step regards the generation of $2^{d}$ different transformations of
the input data in superposition, each entangled with a possible state of the
$control$ register. The single quantum state of the superposition encodes a
specific transformation of the data and it depends on a set of parameters
$\Theta_{k}$. To this end, a unitary $G(\theta_{1},\dots,\theta_{2^{d}})$ that
performs the following operation is assumed222Notice that the definition of
$G(\theta_{1},\dots,\theta_{2^{d}})$ unitary in terms of quantum gates depend
on the specific algorithm in use (Sections 4.2 and 4.1).:

$\displaystyle\ket{\Phi_{1}}$
$\displaystyle=G\left(\theta_{1},\dots,\theta_{2^{d}}\right)\ket{\Phi_{0}}=\frac{1}{\sqrt{2^{d}}}\left(\sum_{h=1}^{2^{d}}\beta^{*}_{h}\ket{h}\ket*{l(x;\Theta_{h})}\right).$
(6)

where the implementation of $G(\theta_{1},\dots,\theta_{2^{d}})$ can be
accomplished in only $d$ steps. Each step consists in the entanglement of the
$i^{th}$ $(i=1,\dots d)$ control qubit with two transformations
$G\left(\theta_{i,1}\right)$ and $G\left(\theta_{i,2}\right)$ of $\ket{x}$
based on two sets of parameters, $\theta_{i,1}$ and $\theta_{i,2}$. Let us
consider a unitary $G\left(\theta_{i,j}\right)$ that implements the
transformation $l\left(x;\theta_{i,j}\right)$. The most straightforward way to
obtain the quantum state in Eq. (6) is to apply $G\left(\theta_{i,j}\right)$
through controlled operations, using as control state the two basis states of
the current control qubit. In particular, the generic $i^{th}$ step involves
the following two transformations:

First, the controlled-unitary $C^{(1)}G\left(\theta_{i,1}\right)$ is executed
to entangle the transformation $G\left(\theta_{i,1}\right)\ket{x}$ with the
excited state $\ket{1}$ of the $i^{th}$ control qubit:

$\displaystyle\ket{\Phi_{i,1}}=$
$\displaystyle\left(C^{(1)}G\left(\theta_{i,1}\right)\right)\ket{c_{i}}\otimes\ket{x}=\left(C^{(1)}G\left(\theta_{i,1}\right)\right)\left(a_{i}\ket{0}+b_{i}\ket{1}\right)\otimes\ket{x}$
$\displaystyle=$
$\displaystyle\left(a_{i}\ket{0}\ket{x}+b_{i}\ket{1}G\left(\theta_{i,1}\right)\ket{x}\right),$
(7)

where $a_{i}$ and $b_{i}$ are the amplitudes of the $i^{th}$ control qubit and
$C^{(1)}G(\theta_{i,1})$ is a controlled operation that entangles the exited
state of the control qubit $\ket{c_{i}}$ to transform the $data$ register
according to the unitary $G(\theta_{i,1})$.

Then, a second controlled-unitary $C^{(0)}G\left(\theta_{i,2}\right)$ is
executed. This time the control state is the $\ket{0}$ basis state:

$\displaystyle\ket{\Phi_{i}}=$
$\displaystyle\left(C^{(0)}G\left(\theta_{i,2}\right)\right)\ket{\Phi_{i,1}}=a_{i}\ket{0}G\left(\theta_{i,2}\right)\ket{x}+b_{i}\ket{1}G\left(\theta_{i,1}\right)\ket{x}.$
(8)

These two transformations are repeated for each qubit in the control register
and two different unitaries $G\left(\theta_{i,1}\right)$ and
$G\left(\theta_{i,2}\right)$ are applied, at each iteration. After $d$ steps,
the control and data registers are fully entangled and $2^{d}$ different
quantum trajectories in superposition are generated. The output of this
procedure can be expressed as follows:

$\displaystyle\ket{\Phi_{d}}=$
$\displaystyle\frac{1}{\sqrt{2^{d}}}\sum_{h=1}^{2^{d}}\beta^{*}_{h}\ket{h}G\left(\Theta_{h}\right)\ket{x}=\frac{1}{\sqrt{2^{d}}}\sum_{h=1}^{2^{d}}\beta^{*}_{h}\ket{h}\ket{l(x;\Theta_{h})}$
(9)

where $G\left(\Theta_{h}\right)$ results from the product of $d$ unitary
matrices $G\left(\theta_{i,j}\right)$ and it represents a single quantum
trajectory. Each trajectory differs from the others for, at least, one unitary
$G\left(\theta_{i,j}\right)$333A detailed example with $d=3$ is described in
Appendix Appendix.

When discussing a specific implementation of QML algorithms (Sections 4.2 and
4.1), we will see that, from a computational point of view, the possibility to
generate $2^{d}$ different transformations in only $d$ steps potentially leads
to scaling exponentially the number of component functions with respect to
classical methods, assuming an efficient implementation of the
$C^{(j)}G(\theta_{i,j})$.

#### (Step 3) Transformation via Interference

Once we generated multiple transformations $l(x;\Theta_{h})$ of the input in
superposition, the third step consists of transforming the $data$ register
through a generic quantum gate $F$ that works via interference:

$\displaystyle\ket{\Phi_{f}}$ $\displaystyle=\left(\mathbb{1}^{\otimes
d}\otimes F\right)\ket{\Phi_{d}}=\left(\mathbb{1}^{\otimes d}\otimes
F\right)\left[\frac{1}{\sqrt{2^{d}}}\sum_{h=1}^{2^{d}}\beta^{*}_{h}\ket{h}\ket{l(x;\Theta_{h})}\right]$
$\displaystyle=\frac{1}{\sqrt{2^{d}}}\sum_{h=1}^{2^{d}}\beta^{*}_{h}\ket{h}\ket*{g^{*}\left(x;\Theta_{h}\right)}=\frac{1}{\sqrt{H}}\sum_{h=1}^{H}\beta^{*}_{h}\ket{h}\ket*{g^{*}_{h}},$
(10)

where $H=2^{d}$. In Eq. (3) the assumption is that the sequential application
of $G(\Theta_{h})$ and $F$ on the quantum state $\ket{x}$ is equivalent to
calculate the function $g_{h}^{*}$ to an input $x$. At this point, different
values of the function $g_{h}^{*}$ are entangled with different states of the
$control$ register.

It is important to notice that a single execution of $F$ allows the
computation of the function $g_{h}^{*}$ for all the quantum trajectories in
superposition. This is extremely useful when, during the computation, the same
operations need to be applied to multiple inputs (e.g., when the activation
function is applied to a huge number of neurons or in the case of ensemble
learning, where the same classifier has to be executed to different sub-
samples of the training set).

#### (Step 4) Measurement

The last step consists of measuring the $data$ register, leaving untouched the
$control$ register:

$\displaystyle\left\langle M\right\rangle$
$\displaystyle=\innerproduct{\Phi_{f}|\mathbb{1}^{\otimes d}\otimes
M}{\Phi_{f}}=\frac{1}{H}\sum_{h=1}^{H}\beta_{h}\innerproduct{h}{h}\otimes\innerproduct*{g^{*}_{h}|M}{g^{*}_{h}}$
$\displaystyle=\frac{1}{H}\sum_{h=1}^{H}\beta_{h}\innerproduct*{g^{*}_{h}|M}{g^{*}_{h}}=\frac{1}{H}\sum_{h=1}^{H}\beta_{h}\left\langle
M_{h}\right\rangle$
$\displaystyle=\frac{1}{H}\sum_{h=1}^{H}\beta_{h}g\left(x;\Theta_{h}\right)=f_{\text{agg}},$
(11)

where $g\left(x;\Theta_{h}\right)=\innerproduct*{g^{*}_{h}|M}{g^{*}_{h}}$ for
a measurement operator $M$, $\beta_{h}=|\beta^{*}_{h}|^{2}$ with
$\sum_{h}|\beta_{h}|^{2}=1$ and $H=2^{d}$.

The expectation value $\left\langle M\right\rangle$ stores the weighted
average of the $2^{d}$ functions $g\left(x;\Theta_{h}\right)$, which is
accessible by measuring the $data$ register. While extracting the single
contribution $g\left(x;\Theta_{h}\right)$ would require an exponential number
of measurements (since those values are in the superposition of $2^{d}$
possible basis states), in a classical supervised learning scenario the
measure of interest is the weighted average of all the functions which can be
directly accessed by measuring the data register and leaving intact the
control register.

To summarise, the proposed architecture allows calculating the aggregation of
multiple and diverse functions described in Eq. (2) using a quantum algorithm.
In particular, it is possible to access the final result by measuring only the
$data$ register while obtaining the weighted average of $2^{d}$ different
transformations $g(x;\cdot)$ of the input data $x$, where $d$ is the size of
the control register. Specifying properly $S_{\beta}$, $S_{x}$,
$\\{G\left(\theta_{i,1}\right),G\left(\theta_{i,2}\right)\\}_{i=1,\dots,d}$
and $F$ allows potentially to reproduce the quantum version of all the ML
algorithms discussed in Section 2. Furthermore, the framework is very generic
and can be adopted for hybrid and fault-tolerant quantum computation. The
quantum circuit for implementing MAQA is depicted in Figure 1.

@C=.7em @R=.6em |0⟩_1 & 3S_β @–[]+<1.4em,2em>;[d]+<1.4em,-7em>
@–[]+<2.2em,2em>;[d]+<2.2em,-7em> @–[]+<2.2em,2em>;[d]+<2.2em,-7em>
…@–[]+<1.2em,2em>;[d]+<1.2em,-7em> 4 4 @–[]+<2.5em,2em>;[d]+<2.5em,-7em>
@–[]+<2.4em,2em>;[d]+<2.4em,-7em>
|0⟩_2 S_β …
… S_β 2 2 …
|0⟩_d S_β 1 1 …
|0⟩^⊗n S_x G(θ_1,1) G(θ_1,2) G(θ_2,1) G(θ_2,2) … G(θ_d,1) G(θ_d,2) F

|Φ_0⟩ |Φ_1⟩ |Φ_2⟩ |Φ_d-1⟩ |Φ_d⟩ |Φ_f⟩

Figure 1: Quantum Circuit for the Multiple Aggregator Quantum Algorithm
(MAQA).

##  4 Discussion

As shown in the previous section, the MAQA allows obtaining a quantum state
that reproduces the idea of ML models as aggregators of functions using the
properties of quantum computing. From a classical ML perspective, relevant
computational advantages are introduced. Given $2^{d}$ component functions,
any classical method that leverages the idea of functions aggregation scales
linearly in $2^{d}$ since it is necessary to compute those functions
explicitly to obtain the overall average. Furthermore, in the worst-case
scenario, each component function has to process all available data; this
implies a linear cost in the training set size multiplied by $2^{d}$. Using
big-$\mathcal{O}$ notation, given a dataset $(x_{i},y_{i})$ for $i=1,\dots N$,
where $x_{i}$ is a $p$-dimensional vector, and $y_{i}$ is the target variable
of interest, the overall time complexity of a model based on the aggregation
of $2^{d}$ functions is:

$\displaystyle\mathcal{O}(2^{d}N^{\alpha}p^{\beta})\qquad\alpha,\beta\geq 1.$

In contrast, MAQA generates a superposition of $2^{d}$ different
transformations of the input in only $d$ steps since the single
transformations are not computed directly, but they result from the
combination of different unitaries $G(\theta_{i,j})$. Then, once the quantum
state in Eq. (3) is generated, any operation (unitary $F$) is propagated to
all the quantum trajectories with a single execution. Using big-$\mathcal{O}$
notation, the time complexity of implementing the MAQA is:

$\displaystyle\mathcal{O}\left(d\times 2C_{G}+C_{F}\right),$

where $C_{G}$ is the cost the controlled operation $C^{(j)}G(\theta)$ and
$C_{F}$ is the cost of $F$. Note that the number of different functions grows
exponentially with respect to the parameter $d$, which has a linear impact on
the overall time complexity. This means that it is possible to generate an
exponentially large number of different transformations of the input while
obtaining their average efficiently, at the cost of increasing the depth of
the corresponding quantum circuit linearly by a factor of $2C_{G}$.

However, these advantages come with some compromises. First, the assumption
about the nature of the operator $G(\theta_{i,j})$. In fact, MAQA assumes that
the product of $G(\theta_{i,j})$ for $i=1,\dots,d$ produce a quantum gate
$G\left(\Theta_{k}\right)$:

$\displaystyle
G\left(\Theta_{k}\right)=\prod_{\begin{subarray}{c}i=1,\dots,d\\\
j=1,2\end{subarray}}G\left(\theta_{i,j}\right).$ (12)

In practice, this means that multiple applications of the unitaries that
depend on some set of parameters $\theta_{i,j}$ result in a single
transformation of the same nature that depends on a derived set of parameters
$\Theta_{k}$. Although any quantum circuit can be expressed as the product of
different unitary matrices, the design of these gates in the context of
supervised learning needs to be accomplished such that the final measurement
provides the target variable of interest.

Finally, when comparing classical and quantum algorithms, it is important to
consider that quantum computation introduces a new complexity class, the
Bounded-error Quantum Polynomial time, representing the class of problems
solvable in polynomial time by an innately probabilistic quantum Turing
machine [13]. Nevertheless, quantum algorithms need to be evaluated in terms
of gate complexity. Thus, it is necessary that the exponential scaling
introduced with respect to $d$ is preserved when considering a specific QML
model.

###  4.1 MAQA as Hybrid Quantum-Classical Algorithm

Recently the idea of aggregating two different unitary operators to reproduce
the output of a two-neuron single-layer neural network via quantum circuit has
been proposed (qSLP) [14, 15]. Since multiple aggregations are the basis of
both qSLP and MAQA, the latter can be seen as a natural extension of the
former with an exponentially large number of neurons in the hidden layer. In
fact, the entanglement between the control and data registers implies the
number of linear combinations to be equal to the number of basis states of the
control register. This, in turn, implies that the number of hidden neurons $H$
scales exponentially with the number of states of the control register as a
consequence of each hidden neuron being represented by a quantum trajectory.
This exponential scaling might enable the construction of a qSLP with an
arbitrarily large number of hidden neurons as the amount of available qubits
increases. In other words, by adopting MAQA to generalise the qSLP, we can
build a model with an incredible descriptive power capable of being a
universal approximator.

From a computational point of view, given $H$ hidden neurons and $L$ training
epochs, the training of a classical SLP scales (at least) linearly in $H$ and
$L$ since the output of each hidden neuron has to be calculated explicitly to
obtain the final output. Furthermore, if $H$ is too large (a necessary
condition for an SLP to be a universal approximator [14, 16]), the problem
becomes NP-hard [17]. The adoption of MAQA to generalize the qSLP allows
scaling linearly with respect to $log_{2}(H)=d$, thanks to the entanglement
between the two quantum registers, which allows generating an exponentially
large number of quantum trajectories in superposition.

However, the main challenge to tackle in the near future for qSLP-MAQA is
still the design of a proper activation function – in the sense of the
Universal Approximation Theorem – which is one of the significant issues for
building a complete quantum neural network. Yet, a recent proposal of QSplines
[8] opened the possibility of approximating non-linear activation functions
via a quantum algorithm. Even so, QSplines use the HHL as a subroutine, a
fault-tolerant quantum algorithm that cannot be adopted in hybrid computation
on NISQ devices.

Nevertheless, recently it has been shown that quantum feature maps alongside
functions aggregation is able to achieve universal approximation [18]. Thus, a
possible future work consists of studying the qSLP-MAQA on top of the quantum
feature map to enable it as a universal functions approximator without
implementing a non-linear quantum activation function.

###  4.2 MAQA as Fault-Tolerant Quantum Algorithm

Recently, a quantum algorithm that implements the idea of ensemble methods has
been proposed [19] and further developed [20]. Looking at the specific quantum
circuit in use, it is possible to observe that quantum ensembles can be
considered as a particular instance of MAQA, where the controlled-rotation in
Eq. (8), (3) are implemented using only the basis state $\ket{1}$ as control
state which is transformed through Pauli-X gate at each iteration.
Furthermore, while MAQA allows flexible quantum trajectories in terms of
parametrised quantum gates $S_{\beta}$ and
$\\{G\left(\theta_{i,1}\right),G\left(\theta_{i,2}\right)\\}_{i=1,\dots,d}$,
in the case of quantum ensemble [20] the weights are pre-fixed (uniform
superposition of the control register) and the transformations of the input
data are represented by CSWAP operations. Thus, MAQA potentially extends the
proposed quantum ensemble specifically defined for bagging strategy to other
ensembles such as boosting and randomisation, where the parameters of the
single base model and the correspondent weights are not pre-fixed. Still, the
main drawback of the quantum ensemble remains the underlying assumption to
encode the complete training and test set into two different quantum registers
and use a large number of trajectories in superposition to compute different
subsamples of the training set. This would require an incredibly high number
of qubits in a fault-tolerant quantum computer.

In this respect, the main challenge to tackle to make the ensemble effective
(using MAQA) in the near future is the design of a quantum classifier based on
interference that guarantees a more efficient data encoding strategy (e.g.
amplitude encoding) and can process larger datasets.

##  5 Conclusions

The practical advantages of using quantum resources to solve machine learning
tasks are still to be demonstrated. However, the ground provided by quantum
mechanics is highly appealing since a low number of qubits allows accessing an
exponentially large Hilbert space.

In this work, we tried to take a further step towards the study of how machine
learning can benefit from quantum computation. The proposed quantum framework,
the Multiple Aggregator Quantum Algorithm (MAQA), is capable of reproducing
some of the most important classical machine learning algorithms using quantum
computing resources. MAQA can potentially improve, in terms of time
complexity, all those models that require explicitly computing multiple and
diverse functions to produce a final strong model. In particular, the cost
aggregating $H$ different functions in classical machine learning requires a
computational cost linear in $H$. Instead, the proposed quantum architecture
allows scaling exponentially in $H$, requiring only $log_{2}(H)$ steps under
the assumption that the cost in terms of circuit complexity is unitary for
each step. The advantage comes directly from using superposition and
entanglement as resources for generating different transformations of the
input. Furthermore, quantum interference allows propagating the use of a
specific unitary (gate $F$) to all the quantum trajectories in superposition.
Hence, the application of $F$ impacts additively the overall time complexity,
and the same operation would require a multiplicative cost in classical
computation.

In addition, we discussed how the proposed approach could be adopted as a
fault-tolerant (quantum ensemble) and hybrid quantum-classical (quantum Single
Layer Perceptron) algorithm, though different technical aspects need to be
further investigated for both cases.

We are still in an early stage of QML, and its contribution to solving real-
world in the context of machine learning is yet to be understood. However,
many research findings, including this work, suggest that the potential of
quantum computing is huge, and machine learning will likely benefit from it in
the future.

## Acknowledgment

This work has been partially funded by the German Ministry for Education and
Research (BMB+F) in the project QAI2-QAICO under grant 13N15586.

This version of the article has been accepted for publication, after peer
review but is not the Version of Record and does not reflect post-acceptance
improvements, or any corrections. The Version of Record is available online
at: http://dx.doi.org/10.1007/s11128-023-03901-w.

## Additional Information

This article is based on Chapter four of the author’s PhD thesis [21].

## References

  * [1] A. Peruzzo, J. McClean, P. Shadbolt, M.-H. Yung, X.-Q. Zhou, P. J. Love, A. Aspuru-Guzik, J. L. O’brien, A variational eigenvalue solver on a photonic quantum processor, Nature Communications 5 (2014) 4213.
  * [2] S. M. Venkatesh, A. Macaluso, M. Klusch, Bilp-q: quantum coalition structure generation, in: Proceedings of the 19th ACM International Conference on Computing Frontiers, 2022, pp. 189–192.
  * [3] S. M. Venkatesh, A. Macaluso, M. Klusch, Gcs-q: Quantum graph coalition structure generation, arXiv preprint arXiv:2212.11372 (2022).
  * [4] J. Biamonte, P. Wittek, N. Pancotti, P. Rebentrost, N. Wiebe, S. Lloyd, Quantum machine learning, Nature 549 (7671) (2017) 195.
  * [5] M. Nakahara, T. Ohmi, Quantum computing: from linear algebra to physical realizations, CRC press, 2008.
  * [6] T. A. Gillespie, Spectral theory of linear operators, Proceedings of the Edinburgh Mathematical Society 23 (3) (1980).  doi:10.1017/S0013091500003886.
  * [7] P. Rebentrost, M. Mohseni, S. Lloyd, Quantum support vector machine for big data classification, Physical review letters 113 (13) (2014) 130503.
  * [8] A. Macaluso, L. Clissa, S. Lodi, C. Sartori, Quantum splines for non-linear approximations, in: Proceedings of the 17th ACM International Conference on Computing Frontiers, 2020, pp. 249–252.
  * [9] M. Benedetti, E. Lloyd, S. Sack, M. Fiorentini, Parameterized quantum circuits as machine learning models, Quantum Science and Technology 4 (4) (2019) 043001\.
  * [10] M. Schuld, A. Bocharov, K. Svore, N. Wiebe, Circuit-centric quantum classifiers, arXiv preprint arXiv:1804.00633 (2018).
  * [11] T. J. Hastie, R. J. Tibshirani, Generalized additive models, Vol. 43, CRC press, 1990.
  * [12] T. Hastie, R. Tibshirani, J. Friedman, The Elements of Statistical Learning, Springer Series in Statistics, Springer New York Inc., New York, NY, USA, 2001\.
  * [13] M. A. Nielsen, I. L. Chuang, Quantum Computation and Quantum Information: 10th Anniversary Edition, 10th Edition, Cambridge University Press, USA, 2011.
  * [14] A. Macaluso, L. Clissa, S. Lodi, C. Sartori, A variational algorithm for quantum neural networks, in: International Conference on Computational Science, Springer, 2020, pp. 591–604.
  * [15] A. Macaluso, F. Orazi, M. Klusch, S. Lodi, C. Sartori, A variational algorithm for quantum single layer perceptron, in: Machine Learning, Optimization, and Data Science, Springer Nature Switzerland, Cham, 2023, pp. 341–356.
  * [16] K. Hornik, M. Stinchcombe, H. White, et al., Multilayer feedforward networks are universal approximators., Neural networks 2 (5) (1989) 359–366.
  * [17] J. S. Judd, Neural network design and the complexity of learning, MIT press, 1990\.
  * [18] T. Goto, Q. H. Tran, K. Nakajima, Universal approximation property of quantum feature map, arXiv preprint arXiv:2009.00298 (2020).
  * [19] A. Macaluso, S. Lodi, C. Sartori, Quantum algorithm for ensemble learning, in: Proceedings of the 21st Italian Conference on Theoretical Computer Science, 2020\.
  * [20] A. Macaluso, L. Clissa, S. Lodi, C. Sartori, Quantum ensemble for classification, arXiv preprint arXiv:2007.01028 (2020).
  * [21] A. Macaluso, A novel framework for quantum machine learning, Ph.D. thesis (2021).

## Appendix

We consider the MAQA architecture when $d=3$. To make the notation simpler, we
indicate the parametrised quantum gate $G(\theta_{i,j})$ as $G_{i,j}$.

Without loss of generality, we can express the quantum gate $S_{\beta}$ as the
tensor product of $d=3$ unitary gates $B_{i}$. Then, the state preparation
step can be expressed as follows:

$\displaystyle\ket{\Phi_{0}}$ $\displaystyle=(S_{\beta}\otimes
S_{x})\ket{0,0,0}_{\text{control}}\otimes\ket{0,\dots,0}_{\text{data}}$
$\displaystyle=\left(\overset{d}{\underset{i=1}{\otimes}}B_{i}\otimes
S_{x}\right)\ket{0}^{\otimes 3}\otimes\ket{0}^{\otimes n}$
$\displaystyle=\overset{3}{\underset{i=1}{\otimes}}\left(a_{i}\ket{0}+b_{i}\ket{1}\right)\otimes\ket{x}=\overset{3}{\underset{i=1}{\otimes}}\ket{c_{i}}\otimes\ket{x},$
(13)

where $\ket{c_{i}}$ is the $i^{th}$ control qubit and $a_{i}$ and $b_{i}$ are
the parameters that determine its amplitudes:

$\displaystyle\ket{c_{i}}=a_{i}\ket{0}+b_{i}\ket{1}.$ (14)

Once the two registers are initialised, each qubit in the $control$ register
is entangled with two different random transformations of the $data$ register.
Thus, the first step after state preparation is the following:

#### Step 1 ($i=1$)

First, the controlled-unitary $C^{(1)}G_{1,1}$ is executed to entangle the
transformation $G_{1,1}\ket{x}$ with the excited state of $\ket{c_{3}}$:

$\displaystyle\ket{\Phi_{0,1}}$ $\displaystyle=\left[\mathbb{1}^{\otimes
2}\otimes C^{(1)}G_{1,1}\right]\ket{\Phi_{0}}$
$\displaystyle=\left[\mathbb{1}^{\otimes 2}\otimes
C^{(1)}G_{1,1}\right]\left(a_{3}\ket{0}+b_{3}\ket{1}\right)\otimes\ket{x}$
$\displaystyle=\ket{c_{1}}\otimes\ket{c_{2}}\otimes\left(a_{3}\ket{0}\ket{x}+b_{3}\ket{1}G_{1,1}\ket{x}\right).$
(15)

Then, a second controlled-unitary $C^{(0)}G_{1,2}$ is executed:

$\displaystyle\ket{\Phi_{1}}$ $\displaystyle=\left[\mathbb{1}^{\otimes
2}\otimes C^{(1)}G_{1,2}\right]\ket{\Phi_{1,1}}$
$\displaystyle=\overset{2}{\underset{i=1}{\otimes}}\ket{c_{i}}\otimes\left(a_{3}\ket{0}G_{1,2}\ket{x}+b_{3}\ket{0}G_{1,1}\ket{x}\right).$
(16)

At this point, two different transformations, $G_{1,1}$ and $G_{1,2}$ of the
initial state $\ket{x}$ are generated in superposition and are entangled with
the two basis states of the control qubit $\ket{c_{3}}$.

#### Step 2 ($i=2$)

The same operations are applied using $\ket{c_{2}}$ as control qubit and
different matrices, $G_{2,1}$ and $G_{2,2}$.

First, the controlled-unitary $C^{(1)}G_{2,1}$ is applied to entangle a
transformation of $\ket{x}$ with the excited state of $\ket{c_{2}}$:

$\displaystyle\ket{\Phi_{2,1}}$ $\displaystyle=\left(\mathbb{1}\otimes
C^{(1)}\otimes\mathbb{1}\otimes G_{2,1}\right)\ket{\Phi_{1}}$
$\displaystyle=\ket{c_{1}}\otimes\frac{1}{\sqrt{4}}\left[a_{2}\ket{0}\left(a_{3}\ket{0}G_{1,2}\ket{x}+b_{3}\ket{0}G_{1,1}\ket{x}\right)\right]+$
$\displaystyle\hskip
49.0pt+b_{2}\ket{1}\left(a_{3}\ket{0}G_{2,1}G_{1,2}\ket{x}+b_{3}\ket{1}G_{2,1}G_{1,1}\ket{x}\right),$
(17)

where the position of the gate $C^{(1)}$ indicates the control qubit used for
$G_{2,1}$. Then, a second controlled-unitary $C^{(0)}G_{2,2}$ is executed:

$\displaystyle\ket{\Phi_{2}}$ $\displaystyle=\left(\mathbb{1}\otimes
C^{(1)}\otimes\mathbb{1}\otimes G_{2,2}\right)\ket{\Phi_{2,1}}$
$\displaystyle=\ket{c_{1}}\otimes\frac{1}{\sqrt{4}}\Big{[}a_{2}\ket{0}\big{(}a_{3}\ket{0}G_{2,2}G_{1,2}\ket{x}+b_{3}\ket{1}G_{2,2}G_{1,1}\ket{x}\big{)}$
$\displaystyle\hskip
50.00008pt+b_{2}\ket{0}\big{(}a_{3}\ket{0}G_{2,1}G_{1,1}\ket{x}+b_{3}\ket{1}G_{2,1}G_{1,2}\ket{x}\big{)}\Big{]}.$
(18)

Notice that the entanglement performed in Step 2.1 influences the entanglement
in Step 2.2, and each trajectory describes a different transformation of
$\ket{x}$. Eq. (Step 2 ($i=2$)  ) can be rewritten expressing the four basis
states of the control register using natural numbers:

$\displaystyle\ket{\Phi_{2}}=\ket{c_{1}}\otimes$
$\displaystyle\frac{1}{\sqrt{4}}\Big{[}a_{2}a_{3}\ket{00}G_{2,2}G_{1,2}\ket{x}+a_{2}b_{3}\ket{01}G_{2,2}G_{1,1}\ket{x}$
$\displaystyle\hskip
10.00002pt+b_{2}a_{3}\ket{10}G_{2,1}G_{1,2}\ket{x}+b_{2}b_{3}\ket{11}G_{2,1}G_{1,1}\ket{x}\Big{]}$
$\displaystyle=\ket{c_{1}}\otimes\frac{1}{\sqrt{4}}\sum_{h=1}^{4}\beta^{*}_{h}\ket{h}G(\Theta_{h})\ket{x},$
(19)

where $G(\Theta_{h})$ is the product of $d=2$ unitaries $G_{i,j}$, the
coefficients $\beta^{*}_{h}$ result from the product of two coefficients
$a_{i}$ and $b_{i}$. Thus, using $2$ control qubits $4$ different quantum
trajectories are generated that correspond to $4$ different transformations of
data $\ket{x}$.

#### Step 3 ($i=3$)

Extending the same procedure when $d=3$, the result is the following:

$\displaystyle\ket{\Phi_{3}}=\frac{1}{\sqrt{8}}\Big{[}$
$\displaystyle\beta^{*}_{1}\ket{000}G_{3,2}G_{2,2}G_{1,2}\ket{x}+\beta^{*}_{2}\ket{001}G_{3,2}G_{2,2}G_{1,1}\ket{x}$
$\displaystyle+$
$\displaystyle\beta^{*}_{3}\ket{010}G_{3,2}G_{2,1}G_{1,2}\ket{x}+\beta^{*}_{4}\ket{011}G_{3,2}G_{2,1}G_{1,1}\ket{x}$
$\displaystyle+$
$\displaystyle\beta^{*}_{5}\ket{100}G_{3,1}G_{2,2}G_{1,2}\ket{x}+\beta^{*}_{6}\ket{101}G_{3,1}G_{2,2}G_{1,1}\ket{x}$
$\displaystyle+$
$\displaystyle\beta^{*}_{7}\ket{110}G_{3,1}G_{2,1}G_{1,2}\ket{x}+\beta^{*}_{8}\ket{111}G_{3,1}G_{2,1}G_{1,1}\ket{x}\Big{]}$
$\displaystyle=\frac{1}{\sqrt{8}}\sum_{h=1}^{8}\beta^{*}_{h}\ket{h}G(\Theta_{h})\ket{x},$
(20)

where each $G(\Theta_{h})$ is the product of $3$ unitaries $G_{i,j}$ for
$i=1,2,3$ and $j=1,2$.

Repeating this procedure $d$ times with different control qubits the result is
the following quantum state:

$\displaystyle\ket{\Phi_{d}}$
$\displaystyle=\frac{1}{\sqrt{2^{d}}}\sum_{h=1}^{2^{d}}\beta^{*}_{h}\ket{h}G(\Theta_{h})\ket{x}=\frac{1}{\sqrt{2^{d}}}\sum_{h=1}^{2^{d}}\beta^{*}_{h}\ket{h}\ket{l(x;\Theta_{h})},$
(21)

where each $G(\Theta_{h})$ is the product of $d=3$ unitaries $G_{i,j}$ for
$i=1,\cdots,d$ and $j=1,2$.

Finally, gate $F$ is applied, as shown in Eq. (3) and the measurement of the
data register is performed.