# Quantum gradient evaluation through quantum non-demolition measurements

P. Solinas  e-mail<EMAIL_ADDRESS>Dipartimento di Fisica, Università
di Genova, via Dodecaneso 33, I-16146, Genova, Italy  INFN - Sezione di
Genova, via Dodecaneso 33, I-16146, Genova, Italy  S. Caletti  e-mail:
<EMAIL_ADDRESS>Dipartimento di Fisica, Università di Genova, via
Dodecaneso 33, I-16146, Genova, Italy  INFN - Sezione di Genova, via
Dodecaneso 33, I-16146, Genova, Italy  G. Minuto  e-mail:
<EMAIL_ADDRESS>Dept. of Informatics, Bioengineering, Robotics,
and Systems Engineering (DIBRIS), Polytechnic School of Genoa University,
Genova, Italy  INFN - Sezione di Genova, via Dodecaneso 33, I-16146, Genova,
Italy

###### Abstract

We discuss a Quantum Non-Demolition Measurement (QNDM) protocol to estimate
the derivatives of a cost function with a quantum computer. The cost function,
which is supposed to be classically hard to evaluate, is associated with the
average value of a quantum operator. Then a quantum computer is used to
efficiently extract information about the function and its derivative by
evolving the system with a so-called variational quantum circuit. To this aim,
we propose to use a quantum detector that allows us to directly estimate the
derivatives of an observable, i.e., the derivative of the cost function. With
respect to the standard direct measurement approach, this leads to a reduction
of the number of circuit iterations needed to run the variational quantum
circuits. The advantage increases if we want to estimate the higher-order
derivatives. We also show that the presented approach can lead to a further
advantage in terms of the number of total logical gates needed to run the
variational quantum circuits. These results make the QNDM a valuable
alternative to implementing the variational quantum circuits.

##  I Introduction

The advent of all-purpose quantum computers able to solve hard computational
problems is still decades away. However, it is commonly believed that some
problems intractable with a classical computer could be within reach of
today’s Noisy Intermediate-Scale Quantum (NISQ) computers.

The problems for which a quantum advantage can be reached in the next years
are the ones that require a large working space such as the simulation of
complex physical and chemical systems. The most promising architecture is a
hybrid quantum-classical one. Among these hybrid algorithms the most relevant
ones are the variational quantum circuit, the variational quantum eigensolver
Peruzzo  _et al._ (2014); Kandala  _et al._ (2017) and the quantum approximate
optimization algorithm Farhi  _et al._ (2014). The computational scheme is the
following. Given a certain cost function in a large parameter space to be
minimized, we associate it with a Hermitian operator. We run a quantum circuit
to evaluate the average values of the Hermitian operator, i.e., the value of
the cost function, in a specified point of the parameter space. Since the
quantum measurements are probabilistic, we need to iterate the process to
reach the desired accuracy. This information is then fed into a classical
computer which elaborates it and determines the following steps of the quantum
computer. Typically, classical computation consists in an optimization
algorithm for which we need information about the derivatives of the cost
functions. Therefore, in these schemes, the main quantum computational task is
to obtain the derivatives of the cost function by measuring quantum
observables with accuracy and minimal cost. From the quantum perspective, the
main resource costs come from the iterations needed to have an accurate
average and the logical gates needed to run the quantum circuit at every
iteration. Because of the limitation in quantum hardware and quantum
operations, any reduction of the cost can bring us closer to obtaining a
quantum advantage in problems of practical interest.

While the usual proposals rely on the direct measurement (DM) of a quantum
observable to extract information about the cost function McArdle  _et al._
(2020); Cerezo  _et al._ (2021); Mari  _et al._ (2021), here, we discuss an
alternative method to estimate directly the derivatives of the observable. A
quantum detector is coupled in sequence with the system from which we want to
extract the information. The information about the observable, its gradient,
or the high derivatives is stored in the detector phase which is eventually
measured. This technique is often called Quantum Non-Demolition Measurement
(QNDM) Clerk (2011) or full counting statistics Bednorz  _et al._ (2012) and
it is rooted in the idea of weak values and weak measurements Aharonov  _et
al._ (1988).

The potential advantages of this approach lie in the fact that, with this
unconventional measurement, we can directly estimate the average of the
variation, i.e., the gradient, of quantum observables which cannot be obtained
with direct measurement. Indeed, the same approach has been used to estimate
the variation of charge Bednorz  _et al._ (2012) and energy Solinas and
Gasparinetti (2015, 2016); Solinas  _et al._ (2021, 2022) in quantum systems.
Since these are related to observables at different times, they are not
associated with any hermitian operator Talkner  _et al._ (2007) and the
measurement process suffers from conceptual and practical subtleties Solinas
and Gasparinetti (2015); Solinas  _et al._ (2022). In this paper, we revert
the process and identify in the QNDM the ideal approach to extract information
about the derivative of an observable.

A similar approach was presented in Refs. O’Brien  _et al._ (2019); Banchi and
Crooks (2021). Here, we discuss it in a more general framework, extend its
applications to estimate high derivatives and discuss the main advantages with
respect to the direct measurement approach. As we show below, the QNDM needs
fewer resources (in terms of iterations and logical operators) than the DM
approach. The advantages increase with the order of the derivatives to be
estimated.

The paper is divided as follows. In Sec. II, we introduce the quantities to be
estimated. In Sec. III, we show how the gradient of an observable can be
measured with the QNDM approach and, in Sec. IV we give an explicit example of
the advantages of the QNDM approach when the cost function is related to a
complex operator. In Sec. V, we extend our framework to the estimate of second
and higher derivative and in Sec. VI we discuss the advantages in terms of
resources needed. Finally, Sec. VII contains the conclusions.

##  II Definition of the fundamental quantities

Let’s suppose to have a quantum state $\mathinner{|{\psi_{0}}\rangle}$ and we
let it evolve with a unitary operator $U(\theta)$ which depends on $m$
parameters ${\vec{\theta}}=(\theta_{1},\theta_{2},....,\theta_{m})$. The final
state is
$\mathinner{|{\psi({\vec{\theta}})}\rangle}=U({\vec{\theta}})\mathinner{|{\psi_{0}}\rangle}$.
In a more physical (and less abstract) framework, we can consider the operator
$U({\vec{\theta}})$ generated by the hamiltonian $H({\vec{\theta}})$ depending
on the $m$ parameters which can depend on time
${\vec{\theta}}(t)=(\theta_{1}(t),\theta_{2}(t),....,\theta_{m}(t))$. In this
case, the evolution is given by
$U({\vec{\theta}})=T\exp\\{-i\int_{0}^{\tau}H({\vec{\theta}}(t))dt\\}$ where
$T$ is the time-ordered product.

To make a specific example for quantum computers, we can consider the system
composed of $n$ qubits. The initial state is
$\mathinner{|{\psi_{0}}\rangle}=\mathinner{|{0}\rangle}=\mathinner{|{00....0}\rangle}$,
and the unitary transformation $U(\theta)$ can be obtained as a sequence of
quantum logical gates.

We assume that the $U({\vec{\theta}})$ operator can be written as

$U({\vec{\theta}})=V_{m}U_{m}({\vec{\theta}_{m}})....V_{2}U_{2}({\vec{\theta}_{2}})V_{1}U_{1}({\vec{\theta}_{1}})$
(1)

where $V_{m}$ are arbitrary transformation (independent of $\vec{\theta}$) and
$U_{m}({\vec{\theta}_{j}})=\exp{\\{-i\theta_{j}H_{j}\\}}$ is a single qubit
operator generated by the Hermitian operator $H_{j}$ Mari  _et al._ (2021).
Supposing that $H_{j}^{2}=\mathds{1}$,

$U_{j}({\vec{\theta}_{j}})=e^{-iH_{j}\theta_{j}/2}=\cos\frac{\theta_{j}}{2}\mathds{1}-i\sin\frac{\theta_{j}}{2}H_{j}.$
(2)

A more general situation where the $H_{j}$s do not satisfy the above condition
is discussed in Ref. Banchi and Crooks (2021).

The quantity we are interested in is the average value of another operator
$\hat{M}$ over the final state $\mathinner{|{\psi({\vec{\theta}})}\rangle}$
Mari  _et al._ (2021)

$f({\vec{\theta}})={\mathinner{\langle{\psi({\vec{\theta}})}|{\hat{M}}|{\psi({\vec{\theta}})}\rangle}}={\mathinner{\langle{0}|{U^{\dagger}({\vec{\theta}})\hat{M}U({\vec{\theta}})}|{0}\rangle}}.$
(3)

Measuring this quantity for two points, we can estimate the derivative of
$f({\vec{\theta}})$ as Mari  _et al._ (2021)

$g_{j_{1}}=\frac{\partial
f({\vec{\theta}})}{\partial\theta_{j_{1}}}=\frac{f(\vec{\theta}+s{\bf
e}_{j_{1}})-f(\vec{\theta}-s{\bf e}_{j_{1}})}{2\sin s}.$ (4)

where ${\bf e}_{j_{1}}$ is the versor along the $\theta_{j_{1}}$ direction.
This method can be generalized to calculate the second derivative (see Eq.
($10$) in Mari  _et al._ (2021))

$\displaystyle g_{j_{1},j_{2}}$ $\displaystyle=$
$\displaystyle\frac{\partial^{2}f({\vec{\theta}})}{\partial\theta_{j1}\partial\theta_{j2}}=$
$\displaystyle=$ $\displaystyle\Big{[}f(\vec{\theta}+s({\bf e}_{j_{1}}+{\bf
e}_{j_{2}}))-f(\vec{\theta}+s(-{\bf e}_{j_{1}}+{\bf e}_{j_{2}}))$
$\displaystyle-$ $\displaystyle f(\vec{\theta}+s({\bf e}_{j_{1}}-{\bf
e}_{j_{2}}))+f(\vec{\theta}-s({\bf e}_{j_{1}}+{\bf
e}_{j_{2}}))\Big{]}\Big{[}2\sin^{2}s\Big{]}^{-1}$

or even higher derivatives defined as Mari  _et al._ (2021)

$g_{j_{1},j_{2},...,j_{d}}=\frac{\partial^{d}f(\vec{\theta})}{\partial\theta_{j1}\partial\theta_{j_{2}}...\partial\theta_{j_{d}}}.$
(6)

The information about the derivatives of $f$ can be used to minimize a cost
function. To this aim, different optimizer algorithms can be used (see, for
example, Mari  _et al._ (2021)). The most common one is the gradient descent
(GD). This optimiser starts with the parameter $\vec{\theta}^{(0)}$ which are
sequentially updated to new values $\vec{\theta}^{(1)}$, $\vec{\theta}^{(2)}$,
…, $\vec{\theta}^{(T)}$. The update of parameters is obtained by the rule
$\vec{\theta}^{(t)}=\vec{\theta}^{(t-1)}-\eta\nabla f(\vec{\theta}^{(t-1)})$
where $\eta>0$.

The access to the second derivatives allows us to implement other optimization
algorithms. For example, in the Newton optimizer, the parameter update rule is
$\vec{\theta}^{(t)}=\vec{\theta}^{(t-1)}-\eta[Hf(\vec{\theta}^{(t-1)})]^{-1}\nabla
f(\vec{\theta}^{(t-1)})$ where $[Hf(\vec{\theta}^{(t-1)})]^{-1}$ is the
inverse of the Hessian matrix estimated using Eq. (II).

Since the Newton optimizer requires the inversion of the Hessian matrix, it
might be resource-consuming for a large parameter space. For this reason,
alternative second-order approaches have been proposed such as the diagonal
Newton optimizer and quantum natural gradient optimizer Cheng (2010); Liu  _et
al._ (2019); Mari  _et al._ (2021). The latter uses the information encoded in
the Fubini-Study metric tensor that can be obtained with similar methods Mari
_et al._ (2021).

##  III Quantum non-demolition measurement: First order derivative - Gradient

The key idea of the QNDM approach is to couple the system to a quantum
detector and store the desired information in the phase of the detector that
is eventually measured O’Brien  _et al._ (2019); Banchi and Crooks (2021). By
choosing properly the system detector interaction Solinas and Gasparinetti
(2015, 2016); Solinas  _et al._ (2017), we can store in the detector the
information of the variation of the average $M$ so that a final measurement
gives us direct access to the gradient (4).

Notice that in terms of operators we can write Eq. (4) as the average of a new
operator $U^{\dagger}(\vec{\theta}+s{\bf
e}_{j_{1}})\hat{M}U(\vec{\theta}+s{\bf
e}_{j_{1}})-U^{\dagger}(\vec{\theta}-s{\bf
e}_{j_{1}})\hat{M}U(\vec{\theta}-s{\bf e}_{j_{1}})=\hat{M}(\vec{\theta}+s{\bf
e}_{j_{1}})-\hat{M}(\vec{\theta}-s{\bf e}_{j_{1}})$. However, being this the
difference between observables, it is not a quantum observable and it cannot
be measured directly but the contributions must be measured separately Talkner
_et al._ (2007); Solinas and Gasparinetti (2015, 2016). The QNDM approach
allows us to access the same information about the averages with direct
measurement of the detector phase, therefore, reducing the resources needed.

More precisely, let us consider the same setup discussed in Sec. II but let us
add an additional quantum system that acts as the detector. We suppose that
the detector has no dynamics, i.e., it effectively evolves over times much
longer than the one needed to perform the full protocol.

The system and the detector interact with the Hamiltonian

$H_{I}=h(t)\lambda\hat{p}\otimes\hat{M}.$ (7)

Here, $\lambda$ is the coupling constant between the system and the detector,
and it can be changed and fixed at the beginning of the evolution. The
operators $\hat{p}$ and $\hat{M}$ act on the detector and the system,
respectively, and we denote the eigenstates of $\hat{p}$ with
$\mathinner{|{p}\rangle}$, i.e.,
$\hat{p}\mathinner{|{p}\rangle}=p\mathinner{|{p}\rangle}$.

The time-dependent function $h(t)$ determines when the system and the detector
are coupled. If the full protocol occurs between time $t=0$ and $t=T$, we take
$h(t)=\delta(T-t)-\delta(t_{1})$ where $\delta(t)$ is the Dirac delta and
$0<t_{1}<T$. By using the Dirac delta, we model the fact that the system
detector interaction occurs over much smaller time scales than all the other
evolutions so that during the coupling we can consider the system dynamics
“frozen”. Under this assumption, the system-detector evolution is associated
with the unitary evolution $U_{\pm}=\exp\\{\pm
i\lambda\hat{p}\otimes\hat{M}\\}$.

In the intervals between the couplings with the detector, the system evolves
with unitary evolution $U_{i}$ ($i=1,2$) so that the total (system$+$detector)
evolution is

$U_{tot}=e^{i\lambda\hat{p}\otimes\hat{M}}U_{2}e^{-i\lambda\hat{p}\otimes\hat{M}}U_{1}$
(8)

where the explicit form of the $U_{i}$ will be determined below.

We consider as initial state
$\mathinner{|{\psi_{0}}\rangle}=\mathinner{|{0}\rangle}\left(\frac{1}{\sqrt{N}}\sum\mathinner{|{p}\rangle}_{D}\right)$
where the sum is over $N$ detector states and the index $D$ denotes the states
in the detector space. The final state is
$\mathinner{|{\psi_{f}}\rangle}=U_{tot}\mathinner{|{\psi_{0}}\rangle}$
associated to the density operator
$\rho_{tot}={\mathinner{|{\psi_{f}}\rangle\langle{\psi_{f}}|}}$. The final
density matrix of the detector is $\rho_{D}^{f}={\rm Tr}_{S}[\rho_{tot}]$
where ${\rm Tr}_{S}$ denotes the trace over the system degrees of freedom.

We define the quasi-characteristic function as Solinas and Gasparinetti (2015,
2016)

$\mathcal{G}_{\lambda}=\frac{\tensor{[}_{D}]{{\mathinner{\langle{\bar{p}}|{\rho_{D}^{f}}|{-\bar{p}}\rangle}}}{{}_{D}}}{\tensor{[}_{D}]{{\mathinner{\langle{\bar{p}}|{\rho_{D}^{0}}|{-\bar{p}}\rangle}}}{{}_{D}}}.$
(9)

where $\mathinner{|{\pm\bar{p}}\rangle}_{D}$ is a specific eigenstates of
$\hat{p}$. From a physical point of view, the quasi-characteristic function is
the phase accumulated between the states
$\mathinner{|{\pm\bar{p}}\rangle}_{D}$ of the detector during the evolution
Solinas and Gasparinetti (2015, 2016). This can be directly measured with
interferometric techniques.

More explicitly, we can now calculate the numerator of Eq. (9). Since $U_{i}$
does not act on the detector state and $\mathinner{|{p}\rangle}_{D}$ are
eigenstates of $\hat{p}$, we have

${}_{D}\mathinner{\langle{\bar{p}|\psi_{f}}\rangle}$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{N}}\sum~{}_{D}\mathinner{\langle{\bar{p}}|}e^{i\lambda\hat{p}\otimes\hat{M}}U_{2}e^{-i\lambda\hat{p}\otimes\hat{M}}U_{1}\mathinner{|{0}\rangle}\mathinner{|{p}\rangle}_{D}$
(10) $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{N}}\sum\delta_{\bar{p},p}e^{i\lambda
p\hat{M}}U_{2}e^{-i\lambda p\hat{M}}U_{1}\mathinner{|{0}\rangle}$
$\displaystyle=$
$\displaystyle\frac{1}{\sqrt{N}}e^{i\lambda\bar{p}\hat{M}}U_{2}e^{-i\lambda\bar{p}\hat{M}}U_{1}\mathinner{|{0}\rangle}.$

Analogously,

$\displaystyle\mathinner{\langle{\psi_{f}|-\bar{p}}\rangle}_{D}$
$\displaystyle=$
$\displaystyle\frac{1}{\sqrt{N}}\mathinner{\langle{0}|}U_{1}^{\dagger}e^{-i\lambda\bar{p}\hat{M}}U_{2}^{\dagger}e^{i\lambda\bar{p}\hat{M}}.$
(11)

Noticing that
$\tensor{[}_{D}]{{\mathinner{\langle{\bar{p}}|{\rho_{D}^{f}}|{-\bar{p}}\rangle}}}{{}_{D}}={\rm
Tr}_{S}[\tensor{[}_{D}]{{\mathinner{\langle{\bar{p}}|{\rho_{tot}}|{-\bar{p}}\rangle}}}{{}_{D}}]$,
we have

$\displaystyle{\rm
Tr}_{S}[\tensor{[}_{D}]{{\mathinner{\langle{\bar{p}}|{\rho_{tot}}|{-\bar{p}}\rangle}}}{{}_{D}}]$
$\displaystyle=$ $\displaystyle\frac{1}{N}{\rm
Tr}_{S}[\mathcal{U}_{\lambda}\rho_{S}^{0}\mathcal{U}^{\dagger}_{-\lambda}]$
(12)

where we have defined the operator $\mathcal{U}_{\lambda}$ as

$\mathcal{U}_{\lambda}=e^{i\lambda\bar{p}\hat{M}}U_{2}e^{-i\lambda\bar{p}\hat{M}}U_{1}.$
(13)

Finally, we can rewrite the quasi-characteristic function as Solinas and
Gasparinetti (2015, 2016)

$\mathcal{G}_{\lambda}={\rm
Tr}_{S}[\mathcal{U}_{\lambda}\rho_{S}^{0}\mathcal{U}^{\dagger}_{-\lambda}].$
(14)

In statistics, the knowledge of the characteristic function gives access to
all the moments of the distribution since they are related to its derivative
calculated in $\lambda$. Since we are dealing with a quantum system, we have a
quasi-characteristic function and the situation is more tricky since the
associate probability distribution is not positively defined Clerk (2011);
Solinas and Gasparinetti (2015, 2016); Solinas  _et al._ (2017, 2021, 2022),
i.e., it is a quasi-probability distribution. This makes it difficult to
assign a precise meaning to the high moments. However, it can be shown that
the first moment is not affected by these problems, it has a practical and
operative interpretation in terms of observable measurements Solinas and
Gasparinetti (2015, 2016).

Given this premise, we are interested in the derivative

$-i\partial_{\lambda}\mathcal{G}_{\lambda}\Big{|}_{\lambda=0}=-i\,{\rm
Tr}_{S}[\partial_{\lambda}\mathcal{U}_{\lambda}\rho_{S}^{0}\mathcal{U}^{\dagger}_{-\lambda}+\mathcal{U}_{\lambda}\rho_{S}^{0}\partial_{\lambda}\mathcal{U}^{\dagger}_{-\lambda}]\Big{|}_{\lambda=0}.$
(15)

This can be numerically obtained once the behavior of $\mathcal{G}_{\lambda}$
is known close to $\lambda=0$.

Figure 1:  Scheme to perform a measure of the gradient using an ancilla qubit
as a quantum detector. The final Hadamard operation followed by the
measurement corresponds to the measurement of the off-diagonal element of the
density matrix and gives information about the accumulated phase.

From definition (13) we can calculate $\mathcal{U}^{\dagger}_{-\lambda}$ and
the corresponding derivatives calculated in $\lambda=0$:
$\partial_{\lambda}\mathcal{U}_{\lambda}\Big{|}_{\lambda=0}$ and
$\partial_{\lambda}\mathcal{U}^{\dagger}_{-\lambda}\Big{|}_{\lambda=0}$. We
obtain

$\displaystyle\partial_{\lambda}\mathcal{U}_{\lambda}\Big{|}_{\lambda=0}$
$\displaystyle=$ $\displaystyle i\bar{p}(\hat{M}U_{2}U_{1}-U_{2}\hat{M}U_{1})$
$\displaystyle\partial_{\lambda}\mathcal{U}^{\dagger}_{-\lambda}\Big{|}_{\lambda=0}$
$\displaystyle=$ $\displaystyle
i\bar{p}(U^{\dagger}_{1}U^{\dagger}_{2}\hat{M}-U^{\dagger}_{1}\hat{M}U^{\dagger}_{2})$
(16)

and, after algebraic manipulation using the cyclic property of the trace, the
first derivative of $\mathcal{G}_{\lambda}$ reads

$-i\partial_{\lambda}\mathcal{G}_{\lambda}\Big{|}_{\lambda=0}=2\bar{p}{\rm
Tr}_{S}[U^{\dagger}_{1}U^{\dagger}_{2}\hat{M}U_{2}U_{1}\rho_{S}^{0}-U^{\dagger}_{1}\hat{M}U_{1}\rho_{S}^{0}].$
(17)

The $2\bar{p}$ can be taken equal to $1$ with opportune rescaling of
$\hat{p}$.

Notice that in Eq. (17) we have the difference of the average value of
$\hat{M}$ evolved with $U_{2}U_{1}$ and the one evolved with $U_{1}$. By
taking $U_{1}=U(\vec{\theta}-s{\bf e}_{j_{1}})$ and
$U_{2}U_{1}=U(\vec{\theta}+s{\bf e}_{j_{1}})$, using $f(\vec{\theta})$
definition in Eq. (3), we recover the numerator in Eq. (4). In the
$\vec{\theta}$ parameter space, $U_{1}$ is the transformation or the circuits
that evolves the states to $\vec{\theta}-s{\bf e}_{j_{1}}$ and $U_{2}U_{1}$
generates the evolution from $\vec{\theta}-s{\bf e}_{j_{1}}$ to
$\vec{\theta}+s{\bf e}_{j_{1}}$ (see Appendix A).

###  III.1 Practical Implementation

A practical implementation of the protocol is shown in Fig. 1. We need $n$
logical qubits plus one detector/ancilla qubit. The latter is initialised to
$\mathinner{|{\eta_{D}^{0}}\rangle}=(\mathinner{|{0}\rangle}+\mathinner{|{1}\rangle})/\sqrt{2}$.
This means that in the following $\hat{p}=\sigma_{z}$ (where $\sigma_{z}$ is
the usual Pauli operator). However, in the following, to distinguish the
detector qubit from the logical ones, we keep using the $\hat{p}$ notation.

To calculate the gradient close to a value $\theta$ we need to choose
$\lambda\ll 1$ (see Appendix B).

The protocol to estimate the gradient of $f$ in the point $\vec{\theta}$ is
the following

  1. 1.

Fix $\vec{\theta}$ and $\lambda\ll 1$.

    1. (a)

apply a Hadamard operator to the detector qubit

    2. (b)

run the circuit $U({\vec{\theta}}-s{\bf e}_{j_{1}})$

    3. (c)

couple the system and the detector with the transformation
$e^{-i\lambda\hat{p}\otimes\hat{M}}$

    4. (d)

run the circuits from ${\vec{\theta}}-s{\bf e}_{j_{1}}$ to
${\vec{\theta}}+s{\bf e}_{j_{1}}$: $U({\vec{\theta}}+s{\bf
e}_{j_{1}},{\vec{\theta}}-s{\bf e}_{j_{1}})$

    5. (e)

couple the system and the detector with the transformation
$e^{i\lambda\hat{p}\otimes\hat{M}}$

    6. (f)

apply a Hadamard operator to the detector qubit

    7. (g)

measure the detector qubit

  2. 2.

repeat $m$ times to determine the detector population.

The detector phase can be measured with interferometric techniques. It can be
shown that the measured population of the detector is related to the
accumulated phase $\phi(\lambda)$ by the relations (Appendix B)

$\displaystyle P_{0}$ $\displaystyle=$
$\displaystyle\frac{1}{2}(1+\cos^{2}\phi)=\cos^{2}\frac{\phi}{2}$
$\displaystyle P_{1}$ $\displaystyle=$
$\displaystyle\frac{1}{2}(1+\cos^{2}\phi)=\sin^{2}\frac{\phi}{2}$

The error on the population scales with the number of repetitions $m$ as
$1/\sqrt{m}$. For small error in the phase $\delta\phi$ we obtain

$|\delta\phi|=\frac{2}{\sqrt{m}\sin\left(2\phi(\lambda)\right)}.$ (18)

The gradient of $f$ estimated with the QNDM approach can be used to determine
the updated vector parameter according to the optimizer algorithm.

##  IV Example and resources estimate

To better understand the QNDM approach we present an example that is general
enough to have many practical applications. For the sake of discussion,
suppose we want to use a quantum computer to calculate the minimum average of
a complex operator $\hat{M}$ in a high-dimensional quantum system. The
$\hat{M}$ operator can be written as the sum of weighted contributions McArdle
_et al._ (2020)

$\hat{M}=\sum_{j=1}^{J}h_{j}\hat{P}_{j}=\sum_{j=1}^{J}h_{j}\prod_{i=1}^{I}\sigma_{i}^{j}$
(19)

where $h_{j}$ is a real coefficient, $\sigma_{i}^{j}$ is one of the Pauli
operators where $i$ is the qubit index and $j$ denotes the term in the
Hamiltonian McArdle  _et al._ (2020) . Concurrently, $I$ and $J$ are two
natural numbers such that we can consider $H$ a complex operator and the $i$
($j$) index runs from $1$ to $I$ ($J$). The operators
$\hat{P}_{j}=\prod_{i}\sigma_{i}^{j}$ are usually called Pauli string
operators McArdle  _et al._ (2020).

These kinds of Hamiltonians describe a large number of different physical
systems from condensed matter to quantum chemistry Wecker  _et al._ (2015);
McArdle  _et al._ (2020). Before proceeding, two additional remarks must be
made. First, the length of a Pauli string is bounded by the number of qubits
$n$. Second, in physical systems because of the nature of the interaction
between particles, usually, we have a maximum of four terms composing the
Pauli strings, i.e., $I_{max}=4$ Wecker  _et al._ (2015); McArdle  _et al._
(2020).

Figure 2: The implementation of the exponential operator
$\exp\\{i\lambda\hat{p}\otimes\hat{M}\\}$ in a quantum computer with an
additional ancilla qubit $a$ (see Appendix C for details). The operator is
decomposed in a product of Pauli string operator
$\exp\\{i\lambda\hat{p}\otimes\hat{P}_{j}\\}$. The number of elementary
operations needed for every operator is shown below. The total number of the
elementary operators is $4nJ$ where $n$ is the number of logical qubits and
$J$ is the number of Pauli strings in $\hat{M}$.

If $\hat{M}$ is a complex operator, the implementation of the exponential
coupling $\exp\\{i\lambda\bar{p}\hat{M}\\}$ needed for the QNDM is, in
general, very resource-consuming. This might seem a critical drawback, but the
QNDM approach offers a clear and elegant way to bypass this problem. For the
purpose of measuring the average of the gradient in the QNDM approach, instead
of implementing the operator
$\exp{\\{i\lambda\hat{p}\otimes\sum_{j}h_{j}\hat{P}_{j}\\}}$, we can implement
$\prod_{j}\exp{\\{i\lambda h_{j}\hat{p}\otimes\hat{P}_{j}\\}}$ as shown in
Fig. 2.

To clarify this point, let us consider the operator
$\hat{M}=h_{1}\hat{P}_{1}+h_{2}\hat{P}_{2}+h_{3}\hat{P}_{3}$. If we implement
the operator $U_{full}=\exp{\\{i\lambda\hat{p}\otimes\hat{M}\\}}$ and then
take the derivative with respect to $\lambda$ to obtain
$\mathcal{G}_{\lambda}$, we obtain

$i\partial_{\lambda}U_{full}\Big{|}_{\lambda=0}=\bar{p}(h_{1}\hat{P}_{1}+h_{2}\hat{P}_{2}+h_{3}\hat{P}_{3}).$
(20)

If we implement the operator $U_{simpl}=\prod_{j}U_{j}=\exp{\\{i\lambda
h_{1}\hat{p}\otimes\hat{P}_{1}\\}}\exp{\\{i\lambda
h_{2}\hat{p}\otimes\hat{P}_{2}\\}}\exp{\\{i\lambda
h_{3}\hat{p}\otimes\hat{P}_{3}\\}}$ (with $U_{j}=\exp{\\{i\lambda
h_{j}\hat{p}\otimes\hat{P}_{j}\\}}$), the derivative reads

$\displaystyle i\partial_{\lambda}U_{simpl}\Big{|}_{\lambda=0}$
$\displaystyle=$
$\displaystyle\bar{p}\Big{[}(h_{1}\hat{P}_{1})U_{1}U_{2}U_{3}+$ (21)
$\displaystyle+$ $\displaystyle
U_{1}(h_{2}\hat{P}_{2})U_{2}U_{3}+U_{1}U_{2}U_{3}(h_{3}\hat{P}_{3})\Big{]}\Big{|}_{\lambda=0}$
$\displaystyle=$
$\displaystyle\bar{p}(h_{1}\hat{P}_{1}+h_{2}\hat{P}_{2}+h_{3}\hat{P}_{3}).$

Therefore, despite being different, the unitary operators $U_{full}$ and
$U_{simpl}$ contribute in the same way to the averages. The main advantage of
this approach lies in the fact that the single Pauli string operator
$\hat{P}_{j}$ can be more easily implemented in a quantum computer in terms of
single qubit operators.

In Fig. 2, we show how the product of exponential can be implemented. The
additional ancilla qubit is needed to properly implement the $\exp{\\{i\lambda
h_{j}\hat{p}\otimes\hat{P}_{j}\\}}$ operator. The technical details can be
found in Appendix C where we also discuss the resources needed in terms of
elementary logical quantum gates. It turns out that every
$\exp\\{i\lambda\hat{p}P_{j}\\}$ block needs a maximum of $4n$ logical
operation (Appendix C). If $\hat{M}$ is composed of $J$ Pauli strings we need
a maximum of $4nJ$ logical operations to implement the full
$\exp\\{i\lambda\hat{p}\hat{M}\\}$ operator.

At this point, we can estimate the cost of the circuit for the QNDM in Fig. 1.
Once we have implemented the operator $U_{1}=U(\vec{\theta}-s{\bf e}_{j_{1}})$
with $k$ logical operators, usually, we need a negligible number of logical
gates to arrive at the point $\vec{\theta}+s{\bf e}_{j_{1}}$. In fact, the
shift along the ${\bf e}_{j_{1}}$ direction involves only a single qubit
rotation as in Eq. (1). Therefore, the total cost of $U_{1}$ and $U_{2}$
operators is approximately $k$. The two system-detector couplings $\exp\\{\pm
i\lambda\bar{p}\hat{M}\\}$ are performed with $8nJ$ elementary operations. The
total cost of the QNDM circuit in Fig. 1 is $8nJ+k$.

Since the circuit must be repeated $m$ times to have the desired statistical
accuracy, the total cost to estimate the gradient of a quantum operator is of
$(k+8nJ)m$ elementary logical operators.

As an additional side remark, we want to emphasize that this result does not
rely on any approximation differently from the usual Lie-Trotter-Suzuki
decomposition McArdle  _et al._ (2020). The latter depends on the assumption
that $\lambda ph_{j}\ll 1$. This implies that it might be necessary to iterate
the process to implement the desired operator for any $\lambda ph_{j}$. On the
contrary, the present methods can be used for any value of $\lambda ph_{j}$,
although a small value of $\lambda$ might be helpful for phase evaluation.

##  V Second and Higher derivatives

The QNDM is even more convenient if we want to calculate the second derivative
in Eq. (II). In this case, we just need to couple the system and the detector
four times in a single run. The unitary transformation we implement is

$\mathcal{U}_{\lambda}=e^{i\lambda\bar{p}\hat{M}}U_{4}e^{-i\lambda\bar{p}\hat{M}}U_{3}e^{i\lambda\bar{p}\hat{M}}U_{2}e^{-i\lambda\bar{p}\hat{M}}U_{1}$
(22)

with the transformations (see Appendix A)

$\displaystyle U_{1}$ $\displaystyle=$ $\displaystyle
U\big{(}\vec{\theta}+s({\bf e}_{j_{1}}-{\bf e}_{j_{2}}),0\big{)}$
$\displaystyle U_{2}$ $\displaystyle=$ $\displaystyle
U\big{(}\vec{\theta}-s({\bf e}_{j_{1}}+{\bf e}_{j_{2}}),\vec{\theta}+s({\bf
e}_{j_{1}}-{\bf e}_{j_{2}})\big{)}$ $\displaystyle U_{3}$ $\displaystyle=$
$\displaystyle U\big{(}\vec{\theta}+s(-{\bf e}_{j_{1}}+{\bf
e}_{j_{2}}),\vec{\theta}-s({\bf e}_{j_{1}}+{\bf e}_{j_{2}})\big{)}$
$\displaystyle U_{4}$ $\displaystyle=$ $\displaystyle
U\big{(}\vec{\theta}+s({\bf e}_{j_{1}}+{\bf e}_{j_{2}}),\vec{\theta}+s(-{\bf
e}_{j_{1}}+{\bf e}_{j_{2}})\big{)}$ (23)

where we have denoted with $U(\vec{\theta}_{i},\vec{\theta}_{j})$ the operator
that allows us to pass from $\vec{\theta}_{j}$ to $\vec{\theta}_{i}$ in the
$\vec{\theta}$ parameter space (Appendix A). In Figure 3 the process is shown
in a bidimensional space, i.e., $\vec{\theta}=(\theta_{1},\theta_{2})$.
Formally the $\mathcal{G}_{\lambda}$ and its derivative are written as in Eqs.
(14) and (15). By direct calculation, we obtain (see Appendix A)

$\displaystyle-i\partial_{\lambda}\mathcal{G}_{\lambda}\Big{|}_{\lambda=0}$
$\displaystyle=f(\vec{\theta}+s({\bf e}_{j_{1}}+{\bf
e}_{j_{2}}))-f(\vec{\theta}+s(-{\bf e}_{j_{1}}+{\bf e}_{j_{2}}))$ (24)
$\displaystyle-f(\vec{\theta}+s({\bf e}_{j_{1}}-{\bf
e}_{j_{2}}))+f(\vec{\theta}-s({\bf e}_{j_{1}}+{\bf e}_{j_{2}})).$

By dividing this by $2\sin^{2}s$ we obtain the expression in Eq. (II).

Figure 3:  Path in the bidimensional $(\theta_{1},\theta_{2})$ parameter space
of the unitary operator $U(\vec{\theta})$ to obtain the second derivative.
First, we implement $U\big{(}\vec{\theta}-s({\bf e}_{j_{1}}+{\bf
e}_{j_{2}})\big{)}$ (point $(\theta_{1}-s,\theta_{2}-s)$), the implement the
shift of $2s$ along the $\theta_{1}$ coordinate to arrive to point $2$, i.e.,
$(\theta_{1}+s,\theta_{2}+s)$. Then, we move to point $3$, i.e.,
$(\theta_{1}-s,\theta_{2}+s)$ and, finally, to point $4$,
$(\theta_{1}+s,\theta_{2}-s)$. Notice that the path is not unique. In fact,
one can equivalently go around in the other direction (i.e. points 1, 4, 3 and
2) or start from point 3 (both directions).

Interestingly, this can be evaluated with

  * •

$n$ qubits

  * •

a single circuit of approximate $k$ the logical gates (needed for
$U(\vec{\theta}+s({\bf e}_{j_{1}}+{\bf e}_{j_{2}}))$)

  * •

As above, we are measuring the detector phase and we need only $m$ repetitions
to have an error of $1/\sqrt{m}$.

Notice that in the DM approach the measure of the second derivative needs the
estimate of the average of four observables, i.e., $M(\vec{\theta}+s(\pm{\bf
e}_{j_{1}}\pm{\bf e}_{j_{2}}))$. On the contrary, with the QNDM we obtain the
same information with a single-phase estimation. The same scheme can be used
to extract information about the derivative of arbitrary order. Interestingly,
the increase in the derivative order implies an increase in the resources
(logical gate and measurements) for the DM while only weakly affecting the
complexity of the QNDM. Following the gate counting procedure discussed in
Sec. IV, we estimate that to estimate the gradient of $f$ with the QNDM, we
need a maximum of $(k+16nJ)m$ elementary logical operators.

##  VI Comparison with other approaches

We are now in a position to discuss the advantages of the QNDM with respect to
the usual DM protocol Kandala  _et al._ (2017); Havlíček  _et al._ (2019);
Mari  _et al._ (2021).

We must stress that the cost of the resources (in terms of circuit
repetitions, number of unitary operations used and so on) depends critically
on the specific problem we want to solve (see, for example, Ref. Wecker  _et
al._ (2015)). A second factor to be kept into account in the comparison is the
quantum hardware used. Usually, the approaches that use the quantum phase to
store the information (this article and O’Brien  _et al._ (2019); McArdle  _et
al._ (2020); Abrams and Lloyd (1999); Banchi and Crooks (2021)) need efficient
hardware, error mitigation and correction techniques since they must preserve
quantum coherences for long times. Therefore, for noisy hardware, the quantum
phase approaches can be less efficient than the direct measure ones McArdle
_et al._ (2020); Mari  _et al._ (2021).

Given these premises, we are interested in giving a rough estimate of the
advantages and pointing out where these can be found without entering the
details of the simulated physical system. In addition, we suppose that the
noise can be eliminated so that the dephasing effect can be neglected.

###  VI.1 Gradient and second derivative measurement

We suppose that we need $\approx k$ logical operators to implement both
$U(\vec{\theta}\pm s{\bf e}_{j_{1}})$. The more straightforward way to
evaluate the function $f(\vec{\theta})$ is to run the circuit/dynamics and
then measure the operator $\hat{M}$ Kandala  _et al._ (2017); Havlíček  _et
al._ (2019); Mari  _et al._ (2021).

As discussed above, since the complexity of the system usually increases with
its dimension, for realistic and practical situations this can be an extremely
difficult task McArdle  _et al._ (2020); Verteletskyi  _et al._ (2020); Yen
_et al._ (2020); Yen and Izmaylov (2021). Therefore, we need to separately
measure every single Pauli string $\hat{P}_{J}$.

The procedure is the following. First, we run the circuit with the unitary
operator $U(\vec{\theta}-s{\bf e}_{j_{1}})$ and measure $\hat{P}_{j}$.
Repeating this protocol $m$ times we can obtain the average value
$\langle\hat{P}_{j}(\vec{\theta}-s{\bf e}_{j_{1}})\rangle$ with precision
scaling as $1/\sqrt{m}$. As discussed in the appendix C, since the measure in
a quantum computer are always performed in the standard basis, i.e.,
eigenstates of $\sigma_{z}$, a generic Pauli string we have to rotate a
maximum of $n$ qubit with Hadamard gates. For example, to measure the Pauli
string $\sigma_{x}^{1}\sigma_{x}^{2}\sigma_{x}^{3}$ we need to apply the
Hadamard operators $H_{had}^{1}H_{had}^{2}H_{had}^{3}$.

We repeat the procedure for all the $\hat{P}_{j}$ for a total of $J$
repetitions. By summing the weighted averages we obtain
$\langle\hat{M}(\vec{\theta}-s{\bf
e}_{j_{1}})\rangle=\sum_{j}h_{j}\langle\hat{P}_{j}(\vec{\theta}-s{\bf
e}_{j_{1}})\rangle=f(\vec{\theta}-s{\bf e}_{j_{1}})$. Then, we repeat the same
procedure $mJ$ times running the circuit with the unitary operator
$U(\vec{\theta}+s{\bf e}_{j_{1}})$ and measure $\hat{M}$ to obtain $\langle
M(\vec{\theta}+s{\bf e}_{j_{1}})\rangle=f(\vec{\theta}+s{\bf e}_{j_{1}})$.

In terms of resources, summing all the contributions we must apply $2mJ(k+n)$.

For the second derivative, the counting is analogous but we have to consider
that we need four averages of the $\hat{M}$ operators and, this, a total of
$4mJ(k+n)$.

| qubits | iterations | gates
---|---|---|---
$\partial_{1}f$–DM | $n$ | $2m$ | $\approx 2mJ(k+2n)$
$\partial_{2}f$–DM | $n$ | $4m$ | $\approx 4mJ(k+2n)$
$\partial_{1}f$–QNDM | $n+2$ | $m$ | $\approx m(k+8nJ)$
$\partial_{2}f$–QNDM | $n+2$ | $m$ | $\approx m(k+16nJ)$
Table 1: Comparison between different methods to measure the gradient and
second derivative of $f(\vec{\theta})$: Direct Measurement (DM) and Quantum
Non-Demolition Measurement (QNDM). The parameter $m$ determines the error on
the estimate of the gradient. The parameter $k$ determines the number of gate
operations needed to implement the operators $U(\vec{\theta})$. Since the
shift $s$ occurs along a single direction ${\bf e}_{j_{1}}$, we suppose that
the number of gates to implement $U(\vec{\theta}\pm s{\bf e}_{j_{1}})$. For
small $s$ is approximatively the same and equal to $k$.

This resource counting must be compared to the one estimate for the QNDM in
Sec. III and V. Table 1 summarizes the resources needed for the two
approaches. As can be seen, the QNDM has a clear advantage in the number of
iterations of the circuit. This is because in a single step we can store
information about two average values of $\hat{M}$, i.e., about the gradient.

We have to stress that in the QNDM approach we usually run a more complex and
deeper quantum circuit. For a comparison of the logical quantum gates, we can
refer to Table 1. We can consider two limits. First, we notice that the number
of logical gates $k$ needed to apply the unitary operator $U(\vec{\theta})$ is
usually greater than $n$ Wecker  _et al._ (2015). In the limit of $k\gg nJ$,
we have that the resource ratio scales as ${\rm DM}/{\rm QNDM}\approx J$.
Therefore, depending on the number of Pauli strings to be measured, we can
have a substantial advantage with the QNDM. This regime is likely to be
interesting for the simulation in quantum chemistry. For example, for the
simulations of medium complex molecules, we could estimate Wecker  _et al._
(2015): $k\approx 10^{9}-10^{10}$ operations to implement $U(\vec{\theta})$,
$n\approx 10^{2}-10^{3}$ (that is the dimension of quantum computer available
in the next years) and $J>10^{3}$ number of Pauli string in the Hamiltonian.
With these numbers, $k\gg nJ$ and we have a reduction of roughly $J$ logical
operations.

In the opposite limit $k\ll nJ$, we have ${\rm DM}/{\rm QNDM}\approx k/n$.
Since the number of operations to implement $U(\vec{\theta})$ is greater than
the number of qubit $n$, also, in this case, the QNDM has an advantage with
respect to the DM.

This advantage increases with the order of the derivative we need. While the
number of repetitions of the DM circuit increases, for the QNDM is constant.
In most common approaches and optimizers, only the first and second
derivatives of the cost function are usually exploited. However, there are
some indications that derivatives of a high order can be used to treat
classically hard optimization problems Ahookhosh and Nesterov (2021a, b). In
addition, in some specific problems in finance, the derivative themselves
carry the relevant information Stamatopoulos  _et al._ (2022). In these cases,
the possibility to efficiently measure the high-order derivatives of the cost
function with the QNDM can give substantial advantages.

##  VII Conclusions

We have discussed an alternative approach to estimate the gradient and higher-
order derivatives of a cost function on quantum hardware by using a quantum
detector. The relevant information is stored in the phase of the detector.
Since the evolution of the system is not affected by the detector, this is
usually referred to as Quantum Non-Demolition Measurement Clerk (2011). At its
roots, the idea is related to the weak values and weak measurements in quantum
mechanics Aharonov  _et al._ (1988); Bednorz  _et al._ (2012).

The main advantage of the QNDM is that a single quantum circuit allows us to
extract the information about the gradient, i.e., the first derivative, and
higher derivative of the average values of a quantum observable. This reduces
the number of circuit iterations needed for the application of the most common
optimization algorithms such as gradient descent and Newton optimization.

In terms of logical operation needed, the QNDM approach implies the use of
more complex quantum circuits since we need to store information in the
detector phase. A quantitative comparison with the Direct Measurement
approaches depends critically on the problem we want to solve. However, by
considering a fairly general operator structure [see Eq. (19)], we were able
to extract some general trends. The QNDM approach is usually advantageous in
terms of logical operator needed. This advantage increases in the cases in
which the quantum operator is written as the sum of many Pauli strings. As
above, for higher derivatives, there is a further advantage in the reduction
of logical operators.

The fast and accurate estimation of the gradient and high derivatives plays a
key role in the variational quantum circuits McArdle  _et al._ (2020); Cerezo
_et al._ (2021). These can be applied to numerous problems in physics and
chemistry as well as optimization problems. The gradient estimation performed
by a quantum computer is fed to a classical computer that processes the
information to reach the minimum of the cost function. This quantum-classical
hybrid scheme has the potential to be one of the first practical cases in
which we have a quantum advantage over classical computers.

In this perspective, any approach that decreases the quantum resources brings
these applications closer to the actual implementation on a quantum computer.
Further studies are needed to fully gauge the impact and advantages of the
QNDM approach. However, especially for the NISQ computers where the resources
are limited, the QNDM approach could offer a valuable alternative to reduce
the resources and efficiently implement some variational quantum circuits.

###### Acknowledgements.

The authors acknowledge financial support from INFN.

## References

  * Peruzzo  _et al._ (2014) A. Peruzzo, J. McClean, P. Shadbolt, M.-H. Yung, X.-Q. Zhou, P. J. Love, A. Aspuru-Guzik, and J. L. O’Brien, Nature Communications 5, 4213 (2014).
  * Kandala  _et al._ (2017) A. Kandala, A. Mezzacapo, K. Temme, M. Takita, M. Brink, J. M. Chow, and J. M. Gambetta, Nature 549, 242 (2017).
  * Farhi  _et al._ (2014) E. Farhi, J. Goldstone, and S. Gutmann, arXiv preprint arXiv:1411.4028 (2014).
  * McArdle  _et al._ (2020) S. McArdle, S. Endo, A. Aspuru-Guzik, S. C. Benjamin, and X. Yuan, Rev. Mod. Phys. 92, 015003 (2020).
  * Cerezo  _et al._ (2021) M. Cerezo, A. Arrasmith, R. Babbush, S. C. Benjamin, S. Endo, K. Fujii, J. R. McClean, K. Mitarai, X. Yuan, L. Cincio, and P. J. Coles, Nature Reviews Physics 3, 625 (2021).
  * Mari  _et al._ (2021) A. Mari, T. R. Bromley, and N. Killoran, Phys. Rev. A 103, 012405 (2021).
  * Clerk (2011) A. A. Clerk, Phys. Rev. A 84, 043824 (2011).
  * Bednorz  _et al._ (2012) A. Bednorz, W. Belzig, and A. Nitzan, New J. Phys. 14, 013009 (2012).
  * Aharonov  _et al._ (1988) Y. Aharonov, D. Z. Albert, and L. Vaidman, Phys. Rev. Lett. 60, 1351 (1988).
  * Solinas and Gasparinetti (2015) P. Solinas and S. Gasparinetti, Phys. Rev. E 92, 042150 (2015).
  * Solinas and Gasparinetti (2016) P. Solinas and S. Gasparinetti, Phys. Rev. A 94, 052103 (2016).
  * Solinas  _et al._ (2021) P. Solinas, M. Amico, and N. Zanghì, Phys. Rev. A 103, L060202 (2021).
  * Solinas  _et al._ (2022) P. Solinas, M. Amico, and N. Zanghì, Phys. Rev. A 105, 032606 (2022).
  * Talkner  _et al._ (2007) P. Talkner, E. Lutz, and P. Hänggi, Phys. Rev. E 75, 050102 (2007).
  * O’Brien  _et al._ (2019) T. E. O’Brien, B. Senjean, R. Sagastizabal, X. Bonet-Monroig, A. Dutkiewicz, F. Buda, L. DiCarlo, and L. Visscher, npj Quantum Information 5, 113 (2019).
  * Banchi and Crooks (2021) L. Banchi and G. E. Crooks, Quantum 5, 386 (2021).
  * Cheng (2010) R. Cheng, arXiv preprint arXiv:1012.1337 (2010).
  * Liu  _et al._ (2019) J. Liu, H. Yuan, X.-M. Lu, and X. Wang, Journal of Physics A: Mathematical and Theoretical 53, 023001 (2019).
  * Solinas  _et al._ (2017) P. Solinas, H. J. D. Miller, and J. Anders, Phys. Rev. A 96, 052115 (2017).
  * Wecker  _et al._ (2015) D. Wecker, M. B. Hastings, and M. Troyer, Phys. Rev. A 92, 042303 (2015).
  * Havlíček  _et al._ (2019) V. Havlíček, A. D. Córcoles, K. Temme, A. W. Harrow, A. Kandala, J. M. Chow, and J. M. Gambetta, Nature 567, 209 (2019).
  * Abrams and Lloyd (1999) D. S. Abrams and S. Lloyd, Phys. Rev. Lett. 83, 5162 (1999).
  * Verteletskyi  _et al._ (2020) V. Verteletskyi, T.-C. Yen, and A. F. Izmaylov, The Journal of Chemical Physics 152, 124114 (2020), https://doi.org/10.1063/1.5141458 .
  * Yen  _et al._ (2020) T.-C. Yen, V. Verteletskyi, and A. F. Izmaylov, Journal of Chemical Theory and Computation 16, 2400 (2020), pMID: 32150412, https://doi.org/10.1021/acs.jctc.0c00008 .
  * Yen and Izmaylov (2021) T.-C. Yen and A. F. Izmaylov, PRX Quantum 2, 040320 (2021).
  * Ahookhosh and Nesterov (2021a) M. Ahookhosh and Y. Nesterov, arXiv preprint arXiv:2107.05958 (2021a).
  * Ahookhosh and Nesterov (2021b) M. Ahookhosh and Y. Nesterov, arXiv preprint arXiv:2109.12303 (2021b).
  * Stamatopoulos  _et al._ (2022) N. Stamatopoulos, G. Mazzola, S. Woerner, and W. J. Zeng, Quantum 6, 770 (2022).
  * McClean  _et al._ (2018) J. R. McClean, S. Boixo, V. N. Smelyanskiy, R. Babbush, and H. Neven, Nature Communications 9, 4812 (2018).

Appendix

##  Appendix A Unitary evolution in the parameter space

In the main text, we have written the total unitary transformation to measure
the gradient as
$U_{tot}=e^{i\lambda\hat{p}\otimes\hat{M}}U_{2}e^{-i\lambda\hat{p}\otimes\hat{M}}U_{1}$.
The $U_{1}$ and $U_{2}$ operation can be thought of as a sequence of non-
commuting logical quantum gates or, alternatively using a more physical
language, generated by a time-dependent Hamiltonian $H(\theta(t))$.

We consider the change in ${\vec{\theta}}$ along a single direction ${\bf
e}_{j_{1}}$ so that the difference between the points ${\vec{\theta}}\pm s{\bf
e}_{j_{1}}$ is parametrized only by $s$. To build the operator
$U({\vec{\theta}})$, we divide the transformation into $N$ steps as usually
done for the evolution of a quantum system or the implementation in a quantum
circuit. Similarly, we suppose that to implement the operators
$U({\vec{\theta}}-s{\bf e}_{j_{1}})$ and $U({\vec{\theta}}+s{\bf e}_{j_{1}})$
we need $N-N_{1}$ and $N+N_{2}$ operations, respectively. Denoting the single
unitary operation a $u_{k}$, we can write

$U_{1}=u_{N-N_{1}}...u_{1}=U({\vec{\theta}}-s{\bf e}_{j_{1}})$ (S1)

and $U_{2}=\Pi_{j=1}^{N+N_{2}}u_{j}=\Pi_{j=N-N_{1}}^{N+N_{2}}u_{j}U_{1}$.
Therefore, we can obtain $U_{2}$ by implementing first $U_{1}$ and then
applying other $N_{2}+N_{1}$ operators.

Using the mapping $U_{1}=U({\vec{\theta}}-s{\bf e}_{j_{1}})$ and
$U_{2}=U({\vec{\theta}}+s{\bf e}_{j_{1}})$, we rewrite the above equation as

$U({\vec{\theta}}+s{\bf e}_{j_{1}})=U({\vec{\theta}}+s{\bf
e}_{j_{1}},{\vec{\theta}}-s{\bf e}_{j_{1}})U({\vec{\theta}}-s{\bf e}_{j_{1}})$
(S2)

where we have implicitly denoted with $U({\vec{\theta}}+s{\bf
e}_{j_{1}},{\vec{\theta}}-s{\bf e}_{j_{1}})$ the operator that allows us to
pass from ${\vec{\theta}}-s{\bf e}_{j_{1}}$ to ${\vec{\theta}}+s{\bf
e}_{j_{1}}$ with $N_{2}+N_{1}$ operators.

In a similar way, we can implement the series of operators
$\mathcal{U}_{\lambda}=e^{i\lambda\bar{p}\hat{M}}U_{4}e^{-i\lambda\bar{p}\hat{M}}U_{3}e^{i\lambda\bar{p}\hat{M}}U_{2}e^{-i\lambda\bar{p}\hat{M}}U_{1}$
to measure the second derivative. The $U_{j}$ operators (with $j=1,..,4$)
correspond to the operators working in the $\vec{\theta}$ space.

We define the quasi-characteristic function as in Eq. (9) and use the above
expression for $\mathcal{U}_{\lambda}$. Taking the first derivative with
respect to $\lambda$, and using the mapping

$\displaystyle U(\vec{\theta}+s({\bf e}_{j_{1}}-{\bf e}_{j_{2}}))$
$\displaystyle=$ $\displaystyle U_{1}$ $\displaystyle U(\vec{\theta}-s({\bf
e}_{j_{1}}+{\bf e}_{j_{2}}))$ $\displaystyle=$ $\displaystyle U_{2}U_{1}$
$\displaystyle U(\vec{\theta}+s(-{\bf e}_{j_{1}}+{\bf e}_{j_{2}}))$
$\displaystyle=$ $\displaystyle U_{3}U_{2}U_{1}$ $\displaystyle
U(\vec{\theta}+s({\bf e}_{j_{1}}+{\bf e}_{j_{2}}))$ $\displaystyle=$
$\displaystyle U_{4}U_{3}U_{2}U_{1},$ (S3)

we obtain

$\displaystyle-i\partial_{\lambda}\mathcal{G}_{\lambda}\Big{|}_{\lambda=0}$
$\displaystyle=$ $\displaystyle{\rm Tr}_{S}[U^{\dagger}(\vec{\theta}+s({\bf
e}_{j_{1}}+{\bf e}_{j_{2}}))MU(\vec{\theta}+s({\bf e}_{j_{1}}+{\bf
e}_{j_{2}}))\rho_{S}^{0}-U^{\dagger}(\vec{\theta}+s(-{\bf e}_{j_{1}}+{\bf
e}_{j_{2}}))MU(\vec{\theta}+s(-{\bf e}_{j_{1}}+{\bf e}_{j_{2}}))\rho_{S}^{0}$
(S4) $\displaystyle-$ $\displaystyle U^{\dagger}(\vec{\theta}+s({\bf
e}_{j_{1}}-{\bf e}_{j_{2}}))MU(\vec{\theta}+s({\bf e}_{j_{1}}-{\bf
e}_{j_{2}}))\rho_{S}^{0}+U^{\dagger}(\vec{\theta}-s({\bf e}_{j_{1}}+{\bf
e}_{j_{2}}))MU(\vec{\theta}-s({\bf e}_{j_{1}}+{\bf e}_{j_{2}}))\rho_{S}^{0}]$
$\displaystyle=$ $\displaystyle f(\vec{\theta}+s({\bf e}_{j_{1}}+{\bf
e}_{j_{2}}))-f(\vec{\theta}+s(-{\bf e}_{j_{1}}+{\bf
e}_{j_{2}}))-f(\vec{\theta}+s({\bf e}_{j_{1}}-{\bf
e}_{j_{2}}))+f(\vec{\theta}-s({\bf e}_{j_{1}}+{\bf e}_{j_{2}}))$

where in the last line we have used the definition of $f(\vec{\theta})$ in Eq.
(3). A part from a $2\sin s$ factor, this is the second derivative as in Eq.
(II). By denoting with $U(\vec{\theta}_{2},\vec{\theta}_{1})$ the operator
that allows us to pass from $\vec{\theta}_{1}$ to $\vec{\theta}_{2}$ , we
obtain the explicit form of the $U_{j}$ operators as in Eq. (23) of the main
text.

##  Appendix B Detection with a qubit

We consider the case in which a qubit, i.e., a two-level system, is used as a
detector. The initial state of the detector is fixed at
$\mathinner{|{\eta_{0}}\rangle}=\cos\alpha/2\mathinner{|{0}\rangle}+\sin\alpha/2\mathinner{|{1}\rangle}$.
The initial contribution to $\mathcal{G}_{\lambda}$ is
$\mathinner{\langle{0|\eta_{0}}\rangle}\mathinner{\langle{\eta_{0}|1}\rangle}=1/2~{}\sin\alpha$.

After the evolution the detector acquires a phase $\phi(\lambda)$ so that its
final state is
$\mathinner{|{\eta_{f}}\rangle}=\cos\alpha/2\mathinner{|{0}\rangle}+\sin\alpha/2e^{-i\phi(\lambda)}\mathinner{|{1}\rangle}$.
The numerator of the $\mathcal{G}_{\lambda}$ is
$\mathinner{\langle{0|\eta_{f}}\rangle}\mathinner{\langle{\eta_{f}|1}\rangle}=1/2~{}\sin\alpha
e^{i\phi(\lambda)}$ and

$\mathcal{G}_{\lambda}=e^{i\phi(\lambda)}.$ (S5)

is exactly the phase accumulated during the evolution.

We can verify that $\mathcal{G}_{\lambda}$ has the properties discussed in the
previous section recalling that for $\lambda=0$ there is no system-detector
interaction and, thus, no phase is accumulated, i.e., $\phi(0)=0$.

In particular,
$-i\partial_{\lambda}\mathcal{G}_{\lambda}|_{\lambda=0}=\partial_{\lambda}\phi(\lambda)|_{\lambda=0}$;
therefore, the average value and the gradient are related to the derivative of
the accumulated phase evaluated in the origin.

In the experiments, we need to calculate
$\partial_{\lambda}\phi(\lambda)|_{\lambda=0}$ numerically. We can approximate
it as

$-i\partial_{\lambda}\mathcal{G}_{\lambda}|_{\lambda=0}=\partial_{\lambda}\phi(\lambda)|_{\lambda=0}\approx\frac{\phi(\Delta\lambda)-\phi(0)}{\Delta\lambda}=\frac{\phi(\Delta\lambda)}{\Delta\lambda}$
(S6)

We can calculate the desired quantity by evaluating the protocol and the phase
at a small $\Delta\lambda\ll 1$.

Let us consider the final detector state
$\mathinner{|{\eta_{f}}\rangle}=\cos\alpha/2\mathinner{|{0}\rangle}+\sin\alpha/2e^{-i\phi(\lambda)}\mathinner{|{1}\rangle}$.
By applying a Hadamard gate and measuring we obtain the probabilities to
measure $0$ and $1$

$\displaystyle P_{0}$ $\displaystyle=$
$\displaystyle\frac{1}{2}(1+\sin\alpha\cos^{2}\phi)$ $\displaystyle P_{1}$
$\displaystyle=$ $\displaystyle\frac{1}{2}(1+\cos\alpha\cos^{2}\phi),$

respectively. For the initial state with $\alpha=\pi/4$, we recover the
results of the main text. The statistical error on, say, $P_{0}$ scales as
$\delta P_{0}=\beta/\sqrt{m}$ where $m$ is the number of repetitions and
$\beta$ a proportionality constant. At the same time, if we consider how
$P_{0}$ changes under a small $\delta\phi$ variation, we have that

$\delta P_{0}=-\cos\phi\sin\phi\delta\phi$ (S7)

and, then,

$|\delta\phi|=\frac{2\beta}{\sqrt{m}\sin 2\phi}.$ (S8)

##  Appendix C Measuring the Pauli string - QNDM and DM circuit implementation

Figure S1: The circuit to implement the operator
$\exp\\{i\lambda\sigma_{z}^{D}\prod_{i}\sigma_{i}^{j}\\}$. We need a maximum
of $2(n+1)$ Hadamard-like gates (to diagonalize the $\sigma_{x}^{k}$ or
$\sigma_{y}^{k}$), $2(n+1)$ CNOT gates and a phase gate (denoted as $ph$) for
a total of $\approx 4n$ elementary gates. Notice that if the $\sigma_{z}^{k}$
operator appears in $\hat{P}_{j}$, we do not need to apply any Hadamard
operator and, thus, the $2(n+1)$ is the estimated maximum number of Hadamard-
like operators needed.

Let us consider the implementation of the exponential operator
$\exp\\{i\lambda\hat{p}\otimes\hat{P}_{j}\\}$ where the Pauli string
$\hat{P}_{j}$ is a product of Pauli operators:
$\hat{P}_{j}=\prod_{i}\sigma_{i}^{j}$.

Having in mind the implementation in a quantum computer, we can substitute the
detector $\hat{p}$ operator with a Pauli operator $\sigma_{z}^{D}$. The
exponential operator reads
$\exp\\{i\lambda\sigma_{z}^{D}\prod_{i}\sigma_{i}^{j}\\}$; therefore, the
presence of the detector simply increases the dimension of the Pauli string of
one qubit.

To implement the exponential operator, we first rotate the Pauli operator to
have a product of $\sigma_{z}$ operators. This can be done by applying
Hadamard operators $H_{had,x}$ for the $\sigma_{x}$ and the corresponding
$H_{had,y}$ for the $\sigma_{y}$. That is $\left(\prod
H_{had,l}\right)\exp\\{i\lambda\sigma_{z}^{D}\hat{P}_{j}\\}\left(\prod
H_{had,l}\right)=\exp\\{i\lambda\sigma_{z}^{D}\prod_{p}\sigma_{z}^{p}\\}$.

The diagonalized exponential operator generates a phase factor $\exp\\{\pm
i\lambda\\}$ depending on the logical string to whom it is applied. More
specifically, taking
$\sigma_{z}\mathinner{|{0}\rangle}=-\mathinner{|{0}\rangle}$ and
$\sigma_{z}\mathinner{|{1}\rangle}=\mathinner{|{1}\rangle}$, the
transformation is
$\mathinner{|{x_{1}x_{2}...x_{n}}\rangle}\rightarrow\exp\\{i\lambda\\}\mathinner{|{x_{1}x_{2}...x_{n}}\rangle}$
if the number of $x_{i}=0$ is even and
$\mathinner{|{x_{1}x_{2}...x_{n}}\rangle}\rightarrow\exp\\{-i\lambda\\}\mathinner{|{x_{1}x_{2}...x_{n}}\rangle}$
if the number of $x_{i}=0$ is odd.

To count the parity of $x_{i}=0$ in a string, we add an ancilla qubit in
position $n+2$ (the $n+1$ qubit is the detector) and apply a sequence of $n$
${\rm C}_{i}{\rm NOT}_{n+2}$ logical gate with $i=1,2,...,n+1$. The ${\rm
C}_{i}{\rm NOT}_{n+2}$ is controlled on the $i$-th logical qubit and acts
always in the ancilla qubit. As a consequence, if we have an even number of
$x_{i}=0$, the ancilla qubit will end up in the $\mathinner{|{0}\rangle}$
state, otherwise, i.e., for an odd number of $x_{i}=0$, in the
$\mathinner{|{1}\rangle}$ state. Then we apply a phase gate on the ancilla
qubit $\exp\\{i\lambda\sigma_{z}^{a}\\}$ and reverse the ${\rm C}_{i}{\rm
NOT}_{n+2}$ operations to reset the ancilla qubit to its original state
$\mathinner{|{0}\rangle}$. The full circuit is shown in Fig. S1.

To implement this operator we need $2(n+1)$ Hadamard gates, $2(n+1)$ CNOT
gates and a phase gate for a total of $4(n+1)+1\approx 4n$ elementary gates
for $n\gg 1$. Notice that this is the worst case where we have to apply $2n$
Hadamard gates. In many cases, the Pauli string operators are the product of a
limited number of Pauli operators so the number of Hadamard gates can be
significantly reduced. For example, in physical problems where we need to find
the energy of the ground state of a Hamiltonian, the Pauli strings are usually
the product of only four Pauli operators McClean  _et al._ (2018); Wecker  _et
al._ (2015).

For the DM approach, we have an analogous situation. In quantum computers, the
measurement is performed in a fixed basis (usually the one in which
$\sigma_{z}$ is diagonal). Since we cannot directly measure an alternative
Pauli observable ($\sigma_{x}$ or $\sigma_{y}$), we have to rotate the qubit
state and perform the measurement in the predetermined basis. As above, this
is done by applying an Hadarmard $H_{had,l}$ operator. Therefore, to measure
the average of a Pauli string, we need a maximum of $n$ $H_{had,l}$ logical
gates.