Efficient Gradient Estimation of Variational Quantum Circuits with Lie Algebraic Symmetries

Mohsen Heidari<EMAIL_ADDRESS>
Department of Computer Sciences, Indiana University, Bloomington, IN, USA
Masih Mozakka<EMAIL_ADDRESS>
Department of Computer Sciences, Indiana University, Bloomington, IN, USA
Wojciech Szpankowski<EMAIL_ADDRESS>
Department of Computer Sciences, Purdue University,    West Lafayette, IN, USA

Hybrid quantum-classical optimization and learning strategies are among the most promising approaches to harnessing quantum information or gaining a quantum advantage over classical methods. However, efficient estimation of the gradient of the objective function in such models remains a challenge due to several factors including the exponential dimensionality of the Hilbert spaces, and information loss of quantum measurements. In this work, we study generic parameterized circuits in the context of variational methods. We develop a framework for gradient estimation that exploits the algebraic symmetries of Hamiltonian characterized through  Lie algebra or group theory. Particularly, we prove that when the dimension of the dynamical Lie algebra is polynomial in the number of qubits, one can estimate the gradient with polynomial classical and quantum resources. This is done by a series of Hadamard tests applied to the output of the ansatz with no change to its circuit. We show that this approach can be equipped with classical shadow tomography to further reduce the measurement shot complexity to scale logarithmically with the number of parameters.

§ INTRODUCTION

VQA are among the hybrid quantum-classical strategies most prominent for quantum optimization and learning  [4, 5].
First introduced as VQE [30], VQA have been studied in a wide range of topics including optimization [10], quantum chemistry [20, 1, 13, 9], and quantum machine learning [12, 33, 28, 24, 15, 16].

VQAs are implemented via an ansatz which is a PQC. Several approaches suggest ways of creating a compact ansatz with shallow circuit depth and high accuracy. The ansatz can be layers of simple parametric gates acting on one or two qubits. More complex structures have been studied for HVA [40, 21].   The unitary variant of coupled cluster theory truncated at single and double excitations (UCCSD) is the first of such an ansatz [30].  Compactness and robustness against barren plateaus are some of the benefits of such models [19, 11, 23].

In general, the ansatz is represented by a generic parameterized operator as
where $\overrightarrow{a}\in \RR^p$ is the vector of adjustable model parameters, and $p$ the number of the parameters. Whichever way one implements the parameters and the ansatz, the objective is to minimize a predefined loss evaluated via a (Hermitian) observable ${O}$  for a specifically prepared input state as

\begin{equation}\label{eq:VQA cost function}
\Loss(\parameter) : = \tr{{O} ~ U(\overrightarrow{a}) \rho  U(\overrightarrow{a})^\dagger},
\end{equation}

where $\rho$ is the input (mixed) state density operator. To ensure computational tractability, it is assumed that the number of parameters $p=\poly(d),$ with $d$ being the number of qubits. With this notion, the objective is to find
\min_{\parameter} \Loss(\parameter).
In the context, of quantum machine learning, $\rho$ represents the expectation of the quantum samples, and $O$ is the incurred prediction loss. In the problem of finding the ground state of a Hamiltonian, $O$ represents the underlying Hamiltonian and the $\rho$ often is the initial state, such as the Hartee-Fock state.

Gradient Estimation. Compared to derivative-free methods, gradient-based optimizers have shown a significant advantage in terms of convergence guarantees [18, 36] and have been widely used in the literature [32, 10, 12, 33, 28].
At each iteration, an unbiased estimate of $\nabla \Loss(\parameter)$ is obtained, and the parameters are updated accordingly.  One of the main challenges is to estimate the gradient efficiently in terms of the classical and quantum resources defined through the following measures:

* The number of unique circuit configurations and gate complexity;

* The number of measurement shots (sample complexity);

* Classical computational complexity.

The number of unique circuit configurations is motivated by the existing NISQ implementations and the fact that quantum circuit reconfiguration tends to be a more costly task than generating measurement samples after executing each unique circuit [41]. The complexity of gradient estimation is negatively affected in quantum settings as the dimension of the associated Hilbert space is exponential in $d$ and the terms in the Hamiltonian $A(\parameter)$ are non-commuting. Moreover, other quantum properties that might contribute negatively are the stochasticity of quantum measurements, associated state collapse, and no-cloning.

The focus of this paper is on the efficient estimation of the gradient with polynomial complexity with respect to the above three measures.   The well-known PSR [32, 28] relies on the Hadamard test with Pauli operators to estimate the partial derivatives without any change to the ansatz. This is done via a simple circuit appended to the output [28].  However, PSR is restricted to gates with two distinct eigenvalues.
Therefore, one may ask the following question: is it possible to estimate the gradient using the same Hadamard tests as in PSR for general ansätze and without any change to their circuits?

Existing works introduce alternative approaches to generalize PSR [3, 41, 37]. Such methods require several changes to the ansatz circuit configuration or have high classical complexity.
In this work, we give an affirmative answer to the above question when the Hamiltonian has certain algebraic symmetries characterized through group theoretic structures defined for Pauli matrices.   We show that the gradient can be estimated in polynomial classical time and a linear number of Hadamard tests for Pauli strings.

§.§ Summary of the main results

We consider a generic parameterized unitary of the form $e^{i  (A(\parameter) +B)}$,  where $B$ is a fixed Hermitian component and $A$ is the parameterized component with $p$ parameters. Our approach is based on a series of Hadamard tests on the output of the ansatz followed by polynomial time classical post-processing.

The Hadamard test  corresponding to a Pauli string $P_j$ measures the following quantity

\begin{equation}\label{eq:D_i}
D_j := i\tr{O[P_j, \rho^{out}]},
\end{equation}

where $\rho^{out}=U(\parameter)\rho U(\parameter)^\dagger$ is the output state of the ansatz. This observable can be implemented using a circuit with a few control rotation gates (see Figure <ref>).
Mitarai 28 used the Hadamard test to estimate the derivative when the ansatz has a simple form $U(\theta) = e^{i\theta P_i}$.

However, it is not clear whether it works when the ansatz is the exponential of a generic Hamiltonian $A(\parameter)$. The main difficulty is due to the non-commuting terms appearing in $A(\parameter)$.

We make use of the powerful theory of Lie algebra that offers insightful perspectives in quantum physics and increasingly becoming important in quantum computing. It enables to capture of the essential features of the underlying symmetries and can be used to analyze the spectrum, eigenstates, and dynamics of quantum systems. Typically, a Hamiltonian is described by a linear combination of terms that correspond to a certain physical interaction. Such terms can be used to generate a Lie algebra, which is called the Hamiltonian algebra or DLA [34, 42]. Often, the DLA is a sub-algebra of $\mathfrak{su}(2^d)$ which is the vector space of all skew-Hermitian $2^d \times 2^d$ matrices with the standard commutator.

We first consider a slightly simpler Hamiltonian of the form $A(\parameter) = \sum_{i=1}^p a_i P_i$, where $a_i\in \RR$ and $P_i$'s are Pauli strings.  We define a group
structure on the Pauli operators $P_i$.  We show that when the Pauli terms in the decomposition of $A(\parameter)$ form a subgroup, then Hadamard tests with proper classical post-processing can be used to estimate the gradient $\nabla \Loss$. Our first result is  abbreviated as follows.

Suppose the Pauli strings appearing in  $A(\parameter)$ are closed under the commutation, that is $[P_i, P_j]=P_k$ for some Pauli string appearing in $A(\parameter)$ up to a constant for all $i,j,k\in \{1,\cdots, p\}$. Then, there is an algorithm that estimates $\nabla \Loss(\parameter)$ with an element wise additive error $\epsilon$ using $\tilde{O}(\frac{p}{\epsilon^2})$ Hadamard tests and $O(p^3+pd)$ classical time[The notation $\tilde{O}$ hides  logarithmic factors.].

For a complete statement of the result see Theorem <ref> and Theorem <ref> in Section <ref> and <ref>.
A Hamiltonian with mutually commuting Pauli strings is a trivial example. However, the Pauli strings in the above theorem do not commute in general. Moreover,  the ansatz is not necessarily local, though, a special example is a $k$-junta Hamiltonian, that is when $A$ acts non-trivially on at most $k$ out of $d$ qubits.

Overview of the techniques. Our theoretical results rely on the interplay among various representations captured by Lie algebra, and group theoretic structures. We first make a connection between the partial derivatives of $\Loss(\parameter)$ to the   adjoint operator which characterizes the derivative of a differentiable operator exponential. Then, we represent the expressions using the stabilizer formulation that has been studied in quantum error correction [6, 14]. This representation enables us to make a connection between the matrix differential in Lie algebra and specific group actions involving the vectors over the Klein four-group and $\ZZ_2\oplus \ZZ_2$. With this framework, the partial derivatives of $\Loss(\parameter)$  can be written as an infinite-length linear combination of expectation values of  Hadamard tests on the output of the ansatz. Then, we show that the number of Hadamard tests in this decomposition is bounded by the dimensionality of the Hamiltonian algebra. Next, we show that such an infinite summation can be written as a classical matrix exponential with a size equal to the dimensionality of the Hamiltonian Lie algebra. Hence, when the DLA dimensionality is polynomial in $d$, the number of qubits, the partial derivative can be computed efficiently with $d$ Hadamard tests and $\poly(d)$ classical time.  The Hadamard tests are similar to those developed for single-gate ansätze.  Moreover, the existing results on statistical estimations of quantum systems, such as shadow tomography can be used to measure the partial derivatives simultaneously and to reduce the costs.

Next, we study a more general Hamiltonian structure where $A(\parameter)$ is decomposed into generic Hermitian terms. We show that when the Hamiltonian DLA has  $\poly(d)$ dimensionality, then the gradient $\nabla\Loss$ can be estimated polynomially. This brings us to the following result:

Suppose the variational terms in the Hamiltonian $A(\parameter)$ belong to a sub-DLA with $\poly(d)$ dimensionality, then there is an algorithm that estimates $\nabla \Loss(\parameter)$ with $\poly(d)$ tests and additional $\poly(d)$ classical time.

The complete statement is given as Theorem <ref> is Section <ref>.
Hamiltonians with polynomial-size DLA are especially important in the context of avoiding the barren plateaus [23, 11] corresponding to the exponential diminish of the gradient in the training of VQA [26, 7].  Recently it was shown that the variance of the gradient is inversely proportional to the dimension of the DLA [11], implying that VQAs with polynomial-size DLAs may not exhibit barren plateaus. This motivates us to restrict or attention to polynomial-size DLAs. Nevertheless, there still is a curiosity to understand the gradient estimation when DLA dimensionality grows exponentially with $d$. A classification of dynamical Lie algebras and bounds on their dimensionality has been thoroughly studied in [42]. Here, we provide an example for more intuition.

The transverse-field Ising model is an example of a Hamiltonian such that the dimensionality of the DLA is $O(d^2)$ with the Hamiltonian characterized as

\begin{align*}
H = \sum_{i=1}^d Z_iZ_{i+1} + X_i,
\end{align*}

where $X_i, Z_i$ are the corresponding Pauli operators.

This paper's results are steps toward a unified framework for efficiently estimating the gradient of an arbitrary parameterized circuit without making changes to its unitary. Moreover, the Hadamard tests allow one to further enhance the estimations using existing methods for estimating the expectation values of observables. We highlight some of these implications below.

Shadow tomography. Turning the gradient estimation to a series of Hadamard tests has another benefit that can further reduce the number of shots to $O(\log p)$. This can be done using the classical shadow tomography [17] — a procedure to estimate several observables with minimal sample complexity.  Suppose the observable $O$ in (<ref>) is $k$-local, meaning that it decomposes into a finite sum of observables acting non-trivially on at most $k$ qubits. Then, we can prove that the number of shots scales logarithmically with $p$.

In Theorem <ref>, suppose additionally that the observable $O$ is $k$-local. Then there exists an algorithm that estimates $\nabla \Loss(\parameter)$ with $O(\frac{1}{\epsilon^2} 4^k \log p)$ ansatz uses and additional $\tilde{O}(p^3 2^{\Theta(k)}+pd)$ classical time.

This result relies on the observation that when $O$ is $k$-local, then so are the Hadamard tests in (<ref>). In that case, one can use classical shadow tomography for local observables. See Corollary <ref> in Section <ref> for detailed statement of the result.

Joint measurability.
One benefit of performing the estimations via the proposed approach is the applicability of measuring groups of the partial derivatives jointly. The Hadamard tests in (<ref>) are observables some of which are jointly measurable. For instance, the Ising model in Example <ref> has several terms that are commuting with each other.  Joint measurability helps reduce the ansatz use and sample complexity for estimating the gradient.  For that, the Pauli strings can be grouped into commuting collections. Then, the strings in each collection can be measured simultaneously with a single reference measurement, as they all can be diagonalized by a single unitary.

Suppose $\mathcal{F}_1,\cdots, \mathcal{F}_m$ form the mutually compatible groups of the Pauli strings appearing in a parameterized Hamiltonian $A(\parameter)=\sum_i a_i P_i$. Then, the number of joint Hadamard tests in Theorem <ref> can be reduced to $\tilde{O}(\frac{m}{\epsilon^2})$.

Jointly measuring Pauli strings has been studied extensively in the literature and several methods have been introduced for grouping the Paulies [22, 39, 8, 2].

§.§ Comparison with related methods

We consider the three complexity measures we mentioned earlier.  Firstly, the number of distinct circuits that need to be evaluated to obtain all the partial derivatives. Secondly,  the overall number of measurements or circuit shots. Thirdly, the classical post-processing time. Although this measure is less restrictive, it is important to ensure it scales polynomially with the number of qubits. Our approach in the context of Theorem <ref> and Corollary <ref> requires no change to the ansatz,  with  $\log(p)$ Hadamard unique circuits, and $\poly(p, d)$ classical time. Below, we highlight some of the most relevant approaches for comparison to our work.

Stochastic PSR:
This method is a generalization of PSR [3, 41], where each partial derivative is written as an integral, and a Monte Carlo strategy is used to estimate it.  For an ansatz with $p$ parameters, it requires $p$ unique circuits each with changes to the structure of the ansatz.
In practice, such modifications require a re-evaluation of the schedule of the underlying quantum-control system and hence are at a disadvantage. Moreover, the stochastic PSR has a high estimation variance. This is because an integral is estimated by sampling values of its integrand with a finite-shot estimate. David [41] presented a neat connection to the Discrete Fourier series and introduced a method that is efficient when the Hamiltonian $A$ is promised to have equidistant eigenvalues. However, for a generic ansatz one first needs to compute the spectral decomposition of $A$ to find the pattern of the parameter shifts. In that case, this process would take an exponential classical time as $A$ is an exponentially large matrix.

Nyquist PSR: Recently  [37] proposed a  shift rule for PQC where only the parameters are shifted without any other modifications of the ansatz. The method relies on a beautiful connection between the Nyquist-Shannon Sampling theorem and the Fourier series that was observed earlier in [41, 38]. The number of unique circuits for this estimation scales with $p$ and the difference between the maximum and minimum eigenvalues of $A$ — a quantity bounded by the operator norm $\norm{A}$.   As the authors reported, this method has low approximation error when the parameter value is large enough. More precisely, the approximation error is $O(\frac{1}{c^2})$ as long as $\theta = (1 - \Omega(1))c$, where $c$ is the maximum magnitude of a parameter value.

Approaches based on the special unitary group.
This is another approach [43] based on Lie algebra and a nice connection to the geometry of $SU(2^d)$ matrices and the adjoint operator. The paper argues that the number of unique parameter shift rules is bounded by the dimensionality of the DLA that contains $A(\parameter)$. However, this approach has a classical run time scaling as $p 2^{\theta(d)}$ as it relies on finding the Jacobian matrix of $U(\parameter)$.

Ansatz Changes     Shots     Classical Complexity

SPSR       $p$     $O(p)$     $O(p)$      [3]

NPSR      $\tilde{O}(p\norm{A})$     $\tilde{O}(p\norm{A})$     $\tilde{O}(p\norm{A})$     [37]

$\mathit{SU}(N)$      $p$     $O(p)$     $p \exp(\Theta(d))$      [43]

Theorem 1     $m_A$     $\tilde{O}(m_A)$     $O(p^3)$     Sub DLA HVA

Corollary 4     $0$     ${O}(4^k\log p)$     $O(p^32^{\Theta(k)})$     Sub DLA HVA + $k$-local $O$

Theorem 2     $m_A$     $\tilde{O}(m_A)$     $\poly(d_{DLA})$     general HVA

Corollary 4     $0$     $O(4^k\log d_{DLA})$     $\poly(d_{DLA}, 2^k)$     general HVA + $k$ local $O$

Comparison of various methods for estimating the gradient of a PQC $e^{i(A(\parameter)+B)}$ with $p$ parameters. $m_A$ in the last two rows refer to the number of compatibility groups of terms appearing in $A$ and  $m_A\leq p$. $d_{DLA}$ is the dimensionality of the sub-DLA spanned by the terms in $A$ and is assumed to be $\poly(n)$. For presentation convenience, we ignored the detailed expression including $\frac{1}{\epsilon^2}$ and $\norm{O}_\infty$.

§.§ Summary of the proof techniques (Move to the main text)

To view how the gradient is estimated first note that the the partial derivative of $\Loss$ in (<ref>) can be written as

\begin{align*}
\pdv{\Loss(\parameter)}{a_i} = \tr{{O} ~ \pdv{a_i}( U(\overrightarrow{a})\rho  U(\overrightarrow{a})^\dagger)}.
\end{align*}

Adjoint operator and the partial derivative:
As the last piece of the puzzle, we note a known result in Lie algebra for the connections between the derivative of a differentiable operator exponential and its adjoint:

\begin{align*}
\dv{\exp{X(\tau)}}{\tau} = \exp{X(\tau)} \frac{1-\exp{-\ad_X}}{\ad_X}\dv{X(\tau)}{\tau}.
\end{align*}

Acceleration via a matrix formulation Next, we show that when the Pauli operators of $A$ form a subgroup of $\CK_4^d$, then one can efficiently compute the operator $\frac{1-\exp{-\ad_X}}{\ad_X}$ in a classical computer.  More precisely, we have the following result
Explain how this derivation is challenging and important.

If the Pauli operators of $A(\parameter)$ form a subgroup, then there exists an algorithm to compute the effect of  $\frac{1-\exp{-\ad_X}}{\ad_X}$ in $O(p^2 d)$ classical time and $O(p^2)$ space.

§.§.§ Exploiting group structures.

We use an effective representation to exploit group theoretic structures in the Pauli decomposition of the Hamiltonian $A(\parameter)$. Particularly we use the stabilizer formulation that has been studied to represent the stabilizer codes in quantum information theory [6, 14].

Stabilizer formulation and isomorphism to $(\ZZ_2\oplus \ZZ_2)^2$.
In this formulation, the Pauli strings on $d$-qubits for a group with operator-product as the group action. An appropriate quotient of such group is known to be isomorphic to the binary group  $(\ZZ_2\oplus \ZZ_2)^d$. As an example, the identity along with the single-qubit Pauli operators have the following  representation in $\ZZ_2\oplus \ZZ_2$:

\begin{align*}
I=(0|0), \qquad X=(0|1), \qquad Y=(1|0), \qquad Z=(1|1).
\end{align*}

Next, we equip this space with a symplectic inner product denoted by $\<\cdot, \cdot\>$. Lastly, we proceed by defining a signed product for   $(\ZZ_2\oplus \ZZ_2)^d$ as

\begin{align*}
\bfs \circledast \bfr : = (\<\bfs , \bfr\> \mod 2)  \varepsilon_{\bfs, \bfr}.
\end{align*}

Klein four-group representation. Note that the Pauli strings can be identified by vectors in $\set{0,1,2,3}^d$. Hence, we consider the structure of  the Klein four-group over the set $\set{0,1,2,3}$ with $\circ$ denoting the group action. Next, we consider its vector version $\CK_4^d$. We are interested in Hamiltonians where its Pauli operator decomposition forms a subgroup of $\CK_4^d$.

Connections to Lie Algebra.
Next, we make a connection to the Lie algebra and write the Lie bracket in the above notation.  With the above formulations we present the connection to the Lie bracket and the adjoint operator.

remove the lemma

For any elements $A,B$ of the Lie algebra $\g$, the corresponding adjoint decomposes as

\begin{align*}
\ad_A(B) = 2i\sum_{\bfr, \bfs \in \CK_4^d} \qas b_\bfr (\bfs \circledast \bfr)  \sigma^{\bfs \circ \bfr},
\end{align*}

where $\qas:=\frac{1}{2^d}\tr{A\qps}$ and $b_\bfr:=\frac{1}{2^d}\tr{B\sigma^\bfr}$ are the Pauli coefficients of $A$ and $B$, respectively.

Give them some idea.

§ PRELIMINARIES AND MODEL

Notation.  For any $d\in \NN$, let $H_d$ be a the Hilbert space of $d$-qubits.
By $\mathcal{B}(H_d)$ denote the space of all bounded (linear) operators acting on $H_d$.  The identity operator is denoted by $\id$.
As usual, a quantum state, in its most general form, is denoted by a density operator; that is a Hermitian, unit-trace, and non-negative linear operator.

A quantum measurement/observable is modeled as a POVM. A POVM $\mathcal{M}$ is represented by a set of operators $\mathcal{M}:=\{M_v, v\in\mathcal{V}\}$, where $\mathcal{V}$ is the (finite) set of possible outcomes. The operators $M_v$ are non-negative and form a resolution of identity, that is $M_v=M_v^\dagger\geq 0,  \sum_v M_v =\id.$

The operator norm (infinity norm) for an operator $A$ is defined as
$$\norm{A}= \sup_{\ket{\phi}}\norm{A\ket{\phi}}.$$

§.§ Variational Quantum Algorithms

VQA are used for solving problems in optimization, learning, and simulation.
They consist of an ansatz as a parameterized quantum operation to be tuned in a quantum-classical hybrid loop. The ansatz is modeled as a (unitary) $U(\overrightarrow{a})$, with $\overrightarrow{a}\in \RR^p$ being the vector of adjustable model parameters, and $p$ the number of the parameters. It can be constructed explicitly by concatenating multiple layers of smaller parametric units. This structure is sometimes referred to as a QNN. To ensure computational tractability, it is assumed that $p=\poly(d),$ with $d$ being the number of qubits.
Recall that the objective of a VQA is to minimize the loss given in (<ref>).

Note that the loss is an expectation value and even the mere act of measuring the loss causes state collapse.  Consequently, one expects that quantum optimization and learning problems are more sample-intensive compared to their classical analogs. This is particularly problematic when the quantum states are expensive to produce.

Making iterative progress in the direction of the steepest descent is one of the most popular optimization techniques in VQA, as it has been in classical problems. Ideally, a gradient descent optimizer applies the following update rule at each iteration $t$:

\begin{align}\label{eq:QSGD ideal}
\overrightarrow{a}^{(t+1)} = \overrightarrow{a}^{(t)} - \eta_t \nabla \Loss(\overrightarrow{a}^{(t)})
\end{align}

where $\eta_t \in \RR$ is the learning rate at iteration $t$.

The above update rule is not realistic as the objective function $\Loss(\overrightarrow{a})$ is an expectation value, and the characteristics of $\rho$ is either unknown or computationally not tractable. Hence, one needs gradient estimates —  an extensively studied topic in the literature. A useful approach is one that is efficient in both computational and sample complexity.   Several approaches have been introduced to efficiently estimate the gradient of the loss under various restrictions on the ansatz.

Product ansätze.
A special variant of the ansatz is a concatenation of single parametric or non-parametric unitary circuits of the form

\begin{align*}
U(\parameter) = U_L(a_L)V_L \cdots U_1(a_1)V_1,
\end{align*}

where $V_l$ is non-parametric and each parametric layer is of the form

\begin{equation}\label{eq:ansatz U_l}
U_l(a_l):= e^{i a_{\bfs_l} \sigma^{\bfs_l}}
\end{equation}

with $\sigma^{\bfs_l}$ being a Pauli string. There are several solutions for VQAs using product anästze.

Parameter Shift Rule.
There have been several approaches to estimating the gradient [12, 28, 36, 16, 18, 32, 27, 43]. The zeroth-order approach (e.g., finite differences) evaluates the objective function in the neighborhood of the parameters.  Although it is a generic approach, recent studies showed their drawbacks in terms of convergence rate [18].  First-order methods (e.g., parameter shift rule)  directly calculate the partial derivatives [32]. When the ansatz is of the form $U(\parameter) = \prod_{j=1}^p e^{-i a_j G_j}$, with $G_j$ the parameter shift rule implies that:

\begin{equation*}
\pdv{\Loss}{a_j} = \<\Loss(\parameter + \frac{\pi}{4}e_j)\> -\<\Loss(\parameter - \frac{\pi}{4}e_j)\>,
\end{equation*}

where $e_j\in \RR^p$ is the $j$th canonical vector, that is $e_{j,j}=1$ and $e_{j,r}=0$ for all $r\neq j$.

$\ket{\phi}$   U_≤l    R__l(+π/2)      R__l(-π/2)        U_>l      [r]   $\hat{y}$

$\ket{0}$    H    -1     X    -1      X         [r]   z

Hadamard test for measuring the partial derivative of product with respect to a parameter $a_{\bfs_l}$ appearing at layer $l$. Here $U_{\leq l}$ corresponds to the first $l$ layers of the ansatz, and $U_{>l}$ to the remaining layers. Here, $X$ is the X-gate and $R_{\bfs_l}$ is the controlled rotation around Pauli $\sigma^{\bfs_l}$.

As shown in [28], the partial derivatives can be directly measured via a Hadamard test (see Fig. <ref>). More formally, the partial derivative of a product ansätze can be written in terms of commutators with associated Pauli operators. This statement is summarized below:

Let ${\rho}_l^{out} =U_{\leq l}\ketbra{\phi} U_{\leq l}^\dagger$ denote the density operator of the output state at layer $l$ when the input is $\ket{\phi}$ with label $y$. Then, the derivative of the loss is given by:

\begin{align}\label{eq:losss derivative}
\frac{\partial  \Loss%( a_{\bfs}, \ket{\phi}, y)
}{\partial a_{\bfs_l}} =   \tr{O U_{>l} \big[ \sigma^{\bfs_l}, {\rho}_l^{out}\big]U^\dagger_{>l}},
\end{align}

where $[\cdot, \cdot]$ is the commutator operation.

The expression above can be implemented using a quantum circuit with measurements at the end. The circuit is shown in Fig. <ref> which is a special example of the generalized Hadamard test [27, 18].  In this procedure, given a sample $\ket{\phi}$ and an ancilla qubit $\ket{0}$, we first apply the first $l$ layers of the ansatz and then apply a special circuit with controlled rotations for taking the derivative of the loss. Then, the rest of the layers of the ansatz are applied and measurement is performed.

§.§.§ Generic PQC

A generic PQC  is an ansätze of the form

\begin{align}\label{eq:PQC nonproduct}
U(\parameter) = e^{i (A(\parameter))},
\end{align}

where  $A(\parameter)$ is the parameterized  (traceless) Hermitian operator. For instance,  $A$ can be decomposed in terms of Pauli strings as
$$A(\parameter)=\sum_\bfs f_{\bfs}(\parameter) \qps,$$
where $\bfs\in \set{0,1,2,3}^d$ and $f_{\bfs}:\RR^p \rightarrow \RR$ is a real-valued function on the space of the parameters. The Pauli decomposition states that

Note that the ansatz can be multi-layered, but, for simplicity, we consider only a single-layer PQC.
Due to the non-commutativity of the Pauli products, Fact <ref> does not apply. A brief summary of some of the relevant approaches is provided in Appendix <ref>.

§ MAIN RESULTS

In what follows we present the main results of the paper. We introduce an approach for estimating the gradient of the loss for a generic ansatz. We show that the Hadamard test followed by classical post-processing leads to an approximation of gradient.
We prove in Theorem <ref> of Section <ref> that the partial derivative is a linear combination of several terms similar to the product ansätze as in Fact <ref>.
The coefficients in the linear combination are calculated through a series of group actions that we define over the stabilizer representation of Pauli strings.
Then, in Theorem <ref> in Section <ref> we show that under certain group structures, the linear combination collapses into $\poly(p)$ terms. In Section  <ref>, we extend this result to Hamiltonians with polynomial size DLA.
Before presenting the To present the main results of the paper we need to introduce a few concepts and notations.

§.§ The Stabilizer Formulation

We first define a few actions on the Pauli strings. Basic definitions of group and Lie algebra are given in Appendix <ref>.

First, the $d$-qubit Pauli group, denoted by  $\mathcal{P}_d$, is defined as the set of all Pauli strings of the form $c\qps$, where $c \in \{\pm, \pm i\}$ and $\bfs\in \{0,1,2,3\}^d$ equipped with the matrix product as the group action. Ignoring the scalar $c$, the group quotient is isomorphic to the Binary vector group  $\ZZ_2^{2d}$. Therefore, the Pauli strings abide a compact representation which we present in the following.

The representation associates the identity and each Pauli operator $X,Y,Z$ with an element of $\ZZ_2\oplus \ZZ_2$ which is the binary vector space with dimension two:

\begin{align*}
I=(0|0), \qquad X=(0|1), \qquad Y=(1|0), \qquad Z=(1|1),
\end{align*}

where we use  $(\cdot |\cdot)$ to denote the elements of $\ZZ_2\oplus \ZZ_2$.

For the $d$-qubit Pauli strings, consider the binary vector space $\ZZ_2^{2d}$, where each element is denoted by $(\bfa | \bfb), \bfa, \bfb\in \ZZ^d_2$. With this formulation, we can represent each Pauli string as

\begin{align*}
\qps \equiv (\bfs^0 | \bfs^1),\qquad \forall \bfs \in \{0,1,2,3\}^d, \quad \forall(\bfs^0 | \bfs^1) \in (\ZZ_2\oplus \ZZ_2)^d,
\end{align*}

where $s_j^0$ and $s^1_j$ give the most significant and least significant bit of the binary representation of $s_j\in \set{0, 1, 2 ,3}$.  Next, we define two operations for $(\ZZ_2\oplus \ZZ_2)^d$: the element-wise binary addition and a symplectic inner product given by

\begin{align*}
\<(\bfa | \bfb), (\bfa'| \bfb')\> = \sum_{j=1}^d a_jb_j'+b_ja_j', \qquad \forall (\bfa | \bfb), (\bfa'| \bfb')\in (\ZZ_2\oplus \ZZ_2)^d.
\end{align*}

We need to define another group action.

Klein four-group representation.  Let $\CK_4$ be the Klein four-group over the set $\set{0,1,2,3}$ with $0$ being the group's identity and $\circ$ denoting the group addition. The group operation for  distinct elements $a,b,c \in \CK_4-\{0\}$ is defined as

\begin{align*}
0\circ a =a \qquad a\circ a = 0,\qquad a\circ b = b \circ a = c.
\end{align*}

It is known that $ \ZZ_2\times \ZZ_2\cong \CK_4$. Hence, we can view the $\circ$ as an operation over $\ZZ_2^2$. Next, by  $\CK_4^d$ denote the usual vector of $\CK_4$  with the element-wise addition $\circ$.

Next, we define a special operation on $\CK_4$. For a pair of $s, r \in \CK_4$ define

\begin{align*}
s \odot r := \mod_2(\<(s^0| s^1), (r^0| r^1)\>) = \mod_2(s^0r^1+s^1r^0),
\end{align*}

where $(s^0,s^1)$ and $(r^0,r^1)$ are, respectively, the binary representation of $s$ and $r$. We  extend this notation for vectors  $\bfs, \bfr \in \CK_4^d$ as

\begin{align*}
\bfs \odot \bfr := \mod_2\qty(\<(\bfs^0| \bfs^1), (\bfr^0| \bfr^1)\>) = (s_1 \odot r_1)\oplus \cdots \oplus (s_d \odot r_d)
\end{align*}

Signed product. Lastly, we proceed by defining a signed product over $\CK_4^d$. For any pair $i, j \in \CK_4$, define

\begin{align*}
\varepsilon_{i,j}:=\begin{cases} -1 &  \text{if}~ (i,j) = (2,1), (3,1), (3,2)\\
1 & \text{otherwise} \end{cases}.
\end{align*}

With these notation, the signed product of any pair $\bfs, \bfr \in \CK_4^d$ is given by

\begin{align*}
\bfs \circledast \bfr : = (\bfs \odot \bfr) \prod_{j=1}^d \varepsilon_{s_j, r_j}.
\end{align*}

As an example, consider the single qubit case, where the elements take values from $\{0,1,2,3\}$.  Table <ref> demonstrates the $\circledast$ operation between pairs of such members.

$s/ r$     0     1     2     3

0     0     0     0     0

1     0     0     1     1

2     0     -1     0     1

3     0     -1     -1     0

Table of the product $s\circledast r$ for the single qubit case.

Now, given the definitions for $\circ, \circledast$ and $\odot$, we are ready to present our main results.

§.§ Partial Derivative

Our first result is the decomposition of the partial derivatives into the expectation of the Hadamard tests with Pauli strings as in the product ansätze.

Consider an ansatz of the form $U(\overrightarrow{a})=\exp{i A(\parameter)}$, where $A(\parameter)=\sum_{\bfs\in \mathcal{S}} \qas \qps$ for some $\mathcal{S}\subseteq \CK_4^d$.  For any $\bft\in \CK_4^d$, define
$$D_\bft:=i \tr{O \big[ \sigma^\bft, {\rho}^{out}\big]},$$
where $\rho^{out}$ is the output state.   Then, the  derivative of the objective function $\Loss$  for this ansatz expands in terms of  $D_\bft$ as:

\begin{align}\label{eq:pdv infinite sum}
\frac{\partial \Loss}{\partial a_{\bfr}} = \sum_{k=0}^\infty \frac{(-2)^k}{(k+1)!} \sum_{\bfs_1 \in \mathcal{S}} \cdots \sum_{\bfs_k \in \mathcal{S}}  \prod_{j=1}^k a_{\bfs_j} \qty( \bfs_j \circledast (\bfs_1 \circ \cdots \circ \bfs_{j-1}\circ \bfr))   D_{\bfs_1 \circ \cdots \circ \bfs_k \circ \bfr},
\end{align}

for any $\bfr\in \mathcal{S}$.

The proof of the theorem is delayed until Section <ref>. With this theorem, one can compute the derivative of the loss for a non-product ansatz by treating it as a single-parameter ansatz as in Fact <ref> but with correction terms. In other words, the derivative of loss can be measured by applying the standard Hadamard test followed by classical corrections. To highlight this result we present an example concerning a single qubit ansatz.

Consider a general single-qubit unitary of the form

\begin{equation*}
U(\overrightarrow{a}) = \exp{ia_1 \sigma^1 + ia_2 \sigma^2 +a_3\sigma^3}.
\end{equation*}

Let $O$ be a generic observable and consider the associated loss $\Loss(\parameter)$ as in (<ref>). As an example, consider taking the partial derivative of $\Loss$ with respect to $a_1$ at the point $a_1=0$ and $a_3=0$. In the context of Theorem <ref>, let $D_j$ be the result of the Hadamard test with Pauli $\sigma^j$, where $j=1,2,3$. Then, (<ref>)  in Theorem <ref> simplifies to the following:

\begin{align*}
\frac{\partial \Loss}{\partial a_1}(\parameter = (0,a_2,0)) = \sum_{k=0}^\infty \frac{(-2)^k}{(k+1)!}  \prod_{j=1}^k a_{2} \qty( 2 \circledast (\underbrace{2 \circ \cdots \circ 2}_{k-1 ~\text{times}}\circ 1))   D_{\underbrace{2 \circ \cdots \circ 2}_{k ~\text{times}} \circ 1},
\end{align*}

where we used the fact that only terms with all $\bfs_j = 2$ are surviving.
Note that for even $k$, $\underbrace{2 \circ \cdots \circ 2}_{k ~\text{times}} \circ 1 = 1$, and for odd $k$ it equals to $2\circ 1 = 3$. Therefore, we have that

\begin{align*}
2 \circledast (\underbrace{2 \circ \cdots \circ 2}_{k-1 ~\text{times}}\circ 1) = \begin{cases}   2 \circledast 3 = 1 & \text{even}~ k\\
2 \circledast 1 =-1 & \text{odd}~ k
\end{cases}
\end{align*}

where we used Table <ref>. Plugging such evaluations in the first equation, we have

\begin{align*}
\prod_{j=1}^k a_{2} \qty( 2 \circledast (\underbrace{2 \circ \cdots \circ 2}_{k-1 ~\text{times}}\circ 1))  = a_2^k (-1)\times 1 \times (-1) \times \cdots
&= \begin{cases} a_2^k  (-1)^{k/2} & \text{even}~ k\\
a_2^k  (-1)^{(k+1)/2}  & \text{odd}~ k
\end{cases}
\end{align*}

\begin{align*}
\frac{\partial \Loss}{\partial a_1}(\parameter = (0,a_2,0)) =
\sum_{p=0}^\infty \frac{(-2)^{2p}}{(2p+1)!} a^{2p}_{2} (-1)^p  D_{1} + \sum_{q=0}^\infty \frac{(-2)^{2q+1}}{(2q+2)!} a^{2q+1}_{2} (-1)^{q+1}  D_{3}
\end{align*}

Hence, one needs to measure $D_1$ and $D_3$ to compute the above partial derivative.  Next, by simplifying the summations, it is not difficult to show that

\begin{align*}
\frac{\partial \Loss}{\partial a_1}(\parameter = (0,a_2,0))  &=   \frac{-1}{2a_2} \qty( \sum_p \frac{(-2a_2)^{2p+1}}{(2p+1)!} (-1)^p )D_1 + \frac{-1}{2 a_2} \qty( \sum_q \frac{(-2a_2)^{2q+2}}{(2q+2)!} (-1)^{q+1} )D_3 \\
&= \frac{-1}{2a_2} \qty\big( \sin(-2a_2) D_1  +  (\cos{(-2a_2)}-1) D_3) \\
&= \frac{1}{2a_2} \qty\big( \sin(2a_2) D_1  +  (1- \cos{(2a_2)}) D_3).
% &= \emph{sinc}(a_2)\qty( \cos{(a_2)}~ D_1 + \sin{(a_2)} ~D_3 ).
\end{align*}

Notice the presence of $D_3$ which relates to the Pauli $\sigma^3$ and not appear in the ansatz expression.
One can verify that this is indeed equal to the analytic gradient of this ansatz (see Appendix <ref> for more details).

In this example, we can see how our method is akin to computing the partial derivative in a closed-form expression. Interestingly, under certain structural properties, one can find the closed-form expression of the summation in Theorem <ref>. The power of this approach lies in the group structure of the Pauli matrices, which produces predictable patterns. We highlight examples of such structures in the next section.  Another point worth mentioning is that the terms $D_i$ only make use of $\rho^{out}$ and hence are easily measurable.

§.§ Ansätze with Pauli Subgroup Structures

We consider special cases of the Hamiltonian for which the gradient can be computed efficiently. Note that when $\mathcal{S}$ is a subgroup of $\CK_4^d$, then the terms appearing in (<ref>)  remain inside the subgroup.  In that case, the partial derivative can be written as follows.

Suppose the Hamiltonian parameter set $\CS$ is a subgroup of $\mathcal{K}_4^d$.  Then, there is an algorithm that computes $\nabla \Loss$  in $O(p^3+pd)$ time with $O(p)$ use of the Hadamard test for Pauli strings as in (<ref>).

Note that here the run time depends on the number of relevant qubits that the ansatz acts. When the ansatz is a $k$-junta then the the run time will be $O(p^3+pk)$.

Such examples might not be interesting as the overall unitary can be efficiently simulated in a classical computer.  There are several subgroups of $\CK_4^d$ for which the Hamiltonian is not $k$-junta. A simple example is
$$A = a_1X^{\tensor d} + a_2 Y^{\tensor d}+ a_3 Z^{\tensor d},$$
where the subgroup is isomorphic to $\CK_4$ inside $\CK_4^d$.

Proof sketch: To prove this theorem, we define a matrix representation of the derivative expression.
Let $\CS\leq \CK_4^d$ be a subgroup of size $p$ with elements denoted by $\bfs_1, \cdots, \bfs_p$. Note that  since $\CS$ is a subgroup, the infinite sum in the theorem reduces to a linear combination of their form

\begin{align*}
\frac{\partial \Loss}{\partial a_{\bfr}} = \sum_{\bfs\in \CS}  g_\bfs(\bfr)  D_{\bfs},
\end{align*}

where $g_\bfs(\bfr)\in \RR$. Therefore, if we consider each $D_{\bfs}$ as a basis in a vector space $\RR^p$, then $\pdv{\Loss}{\qar}$ is a vector in that space with the representation

\begin{align*}
\pdv{\Loss}{\qar} \equiv (g_{\bfs_1}(\bfr), \cdots, g_{\bfs_p}(\bfr))\in \RR^p.
\end{align*}

Moreover, with this representation each $D_{\bfs_j}\equiv \bfe_j$, where $\bfe_j\in \RR^p$ is the $j$th  canonical basis vector. Now define the  $p\times p$ matrix $V$ where the $i$th column is given as

\begin{align}\label{eq:V matrix column}
\bfv_i := 2 \sum_{j = 1}^p   a_{\bfj}  (\bfs_j \circledast \bfs_i)~   \bfe_{j \circ i},
\end{align}

where $j\circ i$ is the shorthand for the index of $\bfs_j\circ \bfs_i$. Note that the above summation can be computed in $O(p d)$ time. Moreover, noting that $\bfv_i = V\bfe_i$ and given (<ref>), the vector representation of the partial derivative of $\Loss$ equals to

\begin{align*}
\pdv{\Loss}{a_{\bfs_i}} \equiv \sum_{k\geq 0} \frac{(-1)^k}{(k+1)!} V^k \bfe_{i}, \qquad \forall i\in [p].
\end{align*}

Next, compute the following matrix
$$B = (I-e^{-V})V^{-1},$$
where $V^{-1}$ is the generalized inverse of $V$. Then, it is not difficult to check that
$B$ gives the vector representation of the partial derivatives. With a proper matrix exponentiation algorithm, this matrix is computed in $O(p^3)$. To compute the gradient from the vector representation of the derivative, one first needs to estimate all $D_{\bfs_i}, i\in [p]$. This can be done with $O(p)$ runs of the derivative-taking circuit of Figure <ref>. Lastly, the gradient is computed by multiplying the matrix with the vectors of $D_{\bfs}$'s as

\begin{align*}
\nabla\Loss = B \overrightarrow{D}_\bfs.
\end{align*}

This procedure is summarized in the following algorithm. An advantage of this method is that all the partial derivatives are estimated via the matrix exponentiation of  $V$.

Note that empirical estimations of  $D_{\bfs}$ can be done via several measurement shots. In a straightforward approach each $D_{\bfs}$ is estimated independent measurement shots. Hence, with the standard concentration analysis (McDiarmid inequality) one can show that all the $D_{\bfs}$ in the subgroup can be estimated with probability at least $(1-\delta)$ up to an additive error $\epsilon$ with
$$O\qty(p \frac{1}{\epsilon^2}\log \frac{p}{\delta})$$
measurement shots.

Subgroup Gradient Estimation
Input: Parameter subgroup $\CS$,

[1]

Estimate the expectation value of Hadamard tests $D_\bfs, \bfs\in \CS$ as in (<ref>).
Compute the matrix $V$ where the column $i$ is computed as in (<ref>).
Compute the matrix $B = (I-e^{-V})V^{-1},$ where $V^{-1}$ is the generalized inverse of $V$.

Return $\nabla\Loss \approx B \overrightarrow{D}_\bfs$.

This figure illustrates the concept behind the subgroup gradient estimation algorithm. The expression in Theorem <ref> points to an alternative representation of the problem. Instead of picturing the gradient as a function in the space of the parameters (the left picture) we can view it as a vector in the landscape of Hadamard tests $D_\bfs$.

As a sanity check. We derive again the the same derivative expression as in Example <ref>, this time using Algorithm <ref>. One can verify that for the general single-qubit unitary we have

\begin{equation*}
V = \begin{bmatrix}
0 & 0 & 2a_2\\
0 & 0 & 0 \\
-2a_2 & 0 & 0
\end{bmatrix}.
\end{equation*}

In Example <ref>, we had $a_1 = a_3 = 0$. Using the singular value decomposition, the matrix $B=(I-e^{-V})V^{-1}$ is computed as

\begin{equation*}
B = \begin{bmatrix}
\frac{\sin(2b)}{2b} & 0 & \frac{1- \cos(2b)}{2b} \\
0 & 0 & 0 \\
\frac{\cos(2b)-1}{2b}  & 0 & \frac{\sin(2b)}{2b}.
\end{bmatrix}
\end{equation*}

Next, having the Hadamard tests $D_1, D_2, D_3$, the gradient is computed as

\begin{equation*}
\nabla \Loss(a_1=0, a_2, a_3=0)  =  B\begin{bmatrix}
D_{1}\\ D_2\\ D_3
\end{bmatrix}
= \begin{bmatrix}
\frac{\sin(2b)}{2b} D_1 +  \frac{1- \cos(2b)}{2b}D_3\\
-\frac{1- \cos(2b)}{2b}D_1 + \frac{\sin(2b)}{2b} D_3
\end{bmatrix}
\end{equation*}

The sum of which gives us the exact analytical expression for the derivative we got before. This example illustrates the connection between the space of $D_\mathbf{s}$ and the group nature of the Pauli operators.

§.§ Proof of Theorem <ref>

We start by taking the partial derivative of $\Loss(\parameter)$. Noting that $\rho^{out}:=U\rho U^\dagger$, the partial derivative of the loss with respect to $\qas$ is

\begin{align*}
\frac{\partial \Loss}{\partial \qas} = \tr{ O \frac{\partial }{\partial \qas}\rho^{out}},
\end{align*}

\begin{align}\label{eq:rhou out derivative}
\pdv{\rho^{out}}{\qas} = \frac{\partial U}{\partial \qas} \qty(\rho U^\dagger) + \qty(U\rho) \pdv{U^\dagger}{\qas},
\end{align}

and we used the fact that $\rho^{out}$ is Fréchet differentiable with respect to $\qas$.

Next, we need to introduce some ingredients in Lie algebra.
We equip the space $\GL(2^d, \CC)$ of  $2^d\times 2^d$ complex matrices with the Lie algebra and the standard Lie bracket defined as
$$[A, B]:= AB- BA,$$
for any $A, B\in \GL(2^d, \CC)$. Such Lie algebra is denoted by $\gl(2^d, \CC)$ or $\g$ for shorthand.   For more details on Lie algebra see [31]. Now, we proceed by the adjoint representation.  For any $X\in \g$, define the mapping $\ad_X: \g\rightarrow \g$ by $\ad_X(Y)= [X,Y]$. With this notation
$$\ad_X^k(Y) = [X, \cdots, [X,Y]\cdots ].$$
The adjoint operator is connected to the derivative of the matrix exponential. In fact, from Theorem 5 of [31], the differential of any exponential operator is given by

\begin{align}\label{eq:exp differential}
\dv{\exp{X(\tau)}}{\tau} = \exp{X(\tau)} \frac{1-\exp{-\ad_X}}{\ad_X}\dv{X(\tau)}{\tau}.
\end{align}

where it is assumed that $X(\tau)$ is a differentiable  (linear) operator with respect to a variable $\tau \in \CC$.

Next, we prove the following lemma which is in a different ordering of the above equation.

Suppose $U=\exp{iH(\tau)}$ for some Hermitian operator valued function $H$. Then,

\begin{align*}
\dv{U}{\tau} =  -\frac{1-\exp{i\ad_H}}{\ad_H}(\dv{H(\tau)}{\tau}) U,
\end{align*}

where it is assumed that $H(\tau)$ is differentiable.

First note that for any operator $X\in \g$

$$(\ad_X)^\dagger = \ad_{X^\dagger}.$$

We can write $U=(e^{-iH(\tau)})^\dagger$.  Therefore, given that $\dv{X^\dagger}{\tau}=(\dv{X}{\tau})^\dagger$, from (<ref>) we can write

\begin{align*}
\dv{U}{\tau} &= \qty(\dv{e^{-iH(\tau)}}{\tau})^\dagger = \qty( \exp{-iH(\tau)} \frac{1-\exp{-\ad_{-iH}}}{\ad_{-iH}}\dv{(-iH(\tau))}{\tau})^\dagger\\
&\qty( \exp{-iH(\tau)} \frac{1-\exp{+i\ad_H}}{\ad_H}\dv{H(\tau)}{\tau})^\dagger\\
&= \qty(  \frac{1-\exp{+i\ad_H}}{\ad_H}\dv{H(\tau)}{\tau})^\dagger U.
\end{align*}

Note that for any $H\in \g$ we have the following equality by its convergent power series:

\begin{align}\label{eq:exp adjoint}
\frac{1-\exp{-\ad_H}}{\ad_H} = \sum_{k=0}^\infty \frac{(-1)^k}{(k+1)!} (\ad_H)^k.
\end{align}

Therefore, applying this equation for the operator $\frac{1-e^{i\ad_H}}{\ad_H}$, the derivative equals to the following

\begin{align*}
\qty(-i\sum_{k=0}^\infty \frac{(i)^k}{(k+1)!} (\ad_H)^k (\dv{H}{\tau}))^\dagger  U =i\sum_{k=0}^\infty \frac{(-i)^k}{(k+1)!} \qty((\ad_H)^k (\dv{H}{\tau}))^\dagger  U.
\end{align*}

Note that for any $X,Y\in \g$ the following identity holds

$$(\ad_X(Y))^\dagger = - \ad_X^\dagger(Y^\dagger).$$

\begin{align*}
\dv{U}{\tau} = i\sum_{k=0}^\infty \frac{(i)^k}{(k+1)!} (\ad_H)^k (\dv{H}{\tau})  U = -\frac{1-\exp{+i\ad_H}}{\ad_H}(\dv{H}{\tau}) U.
\end{align*}

This is the desired expression.

From this lemma, we have that

\begin{align*}
\frac{\partial U}{\partial \qas} = -\frac{1-\exp{i\ad_{A}}}{\ad_{A}}(\pdv{A}{\qas}) U.
\end{align*}

This gives the first part of (<ref>).
Next, from (<ref>), the partial derivative of $U^\dagger$ can be  write as

\begin{align*}
\frac{\partial U^\dagger}{\partial \qas} = U^\dagger  \frac{1-\exp{+i\ad_{A}}}{\ad_{A}}(\pdv{A}{\qas}),
\end{align*}

where we used the fact that $A$ is Hermitian. Therefore, the partial derivative of $\rho^{out}$ equals to

\begin{align*}
\pdv{\rho^{out}}{\qas}  &= -\frac{1-\exp{i\ad_{A}}}{\ad_{A}}(\pdv{A}{\qas}) U \rho U^\dagger + U\rho U^\dagger \frac{1-\exp{i\ad_{A}}}{\ad_{A}}(\pdv{A}{\qas})\\
&= \frac{1-\exp{-i\ad_{A}}}{\ad_{A}}(\pdv{A}{\qas})\rho^{out}  + \rho^{out} \frac{1-\exp{+i\ad_{A}}}{\ad_{A}}(\pdv{A}{\qas}).
\end{align*}

By simplifying the terms in the right-hand side, the partial derivative  can be written as the commutator:

\begin{align*}
\pdv{\rho^{out}}{\qas}  = \qty[\rho^{out}, \frac{1-\exp{i\ad_{A}}}{\ad_{A}}(\pdv{A}{\qas})].
\end{align*}

\begin{align*}
\frac{\partial \Loss}{\partial \qas} = \tr{ O \frac{\partial }{\partial \qas}\rho^{out}} = \tr{ O \qty[\rho^{out}, \frac{1-\exp{i\ad_{A}}}{\ad_{A}}(\pdv{A}{\qas})]}.
\end{align*}

Next, we use the definition of $\frac{1-e^{-\ad_X}}{\ad_X}$ in (<ref>) to decompose the above summation. Setting, $X=-iA$ and from the fact that $\pdv{A}{\qas}=\qps$, the above quantity decomposes as

\begin{align*}
\frac{\partial \Loss}{\partial \qas} = -i\sum_{k=0}^\infty \frac{(i)^k}{(k+1)!}   \tr{ O \qty[\rho^{out}, (\ad_A)^k(\qps)]},
\end{align*}

where the factor $-i$ appears because of the denominator.

It remains to write the above equation in terms of the group actions we defined. For that, we need to establish a few connections between the Lie bracket and the group actions.

The Lie bracket between any pair of Pauli strings $\qps, \sigma^{\bfr}$ is given by

\begin{align*}
[\qps, \sigma^{\bfr} ] = 2i (\bfs \circledast \bfr)  \sigma^{\bfs \circ \bfr},
\end{align*}

The proof of this lemma is built upon the following lemma.

The product of any pair of Pauli strings equals the following

\begin{align*}
\qps\qpr = i^{\bfs \odot \bfr} \sigma^{\bfs \circledast \bfr}.
\end{align*}

We study single qubit products. Let $a\neq b\neq c$ be distinct non-zero  elements of $\CK_4$. Also without loss of generality suppose $a<b$. Then, it is not difficult to verify the following relations

\begin{align*}
&0 \circledast a=0, \qquad 0\circ a = a \Rightarrow \sigma^{0}\sigma^a = \sigma^a\\
&a \circledast a=0, \qquad a\circ a = 0 \Rightarrow \sigma^{a}\sigma^a = \sigma^0\\
&a \circledast b = 1, \qquad a\circ b = c \Rightarrow \sigma^{a}\sigma^a = i \sigma^c\\
&b \circledast a = -1, \qquad a\circ b = c \Rightarrow \sigma^{a}\sigma^a = -i \sigma^c,
\end{align*}

where we used the fact that $i^{-1}=-i$. Hence, we established  the  lemma for $d=1$. For general $d$, with the tensor product, we have that

\begin{align*}
\qps \qpr &= \tensor_j \sigma^{s_j}\sigma^{r_j} = \tensor_j i^{s_j \circledast  r_j} \sigma^{s_j \circ r_j}\\
&= i^{\sum_j s_j \circledast  r_j} \sigma^{\bfs \circ \bfr} = i^{\bfs \circledast \bfr} \sigma^{\bfs \circ \bfr}.
\end{align*}

Hence, with the above lemma, we have that

\begin{equation*}
[\qps, \sigma^{\bfr} ] = i^{\bfs \circledast \bfr} \sigma^{\bfs \circ \bfr} - i^{\bfr \circledast     \bfs} \sigma^{\bfr \circ \bfs}.
\end{equation*}

Hence the proof of Lemma <ref> is complete.

To extend this result to arbitrary operators, first, we need to note the Pauli decomposition.

Any bounded  operator $A$ on the Hilbert space of $d$ qubits can be uniquely written as
$$A = \sum_{\bfs \in \{0,1,2,3\}^d} a_{\bfs} ~\sigma^{\bfs},$$
where $a_\bfs\in \CC$ are the Fourier coefficients of $A$ and are given as

With this notation, we can present the following result on the adjoint operator in the Lie algebra $\g$.

For any $A,B\in \g$, the corresponding adjoint decomposes as

\begin{align*}
\ad_A(B) = 2i\sum_{\bfr, \bfs \in \CK_4^d} \qas b_\bfr (\bfs \circledast \bfr)  \sigma^{\bfs \circ \bfr},
\end{align*}

\begin{align*}
\qas:=\frac{1}{2^d}\tr{A\qps} \qquad b_\bfr:=\frac{1}{2^d}\tr{B\sigma^\bfr}
\end{align*}

are the Pauli coefficients of $A$ and $B$, respectively.

Hence, from Lemma <ref> for $\ad_A(\qps)$, we can write the final expression for the partial derivative of $\Loss$

\begin{align*}
\frac{\partial \Loss(\parameter)}{\partial \qas} = -i\sum_{k=0}^\infty  \frac{(-2)^k}{(k+1)!} \sum_{\bfs_1 \in \mathcal{S}} \cdots \sum_{\bfs_k \in \mathcal{S}}  \prod_{j=1}^k \qty( a_{\bfs_j}  (\bfs_j \circledast \bfr_j))   \tr{ O \qty[\rho^{out}, \sigma^{\bfs_1 \circ \cdots \circ \bfs_k \circ \bfr}]}.
\end{align*}

From the definition of $D_\bft$ in the statement of the Theorem, the proof is complete.

§.§ Extension to Polynomial Size DLA

Next, we extend the result of the previous section to general PQCs for which DLA has $\poly(d)$ dimensionality.

For a PQC as in (<ref>), with the unitary of the form $e^{iA(\parameter)}$ where $A(\parameter) = \sum_i a_i G_i$ with $G_i$ being traceless Hermitian operators, the set $\CG = \{G_1, \cdots, G_p\}$ is called the generators of the PQC.

Given a PQC with generators $\CG$ as in Definition <ref>, the Dynamical Lie Algebra (DLA) $\g_A$ is the subalgebra of
$\su(2^d)$ spanned by the repeated nested commutators of the elements in $\CG$, i.e.,

\begin{align*}
i\g_A := \emph{Span}_\RR\<iG_1, \cdots, iG_p\>_{\text{Lie}} \subseteq \su(2^d),
\end{align*}

where $\<\cdot\>_{\text{Lie}}$ denotes the Lie closure, obtained by repeatedly taking the nested commutators, and the span is over the real numbers. We say that $\g_A$ has polynomial size/dimensionality if $\dim(\g_A) = \poly(d)$ as a vector space.

In general, computing the DLA is computationally expensive. One viable approach involves direct construction starting from a set of generators and iteratively commuting them to discover new elements until a basis of the DLA is attained (for example see Algorithm 1 in [23]). The complexity of this approach in general is $O(2^d)$. However, this is not an issue if the dimension of the DLA to be constructed is promised to by $\poly(d)$.

Suppose that the dimension of $\g_A$  as a vector space is polynomial in $d$.

Let a PQC as in (<ref>) has the DLA $\g_A$. Then, the partial derivative of $\Loss$ decomposes as

\begin{align*}
\frac{\partial \Loss}{\partial a_{\bfr}} = \sum_{G\in \g_A}  \alpha_G(r)  i\tr{O[G, \rho^{out}]},
\end{align*}

where $\alpha_G(r)\in \RR$ are the coefficients. Moreover, if  $\g_A$ has $\poly(d)$ dimensionality and its generators have a $\poly(d)$ size decomposition in Pauli strings, then there is an algorithm that estimates $\nabla \Loss(\parameter)$ with $\poly(d)$ Hadamard tests and additional $\poly(d)$ classical time.

The proof of the Theorem follows from a similar argument as in Theorem <ref>. With such an argument we can show that

\begin{align*}
\frac{\partial \Loss}{\partial a_j} =  i\tr{ O \qty[\rho^{out}, \frac{1-\exp{i\ad_{A}}}{\ad_{A}}(\pdv{A}{a_j})]}.
\end{align*}

From the definition of $\frac{1-e^{i\ad_A}}{\ad_A}$ and the fact that $\pdv{A}{a_j}=G_j$, the above quantity decomposes as

\begin{align*}
\frac{\partial \Loss}{\partial a_j} = -i\sum_{k=0}^\infty \frac{(i)^k}{(k+1)!}   \tr{ O \qty[\rho^{out}, (\ad_A)^k(G_j)]},
\end{align*}

Note that $\ad_A$ maps $\g_A$ unto itself as its commutator ideal. That is
$$\ad_A(G) = \sum_{j=1}^{m} c_j [X_j, Y_j],$$
where $X_j, Y_j\in \g_A$. Since $\g_A$ is closed under the Lie bracket, then  $[X_j, Y_j]\in \g_A$. Hence, $\ad_A(G)\in \g_A$ and it equals to a finite sum of the basis elements $E_l$ of $\g_A$ as

\begin{align*}
\ad_A(G_j)=\sum_{l=1}^{d_\g} \alpha_{l,j} E_l.
\end{align*}

Therefore, the above summation can be rewritten as

\begin{align*}
\frac{\partial \Loss}{\partial a_j} = -i\sum_{k=0}^\infty \frac{(i)^k}{(k+1)!}  \sum_{l=1}^{d_{\g_A}} \alpha_{l,j} \tr{ O \qty[\rho^{out}, E_l]},
\end{align*}

Note that $\ad_A(G_j)$ can be viewed as a vector $(\alpha_{1,j},\cdots, \alpha_{d_{\g_A},j})$ in the $E_l$ basis. Moreover, $\ad_A$ can be viewed as a linear transformation (matrix) $T_A$ in the $E_l$ basis. This matrix is $d_{\g_A}\times d_{\g_A}$ hence is $\poly(d)$ in size. With that taken into account, the partial derivative can be written as

\begin{align*}
\pdv{\Loss}{a_{j}} \equiv \sum_{k\geq 0} \frac{(-i)^k}{(k+1)!} T_A^k G_{j}, \qquad \forall j\in [p].
\end{align*}

Where $G_j$ is understood as a vector in the $E_l$ basis. Now let
$$B_{\g_A} = (1-e^{-iT_A}) T_A^{-1}.$$

This matrix can be computed classically in $\poly(d_{\g_A})=\poly(d)$ time.  To compute the gradient from this vector representation, one first needs to estimate all $R_{l}:=i\tr{ O \qty[\rho^{out}, E_l]}$ and compute the following

\begin{align}\label{eq:B DLA}
\nabla\Loss = B_{\g_A} \overrightarrow{R} .
\end{align}

With that the proof is complete.

The following algorithm summarizes this procedure to estimate the gradient.
DLA Gradient Estimation
Input: DLA generators $E_1, \cdots, E_{d_{\g_A}}$,

[1]

Estimate the expectation value of the tests $R_l:= i\tr{O[E_l, \rho^{out}]}$ for $l\in [d_{\g_A}]$.
Compute the $d_\g\times d_{\g_A}$ matrix $B_{\g_A}$ as in (<ref>).
Return $\nabla\Loss \approx B_{\g_A} \overrightarrow{R}$.

§ APPROXIMATIONS FOR GENERAL CASE

When $\CS\subset \CK_4^d$ is not a subgroup, then the terms appearing in the partial derivative in  Theorem <ref> expand beyond $\CS$. Indeed, they form a subgroup denoted by $\<\CS\>$ called the subgroup generated by $\CS$ and is defined as the smallest subgroup containing $\CS$. In that case, one can compute the gradient from Theorem <ref> and Algorithm <ref> for $\<S\>$. This yields a run-time that scales with $|\<\CS\>|^3$. This is tractable as long as the size of $\<\CS\>$ is polynomial in $d$. Otherwise, one needs a different approach. In what follows, we introduce approximations to the gradient computation for general $\CS$.

§.§ Truncation

In what follows, we propose to approximate the partial derivatives by truncating the infinite sum that appeared as (<ref>) in Theorem <ref>.

Consider an ansatz of the form $U(\overrightarrow{a})=\exp{i A(\parameter)}$, where $A(\parameter) = \sum_{\bfs \in \mathcal{S}} \qas \qps$, and $\mathcal{S}\subseteq\CK_4^d$ is the index of the parameters.   Then, the derivative of $\Loss(\parameter)$ for this ansatz is approximated in terms of  $D_\bft$ as:

\begin{align*}
\frac{\partial \Loss}{\partial a_{\bfr}} = \sum_{k=0}^K \frac{(-2)^k}{(k+1)!} \sum_{\bfs_1 \in \mathcal{S}} \cdots \sum_{\bfs_k \in \mathcal{S}}  \prod_{j=1}^k a_{\bfs_j} \qty( \bfs_j \circledast (\bfs_1 \circ \cdots \circ \bfs_{j-1}\circ \bfr))   D_{\bfs_1 \circ \cdots \circ \bfs_k \circ \bfr} + \epsilon_K,
\end{align*}

where $\epsilon_K = \exp{-K\log \frac{K}{ e \pi p \norm{O}}}$.

In the setting of the above theorem, with bounded $\norm{O}$, setting $K=\Omega(p+\log \frac{1}{\epsilon})$ suffices to approximate $\pdv{\Loss}{\qas}$ upto an $\epsilon$ additive error.

The proof of the theorem follows from the fact that

§.§ Unbiased Estimation via Randomization

The infinite sum can be estimated via a randomization technique. We show that the partial derivative in Theorem <ref> can be written as an expectation value of a function of the Poisson random variable. Consider the Poisson probability mass function with rate $\lambda=2$:

\begin{align*}
\prob{X=k}= \frac{2^k}{k!}e^{-2}, \qquad k=0,1,\cdots
\end{align*}

By sampling $k$ from this distribution we can estimate the partial derivative expression of Theorem <ref>. Note that

\begin{align}\label{eq:pdv infinite sum randomized}
\frac{\partial \Loss}{\partial a_{\bfr}} = \sum_{k=0}^\infty \frac{(2)^{k+1}}{(k+1)!} X(k)
\end{align}

\begin{align*}
X(k):=\frac{(-1)^k}{2} \sum_{\bfs_1 \in \mathcal{S}} \cdots \sum_{\bfs_k \in \mathcal{S}}  \prod_{j=1}^k \qty( a_{\bfs_j}  \qty( \bfs_j \circledast (\bfs_1 \circ \cdots \circ \bfs_{j-1}\circ \bfr)))    D_{\bfs_1 \circ \cdots \circ \bfs_k \circ \bfr}.
\end{align*}

In that case,
$$\frac{\partial \Loss}{\partial a_{\bfr}} =e^2\EE[X(K)].$$
Hence, by sampling from this distribution we can compute an unbiased estimate of the partial derivative.

§.§.§ Short-term Hamiltonian

The above procedure can be made significantly more efficient if one guarantees that the parameter range is small, that is $|\qas|\leq \epsilon$ for small enough $\epsilon>0$. This is the case in short-term Hamiltonian simulation or a Hamiltonian with “weak interactions". In short-term simulation, the unitary is of the form $U=e^{i\Delta t H}$, where $\Delta t$ is the time step. In that case, each partial derivative is approximated in $O(1)$ classical computation with one Hadamard test of Figure <ref>.

If $\max_{\bfs\in \CS}|\qas|\leq \frac{\epsilon}{p (1+2\epsilon)}$, then the partial derivative can be approximated with $\epsilon$ additive error via $\pdv{\Loss}{\qas} \approx D_{\bfs}$.

§ SAMPLE COMPLEXITY OF GRADIENT ESTIMATION

Next, we study more sophisticated ways to estimate the gradient in the setting of Theorem <ref>.

§.§ Measurement with Shadow Tomography

Turning the gradient estimation to a series of Hadamard tests has another benefit that can further reduce the number of shots to $O(\log p)$. This can be done using the classical shadow tomography [17] — a procedure to estimate several observables with minimal sample complexity.  Suppose the observable $O$ in (<ref>) is $k$-local, meaning that it decomposes into a finite sum of observables acting non-trivially on at most $k$ qubits. In that case, we can prove that the number of shots scales logarithmically with $p$.

We note that in general the above approach does not extend to existing works even when $O$ is $k$-local. The reason is that the corresponding observables are not themselves $k$-local.

Shadow tomography with Pauli measurements. In shadow tomography [17], a random unitary $V$ is applied to the  input state followed by a measurement in the computational basis.  Then the measurement outcomes are used to generate a classical matrix called the shadow. Then the shadows are used classically to compute the expectation value of the observables of interest.

Different choices of $V$ have been studied in [17]. Here, we describe this process when $V$ is the tensor product of randomly chosen Pauli operators.  In other words,

\begin{align*}
V = V_1\tensor \cdots \tensor V_d \in CL(2)^{\tensor d},
\end{align*}

where each $V_j$ is chosen randomly and uniformly from the Clifford group $CL(2)$. Given the input state $\rho$, the output after this rotation is given by $V\rho V^{\dagger}$. Measuring this state give a string of bits $\hat{b}_1,\cdots, \hat{b}_d\in \{0,1\}$ as the outcome. This process is equivalent to measuring each qubit of $\rho$ randomly on the $X,Y$ or $Z$ basis.  Given, these bits and with the choice of $V_1,\cdots, V_d$, the shadow matrix is computed as

\begin{align*}
\hat{\rho}:=\bigotimes_{j=1}^d \qty(3V_j^\dagger \ketbra{\hat{b}_j}V_j - I).
\end{align*}

Repeating this procedure for $n$ copies of $\rho$ gives a set of $n$ shadows denoted by $\hat{\rho}_1, \cdots, \hat{\rho}_n$.  Suppose one is interested to estimate the expectation value of the observables $M_1,\cdots, M_m$ with respect to $\rho$.  Then, for each $M_j$, one computes $\tr{M_j\hat{\rho}_i}, i\in [n]$ followed by a meadian of means estimator to estimate $\<M_j\>_{\rho}$.

Suppose the observables $M_j , j\in [m]$ act non-trivially on at most $k$ qubits. Then, the shadow tomography with random Pauli measurements estimates $\<M_j\>_\rho$ for all $j\in [m]$ up to an additive error $\epsilon>0$ provided that
$$n=O\qty(\frac{4^k}{\epsilon^2}\log m \max_j \norm{M_j}^2_\infty),$$

copies of $\rho$.  The algorithm runs in $O\qty(2^{\Theta(k)}\log m)$ classical time.

Gradient Estimation With Shadow Tomography
Input: $n, k$, Parameter subgroup $\CS$,

[1]

$i=1$ to $n$
Construct a random Pauli measurement $V_j$ for the non-trivial qubits of $O$.
Apply the ansatz to generate the output state and measure the nontrivial qubits in the selected basis.
Let $\hat{b}_1,\cdots \hat{b}_k$ be the measurement outcomes. Compute the classical shadow $\hat{\rho}_i$ from the $k$ matrices $3V_j^\dagger \ketbra{\hat{b}_j}V_j - I$.
Partition $\hat{\rho}_i$'s into $m$ groups and let $\tilde{\rho}_l$ be the empirical average of the $l$th group.

For each $\bfs \in \CS$ estimate $D_{\bfs}$ by median of measns for $ \hat{D}_{\bfs, l} = i \tr{O[\qps, \tilde{\rho}_l]}, l\in [m]$.
Compute the matrix $V$ where the column $i$ is computed as in (<ref>).
Compute the matrix $B = (I-e^{-V})V^{-1},$ where $V^{-1}$ is the generalized inverse of $V$.

Return $\nabla\Loss \approx B \overrightarrow{D}_\bfs$.

We will use this procedure to estimate the Hadamard tests with reduced sample complexity. All the steps are summarized as Algorithm <ref>.  Based on this theorem and Theorem <ref>, we have the following result for the case when $O$ in the VQA formulation is $k$-local.

In the setting of Theorem <ref>, suppose the observable $O$ is a finite sum $O=O_1+\cdots +O_m$ of observables that act on at most $k$-qubits. Then there exists an algorithm that estimates $\nabla \Loss(\parameter)$ with $O(\frac{1}{\epsilon^2} 4^k \log (m p) ~\max_j \norm{O_j}_\infty)$ ansatz uses and an additional

$$O\qty(\frac{1}{\epsilon^2}2^{\Theta(k)}pm\log (m p) \max_j \norm{O_j}^2_\infty + p^3+pd)$$

classical time.

We consider a more restricted assumption first. Suppose the observable $O$ in the VQA formulation acts on at most $k$ qubits. Note that being $k$-loval means that $O$ is a finite sum of such observables.  Since $O$ acts on at most $k$-qubits then $O = \tilde{O}\tensor I_{d-k}$ for some $k$-qubit observable $\tilde{O}$.   Then, the Hadamard tests appearing in the partial derivative formulation of Theorem <ref> have locality equal to $k$.  To see this, recall that

\begin{align*}
D_\bfs = i\tr{O[\qps, \rho^{out}]}.
\end{align*}

Moreover, note from [28] the following property of the commutator for any operator $B$:

\begin{align*}
[\qps, B] = i\qty(R_\bfs(\frac{\pi}{2}) B R_\bfs(\frac{\pi}{2})^\dagger - R_\bfs(-\frac{\pi}{2}) B R_\bfs(-\frac{\pi}{2})^\dagger),
\end{align*}

where $R_\bfs(\theta):=e^{-i\frac{\theta}{2} \qps}$. With this observation, we obtain

\begin{align*}
D_\bfs &= -\tr{O \qty(R_\bfs(\frac{\pi}{2}) \rho^{out} R_\bfs(\frac{\pi}{2})^\dagger - R_\bfs(-\frac{\pi}{2}) \rho^{out} R_\bfs(-\frac{\pi}{2})^\dagger)}\\
&=-\tr{R_\bfs(\frac{\pi}{2})^\dagger O R_\bfs(\frac{\pi}{2}) \rho^{out}} + \tr{R_\bfs(-\frac{\pi}{2})^\dagger O R_\bfs(-\frac{\pi}{2}) \rho^{out} },
\end{align*}

where we used the cyclic property of the trace. Without loss of generality assume $\tilde{O}$ acts on the first $k$ qubits. Note that since $R_\bfs$ is in a tensor product form, then we have

\begin{align*}
R_\bfs(\frac{\pi}{2})^\dagger O R_\bfs(\frac{\pi}{2}) &= \qty(\qty\Big(\bigotimes_{j=1}^k R_{s_j}(\frac{\pi}{2})^\dagger) \tilde{O} \qty\Big(\bigotimes_{j=1}^k R_{s_j}(\frac{\pi}{2}))) \bigotimes \qty(\qty\Big(\bigotimes_{j>k}^d R_{s_j}(\frac{\pi}{2})^\dagger) I \qty\Big(\bigotimes_{j>k}^d R_{s_j}(\frac{\pi}{2})))\\
& = \qty\Big(\bigotimes_{j=1}^k R_{s_j}(\frac{\pi}{2})^\dagger) \tilde{O} \qty\Big(\bigotimes_{j=1}^k R_{s_j}(\frac{\pi}{2}))\bigotimes I_{d-k}
\end{align*}

This implies that $R_\bfs(\frac{\pi}{2})^\dagger O R_\bfs(\frac{\pi}{2})$ acts non trivially on $k$ qubits. The same holds for $R_\bfs(-\frac{\pi}{2})^\dagger O R_\bfs(-\frac{\pi}{2})$. Hence, $D_\bfs$ can be written as the expectation value of an observable that acts on $k$ qubits of $\rho^{out}$. As a result, one can use the classical shadow tomography with Pauli measurements and Theorem <ref> to obtain estimates of $D_\bfs$ for all $\bfs\in \CS$. Therefore, we need $O\qty(\frac{4^k}{\epsilon^2}\log p \norm{O}^2_\infty)$ copies of $\rho^{out}$. This estimation takes
$$O\qty(\frac{1}{\epsilon^2}2^{\Theta(k)}\log p \norm{O}^2_\infty)$$
classical time. Lastly, according to Theorem <ref>,  we need an additional $O(p^3 + pd)$ classical time to estimate the gradient from the estimated $D_\bfs, \bfs\in \CS$.

The extension to $k$-local follows from the fact that $O$ is a finite sum $O=O_1+\cdots +O_m$ of observables that act on at most $k$-qubits. In this case $D_\bfs$ can be written as a finite sum $D_\bfs = D_{\bfs, 1} + \cdots + D_{\bfs, m}$. We can use the above procedure to estimate each component. Hence, in total, we will have $m p$ observables acting on at most $k$ qubits. This adds up logarithmically in $m$ to the number of shots and linearly to the time.

§.§ Joint Measurability

One benefit of estimating the gradient via Hadamard tests is that one can potentially group all the required Hadamard tests into mutually compatible collections. Joint measurability refers to the possibility of measuring several observables via one reference measurement.

A set of observables $M_j=\{\Lambda^j_u: u\in \mathcal{U}\}, j=1,2,...,k$ are called stochastically compatible if there exist a measurement $M_{ref}=\set{A_w: w\in \mathcal{W}}$ and  stochastic kernel functions $q_j(w| u)$ on $\mathcal{W}\times \mathcal{U}$ such that
$$\Lambda^j_u = \sum_{w\in \mathcal{W}} q_j(w| u) A_w$$
for all $u\in \mathcal{U}$ and all $j\in [k]$.

Suppose $\CF_1,\cdots, \CF_m$ are groups of the Hadamard tests corresponding to $D_\bfs, \bfs\in \CS$ that are mutually compatible. In that case, for each group, we first perform a fine-grained measurement $M_{j}$ on the sample and then extract the estimations from additional post-measurement classical processes. This approach helps reduce the sample complexity.

Suppose $\mathcal{F}_1,\cdots, \mathcal{F}_m$ form the mutually compatible groups of the Pauli strings appearing in a parameterized Hamiltonian $A(\parameter)=\sum_i a_i P_i$. Then, the number of joint Hadamard tests in Theorem <ref> can be reduced to $\tilde{O}(\frac{m}{\epsilon^2})$.

Jointly measuring Pauli strings has been studied extensively in the literature [22, 39, 8, 2].

§ CONCLUDING REMARKS

This paper provides a framework to estimate the gradient of generic PQCs via Hadamard tests for Pauli operators followed by classical post-processing. It is shown that the proposed approach is polynomial in classical and quantum resources when the DLA of the associated Hamiltonian of the PQC has a dimensionality polynomial in the number of qubits. Moreover, this method does not change the ansatz structure and can be used to reduce the measurement shot complexity to scale logarithmically with the number of parameters. The results would be beneficial in various optimization or learning quantum algorithms that rely on the estimation of the gradient.

As a future work, one might be interested in extending the proposed framework to the estimation of the higher order derivatives and measures such as the Hessian or the Fubini–Study metric defined for a parameterized quantum system. For that, one needs to define a second-order Hadamard test strategy for measuring the double derivatives. Another future work is to extend this strategy to multi-layered PQCs where each layer might be a black box unitary. Finally, deriving lower bounds on the classical and quantum resources needed to estimate the gradient or higher-order derivatives is another important direction.

§ BASIC DEFINITIONS

§.§ Group

A group $G$ is a set equipped with a binary operation (usually denoted by $*$) that satisfies the following properties:

* Closure: For any two elements $a, b \in G$, the result of the operation $a * b$ is also an element of $G$.

* Associativity: For all $a, b, c \in G$, $(a * b) * c = a * (b * c)$.

* Identity Element: There exists an element $e \in G$ such that for all $a \in G$, $a * e = e * a = a$.

* Inverse Element: For every element $a \in G$, there exists an element $a^{-1} \in G$ such that $a * a^{-1} = a^{-1} * a = e$, where $e$ is the identity element.

A subgroup $H$ of a group $G$ is a subset of $G$ that forms a group under the same binary operation as $G$. Formally, $H$ is a subgroup of $G$ if it satisfies the following properties:

* Closure: For any two elements $h_1, h_2 \in H$, the result of the operation $h_1 * h_2$ is also an element of $H$.

* Identity Element: The identity element of $G$ is also in $H$.

* Inverse Element: For every element $h \in H$, its inverse $h^{-1}$ is also in $H$.

§.§ Lie Algebra

A Lie algebra is a vector space $\mathfrak{g}$ equipped with a binary operation called the Lie bracket, denoted by $[\cdot, \cdot]$, which satisfies the following properties for all $X, Y, Z \in \mathfrak{g}$:

* Bilinearity: $[aX + bY, Z] = a[X, Z] + b[Y, Z]$ and $[X, aY + bZ] = a[X, Y] + b[X, Z]$ for all $a, b \in \mathbb{R}$.

* Antisymmetry: $[X, Y] = -[Y, X]$.

* Jacobi Identity: $[X, [Y, Z]] + [Y, [Z, X]] + [Z, [X, Y]] = 0$.

A Lie subalgebra $\mathfrak{s}$ of a Lie algebra $\mathfrak{g}$ is a subspace of $\mathfrak{g}$ that is itself a Lie algebra under the same Lie bracket operation. Formally, $\mathfrak{s}$ is closed under the Lie bracket, meaning that for all $X, Y \in \mathfrak{s}$, $[X, Y] \in \mathfrak{s}$.
For example, in a $d$-qubit quantum system, the subspace

\[
\emph{Span}\{ X^{\tensor d},  Y^{\tensor d}, Z^{\tensor d} \},
\]

is closed under the Lie bracket. Moreover, it has dimensionality equal to $3$ inside the $4^d$ dimensional Lie algebra.

§ ANALYTICAL DERIVATION OF THE GRADIENT IN EXAMPLE <REF>

Note that from (<ref>) the partial derivative of the objective function can be written as

\begin{align*}
\frac{\partial \Loss}{\partial \qas} = \tr{ O \qty(\frac{\partial U}{\partial \qas} \qty(\rho U^\dagger) + \qty(U\rho) \pdv{U^\dagger}{\qas})},
\end{align*}

Given that $\pdv{U^\dagger}{\qas}=(\pdv{U}{\qas})^\dagger$. Then, by denoting $\tilde{U}=\pdv{U}{\qas}$ we have that

\begin{align*}
\pdv{\Loss}{\qas} &= \tr{ O \tilde{U}\rho U^\dagger + U\rho\tilde{U}^\dagger}.
\end{align*}

Next, as $UU^\dagger = I$, we have

\begin{align*}
\pdv{\Loss}{\qas} &= \Tr{O\left( \tilde{U} U^\dagger (U \rho U^\dagger) + (U \rho U^\dagger) U \tilde{U}^{\dagger}\right)} \\
&= \Tr{O\left( \tilde{U} U^\dagger \rho^{out} + \rho^{out} (\tilde{U} U^\dagger)^\dagger\right)},
\end{align*}

where $\rho^{out}$ is the ansatz output.
Note that the single qubit ansatz  can also be written as

\begin{equation} \label{eq:single qubit unitary}
U(\overrightarrow{a}) = I\cos{\theta} + i\left( \sum \hat{a}_s \sigma^s\right) \sin{\theta},
\end{equation}

where $\theta = \sqrt{\sum a_s^2}$ is a normalizing parameter and $\hat{a}_s = \frac{a_s}{\theta}$. Now, we can differentiate $U$ with respect to a single parameter $a_s$ appearing in the sum:

\begin{equation} \label{eq:deriv single unitary}
\begin{split}
\frac{\partial U(\overrightarrow{a})}{\partial a_s} = (-\frac{a_s}{\theta}\sin{\theta}) I &+ i\left( \sum_{s^\prime \neq s}\frac{-a_sa_{s^\prime}}{\theta^3}\sigma^{s^\prime}  + \frac{\theta^2 - a_s^2}{\theta^3}\sigma^s\right)\sin{\theta} \\ &+ i \left( \frac{a_s\cos{\theta}}{\theta}\sum_{s^\prime} \hat{a}_{s^\prime} \sigma^{s^\prime}\right).
\end{split}
\end{equation}

Using (<ref>) to find $U^\dagger$ and (<ref>) for $\tilde{U}$, we can analytically find the derivative of the loss function at each point. To better visualize how we converge to the true derivative, consider the example where

\begin{equation} \label{eq:toy example}
U(\overrightarrow{a}) = \exp{i (a_1\sigma^1 + a_2\sigma^2)}
\end{equation}

We will analytically find the derivative at $a_1 = 0$. The equations yield

\begin{equation*}
\tilde{U} U^\dagger = i \frac{\sin{a_2}}{a_2} \left( \cos{a_2} \sigma^1 + \sin{a_2} \sigma^3\right)
\end{equation*}

Which when plugged into the derivative expression gives

\begin{equation*}
\left.\frac{\partial L}{\partial a_1}\right|_{a_1=0} = \frac{\sin{a_2}}{a_2} \left( \cos{a_2} D_1 + \sin{a_2} D_3 \right).
\end{equation*}

Now, we will derive the same expression using Theorem <ref>. Looking at individual terms in the summation, we have the following:

\begin{equation*}
\left.\frac{\partial L}{\partial a_1}\right|_{a_1=0} = D_1 + a_2 D_3 - \frac{4}{6}a_2^2 D_1 - \frac{8}{24}a_2^3D_3 + \dots
\end{equation*}

Note the repeating nature of the $D_i$s. We then rearrange the terms:

\begin{align*}
\left.\frac{\partial L}{\partial a_1}\right|_{a_1=0} &= D_1 \left( 1 - \frac{4}{6}a_2^2 + \dots \right) + D_3 \left( a_2 - \frac{8}{24}a_2^3 + \dots\right) \\
&= D_1 \frac{1}{a_2} \left( \sum_p \frac{(-1)^p}{(2p+1)!}(2a_2)^{2p+1} \right) + D_3 \frac{1}{a_2} \left( \sum_p \frac{(-1)^p}{(2p+2)!}(2a_2)^{2p+2} \right) \\
&= \frac{1}{2a_2} \left( D_1 \sin{2a_2} + D_3 (1 - \cos{2a_2} ) \right) \\
&= \frac{\sin{a_2}}{a_2} \left( \cos{a_2} D_1 + \sin{a_2} D_3 \right).
\end{align*}

Which is exactly the expression we got analytically. In this example, we can see how our method is akin to approximating a function with a series expansion. The power of this approach lies in the group structure of the Pauli matrices, which produces predictable patterns. Thus, we are dealing with a repeating structure that makes use of a subset of $D_i$s. Another point is that $D_i$ only makes use of $\rho_{out}$. We can calculate these $D_i$ and approximate the derivative to any arbitrary precision.

§ SUMMARY OF RELATED WORKS

In this approach the non-product unitary is approximated via the Suzuki-Trotter transformation [35] which states that for any operators $A_1, A_2, \cdots, A_k$, that do not necessarily commute with each other, the following holds

\begin{equation*}
\lim_{n\rightarrow \infty} \Big(\prod_{j=1}^k \exp{A_j/n }\Big)^n = \exp\big\{\sum_{j=1}^k A_j\big\}.
\end{equation*}

The Trotter formula has been used in literature to derive approximations for generic PQCs [44, 25, 29]. However, implementing the above approximation may lead to a high gate complexity and is not preferable in scenarios where direct implementation is available.

Stochastic PSR:
Stochastic parameter shift rule was proposed in [3] to measure the derivative for unitary operators of the form $e^{i (\theta \qps + B)}$, where $\qps$ is a Pauli string and $B$ is an arbitrary Hermitian operator. The authors showed that the derivative is equal to

\begin{align*}
\pdv{\Loss}{\theta}  = \int_0^1 C_+(\theta, t) - C_-(\theta, t) d t,
\end{align*}

where $C_\pm(\theta, t) : = \tr{O V_\pm(\theta, t) \rho V_\pm^\dagger(\theta, t)}$, with

\begin{align*}
V_\pm(\theta, t) :=  e^{i t (\theta \qps + B)}e^{\pm i \qps} e^{i(1-t)(\theta \qps + B)}.
\end{align*}

Then the authors present a Monte Carlo strategy to estimate this integral. A more general variant of the parameter shift rule has been introduced [41]. The considered general ansatz of the form $U(\theta) = e^{i (\theta A + B)}$ for generic Hermitian $A$ and $B$. Given the spectral decomposition of $A$ the paper provides an explicit formula for the derivative of the objective function. This is done via a Discrete Fourier series approach.

Such methods for generic ansatz $e^{i (\theta A + B)}$ require not only to perturb the parameters but also to change the unitaries involved. As a result,
to modify not just the parameters of the PQC, but also change the unitaries that appear. In practice, such modifications require a re-evaluation of the schedule of the underlying quantum-control system and hence are at a disadvantage. Moreover, the stochastic PSR has a high estimation variance. This is because the above integral is estimated by sampling values of $s$ uniformly in the interval $(0, 1)$ and then
calculating the costs with a finite-shot estimate. In addition, this method leads to a bigger number of unique circuits to compute the derivative, increasing the compilation overhead for both hardware and simulator implementations.

Nyquist PSR: Recently  [37] proposed a “proper" shift rule for PQC of the form $e^{i (\theta A + B)}$ where only the parameters are shifted without any other modifications of the ansatz. The method was called Nyquist parameter shift rule and relies on a beautiful connection between the Nyquist-Shannon Sampling theorem and the Fourier series that was observed earlier in [41, 38]. This paper shows that if $f(x) = tr\{O U \rho U^\dagger \}$ with $U=e^{ix H + B}$ and $K$ being the difference between the maximum and minimum eigenvalues of $A$,  then the Fourier spectrum of $f$ is contained in $[-K, K]$. Hence, the Nyquist-Shannon Sampling theorem can be used to estimate the derivative of the objective function.  The proposed method, however, has a low approximation error when the parameter value is large enough. More precisely, the approximation error is $O(\frac{1}{c^2})$ as long as $\theta = (1 - \Omega(1))c$, where $c$ is the maximum magnitude of a parameter value. Our method is suitable when the parameter value is small.

[1]
Gian-Luca R Anselmetti, David Wierichs, Christian Gogolin, and Robert M
Local, expressive, quantum-number-preserving vqe ansätze for
fermionic systems.
New Journal of Physics, 23(11):113010, November 2021.

[2]
Somshubhro Bandyopadhyay, P. Oscar Boykin, Vwani Roychowdhury, and Farrokh
A new proof for the existence of mutually unbiased bases.
quant-ph/0103162, 2001.

[3]
Leonardo Banchi and Gavin E. Crooks.
Measuring analytic gradients of general quantum evolution with the
stochastic parameter shift rule.
Quantum, 5:386, January 2021.

[4]
Marcello Benedetti, Erika Lloyd, Stefan Sack, and Mattia Fiorentini.
Parameterized quantum circuits as machine learning models.
Quantum Science and Technology, 4(4):043001, 2019.

[5]
M. Cerezo, Andrew Arrasmith, Ryan Babbush, Simon C. Benjamin, Suguru Endo,
Keisuke Fujii, Jarrod R. McClean, Kosuke Mitarai, Xiao Yuan, Lukasz Cincio,
and Patrick J. Coles.
Variational quantum algorithms.
Nature Reviews Physics, 3(9):625–644, aug 2021.

[6]
A. R. Calderbank, E. M Rains, P. W. Shor, and N. J. A. Sloane.
Quantum error correction via codes over gf(4).
August 1996.

[7]
M. Cerezo, Akira Sone, Tyler Volkoff, Lukasz Cincio, and Patrick J. Coles.
Cost function dependent barren plateaus in shallow parametrized
quantum circuits.
Nature Communications, 12(1), March 2021.

[8]
Ophelia Crawford, Barnaby van Straaten, Daochen Wang, Thomas Parks, Earl
Campbell, and Stephen Brierley.
Efficient quantum measurement of pauli operators in the presence of
finite sampling error.
arXiv:1908.06942, 2020.

[9]
Alain Delgado, Juan Miguel Arrazola, Soran Jahangiri, Zeyue Niu, Josh Izaac,
Chase Roberts, and Nathan Killoran.
Variational quantum algorithm for molecular geometry optimization.
Physical Review A, 104(5):052402, November 2021.

[10]
Edward Farhi, Jeffrey Goldstone, and Sam Gutmann.
A quantum approximate optimization algorithm, 2014.

[11]
Enrico Fontana, Dylan Herman, Shouvanik Chakrabarti, Niraj Kumar, Romina
Yalovetzky, Jamie Heredge, Shree Hari Sureshbabu, and Marco Pistoia.
The adjoint is all you need: Characterizing barren plateaus in
quantum ansätze.
September 2023.

[12]
Edward Farhi and Hartmut Neven.
Classification with quantum neural networks on near term processors.
February 2018.

[13]
Harper R. Grimsley, Sophia E. Economou, Edwin Barnes, and Nicholas J. Mayhall.
An adaptive variational algorithm for exact molecular simulations on
a quantum computer.
Nature Communications, 10(1), July 2019.

[14]
Daniel Gottesman.
Stabilizer codes and quantum error correction.
California Institute of Technology, 1997.

[15]
Vojtěch Havlíček, Antonio D. Córcoles, Kristan Temme, Aram W. Harrow,
Abhinav Kandala, Jerry M. Chow, and Jay M. Gambetta.
Supervised learning with quantum-enhanced feature spaces.
Nature, 567(7747):209–212, March 2019.

[16]
Mohsen Heidari, Ananth Y. Grama, and Wojciech Szpankowski.
Toward physically realizable quantum neural networks.
Association for the Advancement of Articial Intelligence
(AAAI), 2022.

[17]
Hsin-Yuan Huang, Richard Kueng, and John Preskill.
Predicting many properties of a quantum system from very few
Nature Physics 16, 1050–1057 (2020), February 2020.

[18]
Aram W. Harrow and John C. Napp.
Low-depth gradient measurements can improve convergence in
variational hybrid quantum-classical algorithms.
Physical Review Letters, 126(14):140502, apr 2021.

[19]
Stuart Hadfield, Zhihui Wang, Bryan O’Gorman, Eleanor Rieffel, Davide
Venturelli, and Rupak Biswas.
From the quantum approximate optimization algorithm to a quantum
alternating operator ansatz.
Algorithms, 12(2):34, February 2019.

[20]
Tyson Jones, Suguru Endo, Sam McArdle, Xiao Yuan, and Simon C. Benjamin.
Variational quantum algorithms for discovering hamiltonian spectra.
Physical Review A, 99(6):062304, June 2019.

[21]
Abhinav Kandala, Antonio Mezzacapo, Kristan Temme, Maika Takita, Markus Brink,
Jerry M. Chow, and Jay M. Gambetta.
Hardware-efficient variational quantum eigensolver for small
molecules and quantum magnets.
Nature, 549(7671):242–246, sep 2017.

[22]
Jay Lawrence, Časlav Brukner, and Anton Zeilinger.
Mutually unbiased binary observable sets onNqubits.
Physical Review A, 65(3), feb 2002.

[23]
Martin Larocca, Piotr Czarnik, Kunal Sharma, Gopikrishnan Muraleedharan,
Patrick J. Coles, and M. Cerezo.
Diagnosing barren plateaus with tools from quantum optimal control.
Quantum, 6:824, September 2022.

[24]
Jin-Guo Liu and Lei Wang.
Differentiable learning of quantum circuit born machines.
Physical Review A, 98(6):062324, December 2018.

[25]
Xuanqing Liu, Tesi Xiao, Si Si, Qin Cao, Sanjiv Kumar, and Cho-Jui Hsieh.
How does noise help robustness? explanation and exploration under the
neural sde framework.
In 2020 IEEE/CVF Conference on Computer Vision and Pattern
Recognition (CVPR), pages 279–287, 2020.

[26]
Jarrod R. McClean, Sergio Boixo, Vadim N. Smelyanskiy, Ryan Babbush, and
Hartmut Neven.
Barren plateaus in quantum neural network training landscapes.
Nature Communications, 9(1), November 2018.

[27]
Kosuke Mitarai, Masahiro Kitagawa, and Keisuke Fujii.
Quantum analog-digital conversion.
Physical Review A, 99(1):012301, January 2019.

[28]
K. Mitarai, M. Negoro, M. Kitagawa, and K. Fujii.
Quantum circuit learning.
Physical Review A, 98(3):032309, sep 2018.

[29]
Alexander Miessen, Pauline J. Ollitrault, and Ivano Tavernelli.
Quantum algorithms for quantum dynamics: A performance study on the
spin-boson model.
Physical Review Research, 3(4):043212, dec 2021.

[30]
Alberto Peruzzo, Jarrod McClean, Peter Shadbolt, Man-Hong Yung, Xiao-Qi Zhou,
Peter J. Love, Alán Aspuru-Guzik, and Jeremy L. O'Brien.
A variational eigenvalue solver on a photonic quantum processor.
Nature Communications, 5(1), jul 2014.

[31]
Wulf Rossmann.
Lie groups: An Introduction through Linear Groups.
Number 5 in Oxford graduate texts in mathematics. Oxford University
Press, Oxford [u.a.], 1. publ. in paperback edition, 2006.
Reprinted 2011.

[32]
Maria Schuld, Ville Bergholm, Christian Gogolin, Josh Izaac, and Nathan
Evaluating analytic gradients on quantum hardware.
Phys. Rev. A 99, 032331 (2019), November 2018.

[33]
Maria Schuld and Nathan Killoran.
Quantum machine learning in feature hilbert spaces.
Physical Review Letters, 122(4):040504, February 2019.

[34]
SG Schirmer, ICH Pullen, and AI Solomon.
Identification of dynamical lie algebras for finite-level quantum
control systems.
Journal of Physics A: Mathematical and General, 35(9):2327,

[35]
Masuo Suzuki.
Generalized trotter's formula and systematic
approximants of exponential operators and inner derivations with applications
to many-body problems.
Communications in Mathematical Physics, 51(2):183–190, jun

[36]
Ryan Sweke, Frederik Wilde, Johannes Meyer, Maria Schuld, Paul K. Faehrmann,
Barthélémy Meynard-Piganeau, and Jens Eisert.
Stochastic gradient descent for hybrid quantum-classical
Quantum, 4:314, aug 2020.

[37]
Dirk Oliver Theis.
“proper” shift rules for derivatives of perturbed-parametric
quantum evolutions.
Quantum, 7:1052, July 2023.

[38]
Javier Gil Vidal and Dirk Oliver Theis.
Calculus on parameterized quantum circuits.
arXiv:1812.06323, December 2018.

[39]
Vladyslav Verteletskyi, Tzu-Ching Yen, and Artur F. Izmaylov.
Measurement optimization in the variational quantum eigensolver using
a minimum clique cover.
The Journal of Chemical Physics, 152(12), March 2020.

[40]
Dave Wecker, Matthew B. Hastings, and Matthias Troyer.
Progress towards practical quantum variational algorithms.
Physical Review A, 92(4):042303, October 2015.

[41]
David Wierichs, Josh Izaac, Cody Wang, and Cedric Yen-Yu Lin.
General parameter-shift rules for quantum gradients.
Quantum, 6:677, March 2022.

[42]
Roeland Wiersema, Efekan Kökcü, Alexander F. Kemper, and Bojko N. Bakalov.
Classification of dynamical lie algebras for translation-invariant
2-local spin systems in one dimension, 2023.

[43]
Roeland Wiersema, Dylan Lewis, David Wierichs, Juan Carrasquilla, and Nathan
Here comes the su(n): multivariate quantum gates and gradients.
Quantum, 8:1275, March 2024.

[44]
Xiaodong Yang, Xinfang Nie, Yunlan Ji, Tao Xin, Dawei Lu, and Jun Li.
Improved quantum computing with higher-order trotter decomposition.
Physical Review A, 106(4):042401, oct 2022.