Measurement optimization of variational quantum simulation by classical shadow and derandomization
Kouhei Nakaji1,4 , Suguru Endo2 , Yuichiro Matsuzaki1 , and Hideaki Hakoshima3
1 Device Technology Research Institute, National Institute of Advanced Industrial Science and Technology (AIST),1-1-1 Umezono,

Tsukuba, Ibaraki 305-8568, Japan.
2 NTT Computer and Data Science laboratories, NTT corporation, Musashino, Tokyo 180-8585, Japan
3 Center for Quantum Information and Quantum Biology,

Osaka University, 1-2 Machikaneyama, Toyonaka, Osaka 560-0043,

arXiv:2208.13934v3 [quant-ph] 20 Apr 2023

Japan.
4 Current address: Department of Computer Science, University of Toronto, Toronto, Ontario, Canada

Simulating large quantum systems is the ultimate goal of quantum computing. Variational
quantum simulation (VQS) gives us a tool to achieve the goal in near-term devices by distributing the computation load to both classical and quantum computers. However, as thesize of the quantum system becomes large, the execution of VQS becomes more and morechallenging. One of the most severe challenges is the drastic increase in the number of measurements; for example, the number of measurements tends to increase by the fourth power ofthe number of qubits in a quantum simulation with a chemical Hamiltonian. This work aims
to optimize measurements in VQS by recently proposed measurement optimization techniquessuch as classical shadow and derandomization. Even though previous literature shows thatthe measurement optimization techniques successfully reduce measurements in the variational
quantum optimization (VQO), how to apply them to VQS was unclear due to the gap between
VQO and VQS in measuring observables. In this paper, we bridge the gap by changing the wayof measuring observables in VQS and propose an algorithm to optimize measurements in VQS.
Our numerical experiment shows the validity of using our algorithm with quantum chemicalsystems. Importantly, we theoretically reveal the advantage of using shadow-based strategies,
e.g., classical shadow and derandomization, not only in VQS but also in VQO by calculating
the variance of an observable. To the best of our knowledge, the advantage of using shadowbased strategies in VQO was given only numerically in previous research; therefore, this paper
also gives significant implications for measurement optimization in general variational quantum
algorithms.

1 Introduction
With the advent of noisy intermediate scale quantum (NISQ) computers [1], the next milestone is to realize
their practical applications. A class of quantum algorithms called hybrid quantum-classical algorithms(see [2, 3]) came to attention because a load of computation can be distributed between quantum computers and classical computers; accordingly, the depth of the quantum circuit may be reduced. Amonghybrid algorithms, variational quantum algorithms (VQAs) [4–9] have been extensively studied for diverse
applications such as quantum computational chemistry [5, 6], quantum machine learning [10], quantum
metrology [11, 12], etc. In VQAs, parametrized quantum states are generated from shallow quantum circuits, and the parameters are iteratively updated with the optimization using classical computers. VQAs
can be broadly classified into two types: variational quantum optimization (VQO) [4–8] and variational
quantum simulation (VQS) [9]. VQO computes and optimizes a problem-specific cost function; for example, in the famous variational quantum eigensolver (VQE) [5, 6], the energy of parametrized quantum
states is evaluated by measuring Pauli operators, which constitutes the problem Hamiltonian in quantum
computers, and the minimization of the energy results in the ground state and its energy. Meanwhile, VQS
simulates the dynamics of the quantum system, e.g., quantum real and imaginary time evolution [9, 13]

Accepted in Quantum 2023-03-28, click title to verify. Published under CC-BY 4.0.

1

and open quantum system dynamics [14] by minimizing the distance between the variationally simulated
evolution and the exact evolution. The update rule of the parameters in VQS is derived from variational
principles for time evolution, such as McLachlan’s variational principle [15]. In both VQO and VQS, updates of the parameters are performed according to the output from quantum computers, more specifically,
expectation values of observables. Since the number of circuit runs is still restricted in the current NISQ
era, it is crucial to reduce the effect of shot noise and save the measurement counts for estimating eachobservable.
However, with the larger scale of the problem, reducing the effect of shot noise becomes more and morechallenging. That is because the number of measurements to realize a given shot noise linearly increases
with the number of the terms in the observable (here, we call ‘term’ by an observable evaluatable with aquantum circuit). For example, given the number of qubits as n, the number of terms scales as O(n4 ) in thechemical Hamiltonian, and therefore, the number of measurements also scales as O(n4 ). To mitigate thechallenge, several measurement optimization methods are proposed in the context of VQO [8, 16–33]. In
general, there are two types of measurement optimization methods: the grouping strategy [8,16,18,22–33]
and the shadow-based strategy [17–21]. Previous results numerically show that the applications of grouping
strategies [26, 34] and shadow-based strategies [17, 19, 35] to VQE of specific problems achieve faster or
better convergence in terms of the number of measurements though the performance of the optimization
depends on which we choose an optimizer. Let us briefly review those two types of strategies in thefollowing.
On the one hand, the grouping strategy is a method to efficiently estimate the cost function of VQO,
such as a molecular Hamiltonian, by grouping commuting terms in the Hamiltonian [8, 16, 18, 21–33]. We
can reduce the number of measurements by simultaneously measuring all the terms in each group. On theother hand, the shadow-based strategies construct a classical representation of a quantum state with asmall number of measurements [17–21]. The classical shadow, which is the original shadow-based strategy,
predicts M observables with O(log M ) measurements up to the additive error; the information-theoretic
arguments show that the order of O(log M ) measurements is asymptotically optimal [17]. However, theclassical shadow is inefficient for estimating particular observables since it typically contains unnecessary
measurement outcomes not used for estimation. Previous studies propose variations of classical shadow to
predict specific observables, such as many-body Hamiltonians, including the sum of Pauli terms [18–21].
One of its variations, derandomization, can efficiently predict such observables [19]. Derandomization isa deterministic protocol to minimize the confidence bound of Pauli observables starting from completely
randomized Pauli measurements, and it typically achieves better performance than the grouping strategy
and original classical shadow. These measurement protocols are experiment-friendly even for NISQ devices,
and a photonic system realized an experiment of these protocols [36].
Even though the above measurement optimization strategies mainly target VQO, measurement optimization is arguably more important in VQS than in VQO. In VQO, the evolution of parameters in eachstep does not need to be correct as far as the optimal value of the observable is found. Meanwhile, in thecase of VQS, we expect that the evolution of parameters in each step correctly describes the time evolution,and therefore, we need to efficiently reduce the shot noise by optimizing the measurements. Still, there
is no previous research that applies the above measurement optimization methods to VQS. As far as ourknowledge, that is because the above measurement optimization methods are mainly constructed for the
problem where their observables are written by the sum of Pauli terms. In VQS, we need to evaluate thereal/imaginary part of the expectation of operators [9, 13, 14], and cannot be straightforwardly written in
the form of estimating the sum of Pauli observables. Approaches to optimize measurements in VQS are
highly demanded.
The aim of this paper is to show the way of applying measurement optimization techniques to VQS.
We particularly focus on shadow-based strategies. Our contribution also exists in clarifying the advantage
of utilizing the shadow-based methods both in VQS and VQO from the theoretical aspect. As far as ourknowledge, the advantage of VQO was shown only numerically in previous research. In this paper, weshow the advantage by evaluating the variance of an observable. Let us specify the contribution of ourwork is summarized as follows:• We successfully simplify the quantum circuit for VQS by a circuit synthesis, so that measurement
optimization techniques (shadow-based strategies and grouping strategies) are applicable,• We theoretically show the advantage of using the shadow-based strategy in the VQS by evaluating

Accepted in Quantum 2023-03-28, click title to verify. Published under CC-BY 4.0.

2

the variance of the observable; the analysis is also applicable for assessing the merit of using the
shadow-based strategy in VQO,
• We show the validity of our algorithm and the theoretical discussion by numerical demonstration.
We show that the derandomization [19], one of the shadow-based strategies, greatly reduces thenumber of measurements compared to the original approach.
We note that even though we focus on the shadow-based methods in this paper, we can also utilize thegrouping strategies for measurement optimization in VQS due to our circuit synthesis. Particularly, therecent result [31] shows that they can outperform derandomization in VQE by optimizing parametersin the grouping so that the variance of the observable is minimized. Also, in [32, 33], authors provide
methods that may outperform the method in [31]. Thus, we expect better performance in VQS by using
those grouping methods than shadow-based methods. It is still worth mentioning that those grouping
methods need a reference quantum state for calculating the variance; the reference state is chosen to be aquantum state close to the ground state of the Hamiltonian in their work since they mainly target VQE.
However, in most cases of VQS, our goal is not to find the ground state of the Hamiltonian, and we need toinvestigate the choice of a reference state. We leave combining our methods with those grouping strategies
for future work.
The rest of the paper is organized as follows. In Section 2, we give a brief review of VQS. We also state
the challenge of VQS in the same section. We review the shadow-based strategy in Section 3; particularly,
we discuss the classical shadow [17] and the derandomization [19] in detail. Section 4 is dedicated toproposing our algorithm that applies the shadow-based strategy to VQS. The advantage of utilizing the
shadow-based strategy is discussed in the same section. In Section 5, we numerically demonstrate theperformance of our algorithm. Finally, we conclude our paper with some discussions in Section 6.

2 Variational quantum simulation without measurements optimization
In this section, we review the variational quantum simulation (VQS) proposed and studied in the previous
literature [9, 13, 14] and state a challenge of the algorithm. In Section 2.1, we give a general descriptionof the VQS based on [37]. In Section 2.2, we show the variational real and imaginary time evolution as
examples. Section 2.3 is dedicated to stating a problem of the previous approach in the VQS.

2.1 General processes
We define the quantum time evolution for a state |v(t)i in this paper as the following differential equation,B(t)d|v(t)i = A(t) |v(t)i .
dt

(2.1)

E
~
where A(t) and B(t) denote matrices. In the VQS, we parameterize |v(t)i as v(θ(t))
by using a set ofparameters. From Mclachlan’s principle [15],δ B(t)~
NP ∂ v(θ(t))
X
k=1

∂θk

E
~
θ˙k − A(t) v(θ(t))

E

= 0,

(2.2)

where |||ψi|| ≡ hψ|ψi, NP is the number of parameters, and
δf ≡

NP
X
∂f
j=1

∂ θ̇j

δ θ˙j .

(2.3)

We note that there are two other conventional variational principles than the Mclachlan’s principle [15]:Dirac and Frenkel’s variational principles [38, 39] and the time-dependent variational principle [40, 41]. A
detailed comparison of different variational principals is made in [14]. The Dirac and Frenkel variational
principle is not suitable for VQS since the equation of the parameters may involve complex solutions,

Accepted in Quantum 2023-03-28, click title to verify. Published under CC-BY 4.0.

3

which contradicts the requirement that parameters are real. Even though the time-dependent variational
principle gives a real solution, the time-dependent variational principle is shown to be more unstable andnot applicable to the evolution of density matrices and imaginary time evolution. In contrast, McLachlan’s
principle generally produces stable solutions, and it is also applicable to all VQS problems.
Equation (2.2) is equivalent to
E
E
D
 D
~
~
~
~
NP
∂ v(θ(t))∂ v(θ(t))∂ v(θ(t))∂ v(θ(t))
X θ̇`B † (t)B(t)+B † (t)B(t)
∂θk
∂θ`
∂θ`
∂θk
`=1
E
 D
~
~
E D
∂ v(θ(t))∂ v(θ(t))~~
 = 0.
A† (t)B(t)+ v(θ(t))(2.4)
B † (t)A(t) v(θ(t))
−
∂θk
∂θk
for all k. We can further simplify (2.4) asNP
X

Mk` θ̇` = Vk ,

(2.5)

`=1

where
E
 D
~
~
∂ v(θ(t))∂ v(θ(t))
,Mk` = < B † (t)B(t)
∂θk
∂θ`
 D

~
E
∂ v(θ(t))~
Vk = < B † (t)A(t) v(θ(t))
,
∂θk

(2.6)

(2.7)

with <[w] as the real part of the complex number w.
Once we compute Mk` and Vk , we can update the parameter by the Euler method with a small
parameter δt as~ + δt) = θ(t)~ + M −1 V · δt,
θ(t
(2.8)
where M (V ) is the matrix (vector) representation of Mk` (Vk ). Next, we show how to compute Mk` and
Vk in the following.
Estimation of Mk` and Vk using quantum circuits
Let R be the unitary operator corresponding to the parameterized quantum circuit (PQC); namelyE~v(θ(t))= R|vref i,(2.9)
with |vref i as a reference state. We often call R the ansatz. We expand R as
R ≡ RNP (θNP )RNP −1 (θNP −1 ) · · · R1 (θ1 ),(2.10)

where Rk is the k-th unitary operator applied to the reference state and θk is the parameter of Rk . Let
us describe
Ng∂Rk (θk ) X=
gk,i Rk (θk )Pk,i ,(2.11)∂θki=1where Ng is the number of Pauli observables when expanding the derivative, gk,i is a coefficient, and Pk,iis a tensor product of Pauli operators and identity: Pk,i ∈ {X, Y, Z, I}⊗n . ThenE
Ng~
∂ v(θ(t))
X=
gk,i Rk,i |vref i,(2.12)∂θki=1

Accepted in Quantum 2023-03-28, click title to verify. Published under CC-BY 4.0.

4whereRk,i = RNP (θNP )RNP −1 (θNP −1 ) · · · Rk (θk )Pk,i Rk−1 (θk−1 ) · · · R1 (θ1 ).

(2.13)

In typical settings, Ng = O(1) is satisfied. For example, if we embed each parameter as exp(−iθk P/2),
where P denotes a tensor products of Pauli operators, Ng = 1 and gk,1 = −i/2. Also, let us describe†

B (t)B(t) =NBBX

β q Pq ,(2.14)

αr Pr ,

(2.15)

q=1†

B (t)A(t) =

N
BA
X
r=1

where Pq and Pr are tensor products of Pauli operators and identity, and where NBB and NBA are thenumbers of terms when decomposing B † (t)B(t) and B † (t)A(t) as the sums of Pauli operators.
By substituting (2.12) and (2.14) into (2.6), we obtainMk` =Ng NBBX
X

hi†∗< gk,ig`,i0 βq hvref | Rk,i
Pq R`,i0 |vref i .

(2.16)

i,i0 q=1

Similarly, by substituting (2.12) and (2.15) into (2.7), we obtainVk =Ng NBAX
X

hi†∗
αr hvref | Rk,i< gk,i
Pr R |vref i .

(2.17)

i=1 r=1hi†
For computing (2.16), we need to evaluate each < eiφ hvref | Rk,i
Pq R`,i0 |vref i with eiφ as a phase, whichis possible by using the Hadamard test. In Fig. 1, we show the quantum circuit to compute the value fork ≤ `. Note that sincehMk` = M`k , we do notineed to consider the case when k > `. To compute (2.17),†
we need to evaluate < eiφ hvref | Rk,i
Pr R |vref i , which is possible by using the quantum circuit in Fig. 2.

√(|0i + eiφ |1i)/ 2•|vref i

•

X

...
Rk−1

R1

•

X

...
Pk,i

R`−1

Rk

H

...
P`,i0

R`

RP

Pq

†
Figure 1: A quantum circuit to compute < eiφ hvref | Rk,i
Pq R`,i0 |vref i for k ≤ `. We use the circuit enclosed by thedotted rectangle
inSection
4.
We
note
that
the
outputquantum
state after applying the gates in the dotted box is√φ2 R`,i0 |0i|vref i + eiφ Rk,i |1i|vref i .
|ψk,i;`,i0 i = 1/√(|0i + eiφ |1i)/ 2



•

X

|vref i

R1

•

X

...

H

...
Rk−1

Pk,i

Rk

RP

Pr

†
Pr R |vref i . As in Figure 1, we use the circuit enclosed by
Figure 2: A quantum circuit to compute < eiφ hvref | Rk,ithe dotted rectangle
inSection
4.
We
note
that
the
outputquantum state after applying the gates in the dotted box is
√
φ
|ψk,i
i = 1/ 2(|0iRk,i |vref i + |1iR|vref i).



Accepted in Quantum 2023-03-28, click title to verify. Published under CC-BY 4.0.

5

2.2 Variational real and imaginary time evolution
Here we show Mk` and Vk in real-time evolution (RTE) and imaginary time evolution (ITE) as examplesof the VQS. The real-time evolution is an algorithm to simulate the Schrödinger equation under a timeevolution Hamiltonian H:∂|v(t)i= −iH|v(t)i.
(2.18)∂tThus, A(t) = −iH and B(t) = 1 and we obtain
Mk` =

Ng
X

h
i
†
∗
< gk,ig`,i0 hvref |Rk,iR`,i0 |vref i ,(2.19)

i,i0

Vk =Ng NBAX
X

hi†∗< −igk,i
αr hvref | Rk,i
Pr R |vref i ,(2.20)

i=1 r=1

PNBA
where we expand the Hamiltonian as H = j=1
αr Pr with αr as a real coefficient. Note that we see
NBB = 1 from (2.19).
The imaginary time evolution is the operation obtained by replacing it to τ and H to H − hHi in(2.18), where hHi = hv(t)|H|v(t)i. Namely,∂ |v(τ )i= −(H − hHi) |v(τ )i .
∂τ

(2.21)

Note that the term hHi is added in order to preserve the norm of the quantum state. Importantly, we canshow that, in a large τ limit, |v(τ )i becomes the ground state of the Hamiltonian H. Thus, the imaginary
time evolution is often used as a way of searching for the ground state. The functions A(t) and B(t) are
given by A(t) = −(H − hHi) and B(t) = 1. The term hHi does not affect the computation of Vk , which
can be checked by substituting the values of A(t) and B(t) to (2.7) asE~D∂ v(θ(t))~
,Vk = <  v(θ(t))− (H − hHi)
(2.22)
∂θk
E
D
E

~
~
~
D
∂ v(θ(t))
∂
v(
θ(t))
v(
θ(t))
~
 + 1 hHi
= −<  v(θ(t))
H.
(2.23)∂θk2∂t
E~
The second term, which corresponds to hHi, vanishes since the norm of v(θ(t))
preserves under the~transformation of θ(t),and therefore, we obtain Mk` and Vk by setting B(t) = 1 and A(t) = −H when
computing (2.16) and (2.17). Thus, in the variational imaginary time evolution, we obtain
Mk` =

Ng
X

h
i
†
∗
< gk,ig`,i0 hvref |Rk,iR`,i0 |vref i ,(2.24)

i,i0

Vk = −

Ng NBAX
X

hi†∗< gk,iαr hvref | Rk,i
Pr R |vref i .

(2.25)

i=1 r=1

The circuit used for evaluating each term of Vk in real and imaginary time evolution is Fig. 2, which we
introduced in Section 2.1. On the other hand, a simpler circuit is used for evaluating Mk` in real andimaginary time evolution;
we show the icircuit in Fig. 3. We can readily check that the output of theh†iφcircuit returns < e hvref | Rk,iR`,i0 |vref i .

2.3 Problems of the current approach
The approach mentioned above to performing the VQS is challenging to implement because the numberof measurements linearly increases as the number of terms denoted by NBB and NBA increases. That

Accepted in Quantum 2023-03-28, click title to verify. Published under CC-BY 4.0.

6√(|0i + eiφ |1i)/ 2•

X|vref i

R1

•

X

...

H

...
Rk−1

Pk,i

R`−1

Rk

P`,i0

†
Figure 3: A quantum circuit to compute < eiφ hvref | Rk,iR`,i0 |vref i for k ≤ `.

is because we need to evaluate the circuit corresponding to each term inside the summations of (2.16)
and (2.17) one by one (the number of measurements also linearly increases according to the value of Ng ,but it is not problematic since Ng is typically taken to be small). For example, in RTE and ITE, thenumber of measurements for evaluating Vk linearly scales with the number of terms in H, which could be
impractical when H has many terms. Suppose that we use a molecular Hamiltonian as the time-evolution
Hamiltonian. In that case, the number of terms in the Hamiltonian is O(n4 ) with n as the number of
qubits, and therefore, the number of measurements also increases in proportion to O(n4 ).
The challenge of a large number of Pauli observables in Hamiltonian is also addressed in the variational
quantum eigensolver (VQE), and many solutions have been proposed [6, 16–22, 26, 27, 30, 42–47]. However,to our best knowledge, no previous work has explored this issue in the context of the VQS. Importantly,
there is a crucial difference between VQS and VQE; namely, to evaluate each term in the summation in(2.14) and (2.15), we need the controlled operations associated with each Pauli observable to evaluateeach term in the summation of (2.14) and (2.15) in VQS, while the VQE only requires estimating the
expectation value of each Pauli observable. Thus, it is not straightforward to apply the methods discussed
in [6, 16–22, 26, 27, 30, 42–47] to VQS. In the rest of the paper, we show a way to overcome the difficulties
and discuss how to apply the proposed methods for the VQE to the VQS. Before discussing the main
point, we will review the proposed methods for the VQE in the next section.
It should be noted that in the case of ITE, the parameter shift-rule [48] is applicable to calculate Vk
under a suitable condition (e.g., Ng = 1 and gk,1 6=P
0 is real), meaning that the evaluation of Vk boils down
to the evaluation of the expectation value of H = r αr Pr with respect to a quantum state generated by
a PQC. In that case, we can reduce the number of measurements by utilizing the techniques proposed in
VQE. Still, in most of the applications of VQS, the parameter-shift rule is not applicable; therefore, weneed to construct a measurement optimization procedure tailored for VQS.

3 Measurements optimization for sum of Pauli observables
In this section, we review strategies toPreduce the number of measurements in the VQE, when the HamilK
tonian has many terms such as H = j=1 aj Pj for K  1. In the VQE, we need to estimate the valuePK
hHi = j=1 aj Tr(Pj ρ) where {aj } are real coefficients, {Pj } are Pauli observables (Pj ∈ {X, Y, Z, I}⊗n ),and ρ is a n-qubit quantum state generated by a quantum circuit. When K increases, the number ofmeasurements to achieve a given precision linearly increases as far as we naively evaluate each term.
For improving the scaling of the number of measurements, several solutions have been proposed [6,
16–22, 26, 27, 30, 42–47]. The most promising solutions are arguably shadow-based strategies [17–21],
which we review in this section. In Section 3.1, we first give a general discussion about shadow-based
strategies [17–21]. Next, in Section 3.2, we show the detail of some of the shadow-based strategies that
we will examine in our numerical experiment in Section 5.

3.1 General discussion for shadow-based strategies
Shadow-based strategies for estimating H [17–21] are mostly composed of three steps. The first step ofK
the estimation process takes the number of measurements Nshot , {aj }K
j=1 and {Pj }j=1 as its inputs. Let⊗nus denote a measurement basis by one of the Pauli observables {X, Y, Z} , where each measurement isperformed on the basis where the Pauli observable is diagonalized. For example, when n = 3, the measurement basis can be expressed as XZX, meaning that the quantum state ρ is measured in the computational

Accepted in Quantum 2023-03-28, click title to verify. Published under CC-BY 4.0.

7basis after we operate the Hadamard gates to the first and the third qubits. The first step returns the arrayshotof measurement basis as {Mr }N
r=1 ; how to select those basis depends on the algorithm. We define thefunction corresponding to the first step as buildMeasurements({aj }, {Pj }, Nshot ). The second step takes
Nshot
⊗n
shot
{Mr }r=1
and ρ as its inputs and returns the measurement results of ρ as {br }N.
r=1 , where br ∈ {1, −1}
shot
,
ρ).
We denote the function corresponding to the second step as getMeasurementResults({Mr }N
r=1
Nshot
K
K
shot
The third step takes {Mr }N
r=1 , {br }r=1 , {Pj }j=1 , and {aj }j=1 as its inputs, and returns the following
ν:
Nshot
1 X
ν=
νr ,
(3.1)
Nshot r=1
where
νr =

K
X

aj

j=1

1
f (Pj , Mr )µ(Pj , br ).
q(Pj )(3.2)

Here, f (P, M ) denotes a covering function, q(Pj ) denotes a covering probability function, and µ(P, b)
denotes an estimation function. These are defined so that ν should be an unbiased estimator of Tr(Hρ).
We show the detailed definitions as follows.
Covering function f (P, M )The covering function f (P, M ) takes a Pauli obserable P and a measurement basis M as its inputs, andreturns outputs as follows
f (P, M ) =nYflocal (P [i], M [i]),

(3.3)

i=1

(
flocal (P [i], M [i]) =1
0

if P [i] = M [i] or P [i] = I,otherwise(3.4)

where P [i] and M [i] is the Pauli observable corresponds to the i-th qubit of P and M . For example,
f (XY I, XY X) = 1 and f (ZIX, ZXY ) = 0. We call M covers P if and only if f (P, M ) = 1 as in [19].
Covering probability function
We define the covering probability function q(P ) depending on whether we generate the measurement basis
probabilistically or not. When they are probabilistically generated according to a probability distribution
qm (M ), q(P ) is the probability that the Pauli observable P is covered, namelyXq(P ) =
qm (M )f (P, M ).
(3.5)
M ∈{X,Y,Z}⊗n

When measurements are chosen deterministically, we define the covering probability function by the ratio
shotthat P is covered when measurements are {Mr0 }N, namely,r0q(P ) =

1
Nshot

N
shot
X

f (P, Mr0 ).

(3.6)

r 0 =1shot
Note that in this deterministic case, q(P ) depends on the chosen measurement basis {Mr0 }N
r 0 =1 , but weomit the corresponding index from q(P ) for simplicity.

Accepted in Quantum 2023-03-28, click title to verify. Published under CC-BY 4.0.

8

Estimation function
The second function is the estimation function µ(P, b) that takes the Pauli observable and the bit arrayb ∈ {1, −1}⊗n as its inputs and operates as follows:µ(P, b) =nYµlocal (P [i], b[i]),

(3.7)

i=1

(
µlocal (P [i], b[i]) =1−1

P [i] = I or b[i] = 1,
otherwise(3.8)

where b[i] is the i-th bit value. Let p(ρ, M ; b) be the probability that n-bit bitarray b is obtained as aresult of measuring ρ in a basis M covering P . Then, we can readily showX
p(ρ, M ; b)µ(P, b) = Tr(P ρ).
(3.9)
b

As we show in Appendix A, ν becomes an unbiased estimator of Tr(Hρ), i.e., E[ν] = Tr(Hρ). Also,
the variance of ν is given by


K
K X
2
X
X
a
1
j

+
aj a` g(Pj , P` )Tr(Pj P` ρ) − [Tr(Hρ)]2  , (3.10)
Var(ν) = E[ν 2 ] − (E[ν])2 =Nshot j=1 q(Pj ) j=1
k6=j

where if Mr is sampled from qm (Mr ),X

g(Pj , P` ) =Mr ∈{X,Y,Z}f (Pj , Mr ) f (P` , Mr )qm (Mr ),q(Pj )q(P` )⊗n(3.11)

shotand if {Mr }N
r=1 is deterministically chosen,g(Pj , P` ) =

1

N
shot
X

Nshot r=1

f (Pj , Mr ) f (P` , Mr ).
q(Pj )q(P` )(3.12)

K
We write the three steps of obtaining ν as estimateNu(ρ, {aj }K
j=1 , {Pj }j=1 , Nshot ) and summarize it in
Algorithm 1. It should be noted that the variance Var(ν) determines the number of measurements Nshot
to achieve a given precision ; namely, Nshot scales as Nshot = O(Var(ν)/2 ).

Algorithm 1 estimateNu(ρ, {aj }, {Pj }, Nshot )Input: ρ, {aj }, {Pj }, Nshot
Nshot
1: {Mr }r=1
= buildMeasurements({aj }, {Pj }, Nshot )Nshot
2: Set {br }r=1 = getMeasurementResults({Mr }, ρ)
3: for r = 1...Nshot do
4:

νr =

K
X

aj

j=1

1
f (Pj , Mr )µ(Pj , br ).
q(Pj )

5: end for
PNshot1
return ν = Nshot
r=1 νr .

3.2 Examples of shadow-based strategies
The shadow-based strategy is originally proposed by the literature [17] and various methods motivated by
the technique have been proposed [17–21], which mainly target to reduce the value of Var(ν) in (3.10).
Among them, we review the classical shadow and the derandomization in the following discussion.

Accepted in Quantum 2023-03-28, click title to verify. Published under CC-BY 4.0.

9

Algorithm 2 buildMeasurements({aj }, {Pj }, Nshot ) in Classical shadow
Input: {aj }, {Pj }, Nshot
1: SET n as the number of qubits of the operators {Pj }.
2: SET results = [].
3: for r = 1...Nshot do
4:
SET M = [].
5:
for i = 1...n do
6:Choose one of {X, Y, Z} with equal probability and set it to M [i].
7:
end for
8:
SET results[r] = M .
9: end for
return results
Classical shadow
In the classical shadow, we randomly choose the measurement basis. Namely, each Mr [i] (i-th Pauli matrix
in r-th measurement basis) is chosen from {X, Y, Z} with equal probability, i.e., qm (M ) = 1/3n . We show
the pseudo-code for buildMeasurements(Nshot ) in classical shadow in Algorithm 2. We can readily
show that the covering probability is given byq(Pj ) =X
M ∈{X,Y,Z}

1
f (Pj , M )= locality(P ) ,
n
j
3
3
⊗n

(3.13)

where locality(O) is the number of non-identity Pauli matrix in the Pauli observable O. As we show in
Appendix B, we obtaing(Pj , P` ) = gcs (Pj , P` ) ≡nYglocal (Pj [i], P` [i]),(3.14)

Pj [i] = I or P` [i] = I.
Pj [i] = P` [i] 6= Iotherwise(3.15)

i=1

where



 1
glocal (Pj [i], P` [i]) =30Therefore,
Var(ν) =

1
Nshot



K
X
X

a2j 3locality(Pj ) +aj a` gcs (Pj , P` )Tr(Pj P` ρ) − [Tr(Hρ)]2  .
j=1(3.16)

j6=`

Derandomization
As we see in (3.10), the value of Var(ν) tends to be small if q(Pj ) is large. From this point of view,the random choice of measurements in the classical shadow is not always efficient; particularly, when thelocality of a Pauli term Pj is large, the value of q(Pj ) is exponentially suppressed, as we see in (3.13).
Thus, strategies to choose the measurement basis set that more efficiently cover Pauli observables areproposed as consecutive works [18–21].
One of such strategies is the derandomization [19]. Unlike the classical shadow, a pre-defined measurement basis set is used in the derandomization (thus, it corresponds to the deterministic case). For
shotobtaining the pre-defined measurement basis {Mr }N
r=1 , the procedure of the derandomization iteratively
determines them in the following order:
M1 [1], M1 [2], · · · M1 [n], M2 [1], M2 [2], · · · M2 [n], · · · MNshot [1], MNshot [2], · · · MNshot [n].

(3.17)

000
Suppose that {Mr }rr=1
and {Mr0 +1 [i]}ii=1
are already determined and let M be the concatinate of {Mr }rr=1i0
and {Mr0 +1 [i]}i=1 . Then, Mr0 +1 [i0 +1] ∈ {X, Y, Z} is determined so that a score function C(M, Mr0 +1 [i0 +Accepted in Quantum 2023-03-28, click title to verify. Published under CC-BY 4.0.

101]) is maximized (here we assume i0 < n for simplicity, but we can give the same argument by replacing
Mr0 +1 [i0 + 1] with Mr0 +2 [1] even when i0 = n).
The score function C is constructed so that the increase of C leads to the increase of {q(Pj )}K
j=1 in
our notation. Let
!
r0
iY
0 +1
ηX
f (Pj , Mr ) − log 1 − βi0 +1
flocal (Pj [k], Mr0 +1 [k])V (Pj , M, Mr0 +1 [i0 + 1]) =(3.18)
2 r=1
k=1

with
βi0 +1 = γ

K
Y

1
,
locality(Pj )3
j=i +2(3.19)

0where γ, η are hyperparameters. In [19], they give some options for the choice of the score function C;
among them, the one used in their numerical experiment isC(M, Mr0 +1 [i0 + 1]) =KX

exp (−V (Pj , M, Mr0 +1 [i0 + 1]) /wj ) ,(3.20)

j=1

where
wj =

|aj |
.
maxj |aj |(3.21)For example, suppose that the observable is
H = XXXZ + XXII + IIXZ + Y Y ZX + Y Y II + IIZX

(3.22)and Nshot = 10. By minimizing (3.20) step by step, {Mr }10
r=1 includes five XXXZ and Y Y ZX, whichleads to1
= 2.
(3.23)q(Pr )In the classical shadow,1≥ 32 .
q(Pr )(3.24)

Thus, we see that the covering probability improves in the derandomization compared to the classicalshadow.

4 Measurements optimization in the VQS
In this section, we propose the VQS with shadow-based strategies reviewed in the previous section. In
Section 4.1, we synthesize the circuits in Fig. 1 and Fig. 2 and remove unnecessary two-qubit operations.
Based on the result of the circuit synthesis, we propose an algorithm applying shadow-based strategies
to the VQS in Section 4.2. In Section 4.3, we discuss the advantage of using our proposed algorithm.
Note that we generally discuss our proposed algorithm in this section and demonstrate its performanceby using RTE and ITE in the numerical demonstration in Section 5. Section 4.4 is dedicated to apply our
implication obtained in Section 4.3 to the general shadow-based strategies discussed in Section 3

4.1 Circuit synthesis
Let us define ρφk,i;`,i0 as the quantum state generated by operating the partial circuits enclosed by the√dotted rectangle in Fig. 1 to the reference state: (|0i + eiφ |1i)/ 2 ⊗ |vref i. The specific form of ρφk,i;`,i0 is
given byφφ
ρφk,i;`,i0 = |ψk,i,`;i
0 ihψk,i;`,i0 |,1φR`,i0 |0i|vref i + eiφ Rk,i |1i|vref i .
|ψk,i;`,i0i = √
2

Accepted in Quantum 2023-03-28, click title to verify. Published under CC-BY 4.0.

(4.1)
(4.2)

11

Let CPq also be the controlled-Pq operation and In be the n-qubit identity operator. Then, for each q, k,and `, the output of the circuit shown in Fig. 1 has the following form,
h
i
h
i
†
< eiφ hvref | Rk,i
Pq R`,i0 |vref i = Tr (X ⊗ In )(CPq )ρφk,i;`,i0 (CP†q )
(4.3)


= Tr Oq ρφk,i;`,i0 ,
(4.4)
with
Oq ≡ (CP†q )(X ⊗ In )(CPq ).

(4.5)

Similarly, let ρφk,i be the quantum state generated by operating the partial circuits enclosed by the dottedrectangle in Fig. 2 to the reference state. The specific form of ρφk,i is given byφφ
ρφk,i = |ψk,iihψk,i|1
φ|ψk,ii = √ (|0iRk,i |vref i + |1iR|vref i).
2

(4.6)
(4.7)

The output of the circuit shown in Fig. 2 has the following form,
h
i
h
i
†
< eiφ hvref | Rk,i
Pr R |vref i = Tr (X ⊗ In )(CPr )ρφk,i (CP†r )= Tr Or ρφk,i .

(4.8)
(4.9)

In general, given an n-qubit unitary operator as U , estimating an observable written by Õ ≡ U (X ⊗In )U † requires many two-qubit gates. In fact, previous literature utilizes many two-qubit gates to estimate
Oq in (4.5) as we see in Fig. 1 and Fig. 2. However, in our case, we can write Oq in (4.5) by only one
Pauli observable. That is because CPq is a Clifford operator, and X ⊗ In is a Pauli observable; a Clifford
operator generally maps one Pauli observable to another Pauli observable. Direct calculation gives the
explicit form of Oq as
Oq = CP†q (X ⊗ In )CPq = (|0ih0| ⊗ In + |1ih1| ⊗ Pq )(X ⊗ In )(|0ih0| ⊗ In + |1ih1| ⊗ Pq )= (|0ih1| ⊗ In + |1ih0| ⊗ Pq )(|0ih0| ⊗ In + |1ih1| ⊗ Pq )= |0ih1| ⊗ Pq + |1ih0| ⊗ Pq
= X ⊗ Pq .

(4.10)

To summarize, we do not need the two-qubit operations in Fig. 1 and Fig. 2 and we can compute eachsummand of (2.16) and (2.17) as
h
i
h
i
†
< eiφ hvref | Rk,i
Pq R`,i0 |vref i = Tr (X ⊗ Pq )ρφk,i;`,i0 ,(4.11)hihi†< eiφ hvref | Rk,i
Pr R |vref i = Tr (X ⊗ Pr )ρφk,i .
(4.12)

4.2 VQS with shadow-based strategies
The benefit of writing each Oq and Or as a Pauli observable is not limited to the reduction of the two-qubit
gates. Although we consider a more general case later, we begin by describing a simplified scenario thateach phase of the summand in (2.16) does not depend on the index q and that in (2.17) does not dependon the index r; namely∗gk,ig`,i0 βq = Gk,i,`,i0 ,q eiφk,i,`,i0 (∀q),∗
gk,iαr = Gk,i,r eiφk,i(4.13)(∀r),(4.14)

where Gk,i,`,i0 ,q and Gk,i,r are real values, and eiφk,i,`,i0 and eiφk,i are phases independent of indices q andr. Then, we obtain
Ng
Ng
X
X
i,i0
Mk` =
Mk`
, Vk =
Vki ,(4.15)
i,i0

i=1

Accepted in Quantum 2023-03-28, click title to verify. Published under CC-BY 4.0.

12

0

i,i
with Mk`≡φPNBB

q=1 Gk,i,`,i ,q Tr0hφk,i,`,i(X ⊗ Pq )ρk,i;`,i00iand Vki =PNBA

r=1 Gk,i,r Trhiφ(X ⊗ Pr )ρk,ik,i . Note that

φk,i,`,iρk,i;`,iis the density operator ρφk,i;`,i0 with φ = φk,i,`,i0 and ρk,ik,i is the density operator ρφk,i with φ = φk,i .
0
PK
i,i0
The notable point here is that Mk`
and Vki have the form of hHi = j=1 aj Tr(Pj ρ) with {aj } as realcoefficients, {Pj } as Pauli observables, and ρ as a quantum state. Therefore, we can use shadow-based
i,i0
strategies to efficiently compute the summation inside the brackets. Namely, by our estimating Mk`
andVki as
0

φ

0

NBB
NBB
k,i,`,iνk,i;`,i0 = estimateNu(ρk,i;`,i0 , {Gk,i,`,i0 ,q }q=1 , {X ⊗ Pq }q=1 , Nshot ),(4.16)φNBABA
νk,i = estimateNu(ρk,ik,i , {Gk,i,q }N
q=1 , {X ⊗ Pq }q=1 , Nshot ),the following valuesM̃k` =NgX

νk,i;`,i0 , Ṽk =

i,i0

Ng
X

νk,i ,(4.17)

i=1

become unbiased estimator of Mk` and Vk . We summarize our proposed estimation process of Mk` and Vk
in Algorithm 3 and Algorithm 4. It should be noted that our circuit synthesis is also useful for grouping
i,i0
strategies [8, 16, 18, 22–33], mainly targeting VQO problems, where we need to calculate Mk`
and Vki .
The variances of M̃k` and Ṽk are given by
Var(M̃k` ) =Ng
X

Var (νk,i;`,i0 ) , Var(Ṽk ) =

i,i0

Ng
X

Var(νk,i ).

(4.18)

i=1
φ

0

NBB
k,i,`,iBy using (3.10) with {aj }K
j=1 = {Gk,i,`,i0 ,q }q=1 and ρ = ρk,i;`,i0 , we obtain(∆Mk` )2,Nv shot

uN
N
N
g
uX
BB X
BB
2
X
X
0
(G
)
u
φ
k,i,`,i ,qk,i,`,i0 
2 ,
∆Mk` ≡ t+
Gk,i,`,i0 ,q Gk,i,`,i0 ,q0 g(P̃q , P̃q0 )Tr(P̃q P̃q0 ρk,i;`,i
− Mk`
0 )
q(P
)
q
0
q=1
q=1 0

Var(M̃k` ) =

i,i

q 6=q

(4.19)
φNBAk,iwhere P̃q = X ⊗ Pq , and by using (3.10) with {aj }K
j=1 = {Gk,i,r }q=1 and ρ = ρk,i , we obtain(∆Vk )2,Nv shot

uN
N
N
g
uX
BA
BA X
2
X
X
(G
)
u
φ
k,i,r

Gk,i,r Gk,i,r0 g(P̃r , P̃r0 )Tr(P̃r P̃r0 ρk,ik,i ) − Vk2 ,∆Vk ≡ t
+
q(P
)
r
0
r=1
r=1
i

Var(Ṽk ) =(4.20)

r 6=r

where P̃r = X ⊗ Pr .
The above computed variances when using the shadow-based strategy
fromh are different
i the ones
PNBB
φk,i,`,i0
i,i00
when we naively estimate each summand of Mk` ≡and Vki =q=1 Gk,i,`,i ,q Tr (X ⊗ Pq )ρk,i;`,i0
h
i
PNBA
φk,i
r=1 Gk,i,r Tr (X ⊗ Pr )ρk,i , which is done by the previous literature (we simply call the strategy as
the naive strategy in the rest of the paper). Let Nnaive be the number of measurements for each circuit innaive
the naive strategy, and let the corresponding estimators of Mk` and Vk be M̃k`and Ṽknaive . Then the

Accepted in Quantum 2023-03-28, click title to verify. Published under CC-BY 4.0.

13

naivevariances of M̃k`and Ṽknaive are given by
naive 2(∆Vknaive )2(∆Mk`), Var(Ṽknaive ) =,NnaiveNnaivev
u Ng NBB 


2 
uX X
φk,i,`,i0naive
t2∆Mk` ≡(Gk,i,`,i0 ,q ) − Gk,i,`,i0 ,q Tr P̃q ρk,i;`,i0,

naive
Var(M̃k`)=(4.21)

i,i0 q=1v
u Ng NBA 


2 
uX X
φ
naive
∆Vk
≡t
(Gk,i,r )2 − Gk,i,r Tr P̃r ρk,ik,i
.
i

r=1

The difference in the variances between the naive strategy and the shadow-based strategy results in the
difference in the infidelity due to the shot noise for a given number of measurements. In the following
subsection, we quantitatively evaluate how the variance affects the infidelity and clarifies the benefit ofusing the shadow-based strategy in the VQS.
Note that in the above introduction of our proposed algorithm, we require that each phase of thesummand in (2.16) does not depend on the index q and that in (2.17) does not depend on the index r.
However, we can remove the requirement as follows. Instead of (2.14) , we can expand B † (t)B(t) as
B † (t)B(t) =NBB
X<(βq )Pq +

q=1

N
BB
X

=(βq )Pq ,(4.22)

q=1

PNBB
PNBB
where the symbol =(z) is the imaginary part of z. Then, we can evaluate q=1<(βq )Pq and q=1=(βq )Pq
independently by the shadow-based strategy though the number of measurements increase. That is also
true for the computation of B † (t)A(t).
Hybrid strategy in RTE and ITE
In RTE and ITE, NBB = 1 is set to calculate Mk` and there are no advantages for using a shadow-based
strategy for the computation of Mk` . Therefore, in RTE and ITE, we take the hybrid strategy, meaningthat we estimate Vk by the shadow-based strategy and Mk` by the naive strategy. At the end of Section 4.3,
we discuss the advantage of the hybrid strategy.
Algorithm 3 Estimation of Mk` with a shadow-based strategy.
φ

0

NBB
NBB
k,i,`,iInput: ρk,i;`,i0 , {Gk,i,`,i0 ,q }q=1 , {Pq }q=1 , Nshot , Ng
1: SET result = 0.
2: for i = 1...Ng do
3:for i0 = 1...Ng do
φk,i,`,i0
NBB
NBB
4:
SET νk,i;`,i0 = estimateNu(ρk,i;`,i0 , {Gk,i,`,i0 ,q }q=1 , {X ⊗ Pq }q=1 , Nshot ).
5:
SET result = result + νk,i;`,i0 .
6:
end for
7: end for
return result

4.3 Advantage of using shadow-based strategy
To clarify the advantage of using the shadow-based strategy over the naive strategy, we compare theinfidelity due to the shot noise after the duration T . The time derivative of the parameter vector can beideally computed as∂ θ~ideal (t)~.
= M −1 V(4.23)∂t

Accepted in Quantum 2023-03-28, click title to verify. Published under CC-BY 4.0.

14

Algorithm 4 Estimation of Vk with a shadow-based strategy.
φ

NBA
BA
Input: ρk,ik,i , {Gk,i,q }N
q=1 , {Pq }q=1 , Nshot , Ng
1: SET result = 0.
2: for i = 1...Ng do
3:for i0 = 1...Ng do
φ
NBA
BA
4:
SET νk,i = estimateNu(ρk,ik,i , {Gk,i,q }N
q=1 , {X ⊗ Pq }q=1 , Nshot ).
5:
SET result = result + νk,i .
6:
end for
7: end for
return result~ with a finite number of measurements, and therefore, the estimated time
However, we estimate M and V
derivative includes the shot noise:~∂ θ(t)~ + δV~ ),= (M + δM )−1 (V∂t(4.24)~ are the error to M and V~ due to the shot noise. As in literature [9], we use the following
where δM and δ V
infidelity as a metric for evaluating the effect of error:
D
E2

E
E r
~ )) ,~~(4.25)DI v(θideal (T )) , v(θ(T )) = 1 − v(θ~ideal (T )) v(θ(TEE~ )) evolved by the distorted
between v(θ~ideal (T )) evolved by the ideal time derivative (4.23) and v(θ(Ttime derivative (4.24) at time T .
Suppose that DInaive is an upper bound of DI when using the naive strategy and DIshadow is the one
when using the shadow-based strategy. As we show in (C.19) in Appendix C, DInaive = DIshadow holdswhen max 2
Nshot
∆shadow',(4.26)
Nnaive
∆max
naive
where
∆naive = ||M

−1

sX
sX
−1 2 ~
naive
naive )2 ,2(∆Vk
) + ||M ||F ||V ||2
(∆Mk`
||F
k

k`

(4.27)sX
sX~ ||2∆shadow = ||M −1 ||F(∆Vk )2 + ||M −1 ||2F ||V(∆Mk` )2 ,
k

k`

max
∆max
naive ≡ max0≤t≤T (∆naive ), ∆shadow ≡ max0≤t≤T (∆shadow ), || · ||F is the Frobenius norm, and || · ||2 isthe two-norm.
Note that the total number of measurements for each step in the case of the shadow-based strategy
total
total
is Nshadow= Nshot (NP Ng + NP2 Ng2 ) while that in the naive strategy is Nnaive= Nnaive (NP Ng NBA +
2
2
total
total
NP Ng NBB ). Thus, the advantage of using a shadow-based strategy arises if Nshadow/Nnaive< 1 when(4.26) is satisfied, i.e., max 2
total
Nshadow
1 + NP Ng∆shadow=< 1.
(4.28)
total
NBA + NP Ng NBB∆max
Nnaivenaive

In the rest of this section, we give a rough estimate of the leftest-hand side of the inequality.
total
total
Estimation of Nshadow/Nnaive~ ||2 when using the naive strategy is almost the same
We assume that the values of ||M −1 ||F and ||M ||2F ||V
as those when using the shadow-based strategy. The assumption is valid since the time evolutions withAccepted in Quantum 2023-03-28, click title to verify. Published under CC-BY 4.0.

15

both strategies are close to each other (recall that we consider the regime where DInaive = DIshadow ). Under
the assumption, we obtain the inequality:s P
s P(∆Vk )2(∆Mk` )2
∆max
shadowk+ P k`(4.29)< P
naive
naive )2 .
max
2)∆naivek (∆Vkk` (∆Mk`We evaluate each term on the right-hand side in the following.
For the numerator in the square root of the first term, it holds(∆Vk )2 ≤NgXφVki (ρk,ik,i ),(4.30)i=1

where
Vki (ρ) =NBA X(Gk,i,r )2 X+
Gk,i,r Gk,i,r0 g(P̃r , P̃r0 )Tr(P̃r P̃r0 ρ).
q(Pr )
r=1 0
r=1

N
BA
X

(4.31)

r 6=r

Since Vki (ρ) depends on the quantum state ρ, it is difficult to directly estimate the value. Instead ofevaluating Vki (ρ), let us estimate the typical size of Vki (ρ) in the overall Hilbert space by computingZZ
2hVki iHaar =dU Vki (U |0ih0|U † ), and Var(Vki )Haar =dU Vki (U |0ih0|U † ) − (hVki iHaar )2 , (4.32)
Haar

Haar

R

where Haar is the average over the Haar distribution. By calculating those statistical values, we can see
how Vki (ρ) distributes when ρ is uniformly sampled. As we show in Appendix D, we can readily show
hVki iHaar =Var(Vki )Haar =NBA
X(Gk,i,r )2,q(Pr )
r=1
1

N
BA X N
BA X
X
X

2n+1 (2n+1 + 1) r=1 0

Gk,i,r Gk,i,r0 Gk,i,u Gk,i,u0 g(P̃r , P̃r0 )g(P̃u , P̃u0 )Tr(P̃r P̃r0 P̃u P̃u0 ),r 6=r u=1 u0 6=u(4.33)
where n is the number of qubits in PQC. By using Tr(P̃r P̃r0 P̃u P̃u0 ) ≤ 2n+1 , we can show Var(Vki ) ∝O(1/2n ), which becomes asymptotically zero as n becomes large. Since the variance is exponentiallyφsuppressed, Vki (ρk,ik,i ) takes the value almost equal to the average with high probability when uniformly
φchoosing ρk,ik,i . In the rest of this section, given a function of a density operator as ξ(ρ), if hξ(ρ)iHaar = ξ0holds with ξ0 as a constant and Var(ξ)Haar is exponentially suppressed, we denote ξ(ρ) ' ξ0 . Also, for
given another function of a density operator as ξ 0 (ρ), if ξ 0 (ρ) ≤ ξ(ρ) ' ξ0 , we denote ξ 0 (ρ) . ξ0 . Thus,
with a large n,NBAX(Gk,i,r )2
φ(4.34)Vki (ρk,ik,i ) 'q(Pr )
r=1
holds and(∆Vk )2 .

Ng NBAX
X (Gk,i,r )2
i=1 r=1q(Pr ).

(4.35)In the Eq. (4.29), for the denominator of the square root in the first term,(∆Vknaive )2 =Ng NBAX
X

φGki,r (ρk,ik,i ),(4.36)

i=1 r=1

where

h
2 i
Gki,r (ρ) = (Gk,i,r )2 1 − Tr(P̃r ρ).

Accepted in Quantum 2023-03-28, click title to verify. Published under CC-BY 4.0.

(4.37)

16

As in the case of Vki (ρ), we can compute the average in the overall Hilbert space (See Appendix D) as1,hGki,r (ρ)iHaar = (Gk,i,r )2 1 − n+12+1(4.38)
(Gk,i,r )4 (2 · 22(n+1) + 8 · 22(n+1) )Var Gki,r (ρ)=.
(2n+1 + 1)2 (2n+1 + 3)
HaarThus, for a large enough n, Gki,r (ρ) ' (Gk,i,r )2 and, we obtain
Ng NBAX
X
naive 2'
(Gk,i,r )2 .
∆Vki=1 r=1(4.39)

By combining (4.34) and (4.39), the first term in the right hand side of (4.29) is bounded ass Ps(∆Vk )2
1
kP.
,
naive )2q(P̃r ) Gk,i,rk (∆Vkwhere h·iGk,i,r is a weighted average defined by
PNP PNg PNBA
2
i=1
k=1
r=1 (Gk,i,r ) Fr.
hFr iG = P
P
P
Ng
NBA
NP
2
i=1
k=1
r=1 (Gk,i,r )(4.40)(4.41)

Similarly,
v

u
2
u
(∆M
)
1
k`
t
P k`
naive )2 .
q(P̃q ) Gk` (∆Mk`s P,(4.42)k,i,`,i0 ,qwhere h·iGk,i,r is another weighted average defined by
PNP PNg PNBB
hFq iG

2
0
i,i0
k,`
q=1 (Gk,i,`,i ,q ) Fq
=
.
P
P
P
Ng
NP
NBB
k,i,`,i0 ,q
2
i,i0
k,`
q=1 (Gk,i,`,i0 ,q )By substituting (4.40) and (4.42) to (4.29), we obtain
v
s


u
u
1
1
∆max
shadowt.
+∆max
q(P̃r ) Gk,i,rq(P̃q ) G
naive

(4.43).

(4.44)k,i,`,i0 ,qCombining with (4.28), we obtain
s
v



u
total
u
1
1
Nshadow
1 + NP Ng

t
.
+
total
NBA + NP Ng NBB
Nnaiveq(P̃r ) Gk,i,rq(P̃q ) G

2
 .

(4.45)k,i,`,i0 ,qThe factor (1 + NP Ng )/(NBA + NP Ng NBB ) < 1 corresponds to the benefit of simaltaneously evaluating the NBA terms in the calcualtionof Vk and theq
q NBB terms in the calculation of Mk` in the shadow-based
strategy. The other factor
1/q(P̃r ) G+
1/q(P̃q ) G
corresponds to how well measurements
0k,i,rk,i,`,i ,qin the shadow-based strategy cover the observables; if observables are not covered efficiently, the benefitof the simaltaneous evaluation may disappear due to this factor. Particularly, when we use the classicalshadow as the shadow-based strategy, we can explicitly write the latter factor as
2
s
v



u
q
2
q
u
1
1
 =

+t
3locality (P̃r ) G+
3locality (P̃q ) G,k,i,rk,i,`,i0 ,qq(P̃r ) Gk,i,rq(P̃q ) G
0k,i,`,i ,q(4.46)
where we use (3.13). Thus, the locality of the observables should be small in order for enjoying the benefitof the classical shadow in the VQS. It should be noted that
locality(P̃r ) = locality(Pr ) + 1,

(4.47)

because P̃r = X ⊗ Pr .

Accepted in Quantum 2023-03-28, click title to verify. Published under CC-BY 4.0.

17Advantages of using the shadow-based strategy in RTE and ITE
Finally, we discuss the advantage of using the shadow-based strategy in RTE and ITE. As we noted at
the end of Section 4.2, we take the hybrid strategy that utilizes the shadow-based strategy for computing
Vk and the naive strategy for computing Mk` . Unlike the general case (where we measure both Vk andMk` by using the shadow-based strategy), we can reduce the number of measurements for estimating eachMk` to smaller than that for estimating Vk (=Nshot ); we set the number of measurements as Nshot α with
0 < α < 1.
Let
rP
sX
naive )2−1
2
−1k` (∆Mk`~ ||2(∆Vk )2 + ||M ||F ||V.
(4.48)∆hybrid = ||M ||FαkThen, as shown in (C.24) in Appendix C, we can show that
Nshot'
Nnaive∆maxhybrid
∆max
naive

2
,

(4.49)

if the error due to the shot noise in the naive strategy is almost the same as that in the hybrid strategy,
total
where we define ∆max
hybrid ≡ max0≤t≤T (∆hybrid ). Given Nhybrid as the total number of measurements for
totaleach step in the hybrid strategy, it holds Nhybrid = Nshot (NP Ng +α(NP Ng )2 ). Under the same assumption
to obtain (4.29), it holds
total
Nhybrid

1 + NP Ng α.
total
NBA + NP Ng
NnaiveBy setting α = 4/

D

1
q(P̃r )sE

1
q(P̃r )

r !2
1
1 + NP Ng α
+
<
α
NBA
Gk,i,rs

1q(P̃r )



r !2
1
+
.
α
Gk,i,r(4.50),

Gk,i,rtotal
Nhybrid
total
Nnaive.

1 + Nα×NBA1q(P̃r )

withNα = Np Ng α = D



9
× ,
4
Gk,i,r

4NP NgE

1q(P̃r )

(4.51).

(4.52)

Gk,i,r

Note that the literature [49] states that given a circuit structure, the number of measurements forNBAPestimating {Mk` }Nk,` asymptotically becomes negligible compared to that for estimating {Vk }k=1 as thenumber of qubits n increases and NBA increases faster than quadratically with n. Though the literatureBA
assumes the naive strategy to estimate {Vk }N
k=1 , the asymptotic behavior would be the same in our hybridP
strategy. The condition that the number of measurements for {Mk` }N
k,` is negligible, corresponds to thecase when Nα  1. In that case, it holds


total
Nhybrid
9
1
1
.
×
× .
(4.53)
total
NBA
Nnaiveq(P̃r ) Gk,i,r 4
For example, let us calculate the right-hand side value of the right-hand side of (4.53) when the time
evolution Hamiltonian is given by (3.22). We can write
X ⊗ H = XXXXZ + XXXII + XIIXZ + XY Y ZX + XY Y II + XIIZX,(4.54)and therefore, the number of terms is given by NBA = 6. Recall that each P̃r in (4.53) corresponds toeach term in (4.54). Let us first calculate when we use the classical shadow. Since the locality of eachterm in (4.54) is equal to or larger than three,1≥ 33 ,(4.55)q(P̃r ) Gk,i,r

Accepted in Quantum 2023-03-28, click title to verify. Published under CC-BY 4.0.

18

where we use (3.13). Thus, in the case of the classical shadow, the right-hand side of (4.53) is


1
1
9
×
× ≥ 10,NBA
q(P̃r ) Gk,i,r 4(4.56)and there is no advantage in using the classical shadow with the Hamiltonian (3.22). On the contrary, whenusing derandomization, half of the chosen measurements are XXXXZ, and the other half is XY Y ZX
(see the discussion around (3.22)). Therefore,


1
1
= ,
(4.57)
2q(P̃r ) Gk,i,r
and the right-hand side of (4.53) is
1
×
NBA



1
q(P̃r )

×
Gk,i,r

9
3
= .
4
4

(4.58)

Thus, there is merit in using derandomization with the Hamiltonian (3.22). We note that we utilize the
Hamiltonian with just six terms in this example, and the advantage of using the shadow-based methods islimited. In Section 5, we will demonstrate the result of VQS when using the molecular Hamiltonian andsee the significant advantages of using the shadow-based methods.

4.4 Remark on the general shadow-based strategies
In Section 4.3, we showed that(∆Vk )2 .

Ng NBAX
X Gk,i,rq(Pr )
i=1 r=1, and (∆Vknaive )2 'Ng NBA
X
X(G2k,i,r )i=1 r=1

in the context of VQS. The same argument holds even in the general estimation problem of hHi =PK
j=1 aj Tr(Pj ρ) discussed in Section 3. Let νnaive be the estimator of hHi when estimating each termof hHi one by one with Nnaive measurements. Also, recall that ν is the estimator of hHi when using ashadow-based strategy with Nshot measurements. Then, with the same arguments in Section 4.3, it canbe shown that
Var(ν) .

K
X
a2j,
Nshot j=1 q(Pj )Var(νnaive ) '

1

1

K
X

Nnaive j=1(4.59)

a2j .

(4.60)

total
The total number of measurements is Nshadow= Nshot in the case of using the shadow-based strategy
total
and Nnaive = KNnaive in the case of naively measuring each term one by one. For achieving Var(ν) =Var(νnaive ), the following relation is necessary

PK
a2j
total
1
Nshadow
j=1 q(Pj ).
PK 2 .
total
K
Nnaivej=1 aj(4.61)

Notably, as long as the covering probability q(Pj ) is O(1), we can reduce the number of measurements by
the factor O(1/K) in the shadow-based strategy.
The above observations consolidate the validity of using the shadow-based strategies [17–21] for estimating the sum of Pauli observables. In Eq. (18) of [18], they also give the upper bound for the varianceof ν. However, it is as loose asK
K X a2jVar(ν) ≤,(4.62)
Nshot j=1 q(Pj )Accepted in Quantum 2023-03-28, click title to verify. Published under CC-BY 4.0.

19

in our notation (in [18], they write the number of terms as |H0 | instead of K). We can easily check that
total
total
if Var(ν) takes the value of the upper-bound (4.62), Nshadow> Nnaive, and there are no advantages ofusing shadow-based strategies. Conversely, previous research does not theoretically show the advantagesof using a shadow-based strategy to estimate the sum of Pauli observables. Our upper bound, however,successfully removes the effect from the second term in (3.10) by using the average over the whole quantum
state and shows that in most of the quantum states ρ, there is an advantage to using the shadow-based
strategy.

5 Numerical Demonstration
Here we demonstrate our shadow-based measurement optimization strategy for the VQS. In Section 4,
we generally discussed our strategy, but in this demonstration, we focus on RTE and ITE, which is thealgorithm most actively studied in the VQS algorithms. As we noted at the end of Section 4.2, we take thehybrid strategy, meaning that we obtain Mk` by the naive strategy and Vk by the shadow-based strategy.
We perform two types of demonstrations. The first one is the evaluation of the infidelity due to theshot noise,
E2
D
E r
E

~ )) ,~~DI v(θideal (T )) , v(θ(T )) = 1 − v(θ~ideal (T )) v(θ(Twhich is defined in (4.25), in various Hamiltonians, and in various measurement strategies. The second
one is the evaluation of the variances in each strategy. Particularly, as a result of the Haar average in theHilbert space, we obtained(∆Vk )2 .

Ng NBAX
X (Gk,i,r )2
i=1 r=1q(Pr ),∆Vknaive

2'

Ng NBA
X
X(Gk,i,r )2 ,i=1 r=1

in (4.34) and (4.39) in Section 4.3. We check if it holds in our numerical experiments.
Before going into each experiment, let us give the setting of numerical experiments in Section 5.1.
Then, in the subsequent sections, we demonstrate the results of the numerical experiments.

5.1 Setting of numerical experiments
Here we summarize Hamiltonians, the reference state, and the ansatz used in the experiments. Also, we
describe strategies to estimate Mk` and Vk .
Hamiltonians
We use three types of problem Hamiltonians H:• H2 molecule in the 6-31g basis with the Bravyi-Kitaev transformation (8-qubit, NBA = 184);• LiH molecule in the STO-3g basis with the Bravyi-Kitaev transformation (12-qubit, NBA = 630);• Heisenberg model Hamiltonian (6-qubit, NBA = 24).
We write the number of qubits in the Hamiltonian as n. Here we define our Heisenberg model Hamiltonian
by
6
6
X
X
X
H = 0.1
σaj σaj+1 + 0.1
σzj ,
(5.1)
j=1 a=x,y,z

j=1

with σaj as a Pauli operator operating on the j-th qubit, and σa7 = σa1 .

Accepted in Quantum 2023-03-28, click title to verify. Published under CC-BY 4.0.

20

Reference state and ansatz
As the reference state |vref i, we use the |vref i = H⊗n |0i⊗n . As the choice of the ansatz R (defined in(2.10)), we use a hardware efficient ansatz, which contains multiple layers of single-qubit rotation gatesand entanglers. Each parameter is embedded in each single-qubit rotation gate as exp(−iσa θj /2), where
σa is a Pauli matrix with a ∈ {x, y, z} (therefore, Ng = 1 and gk,1 = −i/2). The type of each single-qubit
rotation gate ‘a’ is randomly chosen at the initialization. As the two-qubit gates in entanglers, we use
the controlled-Z gates. At each layer, we first perform the single qubit gates and, secondly, perform thecontrolled-Z gates. The way of operating the controlled-Z gates is different in the odd layers and the evenlayers. Let us define a set of integers N1 by N1 = {m1 |m1 ∈ N and 1 < 2m1 ≤ n} and another set of
integers N2 by N2 = {m2 |m2 ∈ N and 1 < 2m2 + 1 ≤ n}. For example, when n = 8, N1 = {1, 2, 3, 4}
and N2 = {1, 2, 3}. Also, we denote by (j, k) the pair of control and target qubits on which we perform a
two-qubit operation. Then, in each odd layer, pairs of (2m1 − 1, 2m1 ) for all m1 ∈ N1 are firstly entangled
by the controlled-Z gates, and next, pairs of (2m2 , 2m2 +1) for all m2 ∈ N2 are entangled. On the contrary,in each even layer, pairs (2m2 , 2m2 + 1) for all m2 ∈ N2 are firstly entangled by the controlled-Z gates,
and next, pairs (2m1 − 1, 2m1 ) for all m1 ∈ N1 are entangled. In Fig. 4, we show our ansatz with thereference state |vref i when n = 8 as an example. For all experiments, we fix the number of layers as four,and therefore, NP = 4n. All parameters are initialized as zero, and therefore, the initial state generated
by the ansatz is |vref i.

Figure 4: The reference state and the ansatz used in the numerical experiment. As an example, we shot the case whenn = 8.

Strategies to estimate Mk` and Vk
For the estimation of Vk , we compare following four strategies including two shadow-based strategies:• Classical shadow (=a shadow-based strategy);• Derandomization (=another shadow-based strategy);• Naively estimate each term of Vk one by one;• Largest degree first (LDF) grouping [16, 18].
The LDF grouping is a method to reduce the number of measurements by grouping the simultaneously
measurable Pauli observables. In the LDF grouping, the measurement basis is chosen so that all of theobservables in one of the groups should be covered by each measurement. The group to be covered by eachmeasurement is chosen according to the probability distribution weighted by the importance of each group,

Accepted in Quantum 2023-03-28, click title to verify. Published under CC-BY 4.0.

21

i.e., the larger the norm of the coefficients of the observables in the group is, the larger the probability
becomes. As shown in previous literature [18–20], the LDF grouping successfully decreases the numberof measurements in the VQE problems without increasing the depth of the circuit. Due to the circuit
synthesis in Section 4.1, Mk` and Vk are written as the sum of Pauli observables, and the LDF groupingis also usable for the evaluation of Vk and Mk` . Thus, it can be a good benchmark for evaluating the
shadow-based strategy. For a more detailed explanation of the LDF grouping, please see [18].
For the estimation of Mk` , there is no room for measurement optimization since NBB = 1. To emphasize the difference in the performances of the above four estimation strategies of Vk , we assume thatMk` can be exactly computed (which corresponds to the case when the number of measurements for estimating {Mk` } is negligible). We can obtain the exact value of Mk` in Eq. (2.16) by directly calculating†
hvref | Rk,i
Pq R`,i0 |vref i with the state vector representation. Finally, throughout our numerical experiment,
for simulating the quantum process, we use Qulacs [50].

5.2 Evaluation of the infidelity
Under the setting we described in Section 5.1, we execute the RTE and ITE for the initial
five stepsE such asE~~
evolves
T = xδt (x ∈ {0, 1, 2, 3, 4, 5}) where we set δt = 0.01, and see how DI v(θideal (T )) , v(θ(T ))
E
E

~ ))
= 0.
depending on the strategies to estimate Vk . For T = 0, we have DI v(θ~ideal (T )) , v(θ(TE
For computing v(θ~ideal (T )) , we exactly compute V by using the state vector representation, while forE~ )) , we evaluate V by using each measurement strategy. Given the total number of
estimating v(θ(Tmeasurements for estimating each Vk as Ntotal , it is set to be Ntotal = 5 × NBA for each strategy (thus,
Nshot = 5 × NBA for the shadow-based strategy and Nnaive = 5 for the naive strategy). For example, inthe case of the H2 molecule Hamiltonian,NBA =E 184, and theE total number of measurements is 920.
~~
The resulting evolutions of DI v(θideal (T )) , v(θ(T ))for each T are shown in Fig. 5. There are
four subfigures in it, which correspond to the followings: (a) ITE with H2 molecule Hamiltonian, (b)RTE with H2 molecule Hamiltonian, (c) ITE with Heisenberg model Hamiltonian and (d) ITE with LiH
molecule Hamiltonian. We perform five trials and write the mean and std error of the mean for each data
point. For each trial, different gate types are randomly chosen in PQC at initialization, but the chosengate types are the same in every strategy as far as they are in the same trial. For example, the gate types
chosen in the first trial of the naive strategy are the same as those in the first trial of the classical shadow
strategy.
In Fig. 5, we see that the classical shadow outperforms the naive method other than the problems with
the Heisenberg model Hamiltonian. The reason why the classical shadow does not work well in the case of
the Heisenberg model Hamiltonian is as follows. The Heisenberg model Hamiltonian (5.1) has NBA = 24
terms, and in the classical shadow, it holds

1' 33 ,(5.2)
q(P̃r ) Gk,i,r
where we use that the locality of most of the terms in (5.1) is two and in the observable to evaluate Vk ,the locality increases by one (see (4.47)). By substituting these into (4.53),
total
Nhybrid

1
.
×
total
N
NnaiveBA1q(P̃r )

×
Gk,i,r9' 2.5,
4(5.3)

meaning that there is not an advantage of using the classical shadow algorithm in the Heisenberg model
Hamiltonian as far as we see the upper bound.
Even though the classical shadow algorithm successfully works in the other problems compared to thenaive strategy, the LDF grouping outperforms the classical shadow in every problem. That is because theclassical shadow contains inefficiency regarding the choice of the measurement basis. The most significant
inefficiency might exist in the choice of the measurement basis for the ancilla qubit. For evaluating Vk , we

Accepted in Quantum 2023-03-28, click title to verify. Published under CC-BY 4.0.

22

(a)

(b)

ITE with H2 molecule Hamiltonian

0.05

naive
classical shadow
ldf grouping
derandomization

0.06

0.04

0.05

DI
0.02

DI

naive
classical shadow
ldf grouping
derandomization

0.03

0.04
0.03
0.02

0.01
0.00

RTE with H2 molecule Hamiltonian

0.07

0.01
0.00

0.01

0.02 0.03
time

0.04

0.00

0.05

0.00

0.01

(c)

DI

DI

0.015
0.010
0.005
0.000

0.00

0.01

0.05

ITE with LiH molecule Hamiltonian

naiveclassical shadow
ldf grouping
derandomization

0.020

0.04

(d)

ITE with Heisenberg Hamiltonian
0.025

0.02 0.03
time

0.02 0.03
time

0.04

0.05

0.040
0.035
0.030
0.025
0.020
0.015
0.010
0.005
0.000

naive
classical shadow
ldf grouping
derandomization

0.00

0.01

0.02 0.03
time

0.04

0.05

~ )) for each T . Four subfigures correspond to the followings: (a)Figure 5: The evolutions of DI v(θ~ideal (T )) , v(θ(TITE with H2 molecule Hamiltonian, (b) RTE with H2 molecule Hamiltonian, (c) ITE with Heisenberg model Hamiltonian,and (d) ITE with LiH molecule Hamiltonian.



need to evaluate the following summation of the observables:N
BA
X

αr X ⊗ Pr .

(5.4)

r

Since all the observables have X in the first qubit, the measurement basis that does not have X in thefirst qubit does not give any information about Vk . Conversely, only the measurement basis that has X in
the first qubit covers the terms in (5.4). However, the measurement basis is chosen with the probability of
1/3, which is obviously inefficient in the computation of Vk . The inefficiency is the reason why the LDF
grouping outperforms the classical shadow.
We enjoy the benefit of the shadow-based strategy with a more intelligent selection of the measurement
basis. The derandomization, which optimizes the measurement basis so that 1/q(Pr ) becomes smaller,
outperforms the LDF grouping other than the problem with the Heisenberg model Hamiltonian. In the
Heisenberg model Hamiltonian, all observables are covered by the only three measurement basis:
XXXXXXX, XY Y Y Y Y Y, XZZZZZZ.

(5.5)

Finding such a measurement basis is easy both for the LDF grouping and the derandomization. In our
numerical experiment, the LDF grouping successfully groups all the observables into three groups, where

Accepted in Quantum 2023-03-28, click title to verify. Published under CC-BY 4.0.

23

each group is covered by one of the above measurement basis. Thus, the measurements in the LDF
grouping are nearly optimal and equal to the derandomization in the Heisenberg model Hamiltonian.
However, in the other problems, the observables are not covered by a few measurement basis. Thus thederandomization, which is a relatively intelligent way of finding the measurement basis, outperforms the
LDF grouping.
EE~ )) at T = 1 with various total numbers of measurements
We also compute DI v(θ~ideal (T )) , v(θ(TNtotal . Similarly to Fig 5, each data point is the mean and the standard error of the mean in five
trials. The results are shown in Fig. 6 (the labels of the subfigures are the same as those of Fig. 5). We
note that both the horizontal and vertical axes are in the logarithmic scale in the figure. We see thatthe same conclusionE the performance of each measurement strategy. Additionally, we see that
E holds for~~almost decreases with the square root of Ntotal , which is consistent withDI v(θideal (T )) , v(θ(T ))(C.23) in Appendix C.
(a)

(b)

RTE with H2 molecule Hamiltonianslope = -0.5
naivecs
ldf grouping
derandomization

10 2

103

2 × 103DI

DI

ITE with H2 molecule Hamiltonian

10 2

103

3 × 103

Ntotalslope = -0.5
naivecs
ldf grouping
derandomization

(c)

10 2

2 × 103

Ntotal3 × 103(d)

ITE with Heisenberg Hamiltonian

ITE with LiH molecule Hamiltonianslope = -0.5
naivecs
ldf grouping
derandomization

6 × 10 3

slope = -0.5
naivecs
ldf grouping
derandomizationDIDI

4 × 10 3
3 × 10 3

10 2

2 × 10 3

3 × 102

4 × 102

6 × 102

Ntotal103

2 × 103
Ntotal3 × 103~ )) at T = 1 for each Ntotal . As in Fig. 5, four subfigures correspond
Figure 6: The value of DI v(θ~ideal (T )) , v(θ(Tto the followings: (a) ITE with H2 molecule Hamiltonian, (b) RTE with H2 molecule Hamiltonian, (c) ITE with
Heisenberg model Hamiltonian, and (d) ITE with LiH molecule Hamiltonian. Both the horizontal and the vertical axes
−0.5
are in the logarithmic scale. As a reference, we plot a slope proportional to Ntotal
by a gray dotted line.



5.3 Evaluation of the variance
Finally, we numerically compare the variance Var(Ṽk ) and Var(Ṽknaive ) for a given total number of measurements. We set the total number of measurements as 5 × NBA . We compute the variance in the two

Accepted in Quantum 2023-03-28, click title to verify. Published under CC-BY 4.0.

24

problems: ITE with H2 molecule Hamiltonian, which is shown in TABLE 1, and ITE with Heisenberg
model Hamiltonian, which is shown in TABLE 2.
The detail of the contents of the tables is as follows. In the ‘Variance’ column, we show the average ofthe variance; namely,NP
1 X
Var(Ṽk ),(5.6)
Variance =NP
k=1

for the shadow-based strategy andN

Variance =P
1 X
Var(Ṽknaive ),
NP

(5.7)

k=1

for the naive strategy. Additionally, checking
2(∆Vk ) .

Ng NBAX
X (Gk,i,r )2
i=1 r=1q(Pr ),

Ng NBAX
X
naive 2'
(Gk,i,r )2 ,∆Vki=1 r=1

in the ‘Approximation’ column, we show the approximate value of ‘Variance’, which is denoted by
Ng NBANPX
X (Gk,i,r )2
1
1 X

,
(5.8)
Approximation =NP
Nshot i=1 r=1 q(Pr )k=1

for the shadow-based strategy, and
Ng NBA
NP
X
X
X
1
 1
Approximation =(Gk,i,r )2  ,NP
Nnaive i=1 r=1

(5.9)

k=1

for the naive strategy. Finally, in the ‘Diff’ row, we show the averege difference between the true value
and the approximate value of the variance, namely,
Diff =

Ng NBANP(Gk,i,r )2
1 XX
1 X
Var(Ṽk ) −,NP
Nshot i=1 r=1 q(Pr )(5.10)

k=1

for the shadow-based strategy, andNg NBA
NP
1 X
1 XX
naive
Diff =
Var(Ṽk)−(Gk,i,r )2 ,NP
Nnaive i=1 r=1(5.11)

k=1

for the naive strategy. The value of ‘Variance’ depends on the parameters of the ansatz. We choose
parameters in two different ways when computing ‘Variance’: (i) all zero, where all parameters are set tobe zero, and (ii) random, where all parameters are randomly selected. For the random parameters choicecase, we generate five patterns of parameters.
There are several implications from the results in the tables. First, both in TABLE 1 and TABLE 2,
we see that all ‘Variance’ are almost equal to the corresponding ‘Approximation’, which shows the validityof the approximation in (4.34) and (4.39) in Section 4.3. Second, the variances are consistent with theresults in Section 5.2. Namely, the smaller variance results in the smaller infidelity in Fig. 5 and Fig. 6.
Finally, variances in the derandomization are much smaller than those in the other measurement methods.
Particularly, the large difference in the variances between the classical shadow and the derandomizationshows the importance of optimizing the covering probability q(Pr ) for decreasing the variance (hence
the infidelity) and implies that other shadow-based strategies [18, 20], which also optimizes the covering
probability, also achieve good performance. Studying the effect of using [18, 20] in VQS is left for future
work.

Accepted in Quantum 2023-03-28, click title to verify. Published under CC-BY 4.0.

25

Table 1: Variance, Approximation, and Diff in the problem of ITE with H2 molecule Hamiltonian.

Parameters Measurement
naive
all zeroclassical shadow
derandomization
naive
randomclassical shadow(pattern1) derandomization
naive
randomclassical shadow(pattern2) derandomization
naive
randomclassical shadow(pattern3) derandomization
naive
randomclassical shadow(pattern4) derandomization
naive
randomclassical shadow(pattern5) derandomization

Variance
0.731
0.524
0.020
0.741
0.525
0.020
0.731
0.528
0.019
0.733
0.526
0.020
0.736
0.528
0.020
0.732
0.510
0.019

Approximation
0.743
0.525
0.020
0.743
0.525
0.020
0.743
0.525
0.020
0.743
0.525
0.020
0.743
0.525
0.020
0.743
0.525
0.020

Diff
0.012
0.001
0.000
0.001
0.000
0.000
0.011
0.002
0.000
0.010
0.000
0.000
0.006
0.002
0.000
0.010
0.014
0.000

Table 2: Variance, Approximation, and Diff in the problem of ITE with Heisenberg model Hamiltonian.

Parameters Measurement
naive
all zeroclassical shadow
derandomization
naive
randomclassical shadow(pattern1) derandomization
naive
randomclassical shadow(pattern2) derandomization
naive
randomclassical shadow(pattern3) derandomization
naive
randomclassical shadow(pattern4) derandomization
naive
randomclassical shadow(pattern5) derandomization

Variance
0.04750.04570.00630.04690.04640.00550.04720.04360.00540.0471
0.04430.00540.04670.04490.00650.04690.04560.0058

Accepted in Quantum 2023-03-28, click title to verify. Published under CC-BY 4.0.

Approximation
0.04800.04500.00590.04800.04500.00590.04800.04500.00590.04800.04500.00590.04800.04500.00590.04800.04500.0059Diff
0.00040.00070.00040.00100.00140.00030.00070.00130.00040.00080.00060.00050.00120.00000.00060.00100.00060.000126

6 Discussion
The variational quantum simulation (VQS) successfully gives a way of simulating quantum systems inthose devices, which is considered as one of the practical applications the near-term quantum devices.
However, when applying VQS to a large quantum system where we have practical interests, the numberof running a quantum device, i.e., the number of measurements, dramatically increases to reduce theeffect of shot noise. Even though the grouping methods and the shadow-based strategies have recently
been proposed to decrease the number of measurements in the variational quantum optimization (VQO),
they are not utilized in VQS because the way of evaluating the observables in VQS is different from theone in VQO and we cannot straightforwardly apply the shadow-based strategies to VQS. In this paper,by changing the circuit and observables in VQS, we make it possible to apply shadow-based strategies
and grouping methods to VQS. We particularly focus on the shadow-based strategies; we theoreticallyshow that with the shadow-based strategy, we can greatly decrease the number of measurements in VQS
compared to the naive strategy, which has been utilized since the original proposal of VQS [9, 13, 14]. Our
numerical demonstration of ITE and RTE using various Hamiltonians, including chemical Hamiltonians,
also implies the validity of using the shadow-based strategy. Particularly, the derandomization [19], oneof the shadow-based strategies, achieves the best performance in most of the Hamiltonians.
Though we examine only the classical shadow and the derandomization in our numerical experiments
as the shadow-based strategies, other shadow-based strategies [20, 21] can also be utilized to reduce thenumber of measurements in VQS. For example, the literature [20] shows that their proposed shadow-based
strategy using decision diagrams effectively reduces the estimator’s variance to smaller than that in the
LDF grouping, which achieves the second-best performance in our numerical demonstration. Additionally,
as we note in Section 1, recent literature [31] reports that their grouping strategies outperform shadowbased strategies, such as derandomization, by optimizing parameters in the grouping so that the varianceof the observable is minimized. Studying the effect of using those shadow-based strategies and the grouping
strategies in VQS is a good direction for future work. We believe that our work takes the first step toward
measurement optimization in VQS and paves the way for building practical applications for simulating
quantum systems.

Acknowledgement
This work was supported by Grant-in-Aid for JSPS Research Fellow 22J01501 and by Leading Initiativefor Excellent Young Researchers MEXT Japan and JST presto (Grant No. JPMJPR1919) Japan. This
work was supported by MEXT Quantum Leap Flagship Program (MEXT Q-LEAP) Grant Number JPMXS0120319794. This paper is partly based on results obtained from a project, JPNP16007, commissionedby the New Energy and Industrial Technology Development Organization (NEDO), Japan.

A Statistical properties of ν
Proof that ν is an unbiased estimator of Tr(Hρ)
In this Section, we show that ν is the unbiased estimator of Tr(Hρ). When measurements are chosenprobabilistically,E[νr ] =

K
X

aj

j=1

=

K
X
j=1

=

K
X

aj

1
q(Pj )
1
q(Pj )X

f (Pj , Mr )qm (Mr )X

Mr ∈{X,Y,Z}⊗nX
Mrp(ρ, Mr ; b)µ(Pj , b)bf (Pj , Mr )qm (Mr )Tr(Pj ρ)(A.1)∈{X,Y,Z}⊗naj Tr(Pj ρ)

j=1= Tr(Hρ),

Accepted in Quantum 2023-03-28, click title to verify. Published under CC-BY 4.0.

27

where we use (3.9) in the second equality, and we use (3.5) in the third equality. Thus, E[ν] =Tr(Hρ).
When measurements are chosen deterministically,E[νr ] =

K
X

aj

j=1

=

r=1

E[νr ]/Nshot =X1f (Pj , Mr )p(ρ, Mr ; b)µ(Pj , b)q(Pj )
b

K
X

aj

1
f (Pj , Mr )Tr(Pj ρ)q(Pj )

aj

Nshot f (Pj , Mr )Tr(Pj ρ),
PNshot0r 0 =1 f (Pj , Mr )j=1=

PNshot

K
X
j=1

(A.2)

where we use (3.6) in the last equality. Therefore,
N
shot
X

1

E[ν] =

E[νr ]
Nshot r=1
PNshotKX
f (Pj , Mr )Tr(Pj ρ)=
aj r=1
PNshot0r 0 =1 f (Pj , Mr )j=1
X
aj Tr(Pj ρ)
=

(A.3)

j=1

= Tr(Hρ).
Variance of ν
The variance of ν is defined by2

12

Var(ν) = E[ν ] − (E[ν]) =

Nshot

1

N
shot
X

Nshot r=1!
E[νr2 ] − [Tr(Hρ)]2

.

(A.4)

When measurements are probabilistically chosen,E[νr2 ] =K
K X
X
j=1 `=1

X

aj a`

Mr ∈{X,Y,Z}f (Pj , Mr ) f (P` , Mr )qm (Mr )q(Pj )q(P` )⊗nX

p(ρ, Mr ; br )µ(Pj , br )µ(P` , br ).

br ∈{−1,1}⊗n(A.5)
As long as Mr covers both Pj and P`X
brp(ρ, Mr ; br )µ(Pj , br )µ(P` , br ) = Tr(Pj P` ρ).

(A.6)∈{−1,1}⊗n

Substituting (A.6) to (A.5), we obtainE[νr2 ] =

K X
K
X

X

aj a`

j=1 `=1

Mr ∈{X,Y,Z}f (Pj , Mr ) f (P` , Mr )qm (Mr )Tr(Pj P` ρ).
q(Pj )q(P` )⊗n(A.7)

By using f (Pj , Mr )2 = f (Pj , Mr ),X
Mr ∈{X,Y,Z}⊗nf (Pj , Mr )q(Pj )2
qm (Mr ) =Accepted in Quantum 2023-03-28, click title to verify. Published under CC-BY 4.0.

1,
q(Pj )(A.8)

28

and therefore,E[νr2 ] =

K
X
X
a2j
+
aj a`
q(Pj )
j=1
j6=`

X
Mr ∈{X,Y,Z}f (Pj , Mr ) f (P` , Mr )qm (Mr )Tr(Pj P` ρ).
q(Pj )q(P` )⊗n(A.9)

As a result,

K
X
X
a2j

Var(ν) =+aj a`
Nshot j=1 q(Pj )
1

j6=k

X
Mr ∈{X,Y,Z}f (Pj , Mr ) f (P` , Mr )qm (Mr )Tr(Pj P` ρ) − [Tr(Hρ)]2  .
q(Pj )q(P` )⊗n(A.10)In the deterministic case, by using (A.6),
1

N
shot
X

Nshot r=1

E[νr2 ] =

K X
K
X

1

Nshot j=1

aj a`

N
shot
X

f (Pj , Mr ) f (P` , Mr )Tr(Pj P` ρ).
q(Pj )q(P` )
r=1`=1(A.11)

By using a calculation similar to (A.8), we can separate the diagonal term and obtain


N
K
shot
X
1 X a2j
1 X
f (Pj , Mr ) f (P` , Mr )Var(ν) =+aj a`Tr(Pj P` ρ) − [Tr(Hρ)]2  . (A.12)
Nshot j=1 q(Pj ) Nshot
q(P
)
q(P
)
j
`
r=1
j6=`

B Variance in classical shadow
Here we derive the variance of ν in the classical shadow [17]. From (3.11), we obtainX1f (Pj , Mr )f (P` , Mr )
g(Pj , P` ) = 3locality(Pj ) 3locality(P` )  n
3
⊗n

(B.1)

Mr ∈{X,Y,Z}= 3locality(Pj ) 3locality(P` )

n
Y
1
i=1

3

X

flocal (Pj [i], Mr [i])flocal (P` [i], Mr [i]).

(B.2)

Mr [i]∈{X,Y,Z}

We can compute the value inside the product as

1



 1/3

Pj [i] = P` [i] = I(Pj [i] = I or P` [i] = I) and Pj [i] 6= P` [i]1flocal (Pj [i], Mr [i])flocal (P` [i], Mr [i]) =
.

3
1/3
Pj [i] = P` [i] 6= IMr [i]∈{X,Y,Z}
0 otherwise(B.3)
X

Also, we can expand
3locality(Pj ) 3locality(P` ) =nYclocal (Pj [i], P` [i]),

(B.4)

i=1

where



 1 Pj [i] = P` [i] = I
clocal (Pj [i], P` [i]) =3 (Pj [i] = I or P` [i] = I) and Pj [i] 6= P` [i] .

9 Pj [i] = P` [i] 6= I

(B.5)

By substituting (B.3) and (B.4) into (B.2), we obtaing(Pj , P` ) = gcs (Pj , P` ) ≡nYglocal (Pj [i], P` [i]),

(B.6)

i=1

Accepted in Quantum 2023-03-28, click title to verify. Published under CC-BY 4.0.

29

where 1glocal (Pj [i], P` [i]) =3
0

Pj [i] = I or P` [i] = I.
Pj [i] = P` [i] 6= Iotherwise(B.7)

Consequently,
Var(ν) =

1
Nshot



K
X
X

a2j 3locality(Pj ) +aj a` gcs (Pj , P` )Tr(Pj P` ρ) − [Tr(Hρ)]2  .
j=1

(B.8)

j6=`

C Derivation of the bound for shot noise
In [9], the following infidelity
DI



E2
E r
D
E
~ ))
~
~
v(θideal (T )) , v(θ(T )) = 1 − v(θ~ideal (T )) v(θ(Tis evaluated when not using the shadow-based strategy. In this section, we firstly show the value of thepreviously computed DI and next, we evaluate DI when using the shadow-based strategy.
Suppose that there are N -steps for the update of the parameters, i.e., N δt = T , and let us writerr
the state evolvedE by (4.24) and the ideal state
E evolved by (4.23) at r-th step as |v i and |videal i. Then,~ )) and v N= v(θ~ideal (T )) . We also assume that |v 0 i = |v 0 i. By using the triangle
v N = v(θ(T
ideal

ideal

inequality, we obtainr−1rDI (|videal
i , |v r i) = DI Ur videal, |v r i
r−1≤ DI Ur videal, Ur v r−1 + DI Ur v r−1 , |v r ir−1= DI videal, v r−1 + DI (|v0r i , |v r i) ,(C.1)

for r ≥ 1, where Ur is the unitary transformation correspoinding to the r-th ideal update of the parametersby (4.23), and |v0r i ≡ Ur v r−1 . By applying (C.1) repeatedly, we obtainNvideal
, vN

DI



≤

N
X

DI (|v0r i , |v r i) .

(C.2)

r=1

Let
δ θ̇~r ≡~∂ θ(t)∂t−t=rδt∂ θ~ideal (t)
∂t

.

(C.3)

t=rδt

Then, the direct calculation shows [9]DI (|v0r i, |v r i) =q

δ(θ~˙r )T W r δ θ~˙r (δt)2 + O((δt)3 ),(C.4)

where each element of the matrix W r is given byr
Wjk=∂ v r−1 ∂ v r−1∂ v r−1
−
∂θj
∂θk
∂θj

r−1

v r−1

v r−1

∂ v r−1
.
∂θk

(C.5)

r−1

∂
Here, we use ∂hv∂t | |v r−1 i = −hv r−1 | ∂|v∂t i , which is derived from ∂thv r−1 |v r−1 i = 0. By ignoring the3
O((δt) ) term inside the square root function, we obtainDIN
videal

, v

N



.

N
X

q
δt δ(θ~˙r )T W r δ θ~˙r
(C.6)

r=1

≤T

p

~
||W ||max ||δ θ̇||max ,Accepted in Quantum 2023-03-28, click title to verify. Published under CC-BY 4.0.

30

~˙
~˙ max = max0≤t≤T ||δ θ(t)||where ||W ||max = maxr ||W r ||F and ||δ θ||
2.
~
The value δ θ̇(t) is related with the amount of the shot noise. By substituting (M + δM )−1 = M −1 −M −1 δM M −1 + O((δM )2 ) in (4.24), we obtain~∂ θ(t)~ + M −1 δ V~ − M −1 δM M −1 V~ ),' (M −1 V
∂t

(C.7)

where we ignore the higher order terms of δM and δV with the assumption that δM and δV are sufficiently
small. Then,~|δ θ̇(t)|=~∂ θ(t)~− M −1 V∂t(C.8)~ − M −1 δM M −1 V~
' M −1 δ V~ ||2 + ||M −1 ||2F ||δM ||F ||V~ ||2 ,
≤ ||M −1 ||F ||δ V
where || · ||F is the Frobenius norm and || · ||2 is the two-norm. Therefore, let~ ||2 + ||M −1 ||2F ||δM ||F ||V~ ||2 ,
∆M V ≡ ||M −1 ||F ||δ V

(C.9)

then
N
, vN
videal

DI



.T

p
||W ||max ||∆M V ||max .

(C.10)~ as we see in the following.
The value ∆M V depends on how to estimate M and V
The infidelity in the case of using naive strategy
In case we naively measure each term of M and V one by one, from the central limit theorem,
pP
pP
naive )2
naive )2k` (∆Mk`k (∆Vk~ ||2 '√√||δM ||F ', ||δ V,
NnaiveNnaive(C.11)

where we use (4.21). Therefore, we obtain∆naive,∆M V ' √Nnaivewhere∆naive = ||M

−1

||F

sX

(C.12)~ ||2(∆Vknaive )2 + ||M −1 ||2F ||VsXknaive )2 .
(∆Mk`

(C.13)

k`Thus, we obtainDIN
videal, vN. DInaive , with DInaive = Tpp||W ||max ∆max
naive / Nnaive ,(C.14)

where ∆naive
max ≡ ||∆naive ||max .
The infidelity in the case of using the shadow-based strategy
Similarly, in case of shadow-based measurements, we obtain
pP
pP
2(∆Vk )2k` (∆Mk` )~ ||2 '√√k, ||δ V,||δM ||F '
Nshot
Nshot

(C.15)

from the central limit theorem, where we use (4.19) and (4.20). Thus,∆shadow∆M V ' √,Nshot

Accepted in Quantum 2023-03-28, click title to verify. Published under CC-BY 4.0.

(C.16)

31

where∆shadow = ||M −1 ||F

sX

~ ||2
(∆Vk )2 + ||M −1 ||2F ||VksX(∆Mk` )2 .

(C.17)k`

As a result,DIN
videal, vN. DIshadow with DIshadow = Tpp||W ||max ∆max
shadow / Nshot ,(C.18)

where ∆shadow≡ max0≤t≤T (∆shadow ). From the equations (C.14) and (C.18), we see that DInaive =maxshadowDI
holds when max 2
Nshot
∆shadow,(C.19)'
Nnaive∆max
naive
where we assume that the value of ||W ||max is almost the same in the naive strategy and the shadow-based
strategy.
The infidelity in the case of using the hybrid strategy
In case of the hybrid strategy, where we use in case NBB = 1, the number of measurements for estimatingeach Mk` can be smaller than that for estimating Vk (=Nshot ); we set the number of measurements as
Nshot α with 0 < α < 1. Then, in the hybrid strategy,
pP
pP
naive )2(∆Vk )2k` (∆Mk`~√√k, ||δ V ||2 ',(C.20)||δM ||F 'Nshot αNshot
from the central limit theorem. We obtain∆hybrid∆M V ' √
,
Nshot
where
∆hybrid = ||M(C.21)rPsX
naive )2
−1 2 ~k` (∆Mk`
2(∆Vk ) + ||M ||F ||V ||2.
||F
α(C.22). DIhybrid with DIhybrid = T(C.23)

−1

k

Then,
DI

N
videal
, vN



p

p
||W ||max ∆max
hybrid / Nshot ,≡ max0≤t≤T (∆hybrid ). From the equations (C.14) and (C.23), we can show that DInaive =where ∆hybridmax
hybridDIholds when max 2
∆hybrid
Nshot',(C.24)
Nnaive∆max
naive
where we assume that the value of ||W ||max is almost the same in the naive strategy and the hybrid
strategy.

D Derivation of the Haar averages
In this section, we derive results of the Haar integration. We use the following itegration formulae [51]:Z
Tr(A)Tr(B)dU Tr(AU BU † ) =
,
(D.1)
N
Haar
Z
Tr(A)Tr(B)Tr(C)Tr(D) + Tr(AC)Tr(BD)dU Tr(AU BU † )Tr(CU DU † ) =,(D.2)
N (N + 1)
HaarZ
4
6Tr(A4 ) + 3[Tr(A2 )]2dU Tr(AU |0ih0|U † ) = 4,(D.3)
N + 6N 3 + 11N 2 + 6N
Haar
where N is the dimension of the Hilbert space of U .

Accepted in Quantum 2023-03-28, click title to verify. Published under CC-BY 4.0.

32

Derivation of (4.33)
By using (D.1) with A = P̃r P̃r0 and B = |0ih0|, we obtainhVki iHaar =NBA
X(Gk,i,r )2+q(Pr )
r=1

Z
dU
Haar

N
BA X
X

Gk,i,r Gk,i,r0 g(P̃r , P̃r0 )Tr(P̃r P̃r0 U |0ih0|U † ),r=1 r 0 6=r

(D.4)

NBA X
Tr(P̃r P̃r0 )(Gk,i,r )2 X.
+
Gk,i,r Gk,i,r0 g(P̃r , P̃r0 )=q(Pr )
2n+1
r=1 0
r=1
N
BA
X

r 6=r

Since P̃r 6= P̃r0 , it holds P̃r P̃r0 6= I, and therefore, Tr(P̃r P̃r0 ) = 0. As a result,hVki iHaar =NBA
X(Gk,i,r )2.
q(Pr )
r=1(D.5)

By using (D.2) with A = P̃r P̃r0 , B = P̃u P̃u0 , and C, D = |0ih0|.


N
BA
X

Var(Vki )Haar = 

2

N
BA X
X

(Gk,i,r )+q(Pr )
r=1 0
r=1

2
† 

Gk,i,r Gk,i,r0 g(P̃r , P̃r0 )Tr(P̃r P̃r0 U |0ih0|U )

r 6=r

=

N
BA X
BA X NX
X

N
BA X NBA X
X
X(Gk,i,r )2q(Pr )
r=1!2

Gk,i,r Gk,i,r0 Gk,i,u Gk,i,u0 g(P̃r , P̃r0 )g(P̃u , P̃u0 )Tr(P̃r P̃r0 )Tr(P̃u P̃u0 ) + Tr(P̃r P̃r0 P̃u P̃u0 )2n+1 (2n+1 + 1)

Gk,i,r Gk,i,r0 Gk,i,u Gk,i,u0 g(P̃r , P̃r0 )g(P̃u , P̃u0 )Tr(P̃r P̃r0 P̃u P̃u0 ),2n+1 (2n+1 + 1)

r=1 r 0 6=r u=1 u0 6=u

=

−

N
BA
X

r=1 r 0 6=r u=1 u0 6=u(D.6)
where in the last equality, we use Tr(P̃r P̃r0 ) = 0.
Derivation of (4.38)
By using (D.2), with A, C = P̃r and B, D = |0ih0|,Z2 hGki,r iHaar = (Gk,i,r )2
dU 1 − Tr(P̃r U |0ih0|U † )Haar!
2(Tr(P̃r ))2 + (Tr(P˜r ))2= (Gk,i,r ) 1 −2n+1 (2n+1 + 1)


1
2
.
= (Gk,i,r ) 1 − n+1
2
+1

(D.7)

By using (D.3) with A = P̃r ,
Z

2 2
i,r
4
Var(Gk ) = (Gk,i,r )dU 1 − Tr(P̃r U |0ih0|U † )− (hGki,r iHaar )2Haar
6Tr(P̃r4 ) + 3[Tr(P̃r2 )]2
2 · 2n+1 − 1= (Gk,i,r )4 4(n+1)+
2+ 6 · 23(n+1) + 11 · 22(n+1) + 6 · 2n+1(22(n+1) + 1)2
=

(D.8)

(Gk,i,r )4 (2 · 22(n+1) + 8 · 22(n+1) ).
(2n+1 + 1)2 (2n+1 + 3)

References[1] John Preskill. “Quantum computing in the NISQ era and beyond”. Quantum 2, 79 (2018).

Accepted in Quantum 2023-03-28, click title to verify. Published under CC-BY 4.0.

33

[2] 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, 032001 (2021).
[3] Marco 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”. Nature Reviews Physics 3, 625–644 (2021).
[4] Edward Farhi, Jeffrey Goldstone, and Sam Gutmann. “A quantum approximate optimization algorithm” (2014). arXiv:1411.4028.
[5] Alberto Peruzzo, Jarrod McClean, Peter Shadbolt, Man-Hong Yung, Xiao-Qi Zhou, Peter J Love,
Alán Aspuru-Guzik, and Jeremy L O’brien. “A variational eigenvalue solver on a photonic quantum
processor”. Nature communications 5, 1–7 (2014).
[6] Abhinav Kandala, Antonio Mezzacapo, Kristan Temme, Maika Takita, Markus Brink, Jerry Chow,and Jay Gambetta. “Hardware-efficient variational quantum eigensolver for small molecules andquantum magnets”. Nature 549, 242–246 (2017).
[7] Nikolaj Moll, Panagiotis Barkoutsos, Lev S Bishop, Jerry M Chow, Andrew Cross, Daniel J Egger,
Stefan Filipp, Andreas Fuhrer, Jay M Gambetta, Marc Ganzhorn, et al. “Quantum optimization
using variational algorithms on near-term quantum devices”. Quantum Science and Technology 3,030503 (2018).
[8] Jarrod R McClean, Jonathan Romero, Ryan Babbush, and Alán Aspuru-Guzik. “The theory ofvariational hybrid quantum-classical algorithms”. New Journal of Physics 18, 023023 (2016).
[9] Ying Li and Simon C Benjamin. “Efficient variational quantum simulator incorporating active error
minimization”. Physical Review X 7, 021050 (2017).
[10] Kosuke Mitarai, Makoto Negoro, Masahiro Kitagawa, and Keisuke Fujii. “Quantum circuit learning”.
Physical Review A 98, 032309 (2018).
[11] Raphael Kaubruegger, Pietro Silvi, Christian Kokail, Rick van Bijnen, Ana Maria Rey, Jun Ye,
Adam M Kaufman, and Peter Zoller. “Variational spin-squeezing algorithms on programmable quantum sensors”. Physical review letters 123, 260505 (2019).
[12] Bálint Koczor, Suguru Endo, Tyson Jones, Yuichiro Matsuzaki, and Simon C Benjamin. “Variationalstate quantum metrology”. New Journal of Physics 22, 083038 (2020).
[13] Sam McArdle, Tyson Jones, Suguru Endo, Ying Li, Simon C Benjamin, and Xiao Yuan. “Variational
ansatz-based quantum simulation of imaginary time evolution”. npj Quantum Information 5, 1–
6 (2019).
[14] Xiao Yuan, Suguru Endo, Qi Zhao, Ying Li, and Simon C Benjamin. “Theory of variational quantum
simulation”. Quantum 3, 191 (2019).
[15] AD McLachlan. “A variational solution of the time-dependent Schrodinger equation”. Molecular
Physics 8, 39–44 (1964).
[16] Vladyslav Verteletskyi, Tzu-Ching Yen, and Artur F Izmaylov. “Measurement optimization in thevariational quantum eigensolver using a minimum clique cover”. The Journal of chemical physics
152, 124114 (2020).
[17] Hsin-Yuan Huang, Richard Kueng, and John Preskill. “Predicting many properties of a quantum
system from very few measurements”. Nature Physics 16, 1050–1057 (2020).
[18] Charles Hadfield, Sergey Bravyi, Rudy Raymond, and Antonio Mezzacapo. “Measurements of quantum hamiltonians with locally-biased classical shadows”. Communications in Mathematical Physics
391, 951–967 (2022).
[19] Hsin-Yuan Huang, Richard Kueng, and John Preskill. “Efficient estimation of pauli observables byderandomization”. Physical review letters 127, 030503 (2021).
[20] Stefan Hillmich, Charles Hadfield, Rudy Raymond, Antonio Mezzacapo, and Robert Wille. “Decision
diagrams for quantum measurements with shallow circuits”. In 2021 IEEE International Conferenceon Quantum Computing and Engineering (QCE). Pages 24–34. IEEE (2021).

Accepted in Quantum 2023-03-28, click title to verify. Published under CC-BY 4.0.

34

[21] Bujiao Wu, Jinzhao Sun, Qi Huang, and Xiao Yuan. “Overlapped grouping measurement: A unified
framework for measuring quantum states” (2021).
[22] Pranav Gokhale, Olivia Angiuli, Yongshan Ding, Kaiwen Gui, Teague Tomesh, Martin Suchara,
Margaret Martonosi, and Frederic T Chong. “Minimizing state preparations in variational quantum
eigensolver by partitioning into commuting families” (2019). arXiv:1907.13623.
[23] Artur F Izmaylov, Tzu-Ching Yen, Robert A Lang, and Vladyslav Verteletskyi. “Unitary partitioning
approach to the measurement problem in the variational quantum eigensolver method”. Journal ofchemical theory and computation 16, 190–195 (2019).
[24] Ophelia Crawford, Barnaby van Straaten, Daochen Wang, Thomas Parks, Earl Campbell, andStephen Brierley. “Efficient quantum measurement of Pauli operators in the presence of finite sampling error”. Quantum 5, 385 (2021).
[25] William J Huggins, Jarrod R McClean, Nicholas C Rubin, Zhang Jiang, Nathan Wiebe, K Birgitta
Whaley, and Ryan Babbush. “Efficient and noise resilient measurements for quantum chemistry on
near-term quantum computers”. npj Quantum Information 7, 1–9 (2021).
[26] Ikko Hamamura and Takashi Imamichi. “Efficient evaluation of quantum observables using entangled
measurements”. npj Quantum Information 6, 1–8 (2020).
[27] Sergey Bravyi, Jay M Gambetta, Antonio Mezzacapo, and Kristan Temme. “Tapering off qubits tosimulate fermionic Hamiltonians” (2017). arXiv:1701.08213.
[28] 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, 062322 (2020).
[29] Tzu-Ching Yen, Vladyslav Verteletskyi, and Artur F Izmaylov. “Measuring all compatible operators
in one series of single-qubit measurements using unitary transformations”. Journal of chemical theory
and computation 16, 2400–2409 (2020).
[30] Andrew Jena, Scott Genin, and Michele Mosca. “Pauli partitioning with respect to gate sets” (2019).
arXiv:1907.07859.
[31] Tzu-Ching Yen, Aadithya Ganeshram, and Artur F Izmaylov. “Deterministic improvements of quantum measurements with grouping of compatible operators, non-local transformations, and covarianceestimates” (2022).
[32] Seonghoon Choi, Tzu-Ching Yen, and Artur F Izmaylov. “Improving Quantum Measurements by
Introducing “Ghost” Pauli Products”. Journal of Chemical Theory and Computation 18, 7394–
7402 (2022).
[33] Seonghoon Choi, Ignacio Loaiza, and Artur F Izmaylov. “Fluid fermionic fragments for optimizing
quantum measurements of electronic hamiltonians in the variational quantum eigensolver”. Quantum
7, 889 (2023).
[34] Andrew Arrasmith, Lukasz Cincio, Rolando D Somma, and Patrick J Coles. “Operator sampling forshot-frugal optimization in variational algorithms” (2020). arXiv:2004.06252.
[35] Gregory Boyd and Bálint Koczor. “Training variational quantum circuits with CoVaR: covarianceroot finding with classical shadows” (2022).
[36] Ting Zhang, Jinzhao Sun, Xiao-Xu Fang, Xiao-Ming Zhang, Xiao Yuan, and He Lu. “Experimental
quantum state measurement with classical shadows”. Physical Review Letters 127, 200501 (2021).
[37] Suguru Endo, Jinzhao Sun, Ying Li, Simon C Benjamin, and Xiao Yuan. “Variational quantum
simulation of general processes”. Physical Review Letters 125, 010501 (2020).
[38] P. A. M. Dirac. “Note on exchange phenomena in the thomas atom”. Mathematical Proceedings of
the Cambridge Philosophical Society 26, 376–385 (1930).
[39] Âkov Il’ič Frenkel. “Wave mechanics; advanced general theory”. Bull. Amer. Math. Soc 41, 776 (1935).
[40] Peter Kramer and Marcos Saraceno. “Geometry of the time-dependent variational principle in quantum mechanics”. In Group Theoretical Methods in Physics: Proceedings of the IX International
Colloquium Held at Cocoyoc, México, June 23–27, 1980. Pages 112–121. Springer (2005).

Accepted in Quantum 2023-03-28, click title to verify. Published under CC-BY 4.0.

35

[41] J Broeckhove, L Lathouwers, E Kesteloot, and P Van Leuven. “On the equivalence of time-dependent
variational principles”. Chem. Phys. Lett. 149, 547–550 (1988).
[42] Artur F Izmaylov, Tzu-Ching Yen, Robert A Lang, and Vladyslav Verteletskyi. “Unitary partitioning
approach to the measurement problem in the variational quantum eigensolver method”. Journal ofchemical theory and computation 16, 190–195 (2019).
[43] 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, 062322 (2020).
[44] Tzu-Ching Yen, Vladyslav Verteletskyi, and Artur F Izmaylov. “Measuring all compatible operators
in one series of single-qubit measurements using unitary transformations”. Journal of chemical theory
and computation 16, 2400–2409 (2020).
[45] Giacomo Torlai, Guglielmo Mazzola, Giuseppe Carleo, and Antonio Mezzacapo. “Precise measurement of quantum observables with neural-network estimators”. Physical Review Research 2,
022060 (2020).
[46] Ophelia Crawford, Barnaby van Straaten, Daochen Wang, Thomas Parks, Earl Campbell, andStephen Brierley. “Efficient quantum measurement of Pauli operators in the presence of finite sampling error”. Quantum 5, 385 (2021).
[47] William J Huggins, Jarrod R McClean, Nicholas C Rubin, Zhang Jiang, Nathan Wiebe, K Birgitta
Whaley, and Ryan Babbush. “Efficient and noise resilient measurements for quantum chemistry on
near-term quantum computers”. npj Quantum Information 7, 1–9 (2021).
[48] Maria Schuld, Ville Bergholm, Christian Gogolin, Josh Izaac, and Nathan Killoran. “Evaluating
analytic gradients on quantum hardware”. Physical Review A 99, 032331 (2019).
[49] Barnaby van Straaten and Bálint Koczor. “Measurement cost of metric-aware variational quantum
algorithms”. PRX Quantum 2, 030324 (2021).
[50] Yasunari Suzuki, Yoshiaki Kawase, Yuya Masumura, Yuria Hiraga, Masahiro Nakadai, Jiabao Chen,
Ken M Nakanishi, Kosuke Mitarai, Ryosuke Imai, Shiro Tamiya, et al. “Qulacs: a fast and versatilequantum circuit simulator for research purpose”. Quantum 5, 559 (2021).
[51] Benoı̂t Collins and Piotr Śniady. “Integration with respect to the Haar measure on unitary, orthogonal
and symplectic group”. Communications in Mathematical Physics 264, 773–795 (2006).

Accepted in Quantum 2023-03-28, click title to verify. Published under CC-BY 4.0.

36