††thanks: Corresponding author<EMAIL_ADDRESS>This work was
done when he was a research intern at JD Explore Academy.

# Accelerating variational quantum algorithms with multiple quantum processors

Yuxuan Du  JD Explore Academy  Yang Qian  School of Computer Science, The
University of Sydney  JD Explore Academy  Dacheng Tao  JD Explore Academy

###### Abstract

Variational quantum algorithms (VQAs) have the potential of utilizing near-
term quantum machines to gain certain computational advantages over classical
methods. Nevertheless, modern VQAs suffer from cumbersome computational
overhead, hampered by the tradition of employing a solitary quantum processor
to handle large-volume data. As such, to better exert the superiority of VQAs,
it is of great significance to improve their runtime efficiency. Here we
devise an efficient distributed optimization scheme, called QUDIO, to address
this issue. Specifically, in QUDIO, a classical central server partitions the
learning problem into multiple subproblems and allocate them to multiple local
nodes where each of them consists of a quantum processor and a classical
optimizer. During the training procedure, all local nodes proceed parallel
optimization and the classical server synchronizes optimization information
among local nodes timely. In doing so, we prove a sublinear convergence rate
of QUDIO in terms of the number of global iteration under the ideal scenario,
while the system imperfection may incur divergent optimization. Numerical
results on standard benchmarks demonstrate that QUDIO can surprisingly achieve
a superlinear runtime speedup with respect to the number of local nodes. Our
proposal can be readily mixed with other advanced VQAs-based techniques to
narrow the gap between the state of the art and applications with quantum
advantage.

##  I Introduction

Deep learning techniques have penetrated the world during past decades
goodfellow2016deep . Prototypical applications comprise using deep neural
networks (DNNs) to facilitate online shopping zhang2019deep , to accelerate
molecule design senior2020improved , and to enhance language translation
devlin2019bert . Notably, the success of deep learning heavily relies on
distributed hardware and distributed optimization techniques Jeffrey2012  in
the sense that multiple GPU cards are employed to collaboratively process the
same learning task. For instance, to ensure that the training of deep
bidirectional transformers (BERT) can be completed within a reasonable time
(e.g., $4$ days), in total $64$ GPU cards are used to conduct the distributed
optimization devlin2019bert . However, the price to pay is ten thousands
dollars and emitting $1438$ lbs of carbon dioxide strubell2019energy . This
considerable resource-consumption signifies that deep learning models will
become hard to develop and optimize when the problem size is continuously
enlarged. Therefore, it is highly demanded to seek novel techniques to resolve
complicated problems in a fast, economic, and environmentally friendly way.

The oath of quantum computing is to accomplish certain tasks beyond the reach
of classical computers biamonte2017quantum ; feynman1982simulating ;
harrow2017quantum . During past years, big breakthrough has been achieved
towards this goal, e.g., a demonstration of the quantum supremacy experiment
in the task of sampling the output of a pseudo-random quantum circuit
arute2019quantum . Among various quantum computational models, variational
quantum algorithms (VQAs) benedetti2019parameterized ; bharti2021noisy ;
cerezo2020variational2 ; du2018expressive ; endo2021hybrid , which are formed
by parameterized quantum circuits (PQCs) and a classical optimizer as shown in
the left panel of Fig. 1, have attracted great attention from industry and
academia. The popularity of VQAs origins from the versatility of PQCs, which
guarantees their efficient implementations on the noisy intermediate-scale
quantum (NISQ) machines preskill2018quantum , as well as theoretical evidence
of quantum superiority abbas2020power ; banchi2021generalization ;
bu2021effects ; caro2021encoding ; du2021efficient ; huang2021information ;
huang2021power ; wang2021towards ; wu2021expressivity . Moreover, prior
studies have exhibited the progress of VQAs of accomplishing diverse learning
tasks, e.g., machine learning issues such as data classification
Du_2021_grover ; havlivcek2019supervised ; mitarai2018quantum ;
schuld2019quantum ; wang2020quantum  and image generation
huang2020experimental ; rudolph2020generation ; zhu2019training ,
combinatorial optimization crooks2018performance ; farhi2014quantum ;
hadfield2019quantum ; harrigan2021quantum , finance alcazar2020classical ;
coyle2021quantum ; hodson2019portfolio , quantum information processing
beckey2020variational ; cerezo2020variational ; larose2019variational ,
quantum chemistry and material sciences bauer2016hybrid ; google2020hartree ;
o2019calculating ; peruzzo2014variational , and particle physics
avkhadiev2020accelerating ; kokail2019self .

Despite the tantalizing achievements, modern VQAs generally are plagued by the
expensive or even unaffordable execution time for large-volume data, hampered
by the regulation such that only a single quantum chip is employed in
optimization. For concreteness, let us recall the machinery of quantum neural
networks (QNNs), as a crucial subclass of VQAs, when dealing with
classification tasks beer2020training ; du2020learnability ;
farhi2018classification ; havlivcek2019supervised ; mitarai2018quantum ;
schuld2019quantum . To correctly classify $n$ training examples, the classical
optimizer iteratively feeds each example to PQCs and then leverages the
discrepancy between the obtained $n$ predictions of QNN and $n$ labels to
conduct the gradient-based optimization. In this way, the classical optimizer
needs to query the single quantum processor at least $O(n)$ times to complete
one iteration. Such overhead prohibits the applicability of QNNs for large $n$
qian2021dilemma . With this regard, it is of great importance to enhance the
computational efficiency of VQAs, as the necessary condition to pursue quantum
advantages.

Figure 1: The scheme of VQAs and quantum distributed optimization scheme
(QUDIO). The left panel illustrates the workflow of conventional VQAs. The
right panel presents the workflow of QUDIO, which consists of multiple local
nodes and a central server. Each local node only manipulates a subgroup of the
given problem in parallel, while the central server communicates with all
local nodes and synchronizes their results. In this way, QUDIO can accelerate
various VQAs.

In light of the key role of distributed techniques in deep learning and the
rapid growth of the amount of available quantum chips, a natural way to
accelerate the training of VQAs is involving multiple quantum processors to
fulfill the joint optimization. Noticeably, designing such a scheme is
extremely challenged and sharply inconsistent with classical distributed
optimization methods boyd2011distributed ; Jeffrey2012  or quantum distributed
computation beals2013efficient , because of the following three aspects.
First, the gradients information operated in VQAs is biased, induced by the
system imperfection such as sample error and gate noise, whereas most
classical distributed methods assume the unbiased gradients information. Next,
unlike distributed methods in DNNs, VQAs are immune to the communication
bottleneck, due to the small number of trainable parameters. Last, quantum
distributed computation beals2013efficient  focuses on using multiple less
powerful quantum circuits to simulate a standard quantum circuit, which
differs from our aim. Although some quantum software bergholm2018pennylane ;
Qiskit  showcases parallelling variational quantum eigen-solvers (VQEs)
peruzzo2014variational , how to devise general distributed-VQAs optimization
schemes with both runtime boost and convergence guarantee remains largely
unknown.

To conquer the above issues, here we devise an efficient QUantum DIstributed
Optimization scheme (abbreviated as QUDIO). An attractive property of QUDIO is
adequately utilizing the accessible quantum resources to accelerate VQAs,
owing to its compatibility. Namely, the deployed quantum processors are
allowed to be any type of quantum hardware such as linear optical, ion-trap,
and superconducting quantum chips. Such a compatibility contributes to apply a
wide class of VQAs to manipulate varied large-scale computational problems and
seek potential quantum advantages by unifying quantum powers in a maximum
extent.

Our second contribution is analyzing convergence of QUDIO under both the ideal
and NISQ scenarios. Particularly, we prove an asymptotic convergence between
QUDIO and conventional VQAs under the ideal setting. By contrast, in the NISQ
case, the convergence rate of QUDIO becomes degraded with respect to the
amplified system noise and the decreased number of quantum measurements. These
results can not only be employed as theoretical guidance to assure good
performance of QUDIO, but also motive us to devise more advanced distributed-
VQAs schemes. To our best knowledge, this is the first proposal towards
distributed-VQAs techniques with theoretical guarantee.

Last, we conduct extensive numerical simulations to validate the computational
efficiency of our scheme. In particular, QUDIO is exploited to accomplish
image classification and the ground energy of hydrogen molecule estimation
tasks under both the ideal and noisy scenarios. The achieved simulation
results validate that QUDIO realizes sublinear and even superlinear speedups
compared with conventional QNNs and VQEs.

The remainder of this study is organized as follows. In Section II, we present
the workflow of QUDIO. Subsequently, in Section III, we exhibit how to use
QUDIO to accelerate QNNs with both convergence analysis and numerical
simulations. Next, we explain the usage of QUDIO to accelerate VQEs with
numerical simulations in Section IV. Last, in Section V, we conclude this
study and discuss future research directions.

##  II Quantum distributed optimization scheme

We depict the paradigm of QUDIO in the right panel of Fig. 1 and present the
corresponding Pseudocode in Alg. 1. Conceptually, QUDIO is constituted by a
classical central server and $Q$ local nodes
$\\{\mathcal{Q}_{i}\\}_{i=1}^{Q}$, where each local node consists of a
classical optimizer and a single quantum processor. The algorithmic
implementation of QUDIO include three steps.

  1. 1.

At the preprocessing stage, the central server partitions the given problem
into $Q$ subproblems and allocates each subproblem to each local node (see
Sections III and IV for concrete explanations).

  2. 2.

The training procedure of QUDIO follows an iterative manner. Set the total
number of global and local steps as $T$ and $W$, respectively. At the $t$-th
global step with $t\in[T]$, the central server first dispatches the
synchronized parameters $\bm{\theta}^{(t)}$ to $Q$ local nodes, which
correspond to initial parameters for local updates (Line 5). With a slight
abuse of notation, we denote $\bm{\theta}^{(t,w)}_{i}$ as the trainable
parameters for the $i$-th local node at the $w$-th local step. The initial
parameters for all local nodes satisfy

$\bm{\theta}^{(t)}=\bm{\theta}^{(t,w=0)}_{i},~{}\forall i\in[Q].$

After initialization, all local nodes $\\{\mathcal{Q}_{i}\\}$ proceed $W$
iterations independently (highlighted by the red region in Lines 7-8). Once
all local updates are fulfilled, the central server collects parameters
$\\{\bm{\theta}^{(t,W)}_{i}\\}_{i=1}^{Q}$ from all local nodes to execute
synchronization of trainable parameters (Line 11), i.e.,

$\bm{\theta}^{(t+1)}=\frac{1}{Q}\sum_{i=1}^{Q}\bm{\theta}_{i}^{(t,W)},$

which completes the $t$-th global step.

  3. 3.

Through repeating the above procedure with $T$ global steps, the central
server outputs the synchronized $\bm{\theta}^{(T)}$ as the optimized
parameters.

Algorithm 1 The Pseudocode of QUDIO. The codes highlighted by the yellow and
pink shadows refer to execute them on the local nodes and the central server,
respectively.

1:Input: The initialized parameters $\bm{\theta}^{(0)}\in[0,2\pi)^{d_{Q}}$,
the employed loss function $\mathcal{L}$, the given dataset/Hamilton, the
hyper-parameters $\\{Q,\eta,W,T\\}$

2:The central server partitions the given problem into $Q$ parts and allocates
them to $Q$ local nodes

3:for $t=0,\cdots,T-1$ do

4: for Quantum processor $\mathcal{Q}_{i}$, $\forall i\in[Q]$ in parallel do

5: $\bm{\theta}_{i}^{(t,0)}=\bm{\theta}^{(t)}$

6: for $w=0,\cdots,W-1$ do

7: Compute the estimate gradients $g_{i}^{(t,w)}$

8: Update $\bm{\theta}_{i}^{(t,w+1)}=\bm{\theta}_{i}^{(t,w)}-\eta
g_{i}^{(t,w)}$

9: end for

10: end for

11: Synchronize
$\bm{\theta}^{(t+1)}=\frac{1}{Q}\sum_{i=1}^{Q}\bm{\theta}_{i}^{(t,W)}$

12:end for

13:Output: $\bm{\theta}^{(T)}$

In principle, compared with original VQAs with single quantum processor, this
parallel optimization mechanism enables QUDIO to reduce the computing time by
a constant factor equal to the number of local nodes $Q$. This linear speedup
is warranted by the small amount of trainable parameters for most VQAs.
Notably, this property differs the algorithmic design between distributed DNNs
and distributed VQAs, where the former concerns a considerable communication
bottleneck caused by billions of trainable parameters of DNNs. Furthermore,
QUDIO is highly compatible and can be seamlessly embedded into cloud
computing, since it supports various types of quantum processors to set up
local nodes and its central server is purely classical.

A core component in QUDIO is the approach of decomposing the given problem
into $Q$ parts, which in turn results in the varied forms of the estimated
gradients $\\{g_{i}^{(t,w)}\\}$ in Line 7 of Alg. 1. For the purpose of
elucidating, in the later context, we separately elaborate how to decompose
the given problem and calculate the estimated gradients when applying QUDIO to
speed up the training of QNNs and VQEs in Section III and Section IV,
respectively.

Remark. Although this study concentrates on QNNs and VQEs, our proposal can be
effectively extended to speed up other VQAs such as quantum approximate
optimization algorithms amaro2021filtering ; farhi2014quantum ;
hadfield2019quantum ; zhang2021neural ; zhou2020quantum .

##  III Accelerate QNN by QUDIO

Let us first formalize the classification task discussed in this section.
Denote the given dataset as $\mathcal{D}=\\{\bm{x}_{i},y_{i}\\}_{i=1}^{n}$,
where $\bm{x}_{i}\in\mathbb{R}^{D_{c}}$ and $y_{i}\in\mathbb{R}$ refer to the
features and label for the $i$-th example. The generic form of the output of
QNNs du2020learnability ; havlivcek2019supervised  is

$h(\bm{\theta},O,\rho_{i})=\operatorname{Tr}(OU(\bm{\theta})\rho_{i}U(\bm{\theta})^{\dagger}),~{}\forall
i\in[n],$ (1)

where $\rho_{i}$, $U(\bm{\theta})=\prod_{l=1}^{L}U_{l}(\bm{\theta})$, and
$O\in\mathbb{C}^{2^{N}\times 2^{N}}$ are the encoded $N$-qubit state
corresponding to $\bm{x}_{i}$, the ansatz with the circuit depth $L$ and
$d_{Q}$ parameters (i.e., $\bm{\theta}\in[0,2\pi)^{d_{Q}}$), and the fixed
quantum observable, respectively. Notably, the versatility of QNNs arises from
diverse data embedding methods (e.g., preparing $\rho_{i}$ by basis,
amplitude, and qubit encoding methods larose2020robust ), flexible
architectures of the ansatz $U(\bm{\theta})$ (e.g., hardware-efficient and
tensor-network based ansatzes), and the agile choice of $O$. The aim of QNNs
is to seek the optimal parameters $\bm{\theta}^{*}$ that minimize a predefined
loss $\mathcal{L}$. Throughout the whole work, we specify $\mathcal{L}$ as the
mean square error with $l_{2}$-norm regularizer, i.e.,

$\mathcal{L}(\bm{\theta},\mathcal{D})=\frac{1}{2n}\sum_{i=1}^{n}(h(\bm{\theta},O,\rho_{i})-y_{i})^{2}+\lambda\|\bm{\theta}\|^{2}_{2},$
(2)

where $\lambda\geq 0$ refers to the regularizer coefficient. The optimization
of $\bm{\theta}$ can be accomplished by either using gradient-free or
gradient-based methods.

We now elaborate on how QUDIO in Alg. 1 accelerates the training of QNN
following Eq. (2). Concretely, at the preprocessing stage, the central server
splits the dataset $\mathcal{D}$ into $Q$ subgroups
$\\{\mathcal{D}_{i}\\}_{i=1}^{Q}$ and assigns them into $Q$ local nodes
$\\{\mathcal{Q}_{i}\\}_{i=1}^{Q}$. In the training procedure, QUDIO harnesses
the following iterative strategy to optimize the trainable parameters. At the
$t$-th global step, when $Q$ local nodes receive the synchronized parameters
$\bm{\theta}^{(t)}$ sent by the central server, they proceed $W$ local updates
independently. Let $\bm{\theta}_{i}^{(t,w=0)}=\bm{\theta}^{(t)}$ for $\forall
i\in[Q]$ and $\forall w\in[W]$. The updating rule of $\mathcal{Q}_{i}$
associated with the stochastic gradient descent optimizer goodfellow2016deep
yields

$\bm{\theta}_{i}^{(t,w+1)}=\bm{\theta}_{i}^{(t,w)}-\eta
g_{i}(\bm{\theta}_{i}^{(t,w)},\bm{x}_{i}^{(t,w)})\in[0,2\pi)^{d_{Q}},$ (3)

where $\eta$ is the learning rate, the example
$(\bm{x}_{i}^{(t,w)},y^{(t)}_{i})$ is uniformly sampled from
$\mathcal{D}_{i}$, and $g_{i}(\bm{\theta}_{i}^{(t,w)},\bm{x}_{i}^{(t,w)})$
denotes the estimation of
$\nabla\mathcal{L}(\bm{\theta}_{i}^{(t,w)},\bm{x}_{i}^{(t,w)})$ induced by the
system noise and sample error. Once all local updates are completed, the
central server receives parameters $\\{\bm{\theta}_{i}^{(t,W)}\\}_{i=1}^{Q}$
and synchronizes them to update the global trainable parameters, i.e.,

$\bm{\theta}^{(t+1)}=\frac{1}{Q}\sum_{i=1}^{Q}\bm{\theta}^{(t,W)}_{i}.$ (4)

Through repeating the above process with $T$ times, QUDIO outputs
$\bm{\theta}^{(T)}$ as the trained parameters.

In the remainder of this section, we first explain the acquisition of the
estimated gradients $g_{i}^{(t,w)}$ and analyze the convergence of QUDIO. We
then benchmark performance of QUDIO towards image classification tasks.

###  III.1 The acquisition of the estimated gradients

The explicit form of $g_{i}(\bm{\theta}_{i}^{(t,w)},\bm{x}_{i}^{(t,w)})$ in
Eq. (3) is established on the analytic gradient
$\nabla\mathcal{L}_{i}(\bm{\theta}_{i}^{(t,w)},\bm{x}_{i}^{(t,w)})$. In this
perspective, here we first recap the mathematical expression of
$\nabla\mathcal{L}_{i}(\bm{\theta}_{i}^{(t,w)},\bm{x}_{i}^{(t,w)})$. According
to Eq. (3), the updating rule of $\mathcal{Q}_{i}$ in the ideal scenario
yields

$\bm{\theta}_{i}^{(t,w+1)}=\bm{\theta}_{i}^{(t,w)}-\eta\nabla\mathcal{L}_{i}(\bm{\theta}_{i}^{(t,w)},\bm{x}_{i}^{(t,w)}),$
(5)

where
$\mathcal{L}_{i}(\bm{\theta}^{(t,w)}_{i},\bm{x}_{i}^{(t,w)})=\frac{1}{2}(\hat{y}_{i}^{(t,w)}-y_{i}^{(t)})^{2}+\lambda\|\bm{\theta}^{(t,w)}_{i}\|^{2}_{2}$
and $\hat{y}_{i}^{(t,w)}=h(\bm{\theta}^{(t,w)}_{i},O,\rho_{i}^{(t)})$ refers
to the prediction of QNN for the sampled example $\bm{x}_{i}^{(t,w)}$ as
defined in Eq. (1). The evaluation of the analytic gradient can be achieved
via the parameter shift rule mitarai2018quantum ; schuld2019evaluating , i.e.,
the $j$-th component of
$\nabla_{j}\mathcal{L}_{i}(\bm{\theta}_{i}^{(t,w)},\bm{x}_{i}^{(t,w)})$ for
$\forall j\in[d_{Q}]$ satisfies

$(\hat{y}_{i}^{(t,w)}-y_{i}^{(t)})\frac{\hat{y}_{i}^{(t,w,+_{j})}-\hat{y}_{i}^{(t,w,-_{j})}}{2}+\lambda\bm{\theta}_{i,j}^{(t,w)},$
(6)

where
$\hat{y}_{i}^{(t,w,\pm_{j})}=h(\bm{\theta}^{(t,w,\pm)}_{i},O,\rho_{i}^{(t)})$
denotes the outputs of QNN with shifted parameters
$\bm{\theta}^{(t,w,\pm)}_{i}=\bm{\theta}^{(t,w)}_{i}\pm\frac{\pi}{2}\bm{e}_{j}$.

In the NISQ scenario, the system noise and the sample error forbid the
acquisition of the analytic gradients. Instead, the classical optimizer can
only collect the estimated gradients. More precisely, suppose that the
depolarization noise channel $\mathcal{N}_{p}(\cdot)$ is injected to each
quantum circuit depth, i.e.,

$\mathcal{N}_{p}(\rho)=(1-p)\rho+p\frac{\mathbb{I}}{2^{N}}.$ (7)

The output state before measurements is

$\gamma_{i}^{(t,w)}=(1-\tilde{p})U(\bm{\theta}^{(t,w)})\rho_{i}^{(t)}U(\bm{\theta}^{(t,w)})^{\dagger}+\tilde{p}\frac{\mathbb{I}}{2^{N}},$

where $\tilde{p}=1-(1-p)^{L_{Q}}$ and $L_{Q}$ refers to the total circuit
depth du2020learnability . Suppose $\\{O,\mathbb{I}-O\\}$ in Eq. (1) refers to
a two-outcome positive operator valued measure (POVM) nielsen2010quantum .
Then a quantum measurement on the state $\gamma_{i}^{(t,w)}$ produces the
outcome that can be viewed as a binary random variable with the Bernoulli
distribution, i.e.,
$V_{k}^{(t,w)}\sim\operatorname{Ber}(\operatorname{Tr}(O\gamma_{i}^{(t,w)}))$
. In this way, the sample mean
$\bar{y}_{i}^{(t,w)}=\sum_{k=1}^{K}V_{k}^{(t,w)}/K$ corresponding to
$\operatorname{Tr}(O\gamma_{i}^{(t,w)})$ is obtained after $K$ measurements.
In conjunction with above observations and the analytic gradients in Eq. (6),
the explicit form of the estimated gradients for the $j$-th component
$g_{i,j}(\bm{\theta}_{i}^{(t,w)},\bm{x}_{i}^{(t,w)})$ with $\forall
j\in[d_{Q}]$ satisfies

$(\bar{y}_{i}^{(t,w)}-y_{i}^{(t)})\frac{\bar{y}_{i}^{(t,w,+_{j})}-\bar{y}_{i}^{(t,w,-_{j})}}{2}+\lambda\bm{\theta}_{i,j}^{(t,w)},$
(8)

where $\bar{y}_{i}^{(t,w,\pm_{j})}=\frac{1}{K}\sum_{k=1}^{K}V_{k}^{(\pm_{j})}$
are estimated outputs with shifted parameters, i.e.,
$V_{k}^{(\pm_{j})}\sim\operatorname{Ber}(\operatorname{Tr}(O\gamma_{i}^{(\pm_{j})}))$
and
$\gamma_{i}^{(\pm_{j})}=(1-\tilde{p})U(\bm{\theta}^{(t,w,\pm_{j})})\rho_{i}U(\bm{\theta}^{(t,w,\pm_{j})})^{\dagger}+\tilde{p}\frac{\mathbb{I}}{2^{N}}$.

To better understand the capability of QUDIO, we further analyze its
convergence rate. The convergence is quantified by the utility
$R_{1}=\frac{1}{T}\sum_{t=1}^{T}\mathbb{E}[\|\nabla\mathcal{L}(\bm{\theta}^{(t)})\|^{2}]$,
where the expectation is taken over the randomness of data sample, the
imperfection of quantum system, and the finite quantum measurements.
Intuitively, the metric $R_{1}$ evaluates how far QNN is away to the
stationary points, which is a standard measure in non-convex optimization
theory jain2017non ; sun2019optimization . The following theorem summarizes
the convergence rate of QUDIO, whose proof is provided in Appendix A.

######  Theorem 1.

Assume that the discrepancy between the collected local gradients and the
analytic gradients is bounded, i.e., $\forall i\in[Q]$ and
$\forall\bm{\theta}\in[0,2\pi)^{d_{Q}}$, there exists

$\mathbb{E}_{\bm{x}_{i}^{(t,w)}\sim\mathcal{D}_{i}}\left[\left\|g_{i}^{(t,w)}-\nabla\bar{\mathcal{L}}_{i}\left(\bm{\theta}_{i}^{(t,w)},\bm{x}_{i}^{(t,w)}\right)\right\|^{2}\right]\leq\sigma,.$
(9)

Following notations in Eqs. (1)-(8), when the system noise is modeled by the
depolarization channel in Eq. (7) and the number of measurements to estimate
the expectation value is $K$, under a mild assumption, the convergence of
QUDIO yields

$R_{1}\leq O\left(\lambda
d_{Q}\sqrt{\frac{S}{T}}+\sqrt{\frac{S}{T}}\left(4W^{2}\sigma^{2}+2W^{2}G_{2}^{2}\right)+C_{1}\right),$

where $C_{1}\sim O(W^{2}(\tilde{p}d_{Q}+\tilde{p}d_{Q}/K))$,
$\tilde{p}=(1-(1-p)^{L_{Q}})$, and $L_{Q}$ is the total circuit depth.

Figure 2: Simulation results of QUDIO towards hand-written digits image
classification. (a) A visualization of some training examples sampled from the
MNIST dataset. (b) Scaling behavior of QUDIO in clock-time for increasing
number of local nodes $Q$ for varied number of local steps $W$. The labels
‘$W=a$ (I)’ and ‘$W=a$ (N)’ refer that the total number of local iterations is
$W=a$ under the ideal and NISQ scenarios respectively. The hyper-parameters
settings for the NISQ case are $p=10^{-5}$ and $K=100$. (c) A box plot that
illustrates the achieved test accuracy of QUDIO with varied $W$ and $Q$ in the
NISQ scenario, where the hyper-parameters settings are same with those
described in (b).

The results of Theorem 1 deliver three-fold implications. First, large system
noise and few number of measurements may induce the optimization of QUDIO to
be divergent, since the term $C_{1}$ is independent with $T$ and is amplified
by $p$ and $1/K$. This observation hints the importance of integrating error
mitigation techniques cai2020mitigating ; du2020quantum ; mcclean2020decoding
; strikis2020learning  into QUDIO to enhance its trainability. Second, in the
NISQ scenario, reducing the iteration number of local updating $W$ suggests a
better performance, because $C_{1}$ is proportional to $W$. This phenomenon is
starkly contrast with classical distributed optimization methods, which adopt
large $W$ to alleviate the communication overhead. Last but not least, under
the ideal setting, the convergence rate between conventional QNNs and QUDIO is
identical, i.e., both of them scale with $O(1/\sqrt{T})$ with respect to the
step number $T$ du2020learnability . Celebrated by the joint optimization
strategy, the similar convergence rate warrants that QUDIO promises a linear
runtime speedup with respect to the increased number of local nodes $Q$.

Remark. The developed tools in the proof of Theorem 1 can be extended to
analyze QNNs with other loss functions, quantum noisy models, and optimizers.
Moreover, we would like to emphasize that naively imitating classical
distributed algorithms to design distributed VQAs is suboptimal, since the
inevitable biased gradient information in the quantum scenario may incur a
deficient convergence.

###  III.2 Numerical simulations

We carry out numerical simulations to exhibit how QUDIO accelerates QNNs when
dealing with a standard binary classification task with a large size of
training examples. More precisely, the exploited dataset is distilled from a
hand-written digits images dataset, called MNIST dataset lecun1998mnist ,
which contains $256$ training examples and $500$ test examples labeled with
digits ‘0’ and ‘1’. Fig. 2(a) visualizes some examples in the distilled
dataset. The amplitude encoding method and the hardware-efficient ansatz are
used to set up all $Q$ local nodes. The hyper-parameters settings are as
follows. The number $Q$ ranges from $1$ to $32$. The number of local
iterations has six settings, i.e., $W\in[1,2,4,8,16,32]$. In the NISQ setting,
we set $K\in[5,100]$ and $p\in[10^{-4},10^{-1}]$. Each setting is repeated
with $5$ times to collect the statistical results. See Appendix B for the
omitted implementation details.

To better quantify the performance of QUDIO from different angles, we
introduce two metrics, i.e., the speedup to accuracy and the test accuracy, to
evaluate the achieved results. Namely, the former considers the speedup ratio
of QUDIO compared with the setting $Q=1$, i.e., supposes that the train
accuracy reaches a predefined threshold (e.g., $95\%$) in $T_{1}$ ($T_{2}$)
clock-time for $Q=1$ ($Q=a$), the speedup to accuracy is evaluated by $T1/T2$.
The latter allows us to compare the top test accuracy of QUDIO within a fixed
number of global steps $T$ with varied $Q$ and $W$.

Fig. 2 exhibits our simulation results. As shown in Fig. 2 (b), for both the
ideal and NISQ cases, QUDIO gains the speedup when increasing the number of
local nodes $Q$. Strikingly, QUDIO can even reach a superlinear speedup in the
NISQ scenario when $W=32$, e.g., it achieves $5.355$ times speedup for $Q=4$.
This phenomenon indicates that QUDIO is insensitive to the communication
bottleneck, which differs from distributed-DNNs Jeffrey2012 . Moreover, the
distinct scaling behavior of QUDIO between the ideal and the NISQ cases is
mainly caused by the fact that the evaluation of the analytic gradients in the
ideal case is extremely fast and the communication cost dominates the runtime
cost. Fig. 2(c) pictures the statistical results of test accuracy for QUDIO.
For all settings of $Q$, an increased $W$ generally degrades the performance
of QUDIO in the statistical view. These results partially echo with Theorem 1
such that larger $W$ suggests worse performance. An evidence is when $Q=16$,
QUDIO achieves the best test accuracy in $W=4$. See Appendix B for more
results and comprehensive investigation about the capability of QUDIO.

Figure 3: Simulation results of QUDIO towards the ground state energy
estimation of hydrogen molecule. (a) The speedup ratio with respect to
different number of local nodes $Q$ and local iterations $W$. (b) The
potential energy surface estimated by QUDIO. The black dotted line represents
the exact ground state energy. The inner plot compares the error between the
ground truth and the estimated results of QUDIO with the case of
$0.3\mathrm{\AA}$ inter-atomic distance.

##  IV Accelerate VQE by QUDIO

Variational quantum eigen-solvers (VQEs) Cervera2021meta-variational ;
kandala2017hardware ; peruzzo2014variational ; tang2021qubit  belong to
another pivotal subclass of VQAs and have a broad usage of tackling quantum
chemistry problems such as ground state estimation. The paradigm of VQEs is
analogous to QNNs and other VQAs, which is completed by minimizing a problem-
specific loss via gradient descent methods. Define the input Hamiltonian as

$H=\sum_{i=1}^{n}\alpha_{i}H_{i}\in\mathbb{C}^{2^{N}\times 2^{N}},$

where $H_{i}$ refers to the $i$-th local Hamiltonian term and
$\alpha_{i}\in\mathbb{R}$ is the corresponding coefficient. Without loss of
generality, suppose that
$H_{i}\in\\{\operatorname{X},\operatorname{Y},\operatorname{Z},\operatorname{I}\\}^{\otimes
N}$ is generated by Pauli operators. The loss function of VQEs yields

$\mathcal{L}(\bm{\theta},H)=\operatorname{Tr}(HU(\bm{\theta})\rho_{0}U(\bm{\theta})^{\dagger}),$
(10)

where $\rho_{0}=(\ket{0}\bra{0})^{\otimes N}$ is a fixed $N$-qubit quantum
state and $U(\bm{\theta})$ is the ansatz defined in Eq. (1). Note that to
evaluate $\mathcal{L}(\bm{\theta},H)$, conventional VQEs use a single quantum
chip to sequentially compute the results
$\\{\operatorname{Tr}(H_{i}U(\bm{\theta})\rho_{0}U(\bm{\theta})^{\dagger})\\}_{i=1}^{n}$
followed by a linear combination with $\\{\alpha_{i}\\}$. This implies the
computational hardness of VQEs when the number of terms $n$ becomes large.

When applying QUDIO to accelerate VQEs, the central server splits the set of
local Hamiltonians and their coefficients into $Q$ subgroups
$\\{\mathcal{S}_{i}\\}_{i=1}^{Q}$ and assigns them into $Q$ local nodes
$\\{\mathcal{Q}_{i}\\}_{i=1}^{Q}$. For instance, in the extreme case of $Q=n$,
we have $\mathcal{S}_{i}=\\{(\alpha_{i},H_{i})\\}$. For the case of $Q<n$, $n$
local hamiltonians are divided into $Q$ subgroups such that
$\mathcal{S}_{i}=\cup_{j\in S_{i}}\\{(\alpha_{j},H_{j})\\}$, where $S_{i}$
refers to the $i$-th index set with $\cup_{i=1}^{Q}S_{i}=[n]$ and $S_{i}\cap
S_{j}=\emptyset$ when $i\neq j$. Moreover, the training procedure follows the
same manner with accelerating QNNs. The updating rule of the local node
$\mathcal{Q}_{i}$ for $\forall i\in[Q]$ yields

$\bm{\theta}_{i}^{(t,w+1)}=\bm{\theta}_{i}^{(t,w)}-\eta
g_{i}(\bm{\theta}_{i}^{(t,w)},H_{S_{i}}),$ (11)

where $g_{i}(\bm{\theta}_{i}^{(t,w)},H_{S_{i}})$ refers to the estimated
gradient of
$\nabla\mathcal{L}_{i}(\bm{\theta}_{i}^{(t,w)},H_{s_{i}})=\operatorname{Tr}(H_{S_{i}}U(\bm{\theta}_{i}^{(t,w)})\rho_{0}U(\bm{\theta}_{i}^{(t,w)})^{\dagger})$
with $H_{S_{i}}=\sum_{j\in S_{i}}\alpha_{j}H_{j}$.

In the subsequent subsections, we first exhibit the explicit form of the
estimated gradient $g_{i}(\bm{\theta}_{i}^{(t,w)},H_{S_{i}})$ and then conduct
extensive numerical simulations to validate performance of QUDIO. Remarkably,
to facilitate simulation, here we mainly focus on the scenario in which the
estimation error is caused by the finite number of measurement and the system
is noiseless.

###  IV.1 The acquisition of the estimated gradients

Before diving into deriving the explicit of estimated gradient, let us first
recall the analytic gradient of
$\mathcal{L}_{i}(\bm{\theta}_{i}^{(t,w)},H_{s_{i}})$. Specifically, based on
the parameter shift rule, the $j$-th component of the analytic gradients is

$\nabla_{j}\mathcal{L}_{i}(\bm{\theta}_{i}^{(t,w)},H_{S_{i}})=\frac{\hat{y}_{i}^{(t,w,+_{j})}-\hat{y}_{i}^{(t,w,-_{j})}}{2},$
(12)

where
$\bm{\theta}^{(t,w,\pm)}_{i}=\bm{\theta}^{(t,w)}_{i}\pm\frac{\pi}{2}\bm{e}_{j}$
and $\hat{y}_{i}^{(t,w,\pm_{j})}=\sum_{j\in
S_{i}}\alpha_{j}\operatorname{Tr}(H_{j}U(\bm{\theta}_{i}^{(t,w,\pm)})\rho_{0}U(\bm{\theta}_{i}^{(t,w,\pm)})^{\dagger})$.

When a finite number of measurement is allowable, the trace terms of
$\hat{y}_{i}^{(t,w,\pm_{j})}$ in Eq. (12) can only be acquired with estimation
error. The detailed procedure to estimate
$\operatorname{Tr}(H_{j}U(\bm{\theta}_{i}^{(t,w,\pm)})\rho_{0}U(\bm{\theta}_{i}^{(t,w,\pm)})^{\dagger})$
is as follows. To estimate this result by measuring the quantum state
$\rho=U(\bm{\theta}_{i}^{(t,w)})\rho_{0}U(\bm{\theta}_{i}^{(t,w)})^{\dagger}$
along the computational basis, an alignment operation should be executed.
Mathematically, the quantum state $\rho$ needs to interact with the unitary
operator $R_{i}$ to generate the quantum state $\rho^{\prime}=R_{j}\rho
R_{j}^{\dagger}$, where $R_{j}$ is composed of a sequence of rotational
single-qubit gates whose row vectors are the eigen-basis of $H_{j}$. Define
$\bm{h}_{i}$ as a vector that collects the eigenvalues of $H_{j}$,
$V_{k}\sim\mathrm{Cat}(2^{N},\bm{p}(\rho^{\prime}))$ as a random variable
following the categorical distribution (also called generalized Bernoulli
distribution), and $\bm{p}(\rho^{\prime})\in\mathbb{R}^{2^{N}}$ refers as a
discrete distribution with
$\bm{p}_{j}(\rho^{\prime})=\mathrm{Tr}(\rho^{\prime}\ket{j}\bra{j})$.
Following the above notations, the term
$\operatorname{Tr}(H_{j}\rho^{\prime})$ is estimated by
$\frac{1}{K}\sum_{k=1}^{K}\bm{h}_{i,V_{k}}$, where $\bm{h}_{i,V_{k}}$ is the
$V_{k}$-th eigenvalue of $H_{j}$. Based on Eq. (12), the estimated gradient
satisfies

$g_{i,j}(\bm{\theta}_{i}^{(t,w)},H_{S_{i}})=\frac{\bar{y}_{i}^{(t,w,+_{j})}-\bar{y}_{i}^{(t,w,-_{j})}}{2},\forall
j\in[d_{Q}],$ (13)

where $\bar{y}_{i}^{(t,w,\pm_{j})}=\sum_{i\in
S_{i}}\alpha_{i}\frac{1}{K}\sum_{k=1}^{K}\bm{h}_{i,V_{k}^{(\pm_{j})}}$ and the
definition of
$V_{k}^{(\pm_{j})}\sim\mathrm{Cat}(2^{N},\bm{p}(\rho^{(\pm_{j})}))$ follows
the same manner with Eq. (8).

Considering that the training procedure is exactly identical to the way of
applying QUDIO to accelerate QNNs, Theorem 1 can also describe the convergence
behavior of QUDIO for accelerating VQEs.

###  IV.2 Numerical simulations

We perform numerical simulations to validate the effectiveness of QUDIO to
accelerate conventional VQEs. To do so, we apply QUDIO to estimate the ground
state energy of hydrogen molecule $\mathrm{H}_{2}$ with varied bond distance,
whose Hamiltonian contains $n=15$ local Hamiltonian terms and requires $N=4$
qubits with $H=\sum_{i=1}^{15}\alpha_{i}H_{i}\in\mathbb{C}^{2^{4}\times
2^{4}}$ kandala2017hardware . The implementation of all local nodes mainly
follows the proposal kandala2017hardware  such that the input quantum state is
$\rho_{0}=\ket{1100}\bra{1100}$ and the trainable unitary refers to the
hardware-efficient ansatz. The hyper-parameters settings are as follows. The
number of local nodes and local iterations is set as $Q\in\\{1,2,4,8\\}$ and
$W\in[1,2,4,8]$, respectively. We fix the number of measurements to be
$K=100$. Each setting is repeated with $5$ times to collect the statistical
results. See Appendix C for the omitted implementation details.

The simulation results are shown in Fig. 3. In particular, the left subplot
illustrates the speedup ratio of QUDIO in terms of the factors $Q$ and $W$.
For all settings of $W$, QUDIO gains runtime speedups by involving more local
nodes. For example, we obtain $1.469$ times acceleration when utilizing $8$
local nodes to optimize VQE with $1$ local update. Although QUDIO provides
certain speedups in the task of estimating the ground state energy of hydrogen
molecule, there is a clear gap towards the linear speedup ratio. It is
noteworthy that this gap arises from the simplicity of the manipulated
problem, where the communication overhead dominates the total computational
runtime. We expect that QUDIO has the ability to earn higher speedup ratio for
large-scale tasks.

The potential energy surface estimated by QUDIO is presented in Fig. 3(b). The
outer plot suggests that for all bond distance settings ranging from
$0.3\mathrm{\AA}$ to $2.1\mathrm{\AA}$, QUDIO obtains the best performance
with $W=1$, which is almost the same with the exact values. By contrast, there
exists an apparent separation between the exact values and the estimated
results of QUDIO with $W=8$. The inner plot further evidences this phenomenon.
Specifically, when the bond distance equals to $0.3\mathrm{\AA}$, QUDIO
witnesses the largest approximation error $0.9Ha$ with $W=8$ and $Q=2$, while
the error is reduced to nearly zero with $W=1$ regardless of the number of
local nodes. All of the above observations collaborate with Theorem 1, where
decreasing $W$ warrants a better performance. Refer to Appendix C for deep
comprehension.

##  V Discussion and conclusion

In this study, we devise QUDIO to accelerate VQAs with multiple quantum
processors. We also provide theoretical analysis about how the system noise
and the number of measurements influence the convergence of QUDIO. An
attractive feature is that in the ideal setting, QUDIO obeys the asymptotic
convergence rate with conventional QNNs, which ensures its runtime speedup
with respect to the increased number of local nodes. The achieved numerical
simulation results confirm the effectiveness of our proposal. Particularly, in
the NISQ scenario, QUDIO can achieve superline speedups in the measure of
time-to-accuracy.

We remark that there are three orthogonal research directions towards the
investigation of distributed VQAs. First, instead of employing the
synchronization approach used in QUDIO, it is intrigued to design asynchronous
distributed-VQAs schemes with convergence guarantees, which may further reduce
the communication overhead and maximally utilize quantum processors with
distinct qualities. Second, with the aim of reducing the runtime cost, it is
important to integrate the effective measurement reduction algorithms with
QUDIO and other distributed VQAs. Repressive examples contain grouping
compatible operators kandala2017hardware ; zhao2020measurement  and classical-
shadows based methods huang2020predicting ; struchalin2021experimental . Last,
a promising direction is combining QUDIO with a recent work
zhang2021variational , which splits the input quantum circuits into several
individual quantum circuits with distributed optimization.

The employment of cloud computing to execute QUDIO naturally invokes the
security issue lu2020quantum . For example, how to defend adversarial attack
or prevent the private information leakage during the training procedure.
Initial studies have leveraged some notions such as differential privacy
du2020quantum ; du2021quantum , hypothesis testing weber2021optimal , and
blind quantum computing li2021quantum  to address this issue. However, little
is known about how these strategies effect the convergence rate. A deep
understanding towards this topic is highly desired, which determines the
applicability of distributed VQAs.

For these reasons, QUDIO and its variants, which marry the distributed
techniques with VQAs, could substantially contribute to use NISQ machines to
accomplish real-world problems with quantum advantages.

## References

  * (1) Ian Goodfellow, Yoshua Bengio, Aaron Courville, and Yoshua Bengio.  Deep learning, volume 1.  MIT press Cambridge, 2016.
  * (2) Shuai Zhang, Lina Yao, Aixin Sun, and Yi Tay.  Deep learning based recommender system: A survey and new perspectives.  ACM Computing Surveys (CSUR), 52(1):1–38, 2019.
  * (3) Andrew W Senior, Richard Evans, John Jumper, James Kirkpatrick, Laurent Sifre, Tim Green, Chongli Qin, Augustin Žídek, Alexander WR Nelson, Alex Bridgland, et al.  Improved protein structure prediction using potentials from deep learning.  Nature, 577(7792):706–710, 2020.
  * (4) Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova.  Bert: Pre-training of deep bidirectional transformers for language understanding.  In Proceedings of the 2019 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 1 (Long and Short Papers), pages 4171–4186, 2019.
  * (5) Jeffrey Dean, Greg S. Corrado, Rajat Monga, Kai Chen, Matthieu Devin, Quoc V. Le, Mark Z. Mao, Marc’Aurelio Ranzato, Andrew Senior, Paul Tucker, Ke Yang, and Andrew Y. Ng.  Large scale distributed deep networks.  In Proceedings of the 25th International Conference on Neural Information Processing Systems - Volume 1, page 1223–1231, Red Hook, NY, USA, 2012. Curran Associates Inc.
  * (6) Emma Strubell, Ananya Ganesh, and Andrew McCallum.  Energy and policy considerations for deep learning in nlp.  In Proceedings of the 57th Annual Meeting of the Association for Computational Linguistics, pages 3645–3650, 2019.
  * (7) Jacob Biamonte, Peter Wittek, Nicola Pancotti, Patrick Rebentrost, Nathan Wiebe, and Seth Lloyd.  Quantum machine learning.  Nature, 549(7671):195, 2017.
  * (8) Richard P Feynman.  Simulating physics with computers.  International Journal of Theoretical Physics, 21(6/7), 1982.
  * (9) Aram W Harrow and Ashley Montanaro.  Quantum computational supremacy.  Nature, 549(7671):203, 2017.
  * (10) Frank Arute, Kunal Arya, Ryan Babbush, Dave Bacon, Joseph C Bardin, Rami Barends, Rupak Biswas, Sergio Boixo, Fernando GSL Brandao, David A Buell, et al.  Quantum supremacy using a programmable superconducting processor.  Nature, 574(7779):505–510, 2019.
  * (11) Marcello Benedetti, Erika Lloyd, Stefan Sack, and Mattia Fiorentini.  Parameterized quantum circuits as machine learning models.  Quantum Science and Technology, 4(4):043001, 2019.
  * (12) Kishor Bharti, Alba Cervera-Lierta, Thi Ha Kyaw, Tobias Haug, Sumner Alperin-Lea, Abhinav Anand, Matthias Degroote, Hermanni Heimonen, Jakob S Kottmann, Tim Menke, et al.  Noisy intermediate-scale quantum (nisq) algorithms.  arXiv preprint arXiv:2101.08448, 2021.
  * (13) M Cerezo, Andrew Arrasmith, Ryan Babbush, Simon C Benjamin, Suguru Endo, Keisuke Fujii, Jarrod R McClean, Kosuke Mitarai, Xiao Yuan, Lukasz Cincio, et al.  Variational quantum algorithms.  arXiv preprint arXiv:2012.09265, 2020.
  * (14) Yuxuan Du, Min-Hsiu Hsieh, Tongliang Liu, and Dacheng Tao.  Expressive power of parametrized quantum circuits.  Phys. Rev. Research, 2:033125, Jul 2020.
  * (15) Suguru Endo, Zhenyu Cai, Simon C Benjamin, and Xiao Yuan.  Hybrid quantum-classical algorithms and quantum error mitigation.  Journal of the Physical Society of Japan, 90(3):032001, 2021.
  * (16) John Preskill.  Quantum computing in the nisq era and beyond.  Quantum, 2:79, 2018.
  * (17) Amira Abbas, David Sutter, Christa Zoufal, Aurélien Lucchi, Alessio Figalli, and Stefan Woerner.  The power of quantum neural networks.  arXiv preprint arXiv:2011.00027, 2020.
  * (18) Leonardo Banchi, Jason Pereira, and Stefano Pirandola.  Generalization in quantum machine learning: a quantum information perspective.  arXiv preprint arXiv:2102.08991, 2021.
  * (19) Kaifeng Bu, Dax Enshan Koh, Lu Li, Qingxian Luo, and Yaobo Zhang.  Effects of quantum resources on the statistical complexity of quantum circuits.  arXiv preprint arXiv:2102.03282, 2021.
  * (20) Matthias C Caro, Elies Gil-Fuster, Johannes Jakob Meyer, Jens Eisert, and Ryan Sweke.  Encoding-dependent generalization bounds for parametrized quantum circuits.  arXiv preprint arXiv:2106.03880, 2021.
  * (21) Yuxuan Du, Zhuozhuo Tu, Xiao Yuan, and Dacheng Tao.  An efficient measure for the expressivity of variational quantum algorithms.  arXiv preprint arXiv:2104.09961, 2021.
  * (22) Hsin-Yuan Huang, Richard Kueng, and John Preskill.  Information-theoretic bounds on quantum advantage in machine learning.  Phys. Rev. Lett., 126:190505, May 2021.
  * (23) Hsin-Yuan Huang, Michael Broughton, Masoud Mohseni, Ryan Babbush, Sergio Boixo, Hartmut Neven, and Jarrod R McClean.  Power of data in quantum machine learning.  Nature communications, 12(1):1–9, 2021.
  * (24) Xinbiao Wang, Yuxuan Du, Yong Luo, and Dacheng Tao.  Towards understanding the power of quantum kernels in the nisq era.  arXiv preprint arXiv:2103.16774, 2021.
  * (25) Yadong Wu, Juan Yao, Pengfei Zhang, and Hui Zhai.  Expressivity of quantum neural networks.  arXiv preprint arXiv:2101.04273.
  * (26) Yuxuan Du, Min-Hsiu Hsieh, Tongliang Liu, and Dacheng Tao.  A grover-search based quantum learning scheme for classification.  New Journal of Physics, 23(2):023020, feb 2021.
  * (27) 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, 2019.
  * (28) Kosuke Mitarai, Makoto Negoro, Masahiro Kitagawa, and Keisuke Fujii.  Quantum circuit learning.  Physical Review A, 98(3):032309, 2018.
  * (29) Maria Schuld and Nathan Killoran.  Quantum machine learning in feature hilbert spaces.  Physical review letters, 122(4):040504, 2019.
  * (30) Kunkun Wang, Lei Xiao, Wei Yi, Shi-Ju Ran, and Peng Xue.  Quantum image classifier with single photons.  arXiv preprint arXiv:2003.08551, 2020.
  * (31) He-Liang Huang, Yuxuan Du, Ming Gong, Youwei Zhao, Yulin Wu, Chaoyue Wang, Shaowei Li, Futian Liang, Jin Lin, Yu Xu, et al.  Experimental quantum generative adversarial networks for image generation.  arXiv preprint arXiv:2010.06201, 2020.
  * (32) Manuel S Rudolph, Ntwali Toussaint Bashige, Amara Katabarwa, Sonika Johr, Borja Peropadre, and Alejandro Perdomo-Ortiz.  Generation of high resolution handwritten digits with an ion-trap quantum computer.  arXiv preprint arXiv:2012.03924, 2020.
  * (33) Daiwei Zhu, Norbert M Linke, Marcello Benedetti, Kevin A Landsman, Nhung H Nguyen, C Huerta Alderete, Alejandro Perdomo-Ortiz, Nathan Korda, A Garfoot, Charles Brecque, et al.  Training of quantum circuits on a hybrid quantum computer.  Science advances, 5(10):eaaw9918, 2019.
  * (34) Gavin E Crooks.  Performance of the quantum approximate optimization algorithm on the maximum cut problem.  arXiv preprint arXiv:1811.08419, 2018.
  * (35) Edward Farhi, Jeffrey Goldstone, and Sam Gutmann.  A quantum approximate optimization algorithm.  arXiv preprint arXiv:1411.4028, 2014.
  * (36) Stuart Hadfield, Zhihui Wang, Bryan O’Gorman, Eleanor G Rieffel, Davide Venturelli, and Rupak Biswas.  From the quantum approximate optimization algorithm to a quantum alternating operator ansatz.  Algorithms, 12(2):34, 2019.
  * (37) Matthew P Harrigan, Kevin J Sung, Matthew Neeley, Kevin J Satzinger, Frank Arute, Kunal Arya, Juan Atalaya, Joseph C Bardin, Rami Barends, Sergio Boixo, et al.  Quantum approximate optimization of non-planar graph problems on a planar superconducting processor.  Nature Physics, 17(3):332–336, 2021.
  * (38) Javier Alcazar, Vicente Leyton-Ortega, and Alejandro Perdomo-Ortiz.  Classical versus quantum models in machine learning: insights from a finance application.  Machine Learning: Science and Technology, 1(3):035003, 2020.
  * (39) Brian Coyle, Maxwell Henderson, Justin Chan Jin Le, Niraj Kumar, Marco Paini, and Elham Kashefi.  Quantum versus classical generative modelling in finance.  Quantum Science and Technology, 6(2):024013, 2021.
  * (40) Mark Hodson, Brendan Ruck, Hugh Ong, David Garvin, and Stefan Dulman.  Portfolio rebalancing experiments using the quantum alternating operator ansatz.  arXiv preprint arXiv:1911.05296, 2019.
  * (41) Jacob L Beckey, M Cerezo, Akira Sone, and Patrick J Coles.  Variational quantum algorithm for estimating the quantum fisher information.  arXiv preprint arXiv:2010.10488, 2020.
  * (42) Marco Cerezo, Alexander Poremba, Lukasz Cincio, and Patrick J Coles.  Variational quantum fidelity estimation.  Quantum, 4:248, 2020.
  * (43) Ryan LaRose, Arkin Tikku, Étude O’Neel-Judy, Lukasz Cincio, and Patrick J Coles.  Variational quantum state diagonalization.  npj Quantum Information, 5(1):1–10, 2019.
  * (44) Bela Bauer, Dave Wecker, Andrew J Millis, Matthew B Hastings, and Matthias Troyer.  Hybrid quantum-classical approach to correlated materials.  Physical Review X, 6(3):031045, 2016.
  * (45) Google AI Quantum et al.  Hartree-fock on a superconducting qubit quantum computer.  Science, 369(6507):1084–1089, 2020.
  * (46) Thomas E O’Brien, Bruno Senjean, Ramiro Sagastizabal, Xavier Bonet-Monroig, Alicja Dutkiewicz, Francesco Buda, Leonardo DiCarlo, and Lucas Visscher.  Calculating energy derivatives for quantum chemistry on a quantum computer.  npj Quantum Information, 5(1):1–12, 2019.
  * (47) 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:4213, 2014.
  * (48) A Avkhadiev, PE Shanahan, and RD Young.  Accelerating lattice quantum field theory calculations via interpolator optimization using noisy intermediate-scale quantum computing.  Physical review letters, 124(8):080501, 2020.
  * (49) Christian Kokail, Christine Maier, Rick van Bijnen, Tiff Brydges, Manoj K Joshi, Petar Jurcevic, Christine A Muschik, Pietro Silvi, Rainer Blatt, Christian F Roos, et al.  Self-verifying variational quantum simulation of lattice models.  Nature, 569(7756):355–360, 2019.
  * (50) 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):1–6, 2020.
  * (51) Yuxuan Du, Min-Hsiu Hsieh, Tongliang Liu, Shan You, and Dacheng Tao.  On the learnability of quantum neural networks.  arXiv preprint arXiv:2007.12369, 2020.
  * (52) Edward Farhi and Hartmut Neven.  Classification with quantum neural networks on near term processors.  arXiv preprint arXiv:1802.06002, 2018.
  * (53) Yang Qian, Xinbiao Wang, Yuxuan Du, Xingyao Wu, and Dacheng Tao.  The dilemma of quantum neural networks.  arXiv preprint arXiv:2106.04975, 2021.
  * (54) Stephen Boyd, Neal Parikh, and Eric Chu.  Distributed optimization and statistical learning via the alternating direction method of multipliers.  Now Publishers Inc, 2011.
  * (55) Robert Beals, Stephen Brierley, Oliver Gray, Aram W Harrow, Samuel Kutin, Noah Linden, Dan Shepherd, and Mark Stather.  Efficient distributed quantum computing.  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 469(2153):20120686, 2013.
  * (56) Ville Bergholm, Josh Izaac, Maria Schuld, Christian Gogolin, M Sohaib Alam, Shahnawaz Ahmed, Juan Miguel Arrazola, Carsten Blank, Alain Delgado, Soran Jahangiri, et al.  Pennylane: Automatic differentiation of hybrid quantum-classical computations.  arXiv preprint arXiv:1811.04968, 2018.
  * (57) Héctor Abraham et. al.  Qiskit: An open-source framework for quantum computing, 2019.
  * (58) David Amaro, Carlo Modica, Matthias Rosenkranz, Mattia Fiorentini, Marcello Benedetti, and Michael Lubasch.  Filtering variational quantum algorithms for combinatorial optimization.  arXiv preprint arXiv:2106.10055, 2021.
  * (59) Shi-Xin Zhang, Chang-Yu Hsieh, Shengyu Zhang, and Hong Yao.  Neural predictor based quantum architecture search.  arXiv preprint arXiv:2103.06524, 2021.
  * (60) Leo Zhou, Sheng-Tao Wang, Soonwon Choi, Hannes Pichler, and Mikhail D Lukin.  Quantum approximate optimization algorithm: Performance, mechanism, and implementation on near-term devices.  Physical Review X, 10(2):021067, 2020.
  * (61) Ryan LaRose and Brian Coyle.  Robust data encodings for quantum classifiers.  Physical Review A, 102(3):032420, 2020.
  * (62) Maria Schuld, Ville Bergholm, Christian Gogolin, Josh Izaac, and Nathan Killoran.  Evaluating analytic gradients on quantum hardware.  Physical Review A, 99(3):032331, 2019.
  * (63) Michael A Nielsen and Isaac L Chuang.  Quantum computation and quantum information.  Cambridge University Press, 2010.
  * (64) Prateek Jain and Purushottam Kar.  Non-convex optimization for machine learning.  Foundations and Trends® in Machine Learning, 10(3-4):142–336, 2017.
  * (65) Ruoyu Sun.  Optimization for deep learning: theory and algorithms.  arXiv preprint arXiv:1912.08957, 2019.
  * (66) Zhenyu Cai, Xiaosi Xu, and Simon C Benjamin.  Mitigating coherent noise using pauli conjugation.  npj Quantum Information, 6(1):1–9, 2020.
  * (67) Yuxuan Du, Min-Hsiu Hsieh, Tongliang Liu, Shan You, and Dacheng Tao.  Quantum differentially private sparse regression learning.  arXiv preprint arXiv:2007.11921, 2020.
  * (68) Jarrod R McClean, Zhang Jiang, Nicholas C Rubin, Ryan Babbush, and Hartmut Neven.  Decoding quantum errors with subspace expansions.  Nature communications, 11(1):1–9, 2020.
  * (69) Armands Strikis, Dayue Qin, Yanzhu Chen, Simon C Benjamin, and Ying Li.  Learning-based quantum error mitigation.  arXiv preprint arXiv:2005.07601, 2020.
  * (70) Yann LeCun.  The mnist database of handwritten digits.  http://yann. lecun. com/exdb/mnist/, 1998.
  * (71) Alba Cervera-Lierta, Jakob S. Kottmann, and Alán Aspuru-Guzik.  Meta-variational quantum eigensolver: Learning energy profiles of parameterized hamiltonians for quantum simulation.  PRX Quantum, 2:020329, May 2021.
  * (72) Abhinav Kandala, Antonio Mezzacapo, Kristan Temme, Maika Takita, Markus Brink, Jerry M Chow, and Jay M Gambetta.  Hardware-efficient variational quantum eigensolver for small molecules and quantum magnets.  Nature, 549(7671):242–246, 2017.
  * (73) Ho Lun Tang, VO Shkolnikov, George S Barron, Harper R Grimsley, Nicholas J Mayhall, Edwin Barnes, and Sophia E Economou.  qubit-adapt-vqe: An adaptive algorithm for constructing hardware-efficient ansätze on a quantum processor.  PRX Quantum, 2(2):020310, 2021.
  * (74) Andrew Zhao, Andrew Tranter, William M Kirby, Shu Fay Ung, Akimasa Miyake, and Peter J Love.  Measurement reduction in variational quantum algorithms.  Physical Review A, 101(6):062322, 2020.
  * (75) Hsin-Yuan Huang, Richard Kueng, and John Preskill.  Predicting many properties of a quantum system from very few measurements.  Nature Physics, 16(10):1050–1057, 2020.
  * (76) GI Struchalin, Ya A Zagorovskii, EV Kovlakov, SS Straupe, and SP Kulik.  Experimental estimation of quantum state properties from classical shadows.  PRX Quantum, 2(1):010307, 2021.
  * (77) Yu Zhang, Lukasz Cincio, Christian F. A. Negre, Piotr Czarnik, Patrick Coles, Petr M. Anisimov, Susan M. Mniszewski, Sergei Tretiak, and Pavel A. Dub.  Variational quantum eigensolver with reduced circuit complexity.  arXiv preprint arXiv:2106.07619, 2021.
  * (78) Sirui Lu, Lu-Ming Duan, and Dong-Ling Deng.  Quantum adversarial machine learning.  Physical Review Research, 2(3):033212, 2020.
  * (79) Yuxuan Du, Min-Hsiu Hsieh, Tongliang Liu, Dacheng Tao, and Nana Liu.  Quantum noise protects quantum classifiers against adversaries.  Physical Review Research, 3(2):023153, 2021.
  * (80) Maurice Weber, Nana Liu, Bo Li, Ce Zhang, and Zhikuan Zhao.  Optimal provable robustness of quantum classification via quantum hypothesis testing.  npj Quantum Information, 7(1):1–12, 2021.
  * (81) Weikang Li, Sirui Lu, and Dong-Ling Deng.  Quantum private distributed learning through blind quantum computing.  arXiv preprint arXiv:2103.08403, 2021.
  * (82) Stephen Boyd, Stephen P Boyd, and Lieven Vandenberghe.  Convex optimization.  Cambridge university press, 2004.
  * (83) Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Kopf, Edward Yang, Zachary DeVito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala.  Pytorch: An imperative style, high-performance deep learning library.  In Advances in Neural Information Processing Systems 32, pages 8024–8035. 2019.
  * (84) David E Rumelhart, Geoffrey E Hinton, and Ronald J Williams.  Learning representations by back-propagating errors.  nature, 323(6088):533–536, 1986.
  * (85) Sergey B Bravyi and Alexei Yu Kitaev.  Fermionic quantum computation.  Annals of Physics, 298(1):210–226, 2002.

##  Appendix A The proof of Theorem 1

The outline of this section is as follows. In Appendix A.1, we simplify some
notations and introduce basic concepts in optimization theory for ease of
discussion. Next, in Appendix A.2, we demonstrate the proof details of Theorem
1.

###  A.1 Basic notations and concepts

Notations. Let us first simplify some notations defined in the main text to
facilitate the derivation. Recall that the loss function of the $i$-th local
node is denoted by
$\mathcal{L}_{i}\left(\bm{\theta}_{i}^{(t,w)},\bm{x}_{i}^{(t,w)}\right)$ with
$\bm{x}_{i}^{(t,w)}$ being an example uniformly sampled from the sub-dataset
$\mathcal{D}_{i}$. In this section, when no confusion occurs, we define the
loss for the $i$-th local node $\forall i\in[Q]$ with the whole sub-dataset
$\mathcal{D}_{i}$ as

$\mathcal{L}_{i}\left(\bm{\theta}_{i}^{(t,w)},\mathcal{D}_{i}\right)\equiv\mathcal{L}_{i}\left(\bm{\theta}_{i}^{(t,w)}\right)=\frac{1}{2|\mathcal{D}_{i}|}\sum_{j=1}^{|\mathcal{D}_{i}|}\left(h(\bm{\theta}_{i}^{(t,w)},O,\rho_{j})-y_{j}\right)^{2}+\lambda\|\bm{\theta}_{i}^{(t,w)}\|^{2}_{2},$
(14)

where $\rho_{j}$ and $y_{j}$ refers to the encoded quantum example and
corresponding label with respect to the $j$-th example in $\mathcal{D}_{i}$.
Besides, we simplify the global loss function as

$\mathcal{L}\left(\bm{\theta}_{i}^{(t,w)},\mathcal{D}_{i}\right)\equiv\mathcal{L}\left(\bm{\theta}_{i}^{(t,w)}\right)=\frac{1}{Q}\sum_{i=1}^{Q}\mathcal{L}_{i}\left(\bm{\theta}_{i}^{(t,w)},\mathcal{D}_{i}\right).$
(15)

Following the same routine, the gradients of the loss
$\mathcal{L}_{i}(\bm{\theta}_{i}^{(t,w)})$ (or
$\mathcal{L}\left(\bm{\theta}_{i}^{(t,w)}\right)$) are written as
$\nabla\mathcal{L}_{i}(\bm{\theta}_{i}^{(t,w)},\mathcal{D}_{i})$ and
$\nabla\mathcal{L}_{i}(\bm{\theta}_{i}^{(t,w)})$ (or
$\nabla\mathcal{L}(\bm{\theta}_{i}^{(t,w)},\mathcal{D}_{i})$ and
$\nabla\mathcal{L}(\bm{\theta}_{i}^{(t,w)})$) interchangeably.

In QUDIO, to decrease runtime, the optimizer of the $i$-th local node for
$\forall i\in[Q]$ only requires the gradients of
$\mathcal{L}_{i}\left(\bm{\theta}_{i}^{(t,w)},\bm{x}_{i}^{(t,w)}\right)$
instead of the whole sub-dataset $\mathcal{D}_{i}$. To distinguish with
$\mathcal{L}_{i}(\bm{\theta}_{i}^{(t,w)})$, we denote

$\mathcal{L}_{i}\left(\bm{\theta}_{i}^{(t,w)},\bm{x}_{i}^{(t,w)}\right)=\bar{\mathcal{L}}_{i}\left(\bm{\theta}_{i}^{(t,w)}\right),~{}\text{and}~{}\nabla\mathcal{L}_{i}\left(\bm{\theta}_{i}^{(t,w)},\bm{x}_{i}^{(t,w)}\right)=\nabla\bar{\mathcal{L}}_{i}\left(\bm{\theta}_{i}^{(t,w)}\right).$
(16)

We note that the gradients between
$\nabla\mathcal{L}_{i}\left(\bm{\theta}_{i}^{(t,w)}\right)$ and
$\nabla\bar{\mathcal{L}}_{i}\left(\bm{\theta}_{i}^{(t,w)}\right)$ have the
relationship

$\mathbb{E}_{\bm{x}_{i}^{(t,w)}\sim\mathcal{D}_{i}}\left[\nabla\bar{\mathcal{L}}_{i}\left(\bm{\theta}_{i}^{(t,w)}\right)\right]=\nabla\mathcal{L}_{i}\left(\bm{\theta}_{i}^{(t,w)}\right),$
(17)

where the expectation is taken over the randomness of sampled examples.

Basic concepts in optimization theory. We introduce two definitions, i.e.,
$S$-smooth and $G$-Lipschitz boyd2004convex , which are employed to quantify
properties of loss functions and achieve the proof of Theorem 1.

######  Definition 1.

A function f is $S$-smooth over a set $\mathcal{C}$ if
$\nabla^{2}f(\bm{u})\preceq S\mathbb{I}$ with $S>0$ and
$\forall\bm{u}\in\mathcal{C}$. A function f is $G$-Lipschitz over a set
$\mathcal{C}$ if for all $\bm{u},\bm{w}\in\mathcal{C}$, we have
$|f(\bm{u})-f(\bm{w})|\leq G\|\bm{u}-\bm{w}\|_{2}$.

As proved in the study du2020learnability , the mean square error loss in Eq.
(2) is smooth and Lipschitz.

######  Lemma 1 (Lemma 2, du2020learnability ).

The loss function $\mathcal{L}$ in Eq. (2) is $S$-smooth with
$S=(3/2+\lambda)d^{2}$ and $G_{1}$-Lipschitz with $G_{1}=d(1+3\pi\lambda)$.

###  A.2 Proof details

The proof of Theorem 1 exploits the relation between the analytic and
estimated gradients of QNN, i.e.,
$\nabla_{j}\bar{\mathcal{L}}_{i}(\bm{\theta}_{i}^{(t,w)})$ and
$\nabla_{j}\mathcal{L}_{i}(\bm{\theta}_{i}^{(t,w)})$.

######  Lemma 2.

Denote $\tilde{p}=1-(1-p)^{L_{Q}}$ with $L_{Q}$ being the quantum circuit
depth. The discrepancy between the analytic gradients
$\nabla\mathcal{L}_{i}(\bm{\theta}_{i}^{(t,w)})$ and the estimated gradients
$\nabla\bar{\mathcal{L}}_{i}(\bm{\theta}_{i}^{(t,w)})$ for the $i$-th node
satisfies

$\displaystyle\mathbb{E}_{\bm{\varsigma}_{i,j}^{(t,w)}}\left[\left\|\nabla\bar{\mathcal{L}}_{i}(\bm{\theta}_{i}^{(t,w)})-\nabla\mathcal{L}_{i}\left(\bm{\theta}_{i}^{(t,w)}\right)\right\|^{2}\right]$
$\displaystyle\leq$
$\displaystyle(\tilde{p}-2)^{2}\tilde{p}^{2}\left\|\nabla\mathcal{L}_{i}\left(\bm{\theta}_{i}^{(t,w)}\right)\right\|^{2}+\frac{(1-\tilde{p})^{2}\tilde{p}^{2}d_{Q}}{4}+(2-\tilde{p})^{2}\tilde{p}^{2}G_{1}d_{Q}+\frac{7(1-\frac{\tilde{p}}{2})^{2}+\frac{1}{8}}{K}d_{Q},$
(18)

where the expectation is taking over the randomness of quantum system noise
and the measurement error.

The proof of the above lemma is given in Appendix A.3.

We are now ready to prove Theorem 1.

###### Proof of Theorem 1.

The machinery of QUDIO in Alg. 1 indicates that the difference between the
trainable parameters at the $(t+1)$-th and $t$-th global steps satisfies

$\bm{\theta}^{(t+1)}-\bm{\theta}^{(t)}=-\frac{\eta}{Q}\sum_{i=1}^{Q}\sum_{w=1}^{W}g_{i}^{(t,w)}~{}.$
(19)

Supported by $S$-smooth property of $\mathcal{L}(\bm{\theta})$ in Lemma 1, we
have

$\mathcal{L}\left(\bm{\theta}^{(t+1)}\right)-\mathcal{L}\left(\bm{\theta}^{(t)}\right)\leq\left\langle\nabla\mathcal{L}\left(\bm{\theta}^{(t)}\right),\bm{\theta}^{(t+1)}-\bm{\theta}^{(t)}\right\rangle+\frac{S}{2}\left\|\bm{\theta}^{(t+1)}-\bm{\theta}^{(t)}\right\|^{2}.$
(20)

Combining Eqs. (19) and (20), we obtain

$\displaystyle\mathbb{E}_{\bm{x}_{i}^{(t,w)}\sim\mathcal{D}_{i},\bm{\varsigma}^{(t,w)}_{i,j}}\left[\mathcal{L}\left(\bm{\theta}^{(t+1)}\right)-\mathcal{L}\left(\bm{\theta}^{(t)}\right)\right]$
$\displaystyle\leq$
$\displaystyle\mathbb{E}_{\bm{x}_{i}^{(t,w)}\sim\mathcal{D}_{i},\bm{\varsigma}^{(t,w)}_{i,j}}\left[\left\langle\nabla\mathcal{L}\left(\bm{\theta}^{(t)}\right),\bm{\theta}^{(t+1)}-\bm{\theta}^{(t)}\right\rangle\right]+\frac{S}{2}\mathbb{E}_{\bm{x}_{i}^{(t,w)}\sim\mathcal{D}_{i},\bm{\varsigma}^{(t,w)}_{i,j}}\left[\left\|\bm{\theta}^{(t+1)}-\bm{\theta}^{(t)}\right\|^{2}\right]$
$\displaystyle=$
$\displaystyle\underbrace{\mathbb{E}_{\bm{x}_{i}^{(t,w)}\sim\mathcal{D}_{i},\bm{\varsigma}^{(t,w)}_{i,j}}\left[\left\langle\nabla\mathcal{L}\left(\bm{\theta}^{(t)}\right),-\frac{\eta}{Q}\sum_{i=1}^{Q}\sum_{w=1}^{W}g_{i}^{(t,w)}\right\rangle\right]}_{\text{T2}}+\frac{S\eta^{2}}{2}\underbrace{\mathbb{E}_{\bm{x}_{i}^{(t,w)}\sim\mathcal{D}_{i},\bm{\varsigma}^{(t,w)}_{i,j}}\left[\left\|\frac{1}{Q}\sum_{i=1}^{Q}\sum_{w=1}^{W}g_{i}^{(t,w)}\right\|^{2}\right]}_{\text{T1}}.$
(21)

where the inequality uses the definition of the smooth function and the
equality is supported by Eq. (20).

We next derive the upper bounds of the terms T2 and T1 by first removing the
data sample noise (i.e., calculating
$\mathbb{E}_{\bm{x}_{i}^{(t,w)}\sim\mathcal{D}_{i}}[\cdot]$) and then removing
the system noise (i.e., calculating
$\mathbb{E}_{\bm{\varsigma}_{i,j}^{t,w}}[\cdot]$).

The elimination of the data sample noise. For the term T1, we rewrite
$g_{i}^{(t,w)}$ by $\nabla\bar{\mathcal{L}}_{i}(\bm{\theta}_{i}^{(t,w)})$ to
remove the randomness of sampling training examples. Mathematically, we have

$\displaystyle\mathbb{E}_{\bm{x}_{j}\sim\mathcal{D}_{i}}\left[\left\|\frac{1}{Q}\sum_{i=1}^{Q}\sum_{w=1}^{W}g_{i}^{(t,w)}\right\|^{2}\right]$
$\displaystyle\leq$
$\displaystyle\frac{W}{Q}\sum_{i=1}^{Q}\sum_{w=1}^{W}\mathbb{E}_{\bm{x}_{j}\sim\mathcal{D}_{i}}\left[\left\|g_{i}^{(t,w)}\right\|^{2}\right]$
$\displaystyle=$
$\displaystyle\frac{W}{Q}\sum_{i=1}^{Q}\sum_{w=1}^{W}\mathbb{E}_{\bm{x}_{j}\sim\mathcal{D}_{i}}\left[\left\|g_{i}^{(t,w)}-\nabla\bar{\mathcal{L}}_{i}\left(\bm{\theta}_{i}^{(t,w)}\right)+\nabla\bar{\mathcal{L}}_{i}\left(\bm{\theta}_{i}^{(t,w)}\right)-\nabla\mathcal{L}_{i}\left(\bm{\theta}_{i}^{(t,w)}\right)+\nabla\mathcal{L}_{i}\left(\bm{\theta}_{i}^{(t,w)}\right)\right\|^{2}\right]$
$\displaystyle\leq$
$\displaystyle\frac{2W}{Q}\sum_{i=1}^{Q}\sum_{w=1}^{W}\mathbb{E}_{\bm{x}_{j}\sim\mathcal{D}_{i}}\left[\left\|g_{i}^{(t,w)}-\nabla\bar{\mathcal{L}}_{i}\left(\bm{\theta}_{i}^{(t,w)}\right)+\nabla\bar{\mathcal{L}}_{i}\left(\bm{\theta}_{i}^{(t,w)}\right)-\nabla\mathcal{L}_{i}\left(\bm{\theta}_{i}^{(t,w)}\right)\right\|^{2}\right]+\frac{2W}{Q}\sum_{i=1}^{Q}\sum_{w=1}^{W}\mathbb{E}_{\bm{x}_{j}\sim\mathcal{D}_{i}}\left[\left\|\nabla\mathcal{L}_{i}\left(\bm{\theta}_{i}^{(t,w)}\right)\right\|^{2}\right]$
$\displaystyle=$
$\displaystyle\frac{4W}{Q}\sum_{i=1}^{Q}\sum_{w=1}^{W}\mathbb{E}_{\bm{x}_{j}\sim\mathcal{D}_{i}}\left[\left\|g_{i}^{(t,w)}-\nabla\bar{\mathcal{L}}_{i}\left(\bm{\theta}_{i}^{(t,w)}\right)\right\|^{2}\right]+\frac{4W}{Q}\sum_{i=1}^{Q}\sum_{w=1}^{W}\mathbb{E}_{\bm{x}_{j}\sim\mathcal{D}_{i}}\left[\left\|\nabla\bar{\mathcal{L}}_{i}\left(\bm{\theta}_{i}^{(t,w)}\right)-\nabla\mathcal{L}_{i}\left(\bm{\theta}_{i}^{(t,w)}\right)\right\|^{2}\right]$
$\displaystyle+\frac{2W}{Q}\sum_{i=1}^{Q}\sum_{w=1}^{W}\mathbb{E}_{\bm{x}_{j}\sim\mathcal{D}_{i}}\left[\left\|\nabla\mathcal{L}_{i}\left(\bm{\theta}_{i}^{(t,w)}\right)\right\|^{2}\right]$
$\displaystyle\leq$ $\displaystyle
4W^{2}\sigma^{2}+\frac{4W}{Q}\sum_{i=1}^{Q}\sum_{w=1}^{W}\underbrace{\left\|\nabla\bar{\mathcal{L}}_{i}\left(\bm{\theta}_{i}^{(t,w)}\right)-\nabla\mathcal{L}_{i}\left(\bm{\theta}_{i}^{(t,w)}\right)\right\|^{2}}_{\text{T3}}+2W^{2}G_{1}^{2},$
(22)

where the first inequality employs $\|\sum_{i=1}^{n}\bm{a}_{i}\|^{2}\leq
n\sum_{i=1}^{n}\|\bm{a}_{i}\|^{2}$, the second inequality uses the triangle
inequality with $\|\bm{a}+\bm{b}\|^{2}\leq 2\|\bm{a}\|^{2}+2\|\bm{b}\|^{2}$,
and the last inequality exploits the bounded variance of the estimated
gradients in the assumption.

The derivation of the upper bound of the term T2 is similar to the operations
applied to T1, i.e., the estimated gradient $g_{i}^{(t,w)}$ is substituted
with $\nabla\bar{\mathcal{L}}_{i}(\bm{\theta}_{i}^{(t,w)})$ to remove the
randomness of the sampled data. Mathemetically, we have

$\displaystyle\mathbb{E}_{\bm{x}_{j}\sim\mathcal{D}_{i}}\left[\left\langle\nabla\mathcal{L}\left(\bm{\theta}^{(t)}\right),-\frac{\eta}{Q}\sum_{i=1}^{Q}\sum_{w=1}^{W}g_{i}^{(t,w)}\right\rangle\right]$
$\displaystyle=$
$\displaystyle\left\langle\nabla\mathcal{L}\left(\bm{\theta}^{(t)}\right),\mathbb{E}_{\bm{x}_{i}}\left[-\frac{\eta}{Q}\sum_{i=1}^{Q}\sum_{w=1}^{W}g_{i}^{(t,w)}\right]\right\rangle$
$\displaystyle=$
$\displaystyle\left\langle\nabla\mathcal{L}\left(\bm{\theta}^{(t)}\right),-\frac{\eta}{Q}\sum_{i=1}^{Q}\sum_{w=1}^{W}\nabla\bar{\mathcal{L}}_{i}\left(\bm{\theta}_{i}^{(t,w)}\right)\right\rangle$
$\displaystyle=$
$\displaystyle-\frac{\eta}{2}\left\|\nabla\mathcal{L}\left(\bm{\theta}^{(t)}\right)\right\|^{2}-\frac{\eta}{2}\left\|\frac{1}{Q}\sum_{i=1}^{Q}\sum_{w=1}^{W}\nabla\bar{\mathcal{L}}_{i}\left(\bm{\theta}_{i}^{(t,w)}\right)\right\|^{2}+\frac{\eta}{2}\left\|\nabla\mathcal{L}\left(\bm{\theta}^{(t)}\right)-\frac{1}{Q}\sum_{i=1}^{Q}\sum_{w=1}^{W}\nabla\bar{\mathcal{L}}_{i}\left(\bm{\theta}_{i}^{(t,w)}\right)\right\|^{2}$
$\displaystyle\leq$
$\displaystyle-\frac{\eta}{2}\left\|\nabla\mathcal{L}\left(\bm{\theta}^{(t)}\right)\right\|^{2}+\frac{\eta}{2}\left\|\nabla\mathcal{L}\left(\bm{\theta}^{(t)}\right)-\frac{1}{Q}\sum_{i=1}^{Q}\sum_{w=1}^{W}\nabla\bar{\mathcal{L}}_{i}\left(\bm{\theta}_{i}^{(t,w)}\right)\right\|^{2}$
$\displaystyle\leq$
$\displaystyle-\frac{\eta}{2}\left\|\nabla\mathcal{L}\left(\bm{\theta}^{(t)}\right)\right\|^{2}+\frac{\eta}{2}\frac{W}{Q}\sum_{i=1}^{Q}\sum_{w=1}^{W}\underbrace{\left\|\nabla\mathcal{L}_{i}\left(\bm{\theta}_{i}^{(t,w)}\right)-\nabla\bar{\mathcal{L}}_{i}\left(\bm{\theta}_{i}^{(t,w)}\right)\right\|^{2}}_{T3},$
(23)

where the second equality uses the unbiased estimation property in terms of
the sampled data in Eq. (17), the third equality uses
$\langle\bm{a},\bm{b}\rangle=\frac{1}{2}(\|\bm{a}\|^{2}+\|\bm{b}\|^{2}-\|\bm{a}-\bm{b}\|^{2})$,
the second inequality employs the triangle inequality, and the last inequality
exploits the explicit form of the analytic gradient such that

$\nabla\mathcal{L}\left(\bm{\theta}^{(t)}\right)=\frac{1}{Q}\sum_{i=1}^{Q}\sum_{w=1}^{W}\nabla\mathcal{L}_{i}\left(\bm{\theta}^{(t,w)}\right),$
(24)

and the relation $\|\sum_{i=1}^{n}\bm{a}_{i}\|^{2}\leq
n\sum_{i=1}^{n}\|\bm{a}_{i}\|^{2}$.

In conjunction with Eqs. (A.2), (A.2), and (A.2), we obtain

$\displaystyle\mathbb{E}_{\bm{x}_{i}^{(t,w)}\sim\mathcal{D}_{i},\bm{\varsigma}^{(t,w)}_{i,j}}\left[\mathcal{L}\left(\bm{\theta}^{(t+1)}\right)-\mathcal{L}\left(\bm{\theta}^{(t)}\right)\right]$
$\displaystyle\leq$
$\displaystyle-\frac{\eta}{2}\left\|\nabla\mathcal{L}\left(\bm{\theta}^{(t)}\right)\right\|^{2}+\frac{\eta
W}{2Q}\sum_{i=1}^{Q}\sum_{w=1}^{W}\mathbb{E}_{\bm{\varsigma}^{(t,w)}_{i,j}}\underbrace{\left[\left\|\nabla\bar{\mathcal{L}}_{i}\left(\bm{\theta}_{i}^{(t,w)}\right)-\nabla\mathcal{L}_{i}\left(\bm{\theta}_{i}^{(t,w)}\right)\right\|^{2}\right]}_{T3}$
(25)
$\displaystyle+\frac{S\eta^{2}}{2}\left(4W^{2}\sigma^{2}+2W^{2}G_{1}^{2}\right)+\frac{S\eta^{2}}{2}\frac{4W}{Q}\sum_{i=1}^{Q}\sum_{w=1}^{W}\underbrace{\mathbb{E}_{\bm{\varsigma}^{(t,w)}_{i,j}}\left[\left\|\nabla\bar{\mathcal{L}}_{i}\left(\bm{\theta}_{i}^{(t,w)}\right)-\nabla\mathcal{L}_{i}\left(\bm{\theta}_{i}^{(t,w)}\right)\right\|^{2}\right]}_{\text{T3}}~{}.$

The elimination of the system and measurement noise. By employing Lemma 2, the
term T3 can be upper bounded by

$\displaystyle\mathbb{E}_{\bm{\varsigma}_{i,j}^{(t,w)}}\left[\left\|\nabla\bar{\mathcal{L}}_{i}\left(\bm{\theta}_{i}^{(t,w)}\right)-\nabla\mathcal{L}_{i}\left(\bm{\theta}_{i}^{(t,w)}\right)\right\|^{2}\right]$
$\displaystyle\leq$
$\displaystyle(\tilde{p}-2)^{2}\tilde{p}^{2}\left\|\nabla\mathcal{L}_{i}\left(\bm{\theta}_{i}^{(t,w)}\right)\right\|^{2}+\frac{(1-\tilde{p})^{2}\tilde{p}^{2}d_{Q}}{4}+(2-\tilde{p})^{2}\tilde{p}^{2}G_{1}d_{Q}+\frac{7(1-\frac{\tilde{p}}{2})^{2}+\frac{1}{8}}{K}d_{Q}.$
(26)

Then, supported by Eqs. (A.2) and (A.2), the upper bound in Eq. (20) satisfies

$\displaystyle\mathbb{E}_{\bm{x}_{j}\sim\mathcal{D}_{i},\bm{\varsigma}_{i,j}^{(t,w)}}\left[\mathcal{L}\left(\bm{\theta}^{(t+1)}\right)-\mathcal{L}\left(\bm{\theta}^{(t)}\right)\right]$
$\displaystyle\leq$
$\displaystyle-\frac{\eta}{2}\left\|\nabla\mathcal{L}\left(\bm{\theta}^{(t)}\right)\right\|^{2}+\frac{S\eta^{2}}{2}\left(4W^{2}\sigma^{2}+2W^{2}G_{1}^{2}\right)$
$\displaystyle+\left(\frac{S\eta^{2}}{2}\frac{4W}{Q}+\frac{\eta
W}{Q}\right)\sum_{i=1}^{Q}\sum_{w=1}^{W}\mathbb{E}_{\bm{x}_{j}\sim\mathcal{D}_{i},\bm{\varsigma}_{i,j}^{(t,w)}}\left[\underbrace{\left\|\nabla\bar{\mathcal{L}}_{i}\left(\bm{\theta}_{i}^{(t,w)}\right)-\nabla\mathcal{L}_{i}\left(\bm{\theta}_{i}^{(t,w)}\right)\right\|^{2}}_{\text{T3}}\right]$
$\displaystyle\leq$
$\displaystyle-\frac{\eta}{2}\left\|\nabla\mathcal{L}\left(\bm{\theta}^{(t)}\right)\right\|^{2}+\frac{S\eta^{2}}{2}\left(4W^{2}\sigma^{2}+2W^{2}G_{1}^{2}\right)$
(27) $\displaystyle+\left(\frac{S\eta^{2}4W^{2}}{2}+\eta
W^{2}\right)\left((\tilde{p}-2)^{2}\tilde{p}^{2}G_{1}^{2}+\frac{(1-\tilde{p})^{2}\tilde{p}^{2}d_{Q}}{4}+(2-\tilde{p})^{2}\tilde{p}^{2}G_{1}d_{Q}+\frac{7(1-\frac{\tilde{p}}{2})^{2}+\frac{1}{8}}{K}d_{Q}\right)~{}.$

Rearranging the terms in Eq. (A.2), the norm of the gradients
$\|\nabla\mathcal{L}(\bm{\theta}^{(t)}\|$ is upper bounded by

$\displaystyle\left\|\nabla\mathcal{L}\left(\bm{\theta}^{(t)}\right)\right\|^{2}\leq$
$\displaystyle\frac{2}{\eta}\mathbb{E}_{\bm{x}_{j}\sim\mathcal{D}_{i},\bm{\varsigma_{i}}}\left[\mathcal{L}\left(\bm{\theta}^{(t)}\right)-\mathcal{L}\left(\bm{\theta}^{(t+1)}\right)\right]+S\eta\left(4W^{2}\sigma^{2}+2W^{2}G_{1}^{2}\right)$
$\displaystyle+$ $\displaystyle\left(S\eta
4W^{2}+2W^{2}\right)\left((\tilde{p}-2)^{2}\tilde{p}^{2}G_{1}^{2}+\frac{(1-\tilde{p})^{2}\tilde{p}^{2}d_{Q}}{4}+(2-\tilde{p})^{2}\tilde{p}^{2}G_{1}d_{Q}+\frac{7(1-\frac{\tilde{p}}{2})^{2}+\frac{1}{8}}{K}d_{Q}\right).$
(28)

Summing over $t$ and dividing both sides by $T$ in Eq. (A.2), we obtain

$\displaystyle\frac{1}{T}\sum_{t=1}^{T}\left\|\nabla\mathcal{L}\left(\bm{\theta}^{(t)}\right)\right\|^{2}$
$\displaystyle\leq$
$\displaystyle\frac{2}{\eta}\mathbb{E}_{\bm{x}_{j}\sim\mathcal{D}_{i},\bm{\varsigma_{i}}}\left[\mathcal{L}\left(\bm{\theta}^{(1)}\right)-\mathcal{L}\left(\bm{\theta}^{(T+1)}\right)\right]+S\eta\left(4W^{2}\sigma^{2}+2W^{2}G_{1}^{2}\right)$
$\displaystyle+$ $\displaystyle\left(S\eta
4W^{2}+2W^{2}\right)\left((\tilde{p}-2)^{2}\tilde{p}^{2}G_{1}^{2}+\frac{(1-\tilde{p})^{2}\tilde{p}^{2}d_{Q}}{4}+(2-\tilde{p})^{2}\tilde{p}^{2}G_{1}d_{Q}+\frac{7(1-\frac{\tilde{p}}{2})^{2}+\frac{1}{8}}{K}d_{Q}\right)$
$\displaystyle\leq$ $\displaystyle\frac{2+80\lambda d_{Q}}{\eta
T}+S\eta\left(4W^{2}\sigma^{2}+2W^{2}G_{2}^{2}\right)$ (29)
$\displaystyle+\left(4S\eta
W^{2}+2W^{2}\right)\left((\tilde{p}-2)^{2}\tilde{p}^{2}G_{1}^{2}+\frac{(1-\tilde{p})^{2}\tilde{p}^{2}d_{Q}}{4}+(2-\tilde{p})^{2}\tilde{p}^{2}G_{1}d_{Q}+\frac{7(1-\frac{\tilde{p}}{2})^{2}+\frac{1}{8}}{K}d_{Q}\right),$

where the second inequality exploits the upper bound of the discrepancy of
loss function.

With setting $\eta=\sqrt{1/ST}$, we achieve

$\displaystyle\frac{1}{T}\sum_{t=1}^{T}\left\|\nabla\mathcal{L}\left(\bm{\theta}^{(t)}\right)\right\|^{2}$
$\displaystyle\leq$ $\displaystyle(2+80\lambda
d_{Q})\sqrt{\frac{S}{T}}+\sqrt{\frac{S}{T}}\left(4W^{2}\sigma^{2}+2W^{2}G_{2}^{2}\right)$
(30)
$\displaystyle+\left(4W^{2}\sqrt{\frac{S}{T}}+2W^{2}\right)\left((\tilde{p}-2)^{2}\tilde{p}^{2}G_{1}^{2}+\frac{(1-\tilde{p})^{2}\tilde{p}^{2}d_{Q}}{4}+(2-\tilde{p})^{2}\tilde{p}^{2}G_{1}d_{Q}+\frac{7(1-\frac{\tilde{p}}{2})^{2}+\frac{1}{8}}{K}d_{Q}\right).$

∎

###  A.3 Proof of Lemma 2

As shown in du2020learnability , the estimated gradients, which are caused by
the gates noise and the sample errors, can be explicitly formulated to relate
with its analytic gradients.

######  Lemma 3 (Modified from Theorem 3, du2020learnability ).

Denote $\tilde{p}=1-(1-p)^{L_{Q}}$ with $L_{Q}$ being the quantum circuit
depth. At the $(t,w)$-th iteration, we define five constants with

$C^{(t,w)}_{i,j,a}=\begin{cases}(1-\tilde{p})\tilde{p}(\frac{1}{2}-{y}_{i}^{(t,w)})(\hat{y}_{i}^{(t,w,+_{j})}-\hat{y}_{i}^{(t,w,-_{j})})-(2\tilde{p}-\tilde{p}^{2})\lambda\bm{\theta}_{j}^{(t,w)},&a=1\\\
(1-\tilde{p})\left(\hat{y}_{i}^{(t,w,+_{j})}-\hat{y}_{i}^{(t,w,-_{j})}\right),&a=2\\\
(1-\tilde{p})\hat{y}_{i}^{(t,w)}+\frac{\tilde{p}}{2}-y_{i}^{(t,w)},&a=3\\\
\frac{-(1-\tilde{p})(\hat{y}_{i}^{(t,w)})^{2}+(1-\tilde{p})^{2}\hat{y}_{i}^{(t,w)}+\frac{\tilde{p}}{2}-\frac{\tilde{p}^{2}}{4}}{K},&a=4\\\
\frac{-(1-\tilde{p})((\hat{y}_{i,+_{j}}^{(t,w)})^{2}+(\hat{y}_{i,-_{j}}^{(t,w)})^{2})+(1-\tilde{p})^{2}(\hat{y}_{i,+_{j}}^{(t,w)}+\hat{y}_{i,-_{j}}^{(t,w)})+\tilde{p}-\frac{\tilde{p}^{2}}{2}}{K},&a=5,\end{cases}$

where $\hat{y}_{i}^{(t,w,\pm_{j})}$, $\hat{y}_{i}^{(t,w)}$, and
$y_{i}^{(t,w)}$ are defined in Eq. (8), $K$ refers to the number of quantum
measurements.

The relation between the estimated and analytic gradients of QNN follows

$\nabla_{j}\bar{\mathcal{L}}_{i}(\bm{\theta}^{(t,w)})=(1-\tilde{p})^{2}\nabla_{j}{\mathcal{L}}_{i}(\bm{\theta}^{(t,w)})+C^{(t,w)}_{i,j,1}+\bm{\varsigma}_{i,j}^{(t,w)},\forall
j\in[d]$ (31)

with
$\bm{\varsigma}_{i,j}^{(t,w)}=C^{(t,w)}_{i,j,2}\xi_{i}^{(t,w)}+C^{(t,w)}_{i,j,3}\xi_{i,j}^{(t,w)}+\xi^{(t,w)}\xi_{i,j}^{(t,w)}$,
where $\xi_{i}^{(t,w)}$ and $\xi_{i,j}^{(t,w)}$ are two random variables with
zero mean and variances $C^{(t,w)}_{i,j,4}$ and $C^{(t,w)}_{i,j,5}$,
respectively.

This results can be used to obtain Lemma 2.

###### Proof of Lemma 2.

We first derive the term

$\displaystyle\mathbb{E}_{\bm{\varsigma}_{i,j}^{(t,w)}}\left[\left\|\nabla\bar{\mathcal{L}}_{i}\left(\bm{\theta}_{i}^{(t,w)}\right)-\nabla\mathcal{L}_{i}\left(\bm{\theta}_{i}^{(t,w)}\right)\right\|^{2}\right]$
$\displaystyle=$
$\displaystyle\sum_{j=1}^{d_{Q}}\mathbb{E}_{\bm{\varsigma}_{i,j}^{(t,w)}}\left[\left(\nabla_{j}\bar{\mathcal{L}}_{i}\left(\bm{\theta}_{i}^{(t,w)}\right)-\nabla_{j}\mathcal{L}_{i}\left(\bm{\theta}_{i}^{(t,w)}\right)\right)^{2}\right]$
$\displaystyle=$
$\displaystyle\sum_{j=1}^{d_{Q}}\mathbb{E}_{\bm{\varsigma}_{i,j}^{(t,w)}}\left[\left((\tilde{p}-2)\tilde{p}\nabla_{j}\mathcal{L}_{i}\left(\bm{\theta}_{i}^{(t,w)}\right)+C^{(t,w)}_{i,j,1}+\bm{\varsigma}_{i,j}^{(t,w)}\right)^{2}\right]$
$\displaystyle=$
$\displaystyle\sum_{j=1}^{d_{Q}}\left((\tilde{p}-2)\tilde{p}\nabla_{j}\mathcal{L}_{i}\left(\bm{\theta}_{i}^{(t,w)}\right)\right)^{2}+\left(C^{(t,w)}_{i,j,1}\right)^{2}+\mathbb{E}_{\bm{x}_{i}\sim\bm{\varsigma}_{i,j}^{(t,w)}}\left[\left(\bm{\varsigma}_{i,j}^{(t,w)}\right)^{2}\right]+2(\tilde{p}-2)\tilde{p}\nabla_{j}\mathcal{L}_{i}\left(\bm{\theta}_{i}^{(t,w)}\right)C^{(t,w)}_{i,j,1}$
$\displaystyle~{}~{}~{}~{}+2\left((\tilde{p}-2)\tilde{p}\nabla_{j}\mathcal{L}_{i}\left(\bm{\theta}_{i}^{(t,w)}\right)+C^{(t,w)}_{i,j,1}\right)\mathbb{E}_{\bm{x}_{i}\sim\bm{\varsigma}_{i,j}^{(t,w)}}\left[\bm{\varsigma}_{i,j}^{(t,w)}\right]$
$\displaystyle\leq$
$\displaystyle(\tilde{p}-2)^{2}\tilde{p}^{2}\left\|\nabla\mathcal{L}_{i}\left(\bm{\theta}_{i}^{(t,w)}\right)\right\|^{2}+\frac{(1-\tilde{p})^{2}\tilde{p}^{2}d_{Q}}{4}+\mathbb{E}_{\bm{x}_{i}\sim\bm{\varsigma}_{i,j}^{(t,w)}}\left[\left(\bm{\varsigma}_{i,j}^{(t,w)}\right)^{2}\right]+(2-\tilde{p})(1-\tilde{p})\tilde{p}^{2}G_{1}d_{Q}$
$\displaystyle=$
$\displaystyle(\tilde{p}-2)^{2}\tilde{p}^{2}\left\|\nabla\mathcal{L}_{i}\left(\bm{\theta}_{i}^{(t,w)}\right)\right\|^{2}+\frac{(1-\tilde{p})^{2}\tilde{p}^{2}d_{Q}}{4}+\sum_{j=1}^{d_{Q}}\left(C^{(t,w)}_{i,j,2}\right)^{2}C^{(t,w)}_{i,j,4}+\left(C^{(t,w)}_{i,j,3}\right)^{2}C^{(t,w)}_{i,j,5}+C^{(t,w)}_{i,j,4}C^{(t,w)}_{i,j,5}$
$\displaystyle+(2-\tilde{p})(1-\tilde{p})\tilde{p}^{2}G_{1}d_{Q}$
$\displaystyle\leq$
$\displaystyle(\tilde{p}-2)^{2}\tilde{p}^{2}\left\|\nabla\mathcal{L}_{i}\left(\bm{\theta}_{i}^{(t,w)}\right)\right\|^{2}+\frac{(1-\tilde{p})^{2}\tilde{p}^{2}d_{Q}}{4}+(1-\tilde{p})^{2}d_{Q}\frac{(1-\tilde{p})^{2}+1/4}{K}+\left(1-\frac{\tilde{p}}{2}\right)^{2}d_{Q}\frac{2(1-\tilde{p})^{2}+1/2}{K}$
$\displaystyle+\frac{(2(1-\tilde{p})^{2}+1/2)((1-\tilde{p})^{2}+1/4)d_{Q}}{K^{2}}+(2-\tilde{p})(1-\tilde{p})\tilde{p}^{2}G_{1}d_{Q}$
$\displaystyle\leq$
$\displaystyle(\tilde{p}-2)^{2}\tilde{p}^{2}\left\|\nabla\mathcal{L}_{i}\left(\bm{\theta}_{i}^{(t,w)}\right)\right\|^{2}+\frac{(1-\tilde{p})^{2}\tilde{p}^{2}d_{Q}}{4}+(2-\tilde{p})^{2}\tilde{p}^{2}G_{1}d_{Q}+\frac{7(1-\frac{\tilde{p}}{2})^{2}+\frac{1}{8}}{K}d_{Q}~{},$
(32)

where the first equality comes from the definition of the $l_{2}$ norm, the
second equality uses Eq. (31), the first inequalitu employs
$\mathbb{E}_{\bm{\varsigma}_{i,j}^{(t,w)}}[\bm{\varsigma}_{i,j}^{(t,w)}]=0$
and the Lipschitz property of $\mathcal{L}_{i}$ in Lemma 1, the second
inequality utilizes the upper bounds of $C^{(t,w)}_{i,j,a}$ in Lemma 3, and
the last inequality simplifies the factor $\tilde{p}$.

∎

##  Appendix B Numerical simulation details of QUDIO towards hand-written
digits image classification

This section presents the detailed description about applying QUDIO to
accomplish hand-written digits image classification tasks. First, in Appendix
B.1, the construction of local nodes is explained. Then, the hyper-parameter
configuration of QUDIO is introduced in Appendix B.2. Last, in Appendix B.3,
we provide complementary simulation results of QUDIO together with thorough
discussions.

###  B.1 The setup of local nodes

The distilled dataset. We sample $756$ digit images labeled with ‘0’ and ‘1’
from the MNIST handwritten digit database lecun1998mnist , where $256$ images
(examples) compose the training set and the rest $500$ images (examples) form
the test dataset. Once the training and test sets are collected, the data
preprocessing is applied, i.e., all examples are down-sampled to $8\times 8$
pixels followed by the vectorization and $l_{2}$ normalization.

Given the processed training set
$\mathcal{D}=\\{\bm{x}_{i},y_{i}\\}_{i=1}^{n}$ with $n=256$ and
$\bm{x}_{i}\in\mathbb{R}^{64}$, the central server partitions it into
$\\{\mathcal{D}_{i}\\}_{i=1}^{Q}$ and allocate them to $Q$ local nodes. For
example, when $Q=32$, we have $|\mathcal{D}_{i}|=8$ for $\forall i\in[32]$.

The implementation of local nodes. Due to the identical implementation of all
local nodes, here we mainly consider the setup of the $i$-th local node
$\mathcal{Q}_{i}$. The construction of PQCs of $\mathcal{Q}_{i}$, i.e.,
$h(\bm{\theta},O,\rho_{i})=\operatorname{Tr}(OU(\bm{\theta})\rho_{i}U(\bm{\theta})^{\dagger}),~{}\forall
i\in[n]$ in Eq. (1), is illustrated in Fig. 4. Particularly, the amplitude
encoding method is adopted to encode $\bm{x}_{i}\in\mathcal{D}_{i}$ to the
quantum state $\rho_{i}$, i.e.,

$\rho_{i}=\ket{\bm{x}_{i}}\bra{\bm{x}_{i}},~{}\text{and}~{}\ket{\bm{x}_{i}}=\sum_{j=1}^{64}\bm{x}_{i,j}\ket{j}.$
(33)

The prepared state $\rho_{i}$ is interacted with $U(\bm{\theta})$ implemented
by the hardware-efficient ansatz. Note that we set the block number in
$U(\bm{\theta})$ as $L=4$ and each block contains a single-qubit layer, i.e.,
$\otimes_{i=1}^{N}Rot(\bm{\theta}_{i})=\otimes_{i=1}^{N}\operatorname{R_{Z}}(\bm{\theta}_{i,1})\operatorname{R_{Y}}(\bm{\theta}_{i,2})\operatorname{R_{Z}}(\bm{\theta}_{i,3})$
and an entangled layer formed by CNOT gates, as highlighted in the dashed box
of Fig. 4. The quantum measurement operator $O$ is set as
$O=\mathbb{I}_{2^{5}}\otimes\ket{0}\bra{0}$.

The predicted label of $\bm{x}_{i}$ is assigned as ‘$0$’ if
$h(\bm{\theta},O,\rho_{i})\leq 0.5$; otherwise, $\bm{x}_{i}$ is classified as
‘$1$’. In the NISQ setting, $\bm{x}_{i}$ is predicted as ‘$0$’ if the sample
mean follows $\bar{y}_{i}\leq 0.5$; otherwise, $\bm{x}_{i}$ is classified as
‘$1$’.

Figure 4: The implementation of the local node $\mathcal{Q}_{i}$ in the binary
classification task. The amplitude embedding method is exploited to convert
the classical example $\bm{x}_{i}$ into the quantum state $\ket{\bm{x}_{i}}$.
The hardware-efficient ansatz is adopted to implement the trainable unitary
$U(\bm{\theta})$. The quantum observable is set as
$O=\mathbb{I}_{2^{5}}\otimes\ket{0}\bra{0}$.

###  B.2 Hyper-parameter settings

The source code of QUDIO. The realization of QUDIO is based on Pytorch
NEURIPS2019_9015  and its distributed communication package. To be more
concrete, we select the GLOO backend and ring all-reduce operations to achieve
the communication protocol between processes on CPU. Note that this
distributed optimization framework can be easily extended to coordinate
multiple quantum processors.

The classical optimizers. In QUDIO, the stochastic gradient descent optimizer
is adopted to update parameters. Its hyper-parameters settings are as follows.
The initial learning rate is set as $0.01$, the momentum factor is set as
$0.9$, and the decay-rate is set as $0.1$ every $40$ epochs. For the ideal
case, the gradient of each learnable parameter is calculated by back
propagation rumelhart1986learning . For the NISQ case, the depolarizing noise
and measurement error are introduced and the gradient is estimated by
parameter shift rule.

Hardware parameters. All simulation results in this study are completed by the
classical device with Intel(R) Xeon(R) Gold 6267C CPU @ 2.60GHz and 128 GB
memory.

(a)  (b)

Figure 5: QUDIO accelerates the training of QNNs. Each point refers to the
averaged result with five repetitions. The labels ‘N’ and ‘I’ represent the
noisy and ideal scenarios, respectively. (a) The speedup ratio of QUDIO in the
measure of a fixed number of training iterations $T$. (b) The speedup ratio of
QUDIO in the measure of the time to accuracy.

###  B.3 More simulation results

Here we demonstrate more simulation results about QUDIO that are omitted in
the main text. Specifically, we first evaluate the speedup ratio of QUDIO in a
fine-grained setting. Then, we investigate how the number of local nodes $Q$
effects the learning performance. Finally, we study how the system noise and
the number of measurements influence the performance of QUDIO. These
supplementary results facilitate us to better understand the capability of
QUDIO from different views.

Speedup analysis. The implementation detail of QUDIO in the exploration of the
speedup ratio is as follows. Both the number of local nodes and the number of
local iterations have six varied settings, i.e., $Q,W\in\\{1,2,4,8,16,32\\}$.
In the NISQ scenario, the depolarization rate and the number of measurement
are set as $p=0.0001$ and $K=100$, respectively.

Fig. 5(a) and Fig. 5(b) respectively depict the speedup ratio of QUDIO in the
measure of a fixed number of training iterations $T$ and the time to accuracy,
and Tab. 1 and Tab. 2 record the concrete running time under each setting.
Note that the first metric differs with the second one in the sense that the
latter requires QUDIO to surpass a threshold accuracy instead of a fixed $T$.
For both two metrics, QUDIO attains a sublinear speedup in terms of $Q$.
Besides, a larger number of local iterations $W$ generally promises a higher
speedup ratio. We also notice that there exists a manifest margin between the
noiseless and NISQ settings. This phenomenon is mainly caused by the opposite
role of the communication overhead, i.e., the communication overhead occupies
a large portion of the computational runtime in the noiseless case, while it
becomes negligible in the NISQ case.

Number of local nodes ($Q$) | NISQ | Running time (s)
---|---|---
1 | ✗ | 1742.21
1 | ✓ | 41974.36
2 | ✗ | 1324.12
2 | ✓ | 20992.67
4 | ✗ | 1020.48
4 | ✓ | 11756.49
8 | ✗ | 900.62
8 | ✓ | 6761.86
16 | ✗ | 856.47
16 | ✓ | 4002.11
32 | ✗ | 938.07
32 | ✓ | 2218.29
(a) Time to iteration

Number of local nodes ($Q$) | NISQ | Running time (s)
---|---|---
1 | ✗ | 104.53
1 | ✓ | 2518.46
2 | ✗ | 66.20
2 | ✓ | 1049.63
4 | ✗ | 40.81
4 | ✓ | 470.25
8 | ✗ | 36.02
8 | ✓ | 270.47
16 | ✗ | 34.25
16 | ✓ | 200.10
32 | ✗ | 28.14
32 | ✓ | 93.16
(b) Time to accuracy

Table 1: The runtime of QUDIO with respect to the varied number of local nodes
$Q$ and fixed $W=32$ local steps.

Number of local steps ($W$) | NISQ | Running time (s)
---|---|---
1 | ✗ | 906.21
1 | ✓ | 2833.21
2 | ✗ | 930.84
2 | ✓ | 2589.63
4 | ✗ | 870.69
4 | ✓ | 2253.47
8 | ✗ | 761.29
8 | ✓ | 2215.43
16 | ✗ | 872.54
16 | ✓ | 2354.78
32 | ✗ | 938.07
32 | ✓ | 2218.29
(a) Time to iteration

Number of local steps ($W$) | NISQ | Running time (s)
---|---|---
1 | ✗ | 54.37
1 | ✓ | 169.99
2 | ✗ | 55.85
2 | ✓ | 155.37
4 | ✗ | 52.24
4 | ✓ | 117.18
8 | ✗ | 45.67
8 | ✓ | 115.20
16 | ✗ | 43.62
16 | ✓ | 117.73
32 | ✗ | 28.14
32 | ✓ | 93.16
(b) Time to accuracy

Table 2: The runtime of QUDIO with respect to the varied number of local steps
$W$ and fixed $Q=32$ local nodes.

Accuracy analysis. We next turn to explore how the factors $Q$ and $W$
influence the final test accuracy. Analogous to the speedup analysis, the
number of local nodes and local iterations have six settings, i.e.,
$Q,W\in\\{1,2,4,8,16,32\\}$. Meanwhile, the depolarizing rate and the number
of measurement are set as $p=0.0001$ and $K=100$ respectively.

As shown in Fig. 6, for all settings, the test accuracy achieved by QUDIO is
above $95\%$. This observation reflects the robustness of QUDIO. Moreover,
when the number of local iterations $W$ is kept to be identical, QUDIO gains
the highest test accuracy with the setting $Q=16$. On the contrary, when the
number of local nodes $Q$ is kept to be identical, increasing $W$ subsumes to
a deteriorate test accuracy, which complies with the claim of Theorem 1.

Figure 6: Role of number of quantum processors and local iterations. The box
plot illustrate the final test accuracy corresponding to varied settings of
$Q$ and $W$.

The role of the system noise and the number of measurements. We end this
section by comprehending how the factors $p$ and $K$ influence the performance
of QUDIO. In particular, the depolarization noise rate $p$ scales from
$0.0001$ to $0.0512$, and the number of measurements ranges from $5$ to $100$.
The number of local nodes and local iterations is fixed to be $Q=16$ and
$W=2$, respectively.

Fig. 7 summarizes the simulation results. In particular, when $p<0.0064$, the
performance of QUDIO heavily depends on the number of measurements. For
example, the test accuracy is around $98\%$ with $M=100$, while it drops to
$87\%$ with $M=5$. When $p>0.0064$, both $p$ and $K$ determine the performance
of QUDIO. For example, for the setting $K=5$, the test accuracy of QUDIO with
$p=0.0512$ is reduced by $14\%$ than the setting $p=0.0256$ (i.e., from $77\%$
to $66\%$).

Figure 7: The test accuracy of QUDIO with varied depolarization rate and the
number of measurements. The label ‘$K=a$’ refers that the number of
measurements $K$ is set as $a$.

##  Appendix C Numerical simulation details of QUDIO towards the ground state
energy estimation of hydrogen molecule

This section provides elaboration about applying QUDIO to estimate the ground
state energy of hydrogen molecule tasks. First, the setup of local nodes and
the hyper-parameters settings are shown in Appendix C.1. Then, we provide
complementary simulation results of QUDIO together with thorough discussions
in Appendix C.2.

###  C.1 Implementation of local nodes and hyper-parameters setting

The implementation of QUDIO mainly follows the proposal kandala2017hardware .
Namely, the binary tree encoding method bravyi2002fermionic  is used to map
the hydrogen molecular Hamiltonian into a 4-qubit system, where
$H\in\mathbb{C}^{2^{4}\times 2^{4}}$ consists of $n=15$ local Hamiltonian
terms. In QUDIO, the central server partitions these $15$ local terms into $Q$
subgroups $\\{\mathcal{S}_{i}\\}$ and allocate them to $Q$ local nodes.

The realization of all local nodes follows the same routine. With this regard,
here we only discuss the realization of the node $\mathcal{Q}_{i}$. The input
quantum state is modified to $\rho_{0}=\ket{1100}\bra{1100}$. The
implementation of the ansatz $U(\bm{\theta})$ is shown in Fig. 8(a), which is
formed by $4$ trainable single-qubit gates followed by $3$ CNOT gates. The
prepared state is continuously operated with the observable $H_{S_{i}}$ to
proceed optimization.

Figure 8: The implementation of local nodes in QUDIO and simulation results.
(a) The circuit implementation in the local node. (b)-(c) The estimated
potential energy surface achieved by QUDIO when $Q=2$ and $Q=4$.

###  C.2 More simulation results

We conduct extensive numerical simulations to benchmark how the number of
local nodes $Q$ and the number of local iterations $W$ effect the performance
of QUDIO. Besides, we explore the trainability of QUDIO. Note that for all
settings, we fix $K=100$.

The role of $Q$. Fig. 8(b) and (c) show the potential energy surface estimated
by QUDIO with the different number of local nodes $Q$. Concisely, for both
$Q=2$ and $Q=4$, increasing the number of local steps $W$ incurs an enhanced
estimation error. Moreover, when QUDIO synchronizes trainable parameters at
every local iteration (e.g., $W=1$), the approximation error approaches to be
zero. These outcomes accord with the simulation results obtained in the main
text and the claim of Theorem 1.

Figure 9: Performance of QUDIO with different number of local updates $W$.

The role of $W$. Fig. 9 compares the potential energy surface estimated by
QUDIO with the different number of local iterations $W$ while the factor of
local nodes is set as $Q=\\{2,4,8\\}$. Specifically, when $W=1$, QUDIO
achieves the zero approximation error regardless of the number of local nodes
$Q$. By contrast, when $W=8$, the performance of QDUIO becomes inferior. In
conjunction with the results in Fig. 8, the factor $W$ determines the
performance of QUDIO in the NISQ setting, as shown in Theorem 1.

Figure 10: The estimated energy achieved by QUDIO during the training process.
The inter-atomic distance is fixed to be $0.3\mathrm{\AA}$.

The trainability of QUDIO. Fig. 10 indicates the estimated ground state energy
of QUDIO with respect to the number of global iterations. The three subplots
hint that the number of local iterations $W$ determines the trainability of
QUDIO. Concretely, a smaller number of local updates $W$ assures a faster
convergence.