††thanks: Electronic address<EMAIL_ADDRESS>Electronic
address<EMAIL_ADDRESS>Electronic address:
<EMAIL_ADDRESS>
# Qsun: an open-source platform towards practical quantum machine learning
applications

Chuong Nguyen Quoc  Vietnamese-German University, Ho Chi Minh City, Vietnam,
70000  Le Bin Ho  Ho Chi Minh City Institute of Physics, National Institute of
Mechanics and Informatics, Vietnam Academy of Science and Technology, Ho Chi
Minh City, 700000, Vietnam  Research Institute of Electrical Communication,
Tohoku University, Sendai, Japan, 980-8577  Lan Nguyen Tran  Ho Chi Minh City
Institute of Physics, National Institute of Mechanics and Informatics, Vietnam
Academy of Science and Technology, Ho Chi Minh City, 700000, Vietnam  Hung Q.
Nguyen  Nano and Energy Center, VNU University of Science, Vietnam National
University, Hanoi, Vietnam, 120401

(August 28, 2024)

###### Abstract

Currently, quantum hardware is restrained by noises and qubit numbers. Thus, a
quantum virtual machine that simulates operations of a quantum computer on
classical computers is a vital tool for developing and testing quantum
algorithms before deploying them on real quantum computers. Various
variational quantum algorithms have been proposed and tested on quantum
virtual machines to surpass the limitations of quantum hardware. Our goal is
to exploit further the variational quantum algorithms towards practical
applications of quantum machine learning using state-of-the-art quantum
computers. In this paper, we first introduce a quantum virtual machine named
Qsun, whose operation is underlined by quantum state wavefunctions. The
platform provides native tools supporting variational quantum algorithms.
Especially using the parameter-shift rule, we implement quantum differentiable
programming essential for gradient-based optimization. We then report two
tests representative of quantum machine learning: quantum linear regression
and quantum neural network.

Keywords: quantum virtual machine, quantum machine learning, quantum
differentiable programming, quantum linear regression, quantum neural network

##  I Introduction

The advent of quantum computers has opened a significant turning point for
exponentially speeding up computing tasks that classical computers need
thousand years to execute Arute2019 ; ChinaSupremacy . Although mankind has
witnessed tremendous development in this field theoretically and
experimentally in the last few years, most state-of-the-art quantum computers
still rely on noisy intermediate-scale quantum computers NISQ. Noises and
qubit-number constraints prevent to build high-fidelity quantum computers
capable of substantially implementing quantum algorithms
Preskill2018quantumcomputingin ; 7927034 ; Wang285 ; Kandala2019 . To bypass
these constraints, various hybrid quantum-classical algorithms that use
classical computers to optimize quantum circuits have been proposed Cerezo2021
; PhysRevA.101.010301 ; PhysRevA.103.L030401 . Among these, variational
quantum algorithms (VQAs) may be the most promising ones in the NISQ era.

VQAs generally consist of three essential steps: (i) initializing quantum
states for a given wavefunction ansatz, (ii) measuring a cost function
suitable for problems being considered, and (iii) minimizing the cost function
and updating new parameters. The self-consistency is performed until
convergence. VQAs have been extensively employed to tackle numerous tasks,
including the variational quantum eigensolvers (VQEs) Peruzzo2014 ;
PhysRevResearch.1.033062 ; Kirby2021contextualsubspace ; Ryabinkin2019 ;
PhysRevA.101.052340 ; Gard2020 ; PRXQuantum.2.020337 , quantum dynamics
simulation PhysRevX.7.021050 ; McArdle2019 ; Nishi2021 ; PRXQuantum.2.030307 ;
zhang2020lowdepth ; PhysRevResearch.3.023095 , mathematical applications
PhysRevA.103.052425 ; huang2019nearterm ; PhysRevLett.122.060504 ;
PhysRevA.103.052416 ; bravoprieto2020variational ; xu2019variational ;
PRXQuantum.2.010315 ; Lloyd2014 ; LaRose2019 , quantum machine learning (QML)
PhysRevA.98.032309 ; farhi2018classification ; PhysRevLett.122.040504 ;
Havlicek2019 ; PhysRevA.101.032308 ; farhi2018classification ; Cong2019 ;
pesah2020absence ; zhang2020trainability ; Beer2020 , and new frontiers in
quantum foundations Arrasmith2019 ; wilde_2017 ; Cerezo2021 ;
PhysRevA.101.062310 ; Koczor_2020 ; PhysRevLett.123.260505 ;
beckey2020variational ; Meyer2021 .

Typically, VQAs employ variational quantum circuits to measure the cost
function on a quantum computer and outsource its optimization to a classical
computer. While one can manipulate gradient-free optimizers, such as Nelder-
Mead simplex Mead , to minimize the cost function, using gradient-based
methods like gradient descent can help us speed up and guarantee the
convergence of the optimization. Several quantum algorithms have been proposed
to evaluate the cost function gradient measured on quantum computers
10.1145/3385412.3386011 ; Sch_fer_2020 ; Stokes2020quantumnatural ;
lopatnikova2021quantum ; yamamoto2019natural ; PhysRevLett.126.140502 . Among
those methods, quantum differentiable programming (QDP) has been introduced
and utilized extensively 10.1145/3385412.3386011 ; Sch_fer_2020 ;
guerreschi2017practical ; bergholm2020pennylane ; PhysRevLett.118.150503 ;
PhysRevA.98.032309 ; PhysRevA.99.032331 ; Banchi2021measuringanalytic . It
relies on a technique called the parameter-shift rule that evaluates the
derivative of any differentiable function using quantum circuits
guerreschi2017practical ; bergholm2020pennylane ; PhysRevLett.118.150503 ;
PhysRevA.98.032309 ; PhysRevA.99.032331 ; Banchi2021measuringanalytic .
Therefore, this method is beneficial for developing “on-circuit” gradient-
based optimization techniques, especially for quantum machine learning (QML)
applications where various methods like quantum neural networks (QNNs) demand
the derivative information of the cost function.

While quantum algorithms should be performed on quantum computers, the current
limitation of NISQ computers cause challenges in developing and testing new
quantum algorithms, demanding the use of virtual alternatives called quantum
virtual machines (QVMs). Besides, QVMs are necessary for modeling various
noisy channels to characterize the noises and the efficiency of quantum error
correction. One can classify QVMs into two types according to the way to build
them: (i) the matrix multiplication approach qiskit ;
cirq_developers_2021_5182845 ; smith2017practical ; bergholm2020pennylane ;
7927034 ; 10.1145 ; Villalonga2019 ; roberts2019tensornetwork , and (ii) the
wavefunction approach suzuki2020qulacs ; Guerreschi_2020 ;
Steiger2018projectqopensource ; kelly2018simulating ; Jones2019 ;
DERAEDT2007121 ; DERAEDT201947 . While the former performs matrix
multiplication for all qubits in quantum circuits, the latter represents
quantum circuits by corresponding wavefunctions. On the one hand, the former
can significantly reduce the memory capacity using tensor network contraction
doi:10.1137/050644756 ; PhysRevLett.125.060503 , and the cost to pay is
exponentially increasing the computational time. On the other hand, the
latter, in principle, can require less computational time in a limit of qubits
number and a sufficiently large memory size that can store the total quantum
wavefunction. Therefore, both approaches exist side by side and a possible
hybrid approach Arute2019 ; markov2018quantum  for convenient purposes.

There are several QVM’s libraries developed for QML orientation, such as
TensorFlow Quantum library broughton2020tensorflow  implemented in Cirq
cirq_developers_2021_5182845 , Pennylane bergholm2020pennylane  designed for
photonics devices, and TensorNetwork roberts2019tensornetwork ; Huggins_2019 .
These libraries are all constructed in the matrix-multiplication type of QVMs,
designed to submit quantum tasks to the developed hardware conveniently.

In this work, we develop a QVM platform named Qsun using the wavefunction
approach towards the QML applications. In Qsun, a quantum register is
represented by a wavefunction, and quantum gates are manipulated directly by
updating the amplitude of the wavefunction. Measurement results rely on
probabilities of the wavefunction. Our simple approach yields faster
computation speed for a small number of qubit when compared to other QVMs such
as Qiskit, ProjectQ, or Pennylane. Basing on this generic QVM, we aim to
exploit the advantages of QDP with the parameter-shift rule as the core engine
towards practical applications in quantum machine learning. Two representative
examples of QML are demonstrated: quantum linear regression and quantum neural
network. These algorithms are compared to standard progams: Qiskit, ProjectQ,
and Pennylane, and classical algorithms when applicable. In these comparisons,
Qsun performs slightly better for QDP, QLR, and QNN. All in all, Qsun is an
efficient combination of QVM with QDP features that is oriented toward machine
learning problems. In the following, we introduce the QVM platform Qsun and
its performance compared to others in Sec. II, followed by an introduction to
the QDP implementation within Qsun in Sec. III. We then discuss some QML
applications of the Qsun package in Sec. IV.

##  II Quantum virtual machine implementation in Qsun

In practice, quantum computers work on quantum algorithms by composing a
quantum register, or qubits, operated by a sequence of quantum gates. Results
are then traced out from quantum measurements. We now introduce our quantum
virtual machine (QVM) named Qsun, an open-source platform simulating the
operation of a generic quantum computer QSUN . We aim the platform to the
development of quantum machine learning (QML) and related problems. We develop
it in Python and employ the Numpy library for fast matrix operations and
numerical computations.

###  II.1 Simulating quantum computers using wavefunction basis

Unlike widely-used approaches based on matrix multiplication qiskit ;
cirq_developers_2021_5182845 ; smith2017practical ; bergholm2020pennylane ;
7927034 ; 10.1145 ; Villalonga2019 ; roberts2019tensornetwork , our platform
is developed using the class of “wavefunction” approach suzuki2020qulacs ;
Guerreschi_2020 ; Steiger2018projectqopensource ; kelly2018simulating ;
Jones2019 ; DERAEDT2007121 ; DERAEDT201947 , in which a quantum register is
represented by its wavefunction. The operation of quantum gates is simulated
by updating the wavefunction’s amplitude, and output results are obtained by
measuring wavefunction’s probabilities. We expect that working directly on
wavefunction is beneficial for QML applications, especially for building and
training variational quantum circuits in quantum neural networks (QNNs). As
depicted in Fig. 1, Qsun consists of three main modules Qwave, Qgates, and
Qmeas for quantum register, quantum gates, and quantum measurement,
respectively.

#### Qwave

In general, for a quantum register with $N$ qubits, its quantum states are
represented in the $2^{N}$-dimension Hilbert space as

$\displaystyle|\psi\rangle=\sum_{j=0}^{2^{N}-1}\alpha_{j}|j\rangle,$ (1)

where $\alpha_{j}$ are complex amplitudes obeying a completeness relation
$\sum_{j=0}^{2^{N}-1}|\alpha_{j}|^{2}=1$, and vectors $|j\rangle$ are elements
of the computational basis. We integrate quantum state’s information into the
class Wavefunction as described in Fig. 1 and Table. 1. The class allows us to
access and update amplitudes directly according to the evolution of the
quantum state under the action of the unitary quantum gates. It also measures
probabilities that contain output information.

Figure 1: Qsun quantum virtual machine (QVM) consists of three main modules Qwave, Qgates, and Qmeas, to simulate quantum circuits with a quantum register (qubits), quantum gates, and quantum measurements, respectively. In Qsun, a quantum register is constructed by its wavefunction through the Wavefunction class and its containing objects, as illustrated in the figure. Table 1:  List of methods used in the class Wavefunction. Methods |  Description
---|---
amplitude |  a NumPy array of complex numbers that stores the amplitudes of quantum states.
state |  a NumPy array of strings that labels the quantum states’s basis.
probabilities |  return a list of corresponding probabilities for each basis vector in the superposition.
print_state |  return a string representing a quantum state of the system in bra-ket notations.
visual_circuit |  print a visualization of a quantum circuit.

#### Qgates

To manipulate a single-qubit gate $U$ $=\begin{pmatrix}a&c\\\
b&d\end{pmatrix}$ acting on the $n^{\rm th}$ qubit, the nonzero elements in
amplitude array are updated as Jones2019

$\displaystyle\begin{pmatrix}\alpha_{s_{i}}\\\
\alpha_{s_{i}+2^{n}}\end{pmatrix}\to U\begin{pmatrix}\alpha_{s_{i}}\\\
\alpha_{s_{i}+2^{n}}\end{pmatrix},$ (2)

where $s_{i}={\rm floor}(i/2^{n})2^{n+1}+(i\text{ mod }2^{n})$, for all
$i\in[0,2^{N-1}-1]$. Here, ${\rm floor}(x)$ is the standard floor function
taking the greatest integer less than or equal to the real of $x$. We unify
the implementation of single- and multiple-qubit gates into a common
framework. We outline the operation of single-qubit gates in Algorithm 1 and
an example of the Hadamard gate in Algorithm 2. We emphasize that Qgates only
update nonzero components of wavefunction amplitudes. This way, we can avoid
demanding matrix multiplications that escalate exponentially ($2^{N}\times
2^{N}$) with the number of $N$ qubits. This generic algorithm allows us to
implement arbitrary unitary gates without decomposing them into universal
ones, which may be advantageous to model a general class of neural networks
using quantum circuits.

Result: Wavefunction with new probability amplitudes

w $\leftarrow$ Wavefunction;

states $\leftarrow$ w.state; ampl $\leftarrow$ w.amplitude;

$N$ $\leftarrow$ size(state[$i$][0]); $n$ $\leftarrow$ target qubit;

$n$_ampl $\leftarrow$ [0, …, 0], size($n$_ampl) = size(ampl);

$cut\leftarrow$ $2^{N-n-1}$;

for  _i $\leftarrow$ 0, size(ampl)_ do

if  _state[ $i$][$n$] == 0_ then

$n$_ampl[$i$] $\leftarrow$ $n$_ampl[$i$] + $a$*ampl[$i$];

$n$_ampl[$i+cut$] $\leftarrow$ $n$_ampl[$i+cut$] + $b$*ampl[$i$];

else

$n$_ampl[$i$] $\leftarrow$ $n$_ampl[$i$] + $d$*ampl[$i$];

$n$_ampl[$i-cut$] $\leftarrow$ $n$_ampl[$i-cut$] + $c$*ampl[$i$];

end if

end for

w.amplitude $\leftarrow n$_ampl

Algorithm 1 Operation of a single-qubit gate:

Result: Wavefunction with new probability amplitudes corresponding to Hadamard
state

w $\leftarrow$ Wavefunction;

states $\leftarrow$ w.state; ampl $\leftarrow$ w.amplitude;

$N$ $\leftarrow$ size(state[$i$][0]); $n$ $\leftarrow$ target qubit;

$n$_ampl $\leftarrow$ [0, …, 0], size($n$_ampl) = size(ampl);

$cut\leftarrow$ $2^{N-n-1}$;

for  _i $\leftarrow$ 0, size(ampl)_ do

if  _state[ $i$][$n$] == 0_ then

$n$_ampl[$i$] $\leftarrow$ $n$_ampl[$i$] + $(1/\sqrt{2})$*ampl[$i$];

$n$_ampl[$i+cut$] $\leftarrow$ $n$_ampl[$i+cut$] + $(1/\sqrt{2})$*ampl[$i$];

else

$n$_ampl[$i$] $\leftarrow$ $n$_ampl[$i$] + $(-1/\sqrt{2})$*ampl[$i$];

$n$_ampl[$i-cut$] $\leftarrow$ $n$_ampl[$i-cut$] + $(1/\sqrt{2})$*ampl[$i$];

end if

w.amplitude $\leftarrow n$_ampl

end for

Algorithm 2 Operation of Hadamard gate:

To mimic the actual operation of quantum computers, we introduce noises into
the wavefunctions. In Qsun, the standard quantum depolarizing channel is
implemented as a single-qubit gate $\mathcal{E}$ that is a part of quantum
circuits acting on the wavefunctions. For a given noisy probability $p$,
applying the gate $\mathcal{E}$ on a mixed quantum state $\rho$ will transform
it to nielsen_chuang_2019

$\displaystyle\rho\xrightarrow{\mathcal{E}}(1-p)\rho+\dfrac{p}{2}I,$ (3)

where $I$ is an $2\times 2$ identity matrix. In general, one can decompose the
depolarizing channel into the bit-flip, phase-flip, and phase-bit-flip as

$\displaystyle\rho\xrightarrow{\mathcal{E}}(1-p)\rho+p_{x}X\rho X+p_{y}Y\rho
Y+p_{z}Z\rho Z,$ (4)

where $p_{x},p_{y},p_{z}$ are the probabilities of bit-flip, phase-bit-flip,
and phase-flip, respectively nielsen_chuang_2019 ; 7927034 . In Qsun, we use
$p_{x}=p_{y}=p_{z}=p/3$ and apply to every qubits in the circuit each after
the action of a quantum gate on the circuit.

#### Qmeas

The module Qmeas is designed to execute quantum measurements on a single qubit
or all qubits in the quantum circuit. For a measurement on a single qubit $n$,
the probability for that its outcome is $|0\rangle$ reads

$\displaystyle p(0)=\sum_{i=0}^{2^{N-1}-1}\langle\psi|j_{s_{i}}\rangle\langle
j_{s_{i}}|\psi\rangle=\sum_{i=0}^{2^{N-1}-1}|\alpha_{s_{i}}|^{2},$ (5)

with $|j_{s_{i}}\rangle$ as the basis element and the post-quantum state after
the measurement is given as

$\displaystyle|\psi^{\prime}\rangle=\dfrac{\sum_{i=0}^{2^{N-1}-1}|j_{s_{i}}\rangle\langle
j_{s_{i}}|\psi\rangle}{\sqrt{p(0)}}.$ (6)

Similarly, the probability for getting the outcome $|1\rangle$ reads

$\displaystyle p(1)=1-p(0),$ (7)

and the post-state is

$\displaystyle|\psi^{\prime}\rangle=\dfrac{\sum_{i=0}^{2^{N-1}-1}|j_{s_{i}+2^{n}}\rangle\langle
j_{s_{i}+2^{n}}|\psi\rangle}{\sqrt{p(1)}}.$ (8)

For all-qubit measurement, the post-quantum state will collapse to one of
$\\{|j\rangle\\}$ with the probability of $|\alpha_{j}|^{2}$. In Qsun, we
build these two measurements onto measure_one and measure_all, respectively.

Figure 2: Comparing computational time $t$ (in seconds) of different QVMs as
increasing the number of qubits $N$. Inset: The testing circuit including an
$H$, a $\sqrt{X}$, and a CNOT gates acts on each qubit. The depth is fixed at
10, and the number $N$ of qubits is varied to assess the QVM performance.

###  II.2 Assessing the QVM performance

Let us now assess the performance of Qsun and compare it with three existing
ones in Pennylane, Qiskit, and ProjectQ. Note that Qsun and ProjectQ belong to
the wavefunction class, while the others are in the matrix multiplication
class. We have adopted the testing circuit from Ref. quantest  composed of the
Hadamard, $\sqrt{X}$, and CNOT gates acting on each qubit. We have fixed the
depth of the circuit at 10 and varied the number of qubits $N$. The code for
this test is shown in A.

Fig. 2 represents the change of computational time when the number $N$ of
qubits increases. In general, there are two magnitudes of slope corresponding
to two ways of QVM implementation. While the wavefunction-based approach is
faster than the matrix-multiplication one for small numbers of qubits ($N<8$),
the opposite behavior is observed for larger numbers of qubits. This
observation reflects the basic properties of these two approaches as discussed
above and see also Ref. suzuki2020qulacs . However, there are some available
techniques to improve the performance of the wavefunction approach for larger
numbers of qubits, such as, SIMD (single-instruction, multiple data)
optimization and multi-threading suzuki2020qulacs . We further summarize a
comparison between Qsun and other simulators in terms of practical quantum
algorithms in Table 2.

Table 2: Comparison between Qsun and other simulators Algorithms | Qsun | ProjectQ | Qiskit | Pennylane
---|---|---|---|---
Standard algorithms | ✓ | ✓ | ✓ | ✓
Quantum diferentiable programming | ✓ | ✓ | ✓ | ✓
Quantum Linear Regression | ✓ | ✓ | ✓ | ✓
Quantum Neural Network | ✓ |  | ✓ | ✓

##  III Quantum differentiable programming implementation in Qsun

Given a quantum state $|\psi(\vec{\theta})\rangle$ with $\vec{\theta}$ as
variational parameters and an observable $\widehat{C}$, the task is to seek
the global minimum of the expectation value
$C(\vec{\theta})=\langle\psi(\vec{\theta})|\widehat{C}|\psi(\vec{\theta})\rangle$
with respect to parameters $\vec{\theta}$. For example, if $\widehat{C}$ is a
Hamiltonian, its global minimum is the ground state energy. In general,
$C(\vec{\theta})$ is called as the cost function, and minimizing the cost
function requires its derivative with respect to parameters $\vec{\theta}$,
$\partial C(\vec{\theta})/\partial\theta$. In classical computing, if the
analytical form of $\partial C(\vec{\theta})/\partial\theta$ is unknown,
finite difference methods are often used to evaluate the derivative
approximately. Although this approximation is fast and easy to implement, its
accuracy depends on discretization steps. In contrast to the classical finite
differentiation, quantum differentiable programming (QDP) is an automatic and
exact method to compute the derivative of a function. QDP is thus essential
for accurate gradient computation in multiple VQAs, including QML models.

Result: Derivative of a function

$f\leftarrow$ Function; $c\leftarrow$ Quantum Circuit;

$p\leftarrow$ Params; $s\leftarrow$ Shift;

diff $\leftarrow$ [0, …, 0], size(diff) = size($p$);

for  _i $\leftarrow$ 0, size($diff$)_ do

$p$_plus $\leftarrow$ copy($p$);

$p$_subs $\leftarrow$ copy($p$);

$p$_plus[$i$] $\leftarrow$ $p[i]+s$;

$p$_subs[$i$] $\leftarrow$ $p[i]-s$;

diff[$i$] $\leftarrow(f(c,p$_plus) $-f(c,p$_subs))/(2*sin($s$));

end for

Algorithm 3 Quantum differentiable programming implementation in Qsun

The heart of QDP is the parameter-shift rule that is analytically computed
using quantum circuits. The algorithm is outlined in Algorithm 3. Let us
introduce a parameterized generator $\widehat{V}$ independent of
$\vec{\theta}$ such that
$|\psi(\vec{\theta})\rangle=e^{i\vec{\theta}\widehat{V}}|\psi\rangle$. The
cost function is then rewritten as

$\displaystyle
C(\vec{\theta})=\langle\psi|e^{-i\vec{\theta}\widehat{V}}\widehat{C}e^{i\vec{\theta}\widehat{V}}|\psi\rangle={\rm
Tr}\bigl{(}\widehat{C}e^{\mathcal{Z}}[\rho]\bigr{)},$ (9)

with $\rho=|\psi\rangle\langle\psi|$; $\mathcal{Z}=i\vec{\theta}\widehat{V}$,
and the superoperator $e^{\mathcal{Z}}[\rho]=e^{\mathcal{-Z}}\rho
e^{\mathcal{Z}}$ miller_1990 . The parameter-shift rule for each
$\theta\in\vec{\theta}$ states that

$\displaystyle\dfrac{\partial C(\vec{\theta})}{\partial\theta}={\rm
Tr}\bigl{(}\widehat{C}\partial_{\theta}e^{\mathcal{Z}}[\rho]\bigr{)}=c\bigl{[}C(\vec{\theta}+s)-C(\vec{\theta}-s)\bigr{]},$
(10)

where $c=1/(2\sin(s))$, and $s$ is determined based on the superoperator and
independent of $\vec{\theta}$. The values of the cost function $C$ at
$\vec{\theta}\pm s$ are measured on quantum computers by implementing two
quantum circuits as follows

The derivative is finally obtained by subtracting measurement results from the
two circuits. An advantage of the parameter-shift rule is that it can compute
the derivative of the given function exactly while the shift $s$ can be chosen
arbitrarily large quanshiftrule .  Figure 3: Computational time $t$ in seconds
of the QDP implementation in Qsun, ProjectQ using an engine Simulator() and
Pennylane using a default.qubit device when the number of qubits increases.
Inset: quantum circuits used for the QDP test. We combine the QDP implemented
in Qsun with a classical optimizer like adaptive moment (Adam) optimization
and implement quantum gradient descent, making the optimization procedure
fully quantum. Let us now examine the performance of the QDP in comparison
with Pennylane, another QML-oriented platform. In this test, we measure the
computational time spent for optimizing the ground state of an Ising model.
The code for running this test is given in B, and a tutorial explaining how to
use QDP with sample codes is written in C. Fig 3 presents the change of
computational time when increasing the number of qubits $N$ for three cases of
Qsun, Pennylane, and ProjectQ. As shown in the inset of Fig. 3, for each
qubit, we have two variational parameters $\vec{\theta}$, meaning that the
number of variational parameters is twice of $N$. Qsun is faster than
Pennylane for $N<5$, whereas it is slower for larger numbers of qubits
($N>5$). Among the three cases, ProjectQ consumes the most computational time
for all $N$. Of course, the exponential growth with the number of qubits is
not inevitable.

##  IV Quantum machine learning applications using Qsun

Various quantum machine learning models can be developed with quantum
differentiable programming (QDP) implementations to evaluate the gradient and
employ gradient-based optimization. This section demonstrates QML applications
using Qsun in two well-known models: quantum linear regression (QLR) and
quantum neural network (QNN). Its performances are compared to other standard
tools.

Before digging into detailed examples, let us derive the QDP implementation in
derivative of a standard mean squared error cost function:

$\displaystyle
C(\vec{\theta})=\frac{1}{M}\sum_{i=0}^{M-1}\Big{(}y_{i}-\widehat{y}_{i}\Big{)}^{2},$
(11)

where $\vec{\theta}$ is a tuple of variational parameters, $y_{i}$ represents
the true value and $\widehat{y}_{i}$ stands for prediction value, with
$i\in[0,M-1],M$ the number of samples in the dataset to be trained.
Concretely, we also consider the prediction value is a composition function of
an activation function $f(\vec{\theta})$, i.e.,
$\widehat{y}_{i}(f(\vec{\theta}))$. Then, we derive at a chain rule:

$\displaystyle\frac{\partial
C(\vec{\theta})}{\partial\theta}=\frac{1}{M}\sum_{i=0}^{M-1}\frac{\partial
C(\widehat{y}_{i})}{\partial\widehat{y}_{i}}\frac{\partial\widehat{y}_{i}(f)}{\partial
f}\frac{\partial f(\vec{\theta})}{\partial\theta},$ (12)

where $\theta\in\vec{\theta}$, with

$\displaystyle\frac{\partial
C(\widehat{y}_{i})}{\partial\widehat{y}_{i}}=-2[y_{i}-\widehat{y}_{i}],$ (13)

the derivative $\frac{\partial\widehat{y}_{i}(f)}{\partial f}$ depends on the
particular form of the activation function, and

$\displaystyle\frac{\partial f(\vec{\theta})}{\partial\theta}$
$\displaystyle=\sum_{k}\frac{\partial m_{k}(\vec{\theta})}{\partial\theta}$
$\displaystyle=\sum_{k}c_{k}\Bigl{[}m_{k}(\vec{\theta}+s_{k})-m_{k}(\vec{\theta}-s_{k})\Bigr{]},$
(14)

where
$m_{k}(\vec{\theta})=\langle\psi(\vec{\theta})|\mathcal{\widehat{M}}_{k}|\psi(\vec{\theta})\rangle$
is a measurement outcome of an operator $\mathcal{\widehat{M}}_{k}$ for a
measurement set $\\{k\\}$. To arrive at the second equality in Eq. (IV), we
have applied the parameter-shift rule (10).

###  IV.1 Quantum linear regression

We implement quantum linear regression on a “diabetes” dataset diabetes
available in the scikit-learn package scikit-learn  with 400 samples for the
training set and 10 samples for the testing set. We write the linear
regression model in the form

$\widehat{y}=wx+b,$ (15)

where $w$ and $b$ are the slope and intercept of the linear regression need to
be obtained. To evaluate them on quantum computers, we store their values in
two qubits. The values of $w$ and $b$ now become the expectation values of the
Pauli matrix $\widehat{\sigma}_{z}$. The quantum version of linear regression
model (15) states

$\widehat{y}(\vec{\theta})=k\big{(}\langle\widehat{w}\rangle
x+\langle\widehat{b}\rangle\big{)},$ (16)

where $k$ is the scaling factor that transforms the output data from $[-1,1]$
into $[-k,k]$. Its value should be chosen so that the R.H.S can cover all of
the true values $\\{y_{i}\\}$. Here, $w,b$ are functions of variational
parameters $\vec{\theta}$, for example $w\equiv
w(\vec{\theta})=\langle\psi(\vec{\theta})|\widehat{\sigma}_{z}|\psi(\vec{\theta})\rangle$
with the initial state $|\psi\rangle=|0\rangle$, (see the inset Fig. 4). The
cost function is the mean squared error as described in Eq. (11).

After having the model, we run the circuit in the inset Fig. 4 to find the
optimized values for $w$ and $b$ via updating $\vec{\theta}$. We first
calculate their derivatives using the chain rule as shown in Eq. (12) and the
parameter-shift rule in the QDP as depicted in Algorithm 3, and then update
new parameters using gradient descent

$\displaystyle\theta_{\rm t+1}=\theta_{\rm t}-\eta\frac{\partial
C(\vec{\theta}_{\rm t})}{\partial\theta},$ (17)

where $\eta$ the learning rate, $\vec{\theta}_{\rm t}$ is the set of
parameters at time step t. They are continuously updated until the optimized
values are found. The quantum circuit that used to train the model is
represented in the inset Fig. 4, and the code is given in D.

For comparison, we next provide the analytical solution for minimizing the
cost function (11) with (15). We vectorize:

$\displaystyle\begin{pmatrix}w\\\ b\end{pmatrix}=(X^{T}X)^{-1}X^{T}Y$ (18)

where

$\displaystyle X=\begin{pmatrix}x_{0}&1\\\ \vdots&\vdots\\\
x_{M-1}&1\end{pmatrix},\quad Y=\begin{pmatrix}y_{0}\\\ \vdots\\\
y_{M-1}\end{pmatrix},$ (19)

and $T$ the transpose operator.

The performance of Qsun is compared to those of Pennylane, ProjectQ, and
analytical solution using the same set of parameters. We use the maximum
number of iterations 1000 and $k=10$. For Pennylane, we use
GradientDescentOptimizer() with its default configuration. For ProjectQ, we
use the same optimize algorithm as we have used for Qsun, which is shown in D.
As we can see from Fig. 4, Qsun’s result is closer to the analytical result
than the Pennylane and ProjectQ one, which implies a high performance of Qsun.
In the same case of the “wavefunction” approach, Qsun has a better efficiency
than ProjectQ in the speedup as shown in Fig. 2 and Fig. 3, and in the
optimization process.

Figure 4: Linear regression models trained by quantum programming using Qsun
(blue), Pennylane (green), ProjectQ (orange) and analytical solution (dashed
purple). We use the diabetes dataset diabetes  with 400 samples for the
training set and 10 samples for the testing set. The boundary $k$ = 10 with
1000 iterations. Inset: the quantum circuit that trains the quantum linear
regression models.

###  IV.2 Quantum neural network

In this subsection, we show the ability of Qsun for deep learning by building
and training a quantum neural network (QNN). We model the QNN as a variational
quantum circuit (VQC) parameterized with multiple variables that are optimized
to train the model. Fig. 5 represents the process for building and training
our QNNs. It is a hybrid quantum-classical scheme with five steps summarized
as follows:

  * •

The quantum part:

  * Step 1: Quantizing and encoding dataset into quantum states using the amplitude encoding method PhysRevA.101.032308 ; PhysRevA.102.032420 .

  * Step 2: Building a QNN circuit and measurement.

  * Step 3: Evaluating the derivative of measurement results using QDP.

  * •

The classical part:

  * Step 4: Deriving the derivative of the defined cost (loss) function.

  * Step 5: Running a gradient-based optimization and updating parameters.

Figure 5: A Quantum Neural Network framework. The blue box is the quantum part
of the QNN, and the green box is the classical part of the QNN. At first,
quantify and encode the dataset (training and testing) into quantum states of
qubits $|\psi_{i}\rangle$. Each feature in the dataset encodes into one qubit.
The employed rotation gates will parameterize the quantum circuit, and the
CNOT gates cause entangled in the circuit. These gates repeat L times for L
layers in the quantum neural network. After that, we measure the circuit and
give the corresponding probabilities $p_{k}$. We employ a QDP scheme with the
pentameter-shift rule to calculate the derivative of $p_{k}$ and send the
results to a classical computer to derive the derivative of the loss function.
After that, we implement a gradient-based optimization to obtain new
parameters. When the scheme is not optimal yet, we update the circuit with new
parameters; when it is optimal, we turn it to the testing process.

#### Data quantization and amplitude encoding

Set $x_{j}^{(i)}$ as the $i^{\rm th}$ sample ($i\in[0,M-1]$) of the $j^{\rm
th}$ feature ($j\in[0,N-1]$) in the dataset of $N$ features with $M$ samples
for each feature. Since we only analyze one sample (other samples are treated
similarly), we can omit the indicator $i$ and rewrite $x_{j}^{(i)}$ as
$x_{j}$. We now map $x_{j}$ into a qubit by normalizing its value in range
$[x_{j}^{\rm(min)},x_{j}^{\rm(max)}]$ into $[0,1]$. Using the Min-Max
normalization, we obtain

$\displaystyle\widetilde{x}_{j}=\dfrac{x_{j}-x_{j}^{\rm(min)}}{x_{j}^{\rm(max)}-x_{j}^{\rm(min)}}.$
(20)

Here, the normalized value $\widetilde{x}_{j}\in[0,1]$. We map this value into
the amplitudes of a qubit as
$|\psi_{j}\rangle=\sqrt{\widetilde{x}_{j}}|0\rangle+\sqrt{1-\widetilde{x}_{j}}|1\rangle.$
Then, the quantum state for $N$ features reads

$\displaystyle|\Psi\rangle=\bigotimes^{N-1}_{j=0}\ket{\psi_{j}}=\bigotimes^{N-1}_{j=0}\left(\sqrt{\widetilde{x}_{j}}\ket{0}+\sqrt{1-\widetilde{x}_{j}}\ket{1}\right).$
(21)

The encoding implementation is given in Algorithms 4:

Result: Probablity amplitudes of QNN’s initial quantum states

sample $\leftarrow$ sample_scaled;

$N$ $\leftarrow$ number_of_features;

ampl $\leftarrow$ [0, …, 0], size(ampl) = $N$;

for  _i $\leftarrow$ 0, N-1_ do

ampl[$i$] $\leftarrow$ [sqrt(sample[$i$]), sqrt(1-sample[$i$])];

end for

ampl_vec $\leftarrow$
$\text{ampl[0]}\otimes\text{ampl[1]}\otimes\cdots\otimes\text{ampl[$N$-1]}$;

Algorithm 4 Encoding data in Qsun

#### Building the QNN circuits

The QNN model uses $N$ qubits for $N$ features and includes multiple layers.
It is parameterized through a set of variables
$\vec{\theta}=\\{\theta_{kjj^{\prime}}\\}$ with $k$ as the layer index, $j$ as
qubit (feature) indices, and $j^{\prime}$ indicates the number of rotation
gates implementing on each qubit. Each layer has one $R_{x}$ and $R_{y}$ gates
for each qubit followed by CNOT gates generating all possible entanglements
between them if the number of qubits is more than one.

#### Decoding probabilities into predictions

Here, we map measurement results from the previous step into classical
predictions by using an activation function. We consider the expectation value
of the projection operator $\widehat{\Pi}=|1\rangle\langle 1|$ as

$\displaystyle
p^{(i)}_{j}\equiv\braket{\widehat{\Pi}}=\langle\psi^{(i)}_{j}|1\rangle\langle
1|\psi^{(i)}_{j}\rangle=|\langle 1|\psi^{(i)}_{j}\rangle|^{2}.$ (22)

Note that $|\psi^{(i)}_{j}\rangle$ is the final state of qubit $j$ at sample
$i$. We also emphasize that one can choose the projection operator
$\widehat{\Pi}$ arbitrarily. We use the sigmoid function, $S(x)=1/(1+e^{-x})$
to transform the measurement data into predictions. In our concrete example
below, we have two features represented by two qubits, i.e, $j=0,1$. We use
the prediction rule as

$\displaystyle S_{i}(p)=S\bigl{(}\gamma(p^{(i)}_{0}-p^{(i)}_{1})\bigr{)},$
(23)

where $\gamma$ is a scaling shape for the sigmoid function such that for large
$\gamma$ then $S(x)$ becomes a Heaviside step function. To use a soft
prediction, we introduce a label-conditional probability for a prediction
value as

$\widehat{y}_{i}(l)=\left\\{\begin{array}[]{lr}S_{i}(p),&\text{for }l=0\\\
1-S_{i}(p),&\text{for }l=1\end{array}\right..$ (24)

#### Training the QNN model

In the current version of Qsun, we are using a quantum-classical hybrid scheme
combining QDP and the classical Adam optimization to train our QNN model. The
cost function is defined by

$C(\vec{\theta})=\frac{1}{M}\sum_{i=0}^{M-1}\bigl{[}1-\widehat{y}_{i}(y_{i})\bigr{]}^{2},$
(25)

which is a function of $\vec{\theta}$. Here, $y_{i}\in\\{0,1\\}$ are true
values, $\widehat{y}_{i}(y_{i})$ are predicted probability conditioned on the
label $y_{i}$. We then evaluate the derivative of the cost function using QDP,
followed by a gradient-descent procedure to search for optimal parameters
$\vec{\theta}$.

The derivative of the cost function concerning each
$\theta_{kjj^{\prime}}\in\vec{\theta}$ following the chain rule (12) gives
(hereafter, we omit its indices)

$\displaystyle\frac{\partial
C(\vec{\theta})}{\partial\theta}=\frac{1}{M}\sum_{i=0}^{M-1}\frac{\partial
C(y_{i})}{\partial\widehat{y}_{i}(l)}\frac{\partial\widehat{y}_{i}(l)}{\partial
p}\frac{\partial p(\vec{\theta})}{\partial\theta},$ (26)

with

$\displaystyle\frac{\partial C(y_{i})}{\partial\widehat{y}_{i}(l)}$
$\displaystyle=-2[1-\widehat{y}_{i}(l)],$ (27)
$\displaystyle\frac{\partial\widehat{y}_{i}(l)}{\partial p}$
$\displaystyle=(-1)^{l}\widehat{y}_{i}(l)[1-\widehat{y}_{i}(l)],$ (28)

and

$\displaystyle\frac{\partial p(\vec{\theta})}{\partial\theta}$
$\displaystyle=\gamma\frac{\partial}{\partial\theta}\Bigl{(}p^{(i)}_{0}-p^{(i)}_{1}\Bigr{)}$
$\displaystyle=\gamma\Bigl{\\{}c_{0}\bigl{[}p^{(i)}_{0}(\theta+s_{0})-p^{(i)}_{0}(\theta-
s_{0})\bigr{]}$ $\displaystyle\hskip 14.22636pt-
c_{1}\bigl{[}p^{(i)}_{1}(\theta+s_{1})-p^{(i)}_{1}(\theta-
s_{1})\bigr{]}\Bigr{\\}},$ (29)

here, we have applied the parameter-shift rule (10). We finally update the
model with new parameters using the Adaptive Moment (Adam) optimization
algorithm kingma2017adam  as follows. Let $g_{\rm t+1}$ be the gradient of the
cost function w.r.t. $\theta$ at time step t+1, (t starts from 0):

$\displaystyle g_{\rm t+1}\equiv\frac{\partial C(\vec{\theta}_{\rm
t})}{\partial\theta},$ (30)

where $\vec{\theta}_{\rm t}$ is the set of parameters at time step t. Next, we
estimate the first and second order moments of the gradient as

$\displaystyle v_{\rm t+1}=\beta_{1}v_{t}+(1-\beta_{1})g_{t+1}$ (31)
$\displaystyle w_{\rm t+1}=\beta_{2}w_{t}+(1-\beta_{2})g^{2}_{t+1},$ (32)

and compute

$\displaystyle\widehat{v}_{\rm t+1}=\frac{v_{\rm
t+1}}{1-\beta_{1}^{t+1}},\text{and }\widehat{w}_{\rm t+1}=\frac{w_{\rm
t+1}}{1-\beta_{2}^{t+1}}.$ (33)

Here, $\beta_{1},\beta_{2}$ and $\epsilon$ are nonnegative weights, originally
suggested as $\beta_{1}=0.9,\beta_{2}=0.999$ and $\epsilon=10^{-8}$. Also,
$v_{0}=w_{0}=0$. Finally, the parameter $\theta$ is updated to

$\displaystyle\theta_{\rm t+1}=\theta_{\rm t}-\frac{\eta\widehat{v}_{\rm
t+1}}{\sqrt{\widehat{w}_{\rm t+1}}+\epsilon},$ (34)

where $\eta$ is the learning rate. The implementation of Adam algorthims in a
combination with QDP is given in Algorithm 5.

Result: Updated parameters

$c\leftarrow$ Quantum Circuit;

$p\leftarrow$ Params; $s\leftarrow$ Shift;

$\beta_{1},\beta_{2},\epsilon,\eta\leftarrow$ Beta1, Beta2, Epsilon, Eta;

$v$_adam, $s$_adam $\leftarrow$ $v$_adam, $s$_adam;

$t\leftarrow$ $t^{th}$ iteration;

diff $\leftarrow$ zero_matrix, size(diff) = size($p$);

for  _i $\leftarrow$ 0, size(param)_ do

for  _j $\leftarrow$ 0, size(param[$i$])_ do

for  _k $\leftarrow$ 0, size(param[$i][j$])_ do

$p$_plus $\leftarrow$ copy($p$);

$p$_subs $\leftarrow$ copy($p$);

$p$_plus$[i][j][k]$ $\leftarrow p[i][j][k]+s$;

$p$_subs$[i][j][k]$ $\leftarrow p[i][j][k]-s$;

diff$[i][j][k]\leftarrow$ QDP($p$_plus, $p$_subs, $c$, $s$);

end for

end for

end for

$p$, $v$_adam, $s$_adam $\leftarrow$ Adam(diff, $p$, $v$_adam, $s$_adam,
$\eta$, $\beta_{1}$, $\beta_{2}$, $\epsilon$, t) ;

Algorithm 5 Adam optimization implementation in Qsun

#### Preliminary QNN results

Figure 6: (a) Classical NN model; (b) and (c) QNN models using Pennylane and
Qsun, respectively; (d) and (e) are the zoom-in of (b) and (c), respectively.
Each plot shows from 1 layer (1L) to 5 layers (5Ls). See text for more details
about parameters used for each model.

We apply the QNN model on the Social Network Ads dataset SNA . This is a
categorical dataset to determine whether a user purchases a particular product
or not. We consider two features $Age$ and $EstimatedSalary$ to train the
model with the output $Purchased$. We use 400 samples and split them into 323
samples for the training set and 80 samples for the testing set. Once the data
is normalized, we encode $Age$ and $EstimatedSalary$ into two qubits
$|\psi_{j}\rangle,j={0,1}$, while $y$ = $Purchased$. We evaluate the predicted
value $\widehat{y}$ and minimize the loss function (25) to train the QNN
model. The full code is given in E.

For the comparison of performance, we train the QNN model using Qsun and
Pennylane. We note that while the QDP and Adam algorithm are implemented in
Qsun and Pennylane separately, the QNN circuit, encoding, and decoding
procedures are the same for both approaches. We also compare their results
with those from a classical model to explore advantages of QNN. For the
classical NN model, we use the MLPClassifier function from scikit-learn that
has 100 nodes per layer with the ReLU activation function. For the QNN model,
we have 4 nodes per layer with the sigmoid activation function. We fix
$s_{0}=s_{1}=\pi/20$ in the parameter-shift rule and $\eta=0.1$. For all
optimization, we use the Adam algorithm with $\beta_{1}=0.9,\beta_{2}=0.99$,
and $\epsilon=10^{-6}$.

Fig. 6 represents the loss function versus iteration for different layers as
shown in the figure for three cases: a classical NN model (a) and a QNN model
trained by Pennylane (b) and Qsun (c). In general, the loss function will
reduce and get unchanged after sufficient layers. In our example, the loss
function becomes saturated after 5 layers. As expected, the loss function of
QNN model is much smaller than the the classical one. Interestingly, when
comparing the two quantum approaches (see insets (d, e) for a zoom in), the
Qsun loss is converged faster than the Pennylane one.

We next focus on our QNN model with fixed 5 layers. In Fig. (7)(a), we show
the reduction of cost functions and activation functions during the iteration
for the training and testing processes. They both behave similarly and become
saturated after 100 iterations. Fig. (7)(b, c) represent results for the
training and testing sets. In the figure, both the true value $y$ and the
predicted value $\widehat{y}$ are labeled as ‘0’ and ‘1’ for $no-purchase$ and
$purchase$, respectively. $y\widehat{y}$ = ‘00’ or ‘11’ indicates a correct
prediction, whereas $y\widehat{y}$ = ‘01’ or ‘10’ indicates a wrong
prediction. It can be seen from Fig. 7 that our QNN model has a good
performance, reflected by large values of main diagonal elements
$y\widehat{y}$ = ‘00’ and ‘11’. Importantly, unlike the classical neural
network SNA , it does not require many nodes in each layer, which is one of
the main advantages of QNN.

Figure 7:  Results from our QNN model with five layers and two features on the
Social Network Ads dataset. (a) Plot of loss and accuracy versus iteration for
the training and testing process. (b, c) The confusion matrices for the
training and testing sets. Here we use 400 samples and split them into 320
samples for training and 80 samples for testing with 400 iterations.

##  V Conclusion

We have developed Qsun, an open-source platform of quantum virtual machine,
that simulates the operation of a real quantum computer using the wavefunction
approach. For small qubit numbers, the current version of Qsun runs standard
tasks significantly faster than other platforms do. The quantum differentiable
programming is implemented as a built-in function of Qsun, allowing us to
execute quantum machine learning applications effectively. We have employed
Qsun to implement two standard models in machine learning: linear regression
and neural network. For the former, Qsun yields a quantum regression model
nearly overlap with the classical reference, and somewhat better than the
Pennylane one. For the latter, the QNN model trained using Qsun shows a good
performance with a less number of nodes in each layer than the classical
neural network. Although Qsun is aimed to quantum machine learning problems,
as a generic quantum virtual machine, it can be used for multiple other
purposes, such as developing and testing variational quantum algorithms for
electronic structures or quantum information. It is well-fit to personal use
thanks to its light weight. Qsun is under active development to cover a wide
range of contents in machine learning. Our code is open source and available
on GitHub QSUN .

###### Acknowledgements.

This work was supported by JSPS KAKENHI Grant Number 20F20021.

## References

  * [1] F. Arute et al.  Quantum supremacy using a programmable superconducting processor.  Nature, 574(7779):505–510, Oct 2019.
  * [2] Zhong et al.  Quantum computational advantage using photons.  Science, 370(6523):1460–1463, 2020.
  * [3] John Preskill.  Quantum Computing in the NISQ era and beyond.  Quantum, 2:79, August 2018.
  * [4] N. Khammassi, I. Ashraf, X. Fu, C.G. Almudever, and K. Bertels.  Qx: A high-performance quantum computer simulation platform.  In Design, Automation Test in Europe Conference Exhibition (DATE), 2017, pages 464–469, 2017.
  * [5] Jianwei Wang, Stefano Paesani, Yunhong Ding, Raffaele Santagati, Paul Skrzypczyk, Alexia Salavrakos, Jordi Tura, Remigiusz Augusiak, Laura Mančinska, Davide Bacco, Damien Bonneau, Joshua W. Silverstone, Qihuang Gong, Antonio Acín, Karsten Rottwitt, Leif K. Oxenløwe, Jeremy L. O’Brien, Anthony Laing, and Mark G. Thompson.  Multidimensional quantum entanglement with large-scale integrated optics.  Science, 360(6386):285–291, 2018.
  * [6] Abhinav Kandala, Kristan Temme, Antonio D. Córcoles, Antonio Mezzacapo, Jerry M. Chow, and Jay M. Gambetta.  Error mitigation extends the computational reach of a noisy quantum processor.  Nature, 567(7749):491–495, Mar 2019.
  * [7] M. Cerezo, Andrew Arrasmith, Ryan Babbush, Simon C. Benjamin, Suguru Endo, Keisuke Fujii, Jarrod R. McClean, Kosuke Mitarai, Xiao Yuan, Lukasz Cincio, and Patrick J. Coles.  Variational quantum algorithms.  Nature Reviews Physics, Aug 2021.
  * [8] Michael Lubasch, Jaewoo Joo, Pierre Moinier, Martin Kiffner, and Dieter Jaksch.  Variational quantum algorithms for nonlinear problems.  Phys. Rev. A, 101:010301, Jan 2020.
  * [9] Jacob Biamonte.  Universal variational quantum computation.  Phys. Rev. A, 103:L030401, Mar 2021.
  * [10] Alberto Peruzzo, Jarrod McClean, Peter Shadbolt, Man-Hong Yung, Xiao-Qi Zhou, Peter J. Love, Alán Aspuru-Guzik, and Jeremy L. O’Brien.  A variational eigenvalue solver on a photonic quantum processor.  Nature Communications, 5(1):4213, Jul 2014.
  * [11] Ken M. Nakanishi, Kosuke Mitarai, and Keisuke Fujii.  Subspace-search variational quantum eigensolver for excited states.  Phys. Rev. Research, 1:033062, Oct 2019.
  * [12] William M. Kirby, Andrew Tranter, and Peter J. Love.  Contextual Subspace Variational Quantum Eigensolver.  Quantum, 5:456, May 2021.
  * [13] Ilya G. Ryabinkin, Scott N. Genin, and Artur F. Izmaylov.  Constrained variational quantum eigensolver: Quantum computer search engine in the fock space.  Journal of Chemical Theory and Computation, 15(1):249–255, Jan 2019\.
  * [14] Kazuhiro Seki, Tomonori Shirakawa, and Seiji Yunoki.  Symmetry-adapted variational quantum eigensolver.  Phys. Rev. A, 101:052340, May 2020.
  * [15] Bryan T. Gard, Linghua Zhu, George S. Barron, Nicholas J. Mayhall, Sophia E. Economou, and Edwin Barnes.  Efficient symmetry-preserving state preparation circuits for the variational quantum eigensolver algorithm.  npj Quantum Information, 6(1):10, Jan 2020.
  * [16] Nikolay V. Tkachenko, James Sud, Yu Zhang, Sergei Tretiak, Petr M. Anisimov, Andrew T. Arrasmith, Patrick J. Coles, Lukasz Cincio, and Pavel A. Dub.  Correlation-informed permutation of qubits for reducing ansatz depth in the variational quantum eigensolver.  PRX Quantum, 2:020337, Jun 2021.
  * [17] Ying Li and Simon C. Benjamin.  Efficient variational quantum simulator incorporating active error minimization.  Phys. Rev. X, 7:021050, Jun 2017.
  * [18] Sam McArdle, Tyson Jones, Suguru Endo, Ying Li, Simon C. Benjamin, and Xiao Yuan.  Variational ansatz-based quantum simulation of imaginary time evolution.  npj Quantum Information, 5(1):75, Sep 2019.
  * [19] Hirofumi Nishi, Taichi Kosugi, and Yu-ichiro Matsushita.  Implementation of quantum imaginary-time evolution method on nisq devices by introducing nonlocal approximation.  npj Quantum Information, 7(1):85, Jun 2021.
  * [20] Yong-Xin Yao, Niladri Gomes, Feng Zhang, Cai-Zhuang Wang, Kai-Ming Ho, Thomas Iadecola, and Peter P. Orth.  Adaptive variational quantum dynamics simulations.  PRX Quantum, 2:030307, Jul 2021.
  * [21] Zi-Jian Zhang, Jinzhao Sun, Xiao Yuan, and Man-Hong Yung.  Low-depth hamiltonian simulation by adaptive product formula, 2020.
  * [22] Chee Kong Lee, Pranay Patil, Shengyu Zhang, and Chang Yu Hsieh.  Neural-network variational quantum algorithm for simulating many-body dynamics.  Phys. Rev. Research, 3:023095, May 2021.
  * [23] Kenji Kubo, Yuya O. Nakagawa, Suguru Endo, and Shota Nagayama.  Variational quantum simulations of stochastic differential equations.  Phys. Rev. A, 103:052425, May 2021.
  * [24] Hsin-Yuan Huang, Kishor Bharti, and Patrick Rebentrost.  Near-term quantum algorithms for linear systems of equations, 2019.
  * [25] Yi ğit Subaş ı, Rolando D. Somma, and Davide Orsucci.  Quantum algorithms for systems of linear equations inspired by adiabatic quantum computing.  Phys. Rev. Lett., 122:060504, Feb 2019.
  * [26] Oleksandr Kyriienko, Annie E. Paine, and Vincent E. Elfving.  Solving nonlinear differential equations with differentiable quantum circuits.  Phys. Rev. A, 103:052416, May 2021.
  * [27] Carlos Bravo-Prieto, Ryan LaRose, M. Cerezo, Yigit Subasi, Lukasz Cincio, and Patrick J. Coles.  Variational quantum linear solver, 2020.
  * [28] Xiaosi Xu, Jinzhao Sun, Suguru Endo, Ying Li, Simon C. Benjamin, and Xiao Yuan.  Variational algorithms for linear algebra, 2019.
  * [29] Rolando D. Somma and Yi ğit Subaş ı.  Complexity of quantum state verification in the quantum linear systems problem.  PRX Quantum, 2:010315, Jan 2021.
  * [30] Seth Lloyd, Masoud Mohseni, and Patrick Rebentrost.  Quantum principal component analysis.  Nature Physics, 10(9):631–633, Sep 2014.
  * [31] Ryan LaRose, Arkin Tikku, Étude O’Neel-Judy, Lukasz Cincio, and Patrick J. Coles.  Variational quantum state diagonalization.  npj Quantum Information, 5(1):57, Jun 2019.
  * [32] K. Mitarai, M. Negoro, M. Kitagawa, and K. Fujii.  Quantum circuit learning.  Phys. Rev. A, 98:032309, Sep 2018.
  * [33] Edward Farhi and Hartmut Neven.  Classification with quantum neural networks on near term processors, 2018\.
  * [34] Maria Schuld and Nathan Killoran.  Quantum machine learning in feature hilbert spaces.  Phys. Rev. Lett., 122:040504, Feb 2019.
  * [35] Vojtěch Havlíček, Antonio D. Córcoles, Kristan Temme, Aram W. Harrow, Abhinav Kandala, Jerry M. Chow, and Jay M. Gambetta.  Supervised learning with quantum-enhanced feature spaces.  Nature, 567(7747):209–212, Mar 2019.
  * [36] Maria Schuld, Alex Bocharov, Krysta M. Svore, and Nathan Wiebe.  Circuit-centric quantum classifiers.  Phys. Rev. A, 101:032308, Mar 2020.
  * [37] Iris Cong, Soonwon Choi, and Mikhail D. Lukin.  Quantum convolutional neural networks.  Nature Physics, 15(12):1273–1278, Dec 2019.
  * [38] Arthur Pesah, M. Cerezo, Samson Wang, Tyler Volkoff, Andrew T. Sornborger, and Patrick J. Coles.  Absence of barren plateaus in quantum convolutional neural networks, 2020\.
  * [39] Kaining Zhang, Min-Hsiu Hsieh, Liu Liu, and Dacheng Tao.  Toward trainability of quantum neural networks, 2020.
  * [40] Kerstin Beer, Dmytro Bondarenko, Terry Farrelly, Tobias J. Osborne, Robert Salzmann, Daniel Scheiermann, and Ramona Wolf.  Training deep quantum neural networks.  Nature Communications, 11(1):808, Feb 2020.
  * [41] Andrew Arrasmith, Lukasz Cincio, Andrew T. Sornborger, Wojciech H. Zurek, and Patrick J. Coles.  Variational consistent histories as a hybrid algorithm for quantum foundations.  Nature Communications, 10(1):3438, Jul 2019.
  * [42] Mark Wilde.  Quantum information theory.  Cambridge University Press, 2017.
  * [43] Carlos Bravo-Prieto, Diego García-Martín, and José I. Latorre.  Quantum singular value decomposer.  Phys. Rev. A, 101:062310, Jun 2020.
  * [44] Bálint Koczor, Suguru Endo, Tyson Jones, Yuichiro Matsuzaki, and Simon C Benjamin.  Variational-state quantum metrology.  New Journal of Physics, 22(8):083038, aug 2020.
  * [45] Raphael Kaubruegger, Pietro Silvi, Christian Kokail, Rick van Bijnen, Ana Maria Rey, Jun Ye, Adam M. Kaufman, and Peter Zoller.  Variational spin-squeezing algorithms on programmable quantum sensors.  Phys. Rev. Lett., 123:260505, Dec 2019.
  * [46] Jacob L. Beckey, M. Cerezo, Akira Sone, and Patrick J. Coles.  Variational quantum algorithm for estimating the quantum fisher information, 2020.
  * [47] Johannes Jakob Meyer, Johannes Borregaard, and Jens Eisert.  A variational toolbox for quantum multi-parameter estimation.  npj Quantum Information, 7(1):89, Jun 2021.
  * [48] Nelder A. John and R. Mead.  A simplex method for function minimization.  Computer Journal, 7(4):308–313, 1965.
  * [49] Shaopeng Zhu, Shih-Han Hung, Shouvanik Chakrabarti, and Xiaodi Wu.  On the principles of differentiable quantum programming languages.  In Proceedings of the 41st ACM SIGPLAN Conference on Programming Language Design and Implementation, PLDI 2020, page 272–285, New York, NY, USA, 2020. Association for Computing Machinery.
  * [50] Frank Schäfer, Michal Kloc, Christoph Bruder, and Niels Lörch.  A differentiable programming method for quantum control.  Machine Learning: Science and Technology, 1(3):035009, aug 2020\.
  * [51] James Stokes, Josh Izaac, Nathan Killoran, and Giuseppe Carleo.  Quantum Natural Gradient.  Quantum, 4:269, May 2020.
  * [52] Anna Lopatnikova and Minh-Ngoc Tran.  Quantum natural gradient for variational bayes, 2021.
  * [53] Naoki Yamamoto.  On the natural gradient for variational quantum eigensolver, 2019.
  * [54] Aram W. Harrow and John C. Napp.  Low-depth gradient measurements can improve convergence in variational hybrid quantum-classical algorithms.  Phys. Rev. Lett., 126:140502, Apr 2021.
  * [55] Gian Giacomo Guerreschi and Mikhail Smelyanskiy.  Practical optimization for hybrid quantum-classical algorithms, 2017.
  * [56] Ville Bergholm, Josh Izaac, Maria Schuld, Christian Gogolin, M. Sohaib Alam, Shahnawaz Ahmed, Juan Miguel Arrazola, Carsten Blank, Alain Delgado, Soran Jahangiri, Keri McKiernan, Johannes Jakob Meyer, Zeyue Niu, Antal Száva, and Nathan Killoran.  Pennylane: Automatic differentiation of hybrid quantum-classical computations, 2020.
  * [57] Jun Li, Xiaodong Yang, Xinhua Peng, and Chang-Pu Sun.  Hybrid quantum-classical approach to quantum optimal control.  Phys. Rev. Lett., 118:150503, Apr 2017.
  * [58] Maria Schuld, Ville Bergholm, Christian Gogolin, Josh Izaac, and Nathan Killoran.  Evaluating analytic gradients on quantum hardware.  Phys. Rev. A, 99:032331, Mar 2019.
  * [59] Leonardo Banchi and Gavin E. Crooks.  Measuring Analytic Gradients of General Quantum Evolution with the Stochastic Parameter Shift Rule.  Quantum, 5:386, January 2021.
  * [60] Qiskit: Open-source quantum development.
  * [61] Cirq Developers.  Cirq, August 2021.  See full list of authors on Github: https://github .com/quantumlib/Cirq/graphs/contributors.
  * [62] Robert S. Smith, Michael J. Curtis, and William J. Zeng.  A practical quantum instruction set architecture, 2017.
  * [63] Krysta Svore, Alan Geller, Matthias Troyer, John Azariah, Christopher Granade, Bettina Heim, Vadym Kliuchnikov, Mariia Mykhailova, Andres Paz, and Martin Roetteler.  Q#: Enabling scalable quantum computing and development with a high-level dsl.  In Proceedings of the Real World Domain Specific Languages Workshop 2018, RWDSL2018, New York, NY, USA, 2018. Association for Computing Machinery.
  * [64] Benjamin Villalonga, Sergio Boixo, Bron Nelson, Christopher Henze, Eleanor Rieffel, Rupak Biswas, and Salvatore Mandrà.  A flexible high-performance simulator for verifying and benchmarking quantum circuits implemented on real hardware.  npj Quantum Information, 5(1):86, Oct 2019.
  * [65] Chase Roberts, Ashley Milsted, Martin Ganahl, Adam Zalcman, Bruce Fontaine, Yijian Zou, Jack Hidary, Guifre Vidal, and Stefan Leichenauer.  Tensornetwork: A library for physics and machine learning, 2019.
  * [66] Yasunari Suzuki, Yoshiaki Kawase, Yuya Masumura, Yuria Hiraga, Masahiro Nakadai, Jiabao Chen, Ken M. Nakanishi, Kosuke Mitarai, Ryosuke Imai, Shiro Tamiya, Takahiro Yamamoto, Tennin Yan, Toru Kawakubo, Yuya O. Nakagawa, Yohei Ibe, Youyuan Zhang, Hirotsugu Yamashita, Hikaru Yoshimura, Akihiro Hayashi, and Keisuke Fujii.  Qulacs: a fast and versatile quantum circuit simulator for research purpose, 2020.
  * [67] Gian Giacomo Guerreschi, Justin Hogaboam, Fabio Baruffa, and Nicolas P D Sawaya.  Intel quantum simulator: a cloud-ready high-performance simulator of quantum circuits.  Quantum Science and Technology, 5(3):034007, may 2020.
  * [68] Damian S. Steiger, Thomas Häner, and Matthias Troyer.  ProjectQ: an open source software framework for quantum computing.  Quantum, 2:49, January 2018.
  * [69] Adam Kelly.  Simulating quantum computers using opencl, 2018.
  * [70] Tyson Jones, Anna Brown, Ian Bush, and Simon C. Benjamin.  Quest and high performance simulation of quantum computers.  Scientific Reports, 9(1):10736, Jul 2019.
  * [71] K. De Raedt, K. Michielsen, H. De Raedt, B. Trieu, G. Arnold, M. Richter, Th. Lippert, H. Watanabe, and N. Ito.  Massively parallel quantum computer simulator.  Computer Physics Communications, 176(2):121–136, 2007.
  * [72] Hans De Raedt, Fengping Jin, Dennis Willsch, Madita Willsch, Naoki Yoshioka, Nobuyasu Ito, Shengjun Yuan, and Kristel Michielsen.  Massively parallel quantum computer simulator, eleven years later.  Computer Physics Communications, 237:47–61, 2019.
  * [73] Igor L. Markov and Yaoyun Shi.  Simulating quantum computation by contracting tensor networks.  SIAM Journal on Computing, 38(3):963–981, 2008.
  * [74] Feng Pan, Pengfei Zhou, Sujie Li, and Pan Zhang.  Contracting arbitrary tensor networks: General approximate algorithm and applications in graphical models and quantum circuit simulations.  Phys. Rev. Lett., 125:060503, Aug 2020.
  * [75] Igor L. Markov, Aneeqa Fatima, Sergei V. Isakov, and Sergio Boixo.  Quantum supremacy is both closer and farther than it appears, 2018.
  * [76] Michael Broughton, Guillaume Verdon, Trevor McCourt, Antonio J. Martinez, Jae Hyeon Yoo, Sergei V. Isakov, Philip Massey, Murphy Yuezhen Niu, Ramin Halavati, Evan Peters, Martin Leib, Andrea Skolik, Michael Streif, David Von Dollen, Jarrod R. McClean, Sergio Boixo, Dave Bacon, Alan K. Ho, Hartmut Neven, and Masoud Mohseni.  Tensorflow quantum: A software framework for quantum machine learning, 2020.
  * [77] William Huggins, Piyush Patil, Bradley Mitchell, K Birgitta Whaley, and E Miles Stoudenmire.  Towards quantum machine learning with tensor networks.  Quantum Science and Technology, 4(2):024001, jan 2019.
  * [78] Chuong Nguyen Quoc, Le Bin Ho, Lan Nguyen Tran, and Hung Q. Nguyen.  Qsun, a platform towards quantum machine learning, 2021.  https://github.com/ChuongQuoc1413017/Quantum_Virtual_Machine.
  * [79] Michael A. Nielsen and Isaac L. Chuang.  Quantum computation and quantum information.  Cambridge University Press, 2019.
  * [80] Ryan LaRose.  Overview and Comparison of Gate Level Quantum Software Platforms.  March 2019.  arXiv:1807.02500 [quant-ph].
  * [81] Willard Miller.  Symmetry groups and their applications.  Academic Press, 1990.
  * [82] Gavin E. Crooks.  Gradients of parameterized quantum gates using the parameter-shift rule and gate decomposition.  May 2019.  arXiv:1905.13311 [quant-ph].
  * [83] Bradley Efron, Trevor Hastie, Iain Johnstone, and Robert Tibshirani.  Least angle regression.  Annals of Statistics, page 407, 2004.
  * [84] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay.  Scikit-learn: Machine learning in Python.  Journal of Machine Learning Research, 12:2825–2830, 2011.
  * [85] Ryan LaRose and Brian Coyle.  Robust data encodings for quantum classifiers.  Phys. Rev. A, 102:032420, Sep 2020.
  * [86] Diederik P. Kingma and Jimmy Ba.  Adam: A method for stochastic optimization, 2017.
  * [87] Social network ads knn, lr, neural network.  https://www.kaggle.com/sunilsharanappa/social-network-ads-knn-lr-neural-network.

##  Appendix A Code for benchmarking

In this section, we present the code used to benchmark the performance of Qsun
and compared it to other packages such as Qiskit, ProjectQ, and Pennylane. The
criteria are the time to execute a certain quantum circuit of constant depth.
The circuit contain an H, a $\sqrt{X}$, and a CNOT with $depth=10$ and
measurements on all qubits. With the number of qubits as a variable, the time
needed to execute the circuit is estimated, yielding the package’s
performance. The number of qubits varies from 2 to 15 qubits, which means we
have to implement 14 circuits for each simulation. The implementation’s time
is stored in a Pandas’s Dataframe and is shown in Fig. 2 in the main text.

⬇

1#import module for visualization

2import pandas as pd

3import seaborn as sns

4import time

5import matplotlib.pyplot as plt

6

7#import Qiskit module

8from qiskit import IBMQ, BasicAer

9from qiskit.providers.ibmq import least_busy

10from qiskit import QuantumCircuit, ClassicalRegister, QuantumRegister,
execute

11

12#import ProjectQ module

13from projectq import MainEngine

14import projectq.ops as ops

15from projectq.backends import Simulator

16

17#import Qsun module

18from Qsun.Qcircuit import *

19from Qsun.Qgates import *

20

21#import Pennylane module

22import pennylane as qml

23from pennylane import numpy as np

24

25# Qsun circuit

26def qvm_circuit(n_qubit, depth):

27 circuit = Qubit(n_qubit)

28 for m in range(depth):

29 for i in range(n_qubit):

30 H(circuit, i)

31 Xsquare(circuit, i)

32 for i in range(1, n_qubit):

33 CNOT(circuit, i, 0)

34 return circuit.probabilities()

35

36# Pennylane circuit

37dev = qml.device(’default.qubit’, wires=qubit_max)

38

39@qml.qnode(dev)

40def qml_circuit(n_qubit, depth):

41 for m in range(depth):

42 for i in range(n_qubit):

43 qml.Hadamard(wires=i)

44 qml.SX(wires=i)

45 for i in range(1, n_qubit):

46 qml.CNOT(wires=[i, 0])

47 return qml.probs(wires=range(n_qubit))

48

49# Qiskit circuit

50def qiskit_circuit(n_qubit, depth):

51 qr = QuantumRegister(n_qubit)

52 cr = ClassicalRegister(n_qubit)

53 circuit = QuantumCircuit(qr, cr)

54 for m in range(depth):

55 for i in range(n_qubit):

56 circuit.h(qr[i])

57 circuit.rx(math.pi/2, qr[i])

58 for i in range(1, n_qubit):

59 circuit.cx(qr[i], qr[0])

60 circuit.measure_all()

61 backend = BasicAer.get_backend(’qasm_simulator’)

62 result = execute(circuit, backend=backend, shots=1024).result()

63 return result.get_counts(circuit)

64

65# ProjectQ circuit

66def projectq_circuit(n_qubit, depth):

67 eng = MainEngine(backend=Simulator(gate_fusion=True),engine_list=[])

68 qbits = eng.allocate_qureg(n_qubit)

69 for level in range(depth):

70 for q in qbits:

71 ops.H|q

72 ops.SqrtX|q

73 if q != qbits[0]:

74 ops.CNOT|(q,qbits[0])

75 for q in qbits:

76 ops.Measure|q

77 eng.flush()

78

79#run test for Qsun

80data_qvm = []

81

82for n_qubit in range(qubit_min, qubit_max):

83 start = time.time()

84 qvm_circuit(n_qubit, depth)

85 end = time.time()

86 data_qvm.append(end-start)

87

88#define the test parametes

89qubit_min = 2

90qubit_max = 10

91depth = 10

92

93#run test for Pennylane

94data_qml = []

95

96for n_qubit in range(qubit_min, qubit_max):

97 start = time.time()

98 qml_circuit(n_qubit, depth)

99 end = time.time()

100 data_qml.append(end-start)

101

102#run test for Qiskit

103data_qiskit = []

104

105for n_qubit in range(qubit_min, qubit_max):

106 start = time.time()

107 qiskit_circuit(n_qubit, depth)

108 end = time.time()

109 data_qiskit.append(end-start)

110

111#run test for ProjectQ

112data_projectq = []

113

114for n_qubit in range(qubit_min, qubit_max):

115 start = time.time()

116 projectq_circuit(n_qubit, depth)

117 end = time.time()

118 data_projectq.append(end-start)

119

120df = pd.DataFrame({’QVM’: data_qvm, ’Pennylane’: data_qml,

121 ’Qiskit’: data_qiskit, ’ProjectQ’: data_projectq}, index =
range(qubit_min, qubit_max))

122

123sns.lineplot(data = df)

124plt.xlabel(’Number of Qubits with depth = ’+str(depth))

125plt.ylabel(’Time(s)’)

126plt.grid()

127plt.savefig(’compare.png’)

Code example 1: Implementation of the benchmarking test on four typical
languages: Qsun, Qiskit, ProjectQ, and Pennylane.

##  Appendix B Code for QDP benchmarking

For this code, we run the QDP algorithm on Qsun and measure how long it will
take to update parameters. We implement the circuit described in Fig. 3 and
measure the expected values of output, then record the time. The number of
qubits varies from 1 to 10, and the maximum number of iterations is 1000.
Parameters for QDP, in this case, are $\eta=0.1$ and $s=\pi/20$.

⬇

1from Qsun.Qcircuit import *

2from Qsun.Qgates import *

3import numpy as np

4import time

5

6def circuit(params, n):

7 c = Qubit(n)

8 for j in range(0, 2*n, 2):

9 RX(c, j//2, params[j])

10 RY(c, j//2, params[j+1])

11 return c

12

13def output(params):

14 c = circuit(params, len(params)//2)

15 prob = c.probabilities()

16 return -np.sum([i*prob[i] for i in range(len(prob))])

17

18def cost(params):

19 expval = output(params)

20 return expval

21

22def grad(params, shift, eta):

23 for i in range(len(params)):

24 params_1 = params.copy()

25 params_2 = params.copy()

26 params_1[i] += shift

27 params_2[i] -= shift

28 diff[i] = (cost(params_1)-cost(params_2))/(2*np.sin(shift))

29 for i in range(len(params)):

30 params[i] = params[i] - eta*diff[i]

31 return params

32

33time_qvm = []

34n = 10

35for i in range(1, n+1):

36

37 start = time.time()

38 params = np.random.normal(size=(2*i,))

39 diff = np.random.normal(size=(2*i,))

40

41 for i in range(1000):

42 params = grad(params, np.pi/20, eta=0.01)

43 end = time.time()

44

45 time_qvm.append(end-start)

46 print(cost(params))

Code example 2: Implementation of the quantum differentiable programming test
in Qsun.

##  Appendix C QDP Tutorial Appendix

This appendix demonstrates how to run a gradient descent algorithm using QDP
in Qsun. We use $\eta=0.1$ and $s=\pi/20$ to find the derivative of a circuit.
By doing that, we maximize the expectation values of one qubit. So the
objective function we want to maximize is $|\langle 1|\psi\rangle|^{2}$, where
$\ket{\psi}$ is a quantum state of that qubit.

⬇

1import numpy as np

2from Qsun.Qcircuit import *

3from Qsun.Qgates import *

4

5def circuit(params):

6 c = Qubit(1)

7 RX(c, 0, params[0])

8 RY(c, 0, params[1])

9 return c

10

11def output(params):

12 c = circuit(params)

13 prob = c.probabilities()

14 return 0.*prob[0] + 1*prob[1]

15

16def cost(params):

17 c = circuit(params)

18 prob = c.probabilities()

19 expval = output(params)

20 return np.abs(expval - 1)**2

21

22def grad(params, shift, eta):

23 for i in range(len(params)):

24 params_1 = params.copy()

25 params_2 = params.copy()

26 params_1[i] += shift

27 params_2[i] -= shift

28 diff[i] = (cost(params_1)-cost(params_2))/(2*np.sin(shift))

29 for i in range(len(params)):

30 params[i] = params[i] - eta*diff[i]

31 return params

32

33params = np.random.normal(size=(2,))

34diff = np.random.normal(size=(2,))

35

36for i in range(1000):

37 params = grad(params, np.pi/20, eta=0.01)

38

39print(’Circuit output:’, output(params))

40print(’Final parameter:’, params)

41

42>>> Circuit output: -0.9381264201123851

43>>> Final parameter: [ 0.29017649 -2.93657549]

Code example 3: Implementation of the Gradient Descent algorithm by using QDP
in Qsun.

##  Appendix D Quantum Linear Regression Appendix

Here, a Quantum Linear Regression model is programmed in Qsun, with results
shown in Fig. 4. Its parameters are $k=10$, $\eta=0.1$, and $s=\pi/20$. The
optimization algorithm used in this code is Gradient Descent.

⬇

1from Qsun.Qcircuit import *

2from Qsun.Qgates import *

3from sklearn import datasets

4

5def circuit(params):

6 c = Qubit(1)

7 RX(c, 0, params[0])

8 RY(c, 0, params[1])

9 return c

10

11def output(params):

12 c = circuit(params)

13 prob = c.probabilities()

14 return 1*prob[0] - 1*prob[1]

15

16def predict(x_true, coef_params, intercept_params, boundary=10):

17 coef = boundary*output(coef_params)

18 intercept = boundary*output(intercept_params)

19 return coef*x_true.flatten() + intercept

20

21def errors(x_true, y_true, coef_params, intercept_params, boundary):

22 return mean_squared_error(y_true, predict(x_true, coef_params,
intercept_params, boundary))

23

24def grad(x_true, y_true, coef_params, intercept_params, shift, eta,
boundary):

25

26 coef_diff = np.zeros((2,))

27 intercept_diff = np.zeros((2,))

28

29 for i in range(len(coef_params)):

30 coef_params_1 = coef_params.copy()

31 coef_params_2 = coef_params.copy()

32 coef_params_1[i] += shift

33 coef_params_2[i] -= shift

34 for x, y in zip(x_true, y_true):

35 coef_diff[i] -= x*(y-predict(x, coef_params, intercept_params,
boundary))*(output(coef_params_1)-output(coef_params_2))/(2*np.sin(shift))

36

37

38 for i in range(len(coef_params)):

39 intercept_params_1 = intercept_params.copy()

40 intercept_params_2 = intercept_params.copy()

41 intercept_params_1[i] += shift

42 intercept_params_2[i] -= shift

43 for x, y in zip(x_true, y_true):

44 intercept_diff[i] -= (y-predict(x, coef_params, intercept_params,
boundary))*(output(intercept_params_1)-output(intercept_params_2))/(2*np.sin(shift))

45

46 coef_diff = coef_diff*boundary*2/len(y_true)

47 intercept_diff = intercept_diff*boundary*2/len(y_true)

48

49 for i in range(len(coef_params)):

50 coef_params[i] = coef_params[i] - eta*coef_diff[i]

51

52 for i in range(len(intercept_params)):

53 intercept_params[i] = intercept_params[i] - eta*intercept_diff[i]

54

55 return coef_params, intercept_params

56

57X, y = datasets.load_diabetes(return_X_y=True)

58y = (y - np.min(y)) / (np.max(y) - np.min(y))

59# Use only one feature

60X = X[:, np.newaxis, 2]

61

62# Split the data into training/testing sets

63X_train = X[:400]

64X_test = X[-10:]

65

66# Split the targets into training/testing sets

67y_train = y[:400]

68y_test = y[-10:]

69

70coef_params = np.random.normal(size=(2,))

71intercept_params = np.random.normal(size=(2,))

72

73start = time.time()

74for i in range(1000):

75 coef_params, intercept_params = grad(X_train, y_train, coef_params,
intercept_params, np.pi/20, eta=0.01, boundary=10)

76end = time.time()

77

78y_pred = predict(X_test, coef_params, intercept_params, boundary=10)

Code example 4: Implementation of the linear regression model in qSUN.

##  Appendix E Quantum Neural Network Appendix

The code showing in this appendix implements a process from training a QNN
model with 5 layers and 2 qubits on the "Social_Network_Ads" dataset. From
this dataset, we use a quantum neural network to predict whether a customer
buys a car or not. The function train_test_split() splits the original dataset
into the training and testing datasets. MinMaxScaler() scales the training
dataset to [0, 1] before feeding it into the QNN circuit. The optimization
algorithm in this model is the Adam optimization algorithm, as presented in
the main text. The parameters used in the model are: $\rm test\\_size=0.2;\rm
random\\_state=0;s=\pi/20;\eta=0.1;\beta_{1}=0.9,\beta_{2}=0.99;\epsilon=10^{-6};\rm
number\\_of\\_iterations=150$.

⬇

1# import libraries

2from Qsun.Qcircuit import *

3from Qsun.Qgates import *

4from Qsun.Qmeas import *

5from sklearn.model_selection import train_test_split

6from sklearn.preprocessing import MinMaxScaler

7from sklearn.metrics import confusion_matrix

8import pandas as pd

9import seaborn as sn

10import matplotlib.pyplot as plt

11import time

12

13# one layer with full entanglement

14def layer(circuit, params):

15 circuit_layer = circuit

16 n_qubit = len(params)

17 for i in range(n_qubit):

18 RX(circuit_layer, i, params[i][0])

19 RY(circuit_layer, i, params[i][1])

20 for i in range(n_qubit-1):

21 CNOT(circuit_layer, i, i+1)

22 CNOT(circuit_layer, n_qubit-1, 0)

23 return circuit_layer

24

25# encoding the features

26def initial_state(sample):

27 circuit_initial = Qubit(len(sample))

28 ampli_vec = np.array([np.sqrt(sample[0]), np.sqrt(1-sample[0])])

29 for i in range(1, len(sample)):

30 ampli_vec = np.kron(ampli_vec, np.array([np.sqrt(sample[i]),
np.sqrt(1-sample[i])]))

31 circuit_initial.amplitude = ampli_vec

32 return circuit_initial

33

34# QNN circuit

35def qnn(circuit, params):

36 n_layer = len(params)

37 circuit_qnn = circuit

38 for i in range(n_layer):

39 circuit_qnn = layer(circuit_qnn, params[i])

40 return circuit_qnn

41

42# QNN model

43def qnn_model(sample, params):

44 circuit_model = initial_state(sample)

45 circuit_model = qnn(circuit_model, params)

46 return circuit_model

47

48# activation function

49def sigmoid(x):

50 return 1 / (1 + math.exp(-10*x))

51

52# make a prediction

53def predict(circuit):

54 prob_0 = measure_one(circuit, 0)

55 prob_1 = measure_one(circuit, 1)

56 expval_0 = prob_0[1]

57 expval_1 = prob_1[1]

58 pred = sigmoid(expval_0-expval_1)

59 return [pred, 1-pred]

60

61# make a prediction for exp

62def predict_ex(circuit):

63 prob_0 = measure_one(circuit, 0)

64 prob_1 = measure_one(circuit, 1)

65 expval_0 = prob_0[1]

66 expval_1 = prob_1[1]

67 return [expval_0, expval_1]

68

69# loss function

70def square_loss(labels, predictions):

71 loss = 0

72 for l, p in zip(labels, predictions):

73 loss = loss + (1 - p[l]) ** 2

74 loss = loss / len(labels)

75 return loss

76

77# loss function of QNN

78def cost(params, features, labels):

79 preds = [predict(qnn_model(x, params)) for x in features]

80 return square_loss(labels, preds)

81

82# https://d2l.ai/chapter_optimization/adam.html?highlight=adam

83def adam(X_true, y_true, params, v, s, shift, eta, drop_rate, beta1, beta2,
eps, iter_num):

84 diff = np.zeros(params.shape)

85 for i in range(len(params)):

86 for j in range(len(params[i])):

87 for k in range(len(params[i][j])):

88 dropout = np.random.choice([0, 1], 1, p = [1 - drop_rate, drop_rate])[0]

89 if dropout == 0:

90 params_1 = params.copy()

91 params_2 = params.copy()

92 params_1[i][j][k] += shift

93 params_2[i][j][k] -= shift

94 for x, y in zip(X_true, y_true):

95 circuit = qnn_model(x, params)

96 circuit_1 = qnn_model(x, params_1)

97 circuit_2 = qnn_model(x, params_2)

98 ex_plus = predict_ex(circuit_1)

99 ex_subs = predict_ex(circuit_2)

100 pred = predict(circuit)

101 diff_loss = ((-1)**y)*(-2)*(1-pred[y])*pred[y]*(1-pred[y])

102 diff_ex = 10*((ex_plus[0] - ex_subs[0]) - (ex_plus[1] -
ex_subs[1]))/(2*np.sin(shift))

103 diff[i][j][k] += diff_loss*diff_ex

104

105 diff /= len(y_true)

106 v = beta1 * v + (1 - beta1) * diff

107 s = beta2 * s + (1 - beta2) * (diff**2)

108 v_bias_corr = v / (1 - beta1 ** (iter_num+1))

109 s_bias_corr = s / (1 - beta2 ** (iter_num+1))

110 params -= eta * v_bias_corr / (np.sqrt(s_bias_corr) + eps)

111

112 return params, v, s

113

114# source: https://www.kaggle.com/rakeshrau/social-network-ads

115dataset = pd.read_csv(’Social_Network_Ads.csv’)

116X = dataset.iloc[:, 2:-1].values

117y = dataset.iloc[:, -1].values

118

119# splitting dataset

120X_train, X_test, y_train, y_test = train_test_split(X, y, test_size = 0.2,
random_state = 0)

121

122# scaling feature

123scaler = MinMaxScaler()

124X_train = scaler.fit_transform(X_train)

125X_test = scaler.transform(X_test)

126

127# create parameters

128n_layer = 5

129params = np.random.normal(size=(n_layer, len(X_train[0]), 2,))

130v = np.zeros(params.shape)

131s = np.zeros(params.shape)

132

133# training model

134start = time.time()

135for i in range(150):

136 params, v, s = adam(X_train, y_train, params, v, s,

137 shift=np.pi/20, eta=0.1, drop_rate=0.0, beta1=0.9, beta2=0.999, eps=1e-6,
iter_num=i)

138

139# confusion matrix

140label = y_test

141pred = [predict(qnn_model(x, params)) for x in X_test]

142pred = np.argmax(pred, axis=1)

143con = confusion_matrix(label,pred)

144sn.heatmap(con, annot=True, cmap="Blues")

Code example 5: Implementation of the QNN with two layers and two features in
Qsun.