# Fourier expansion in variational quantum algorithms

Nikita A. Nemkov<EMAIL_ADDRESS>Russian Quantum Center, Skolkovo, Moscow
143026, Russia  National University of Science and Technology “MISIS”, Moscow
119049, Russia  Evgeniy O. Kiktenko  Russian Quantum Center, Skolkovo, Moscow
143026, Russia  National University of Science and Technology “MISIS”, Moscow
119049, Russia  Aleksey K. Fedorov<EMAIL_ADDRESS>Russian Quantum Center,
Skolkovo, Moscow 143026, Russia  National University of Science and Technology
“MISIS”, Moscow 119049, Russia

###### Abstract

The Fourier expansion of the loss function in variational quantum algorithms
(VQA) contains a wealth of information, yet is generally hard to access. We
focus on the class of variational circuits, where constant gates are Clifford
gates and parameterized gates are generated by Pauli operators, which covers
most practical cases while allowing much control thanks to the properties of
stabilizer circuits. We give a classical algorithm that, for an $N$-qubit
circuit and a single Pauli observable, computes coefficients of all
trigonometric monomials up to a degree $m$ in time bounded by
$\mathcal{O}(N2^{m})$. Using the general structure and implementation of the
algorithm we reveal several novel aspects of Fourier expansions in
Clifford+Pauli VQA such as (i) reformulating the problem of computing the
Fourier series as an instance of multivariate boolean quadratic system (ii)
showing that the approximation given by a truncated Fourier expansion can be
quantified by the $L^{2}$ norm and evaluated dynamically (iii) tendency of
Fourier series to be rather sparse and Fourier coefficients to cluster
together (iv) possibility to compute the full Fourier series for circuits of
non-trivial sizes, featuring tens to hundreds of qubits and parametric gates.

###### Contents

  1. I Introduction and results
  2. II Parameterized quantum circuits
    1. II.1 Trigonometric expansion of the unitary
    2. II.2 Fourier expansion of the loss function
    3. II.3 Fourier terms from averages
  3. III Clifford+Pauli variational circuits
    1. III.1 Definition and properties
    2. III.2 Computing averages
  4. IV Classical algorithm
    1. IV.1 Expansion of the dressed Hamiltonian
    2. IV.2 Accounting for the expectation values
    3. IV.3 Truncated Fourier series as an approximation
    4. IV.4 Is there a more efficient algorithm?
  5. V Case studies
    1. V.1 Random circuits
    2. V.2 QAOA
    3. V.3 Hardware-efficient circuits
  6. VI Discussion and outlook
  7. A Structure of a generic Fourier expansion
    1. A.1 Level expansion and number of terms
    2. A.2 Coefficients from averages
  8. B Optimizing Pauli order
  9. C Details of numerical computations
    1. C.1 Random circuits
    2. C.2 QAOA
    3. C.3 Hardware-efficient circuits
    4. C.4 Estimating complexity from Monte-Carlo sampling

##  I Introduction and results

Variational quantum algorithms (VQA) [1] are the leading candidates to make
the most out of current NISQ devices [2, 3]. While the scope of potential VQA
applications is extremely broad, there are also many theoretical and practical
limitations. VQA are hybrid algorithms, using classical optimization to train
parameterized quantum circuits, and in this sense they are similar to the
classical machine learning models. The loss function of a VQA is defined as an
average of some observable in the state, prepared by the parameterized quantum
circuit. The structure of VQA loss landscape is of central importance, because
the efficiency of the classical optimization largely determines the quality of
solutions obtained by VQA. Accessible shapes of the loss function also
determine the expressive power of the quantum machine learning models [4].

Typically, VQA are trained by gradient-based methods and their local
properties are of the most interest. At the same time, the structure of
parametrized quantum circuits makes Fourier series representation a natural
and rich language for description of VQA loss functions. We now briefly survey
some of the relevant results. Under the assumption that generators of
parametric gates have commensurable eigenvalues, the Fourier series in fact
truncates to a trigonometric polynomial. As shown in Ref. [5], accessible
frequencies in this expansion are determined by the spectrum of the
generators, while the coefficients depend on the structure of the circuit and
the observable. In Refs. [5, 6] this observation was used to highlight the
importance of data encoding in quantum machine learning models. Recently,
there has been an interest in quantitatively studying the expressive power of
machine learning models based on the properties of their Fourier expansion [7,
8]. In Refs. [9, 10] the Fourier representation was used to de-quantize a
class of quantum machine learning models. Paper [11] investigates the case
when the Fourier series is rather sparse, so that the loss landscape can be
efficiently recovered with limited experimental data. An interesting proposal
made in Ref. [12] shows that noise in VQA can be detected by observing
inaccessible frequencies and mitigated by filtering them out. We also note
that the fundamental result on the NP-hardness of training general VQA [13]
relies on their loss functions being trigonometric polynomials.

While the Fourier representation can be a very convenient tool to characterize
variational loss functions, it has its limitations. When generators of the
parametric gates square to identity, which is the most common case, the
Fourier series for a VQA with $M$ parameters is a multivariate trigonometric
polynomial containing up to $3^{M}$ terms. Exponential growth of accessible
terms makes the Fourier series an impractical description, unless the actual
distribution of coefficients is very sparse (e.g. the number of non-zero terms
only grows polynomially with $M$). Computing the Fourier coefficients is also
a challenge. For example, evaluating the lowest order, constant Fourier term,
amounts to finding the loss function averaged over all parameter
configurations, and there seems to be no efficient recipe for that in general.

In this paper, we restrict attention to a special class of parametrized
quantum circuits, which we refer to as Clifford+Pauli circuits. Parametric
gates in Clifford+Pauli circuits are exponentials of Pauli strings, while
constant gates are Clifford gates. This is in fact a very large class of
circuits that includes the majority of most studied VQA, such Quantum
Approximate Optimization Algorithm (QAOA) [14], Hardware-Efficient Ansatz
(HEA), [15], Unitary Coupled Cluster Ansatz [16]. Special properties of
stabilizer circuits [17, 18] give an essential technical leverage to study the
Fourier expansion of Clifford+Pauli circuits in quantitative details.

The interplay between properties of stabilizer circuits and VQA have been
explored previously, mainly in the context of quantum chemistry [19, 20, 21].
In particular, initialization methods based on perturbative expansion around
Clifford points [22] or on the discrete search through the space of Clifford
circuits [23, 24] have been developed, ansatz structures [25, 26] and
partitioning schemes [27], based on the properties of Clifford gates have been
proposed. In this work, however, we focus on a different scope of questions.

Our core technical contribution is an efficient classical algorithm computing
all Fourier coefficients in the loss function up to level $m$, with time
complexity bounded by $O(N2^{m})$, where $N$ is the number of qubits. Note
that in general, the complexity of the Fourier series is not limited directly
by the number of qubits $N$, but rather by the total number of parametric
gates $M$. For Clifford+Pauli VQA this observation thus admits a concrete
realization.

The algorithm has both theoretical and practical utility. On the theoretical
side, we show that typical Fourier series are much sparser than anticipated in
the general case. For Clifford+Pauli circuits with a single Pauli observable,
the number of coefficients is upper bounded by $2^{M}$, and for the worst case
behavior expected in practice we find $\left(\frac{3}{2}\right)^{M}$. We also
show, that truncating the Fourier polynomial below the maximum degree gives an
approximation that can be quantified by the $L^{2}$ norm of the loss function
and evaluated dynamically, as the algorithm proceeds. The number of terms
contributing non-negligibly the functional norm is typically an exponential
fraction of all terms, yet still growing exponentially itself.

On the practical side, we perform several case studies to probe the structure
of Fourier expansions in more detail. In all examples we find that Fourier
terms tend to cluster around some mean level, which is exactly $\frac{M}{2}$
for non-local random circuits, but much smaller for the local circuits with
special structure, such as QAOA or HEA, making their Fourier expansion much
sparser and easier to compute.

Our algorithm is based on a simple recursive expansion of the loss function. A
similar approach was described in the context of the Qubit Coupled Cluster
Method [25], and in the context of QAOA in Ref. [28]. However, there are some
important distinctions with the previous work. First, we identify
Clifford+Pauli circuits as the class to which the method is universally
applicable, and treat the problem in general terms, as well as establish its
direct relation to the Fourier series expansion. More importantly, we
demonstrate how to significantly reduce the algorithm cost by pruning some
branches of the recursive expansion early, based on filtering by expectation
values. While this does not change the large $M$ asymptotic, for practical
cases the difference is essential. For example, for random Clifford+Pauli
circuits on $N=50$ qubits, it allows increasing the depth of circuits that can
be handled from $M=30$ to $M=80$ without changing the computational budget.
Finally, we formulate the problem of computing all non-zero coefficients in
the Fourier expansion as an instance of the multivariate boolean quadratic
problem, which allows us to argue that our algorithm is likely not far from
optimal yet points towards potential improvements.

##  II Parameterized quantum circuits

In this section, we establish some notation, describe basic properties of
variational circuits and their loss functions, and discuss how Fourier
expansion arises in this context.

###  II.1 Trigonometric expansion of the unitary

We will assume that a unitary matrix $U(\boldsymbol{\phi})$ of a parameterized
quantum circuit takes the following form:

$\displaystyle U(\boldsymbol{\phi})=C_{M}P_{M}(\phi_{M})\dots
C_{2}P_{2}(\phi_{2})C_{1}P_{1}(\phi_{1})C_{0}\ .$ (1)

Here $C_{m}$ are constant gates, $P_{m}(\phi)=e^{-\imath\frac{\phi}{2}G_{m}}$
are single-parameter rotations, and
$\boldsymbol{\phi}=(\phi_{1},\ldots,\phi_{M})$ is a vector of parameters. We
will assume that Hermitian generators of the parameterized gates square to
identity $G_{m}^{2}=\mathbbm{1}$, so that

$\displaystyle P_{m}(\phi)=\mathbbm{1}\cos\frac{\phi}{2}-\imath
G_{m}\sin\frac{\phi}{2}\ .$ (2)

Applying this relation to each parametric gate in the circuit one obtains the
following formal trigonometric expansion containing $2^{M}$ terms

$\displaystyle
U(\boldsymbol{\phi})=\sum_{I\in\\{0,1\\}^{M}}U_{I}t_{I}\left(\frac{\boldsymbol{\phi}}{2}\right)\
.$ (3)

Here $I=(I_{1},\ldots,I_{M})$ with $I_{m}\in\\{0,1\\}$ is a multi-index,
$t_{I}(\boldsymbol{\phi})$ is a multivariate trigonometric monomial of order
$M$

$\displaystyle t_{I}(\boldsymbol{\phi})=\prod_{m=1}^{M}t_{I_{m}}(\phi_{m})\ ,$
(4)

where each term in the product is defined by

$\displaystyle
t_{i}(\phi)=\cos^{1-i}\phi\sin^{i}\phi=\begin{cases}\cos\phi,\quad i=0\\\
\sin\phi,\quad i=1\end{cases}\ .$ (5)

We note that coefficient matrices $U_{I}$ correspond to the circuit unitary,
evaluated at specific values

$U_{I}=U(\boldsymbol{\phi}=\pi I)=\alpha C_{M}G_{M}^{I_{M}}\ldots
C_{1}G_{M}^{I_{1}}C_{0}\ ,$ (6)

where $\alpha$ is a phase factor $\alpha=(-\imath)^{\sum_{m}I_{m}}$.

###  II.2 Fourier expansion of the loss function

The loss function $F(\boldsymbol{\phi})$ of a variational algorithm is defined
by the average of some Hermitian operator $H$, often referred to as the
Hamiltonian, in the state prepared by the circuit

$\displaystyle
F(\boldsymbol{\phi})=\Braket{0}{U^{\dagger}(\boldsymbol{\phi})HU(\boldsymbol{\phi})}{0}\
.$ (7)

Here and in the following $\ket{0}=\ket{0}^{\otimes N}$ is the all zeros state
of $N$ qubits. Substituting expansion (3) into the loss function gives

$\displaystyle
F(\boldsymbol{\phi})=\sum_{IJ}t_{I}\left(\frac{\boldsymbol{\phi}}{2}\right)t_{J}\left(\frac{\boldsymbol{\phi}}{2}\right)\Braket{0}{U^{\dagger}_{I}HU_{J}}{0}\
.$ (8)

In contrast to the expansion of unitary (3), which is homogeneous, expansion
of the loss function contains trigonometric monomials of various degrees, see
App. A for details. Let us organize the Fourier expansion of the loss function
by level

$\displaystyle F(\boldsymbol{\phi})=\sum_{m=0}^{M}F_{m}(\boldsymbol{\phi})\ .$
(9)

Each level $F_{m}(\boldsymbol{\phi})$ only involves monomials of order $m$.
There are $\binom{M}{m}$ possible parameter subsets at level $m$ each giving
rise to $2^{m}$ trigonometric monomials. Hence, the total number of
independent coefficients in the Fourier expansion is

$\displaystyle\sum_{m=0}^{M}2^{m}\binom{M}{m}=3^{M}\ .$ (10)

###  II.3 Fourier terms from averages

Computing the Fourier series for generic loss functions appears to be a
formidable task. Indeed, let us note that the constant term $F_{0}$ in the
Fourier expansion can be thought of as the loss function, averaged over all
parameters

$\displaystyle
F_{0}=\Braket{F(\boldsymbol{\phi})}_{\boldsymbol{\phi}}:=\frac{1}{(2\pi)^{M}}\int_{0}^{2\pi}\prod_{m=1}^{M}d\phi_{m}\,\,F(\boldsymbol{\phi})\
.$ (11)

This relation holds, because all higher levels $F_{m>0}(\boldsymbol{\phi})$ in
the Fourier series trivially vanish when averaged. Higher level terms can be
obtained similarly, see App. A.

The average in Eq. (11) can be expressed in a succinct form using
orthogonality of trigonometric monomials $t_{I}$ (4)

$\displaystyle\Braket{t_{I}(\boldsymbol{\phi})t_{J}(\boldsymbol{\phi})}_{\boldsymbol{\phi}}=2^{-M}\delta_{IJ}\
.$ (12)

Hence, averaging (8) yields

$\displaystyle
F_{0}=\frac{1}{2^{M}}\sum_{I\in\\{0,1\\}^{M}}\Braket{0}{U_{I}^{\dagger}HU_{I}}{0}\
.$ (13)

Evaluating this expression explicitly seems to be out of reach for generic
circuits. In terms of a classical simulation, computing any single expectation
value in Eq. (13) is difficult on its own for a sufficiently large number of
qubits. Even when the averages can be computed efficiently, either classically
or provided an access to a quantum computer, equation (13) still requires
summing $2^{M}$ terms, infeasible for any significant number of parameters
$M$. As we show in the next section, for Clifford+Pauli quantum circuits
evaluating $F_{0}$ and in fact any particular monomial in the Fourier
expansion is classically efficient.

##  III Clifford+Pauli variational circuits

###  III.1 Definition and properties

Let us first establish some notation relevant for stabilizer circuits. A
single-qubit Pauli operator is simply an $X,Y,Z$ Pauli matrix or an identity,
possibly with a phase $\pm 1,\pm\imath$. An $n$-qubit Pauli operator is a
tensor product of $n$ arbitrary single-qubit Pauli operators. Any two Pauli
operators $P_{1}$ and $P_{2}$ either commute or anti-commute: $P_{1}P_{2}=\pm
P_{2}P_{1}$. Clifford gates $C$ are operators that transform every Pauli gate
$P$ into some Pauli gate $P^{\prime}$: $C^{\dagger}PC=P^{\prime}$. The group
of Clifford gates can be generated by the Hadamard gate H , $\textsf{S
}=\sqrt{Z}$, and controlled NOT gate. Circuits consisting only of the Clifford
gates applied to the stabilizer states (of which $|0\rangle$ is an example)
can be efficiently simulated classically due to the Gottesman-Knill theorem
[18].

We define Clifford+Pauli variational circuits as a subset of variational
circuits (1), where generators of parametric gates are Pauli operators and all
constant gates are Clifford gates. For clarity of exposition, we also assume
that the Hamiltonian $H$ is a Pauli operator. The case that is most relevant
in practice, when the Hamiltonian is a polynomial-sized sum of Pauli
operators, can be handled by linearity. We stress that both the Pauli
generators and the Hamiltonian are allowed to have arbitrary weight, i.e. be
supported on any number of qubits. Note that Pauli rotations with generic
angles are not Clifford gates, and hence Clifford+Pauli circuits can not be
simulated efficiently using the stabilizer formalism.

Clifford+Pauli circuits admit a simple canonical form, where all Clifford
gates are eliminated. First, one uses commutation properties of Clifford and
Pauli operators to drag all the Clifford gates to the very end of the circuit.
Generators of Pauli rotations will generally change during the process. The
Clifford gate $C$ accumulated at the end of the circuit is absorbed into the
Hamiltonian $H\to C^{\dagger}HC$, which remains a Pauli operator. Hence,
without loss of generality, we will assume that Clifford+Pauli variational
circuit takes the following Pauli form

$\displaystyle U(\boldsymbol{\phi})=P_{M}(\phi_{M})P_{M-1}(\phi_{M-1})\dots
P_{1}(\phi_{1})\ ,$ (14)

where $P(\phi)=e^{-\imath\frac{\phi}{2}P}$ for some Pauli string $P$. We will
use notation $(P_{1}\dots P_{M}|H)$ for a Clifford+Pauli circuit in Pauli
form. Note that this form is not unique, as Pauli rotations with commuting
generators can be swapped.

###  III.2 Computing averages

Now let us revisit the computation of the average loss function (13). Eq. (13)
now takes the form

$\displaystyle
F_{0}=\frac{1}{2^{M}}\sum_{I\in\\{0,1\\}^{M}}\Braket{0}{P_{1}^{I_{1}}\dots
P_{M}^{I_{M}}HP_{M}^{I_{M}}\dots P_{1}^{I_{1}}}{0}\ .$ (15)

We claim that this sum vanishes unless $H$ commutes with every $P_{i}$.
Indeed, suppose there is some $P_{m}$ s.t. $P_{m}H=-HP_{m}$. Then, it is
straightforward to see, that any two terms in the sum that only differ in the
value of $I_{m}$ give opposite contributions. Hence, for the circuit
$(P_{1}\dots P_{M}|H)$ we find

$\displaystyle
F_{0}=\begin{cases}\braket{0}{H}{0},\quad[H,P_{i}]=0\quad\forall i\\\
0,\qquad\qquad\text{otherwise}\end{cases}\ .$ (16)

Therefore, computing the average loss function for Clifford+Pauli circuits is
a trivial task for any number of qubits and any number of parameters. In fact,
as we show in the following, this applies to every individual term in the
Fourier series. The difficulty of computing the full Fourier expansion then
stems solely from the fact that the total number of non-vanishing coefficients
can be exponentially large. In the next section, we present an efficient
classical algorithm to compute Fourier expansion level by level.

##  IV Classical algorithm

###  IV.1 Expansion of the dressed Hamiltonian

Figure 1: A sample diagram representing recursive expansion of a dressed
Hamiltonian. Here
$\\{P_{3},H\\}=\\{P_{2},H\\}=\\{P_{1},P_{3}H\\}=\\{P_{1},P_{2}H\\}=[P_{2},P_{3}H]=[P_{1},H]=0$.

Introduce the following notation for an operator conjugated by circuit’s
unitary

$\displaystyle
O[\boldsymbol{\phi}]=U^{\dagger}(\boldsymbol{\phi})OU(\boldsymbol{\phi})\ .$
(17)

Following the quantum chemistry literature, we will call
$H(\boldsymbol{\phi})$ the dressed Hamiltonian. The loss function is the
average of the dressed Hamiltonian in the all-zero state

$\displaystyle F(\boldsymbol{\phi})=\Braket{0}{H[\boldsymbol{\phi}]}{0}\ .$
(18)

Next, we make the following simple observation: for an arbitrary Pauli string
$O$ it holds

$\displaystyle P(\phi)^{\dagger}OP(\phi)=\begin{cases}O,\quad&[P,O]=0\\\
O\cos\phi+\imath PO\sin\phi,\quad&\\{P,O\\}=0\end{cases}\ ,$ (19)

i.e. when the conjugating Pauli rotation $P(\phi)$ commutes with $O$, it
cancels out, while for an anti-commuting Pauli rotation the result can be
written as a sum of two Pauli operators. This gives a recurrence procedure to
expand the dressed Hamiltonian. Indeed, it follows from (19) that for an
arbitrary Pauli string $O$

$\displaystyle
O[\boldsymbol{\phi}^{(m)}]=\begin{cases}O[\boldsymbol{\phi}^{(m-1)}],&[O,P_{m}]=0;\\\
\begin{aligned} &O[\boldsymbol{\phi}^{(m-1)}]\cos\phi_{m}+\\\
&~{}~{}+\imath(P_{m}O)[\boldsymbol{\phi}^{(m-1)}]\sin\phi_{m},\end{aligned}&\\{O,P_{m}\\}=0;\end{cases}$
(20)

Here $\boldsymbol{\phi}^{(m)}$ the subset of the first $m\leq M$ parameters
$\boldsymbol{\phi}^{(m)}:=(\phi_{1},\ldots,\phi_{m})$ and
$O[\boldsymbol{\phi}^{(m)}]$ is defined as in (17) with the conjugating
unitary $U(\boldsymbol{\phi}^{(m)}):=P_{m}(\phi_{m})\ldots P_{1}(\phi_{1})$.

Repeatedly applying (20) to $H[\boldsymbol{\phi}]\equiv
H[\boldsymbol{\phi}^{(M)}]$ represents the dressed Hamiltonian as a sum of
Pauli strings multiplied by trigonometric monomials, i.e. gives the Fourier
expansion of the dressed Hamiltonian with operator coefficients. The recurrent
expansion can be conveniently visualized as a binary tree, see Fig. 1 for an
example. The tree is constructed as follows. The nodes, which we refer to as
the computational nodes, correspond to variational circuits, specified by a
list of Pauli generators and an observable $(P_{1}\dots P_{m}|O)$. If the
observable $O$ anti-commutes with the last Pauli generator $P_{m}$, the node
branches into two $(P_{1}\dots P_{m}|O)\to\cos\phi_{m}(P_{1}\dots
P_{m-1}|O)+\imath\sin\phi_{m}(P_{1}\dots P_{m-1}|P_{m}O)$. For brevity, we
omit coefficients at the diagram. Branching increases the Fourier level by
one. If the last Pauli generator instead commutes with $O$, it is simply
removed $(P_{1}\dots P_{m}|O)\to(P_{1}\dots P_{m-1}|O)$ and the Fourier level
remains unchanged. We depict this by horizontal arrows at the diagram. When
there are no Pauli generators left, the node contains the final observable
encoding a single operator coefficient in the Fourier expansion of the dressed
Hamiltonian.

The graphical representation makes several distinctive features of Fourier
series for Clifford+Pauli circuits manifest. Let $n(m)$ be the number of
resulting Fourier modes at level $m$. Introduce

$\displaystyle\delta(m)=2^{-m}n(m),\quad\Delta(m)=\sum_{k=0}^{m}\delta(k)\ .$
(21)

For any Clifford+Pauli circuit and any Pauli Hamiltonian it holds

$\displaystyle\Delta(M)=\sum_{m=0}^{M}2^{-m}n(m)=1\ ,$ (22)

i.e. the weighted sum of populations at all Fourier levels is an invariant.

This implies certain constraints on the distribution of Fourier terms. For
example, the maximum number of Fourier terms $\sum_{m}n(m)$ is upper bounded
by $2^{M}$ (when the last level is fully populated), cf. the bound for generic
circuits $3^{M}$ (10). Importantly, the presence of any single Fourier term at
level $m<M$ reduces the maximum possible amount of terms at other levels. For
instance, if $n(0)=1$, i.e. $F_{0}\neq 0$, all other Fourier terms vanish.

For $N$-qubit circuits, processing each computational node only involves
multiplying Pauli strings of length $N$, and hence has complexity
$\mathcal{O}(N)$. Since the total number of nodes is at most $2^{M}$, the time
complexity of the algorithm can be bounded by $\mathcal{O}(N2^{M})$. In the
following, we will give more detailed estimates for the expected number of
computational nodes and complexity of the algorithm.

###  IV.2 Accounting for the expectation values

So far, we discussed the expansion of the dressed Hamiltonian. In turn, the
loss function is given by its expectation value in the all-zero state (18). A
likely scenario is that the majority of the final Pauli observables have
vanishing expectations, and hence do not contribute to the loss function. This
observation allows significantly increasing the efficiency of the computation
by pruning unfit branches in advance.

Let $\mathbb{F}_{2}^{n}$ be a vector space of binary strings of length $n$.
For $k=(k_{1},\dots,k_{n})\in\mathbb{F}_{2}^{n}$ denote
$\boldsymbol{Z}(k)=\bigotimes_{i=1}^{n}Z^{k_{i}}$ ($\boldsymbol{X}(k)$ is
defined similarly). For any $n$-qubit Pauli operator $P$ one can define two
vectors $P_{Z},P_{X}\in\mathbb{F}_{2}^{n}$ such that

$\displaystyle P=\alpha\,\,\boldsymbol{Z}(P_{Z})\boldsymbol{X}(P_{X})\ ,$ (23)

where $\alpha$ is a phase factor. One can think of $P_{Z}$ and $P_{X}$ as
coordinates of $P$ in the basis of Pauli $Z$ and $X$ strings, respectively.
With the notation in place, we can explain how expectation values can be taken
into account during the expansion of the dressed Hamiltonian.

A Pauli string $P$ has a non-zero expectation value $\Braket{0}{P}{0}\neq 0$
iff $P_{X}=(0,\ldots,0)$. First assume, for simplicity, that $X$-vectors of
the first $N$ Pauli generators $(P_{1})_{X},\ldots,(P_{N})_{X}$ are linearly
independent and constitute a basis in $\mathbb{F}_{2}^{N}$. This implies that
for every $O$, there is a unique vector $k\in\mathbb{F}_{2}^{N}$ such that

$O_{X}=k_{1}(P_{1})_{X}+\ldots+k_{N}(P_{N})_{X}\ .$ (24)

Therefore, among all $2^{N}$ possible observables of the form
$P_{1}^{t_{1}}\ldots P_{N}^{t_{N}}O$, that can be produced during recursive
expansion of $(P_{1}\dots P_{N}|O)$, only a single one with $t_{i}=k_{i}$ can
yield a non-zero expectation value (all other terms will have a non-zero
$X$-component). Thus, instead of generating the full recursive expansion of
$(P_{1}\dots P_{N}|O)$, which can contain up to $2^{N}$ nodes, we can first
find $k$ from (24), and then check if this operator actually appears in the
dressed Hamiltonian, i.e. it is compatible with the branching rules. This
yields an exponential saving for large $N$.

Now let us lift the restriction of the first $N$ Pauli generators forming the
basis. The necessary condition for $(P_{1}\ldots P_{m}|O)$ to have a non-zero
expectation is that $O_{X}$ is contained in the span of
$(P_{1})_{X},\ldots,(P_{m})_{X}$. Therefore, for each newly generated
computational node, we can test if this condition is satisfied. If it is not,
all final observables stemming from the expansion of this node have zero
expectation values, and the node can be disregarded.

As discussed in App.B, there is some room for further optimization based on
the freedom to permute commuting Pauli generators. Also, at this point, we
would like to spell out explicitly an elementary observation about the Fourier
series of Clifford+Pauli circuits. For a single Pauli observable, all non-zero
coefficients are given by averages of Pauli strings, and hence equal to $\pm
1$.

###  IV.3 Truncated Fourier series as an approximation

Having many terms at low Fourier levels appears to be convenient, because this
partially reduces a proliferation of coefficients at subsequent levels. This
is further reinforced by the observation that each individual term at a lower
level contributes exponentially more to the loss function than a term at a
higher level. Intuitively, this is because the average absolute value of a
trigonometric monomial of order $m$ is $(\pi/2)^{-m}$ and decays exponentially
with degree $m$. At the same time, there can be exponentially more terms at
higher levels. We can quantify this trade-off by evaluating the $L^{2}$ norm
of the loss function. From orthogonality of trigonometric monomials (12) it
follows

$\displaystyle||F||^{2}:=\Braket{}{F(\boldsymbol{\phi})}{{}^{2}}_{\boldsymbol{\phi}}=\sum_{m=0}^{M}2^{-m}l(m)\
,$ (25)

where $l(m)$ is the number of non-zero Fourier terms in the expansion of the
loss function at level $m$, which is upper bounded by the number of non-zero
terms in the expansion of the dressed Hamiltonian ($l(m)\leq n(m)$). Note that
$n(m)-l(m)$ is the number of operators in the dressed Hamiltonian at level $m$
with zero expectation value (i.e. with non-trivial $X$ coordinate). Using
(22), we can then bound the $L^{2}$ norm of the loss function:

$\displaystyle||F||^{2}\leq 1\ .$ (26)

Let $F^{(m)}(\boldsymbol{\phi})$ denote the Fourier series truncated to the
first $m$ levels. Then,

$\displaystyle||F^{(m)}-F||^{2}\leq\sum_{k=m+1}^{M}2^{-k}l(k)\leq 1-\Delta(m)\
.$ (27)

If $\Delta(m)$ is close to 1, i.e. sufficiently many terms are concentrated up
to level $m$, the truncated Fourier series $F^{(m)}(\boldsymbol{\phi})$ gives
a good approximation to the full loss function. Note that our recursive
expansion generates $F^{(m)}(\boldsymbol{\phi})$ level by level, so the
quality of the approximation can be gauged dynamically, and the computation
stopped when the necessary accuracy is reached.

We need to mention two caveats related to the approximation result stated.
First, while closeness in $L^{2}$ norm guarantees good approximation for most
parameter configurations, it does not translate directly into point-wise
convergence. Indeed, while the average absolute value of higher-level
monomials (4) is exponentially suppressed, their maximum values are
independent of the order
($\operatorname{max}_{\boldsymbol{\phi}}|{t}_{I}(\boldsymbol{\phi})|=1$).
Second, the bound (27) may be too weak in practice, as it effectively assumes
that all the final observables in the dressed Hamiltonian expansion above
level $m$ have non-zero expectation values. In practice, we expect that only
an exponentially small fraction of observables contributes to the loss
function norm. Properly taking this into account can significantly strengthen
the bound, but requires accounting for the structure of a particular circuit
at hand. We illustrate this in a random circuit model discussed in Sec. V.

###  IV.4 Is there a more efficient algorithm?

((a))

((b))

Figure 2: (a) A sketch of the hierarchy of scales in the problem. (b) A sample
computational tree taking into account different properties of nodes. Dashed
arrows lead to nodes $(P_{1}\dots P_{m}|O)$ for which $O_{X}$ does not lie in
the span of $(P_{i})_{X}$. These nodes do not contribute to the loss function
and are not actually generated by the algorithm. The number of remaining nodes
quantifies the complexity of the algorithm. Only two final observables with
non-vanishing expectation values are present, and the one at the second level
contributes most to the loss function norm.

As sketched in Fig. 2(a), our discussion features four different scales. The
largest scale (I) is set by the number of terms in the expansion of the
dressed Hamiltonian. It depends only on the structure of the circuit and the
Hamiltonian, and can contain up to $2^{M}$ terms. Note that the total number
of nodes in the computational tree Fig. 1 can only exceed the number of final
observables by a constant factor, so the recursive expansion algorithm is
optimal for computing the Fourier series of the dressed Hamiltonian.

Another scale (III) corresponds to the number of non-zero terms in the Fourier
expansion of the loss function. It quantifies the very complexity of
describing the loss function by its Fourier series, and implies a limit to
when such a description can be practical. Also, as discussed in Sec.(IV.3),
the truncated Fourier series can furnish a good approximation to the full loss
function, while containing only a tiny fraction of all non-zero terms. Hence,
we associate a separate scale (IV) to it.

The most relevant scale in practice, however, is set the by number of
computation nodes (II). It quantifies the complexity of the algorithm. Without
accounting for the expectation values, it simply coincides with the number of
terms in the dressed Hamiltonian. In Sec. IV.2 we explained how to prune the
branches of the dressed Hamiltonian expansion, with the expected saving being
exponential in the number of qubits. Still, in general, this leaves a large
gap between the number of computational nodes and the number of non-zero terms
in the loss function. Indeed, assume for simplicity that the first $N$ Pauli
generators span an $X$-basis, and consider a computational node
$(P_{1},\dots,P_{N}|O)$. For generic $P_{i}$ and $O$ this node is
exponentially unlikely to make a non-zero contribution to the loss function.
Indeed, there is a unique combination $P_{1}^{k_{1}}\dots P_{N}^{k_{N}}O$ that
has a non-zero expectation value in the all-zero state, but this very
combination is unlikely to satisfy all the branching rules, and hence actually
appear in the recursive expansion. Hence, most nodes of the type
$(P_{1},\dots,P_{N}|O)$ do not contribute to the loss function and are not
necessary to generate in the first place. Can a more efficient pruning
algorithm be developed?

Let us formalize the question. For Pauli strings $P_{i},P_{j}$, set $\langle
P_{i},P_{j}\rangle=1$ if they anti-commute and $\langle P_{i},P_{j}\rangle=0$
otherwise. Note that for any three Pauli strings $P_{i},P_{j},P_{k}$ it holds
that $\langle{P_{i},P_{j}P_{k}}\rangle=\langle
P_{i},P_{j}\rangle+\langle{P_{i},P_{k}}\rangle$ (here and in the following
$({\rm mod}~{}2)$ is implied). All possible final observables in the expansion
of the dressed Hamiltonian are of the form

$O(k)=P_{1}^{k_{1}}\dots P_{M}^{k_{M}}H,$ (28)

with some $k\in\mathbb{F}_{2}^{M}$. In order for $O(k)$ to actually appear in
the set of final observables, $k$ has to be consistent with the branching
rules. If $k_{i}=1$, the Pauli string $P_{i}$ must anti-commute with
$P_{i+1}^{k_{i}+1}\dots P_{M}^{k_{M}}H$, while $k_{i}=0$ implies no
constraints. These conditions can be expressed as $M$ equations
($i=1,\dots,M$)

$\displaystyle
k_{i}=k_{i}\braket{P_{i},H}+\sum_{j=i+1}^{M}k_{i}k_{j}\braket{P_{i},P_{j}}\ .$
(29)

The condition for $O(k)$ to have a non-zero expectation value is

$\displaystyle H_{X}=\sum_{i=1}^{M}k_{i}(P_{i})_{X}=0\ .$ (30)

This relation contains $N$ constraints for an $N$-qubit problem.

Together, branching constraints (29) and $X$-constraints (30) present an
instance of a boolean multivariate quadratic problem (Boolean MQ), which is
known to be NP-hard. State-of-the-art algorithms [29, 30] have worst-case time
complexities around $2^{0.69M}$ to find a solution or prove one does not
exist. Due to the recurrent structure of equations (29), our Boolean MQ
instances are significantly simpler than the general case. As shown in the
next section, for random circuits, which are expected to capture the worst
case behavior in practice, time complexity around $2^{0.59M-N}$ is sufficient
to find all solutions. Thus, while generic algorithms for the Boolean MQ
problem are unlikely to be useful directly, there is a possibility that more
efficient pruning techniques can be adopted in our scheme, narrowing the gap
between the number of computational nodes and non-zero coefficients in the
loss function.

##  V Case studies

So far, we discussed general properties of the Fourier expansion for
Clifford+Pauli circuits. In this section, we consider several specific
examples that showcase how the expansions are structured in practice. We will
both make analytic estimates and put the classical algorithm to work in
numeric simulations.

###  V.1 Random circuits

((a))

((b))

Figure 3: (a) Level distribution of the number of terms $l(m)$ and norm
$\nu(m)$ (both normalized) in Fourier expansion of the dressed Hamiltonian for
random circuits. Scatter plots are mean values computed from the simulations,
filled areas quantify standard deviations. Solid curves are theoretical
predictions. (b) Complexity of the algorithm for random Pauli circuits as a
function of the number of qubits $N$ with depth $M=N/\log\frac{3}{2}$. Details
of numerical experiments are discussed in App. C.1.

We first study the case where all the Pauli generators, as well as the
Hamiltonian, are random Pauli operators with the support on all of the qubits.
In this setup, it is simple to give probabilistic estimates for the expected
distribution of Fourier terms in the loss function.

In fact, the assumption that the Pauli generators are random is not necessary
as long as the observable is random. Therefore, we expect this behavior to
capture well the asymptotic limit of most sufficiently deep circuits. Indeed,
even if the original Hamiltonian is local, as we go down the computational
tree (see Fig. 1), the observables at the intermediate computational nodes
become ever more scrambled, and eventually behave as random Pauli operators.
The argument is not rigorous, of course, and circuits that do not confirm to
this pattern can be constructed. Nevertheless, the behavior of random circuits
should give a useful reference point for generic circuits.
Distribution of terms in the dressed Hamiltonian.

We will first look at the coarse-grained characteristics of the dressed
Hamiltonian expansion, such as the number of non-zero terms at level $m$,
denoted by $n_{M}(m)$, and the total number of terms $n_{M}=\sum_{m}n_{M}(m)$
(here we add a subscript $M$ to emphasize dependence on the total number of
parameterized gates). When all $P_{m}$ and $H$ are random, the probability of
branching at each computational node is $\frac{1}{2}$. Therefore, on average,
each iteration of the algorithm increases the total number of nodes $n_{M}$ by
a factor of $\frac{3}{2}$, leading to

$\displaystyle n_{M}=\left(\frac{3}{2}\right)^{M}\ .$ (31)

The same reasoning applies to the number of non-zero computational nodes at
each level $n_{M}(m)$, which hence satisfies the recurrence relation
$n_{M+1}(m)=\frac{1}{2}n_{M}(m)+n_{M}(m-1)$. Solving it yields

$\displaystyle n_{M}(m)=2^{m-M}\binom{M}{m}\ .$ (32)

One can check that $\sum_{m=0}^{M}n_{M}(m)=n_{M}$.
Distribution of terms in the loss function.

So far we discussed the distribution of Fourier terms in the dressed
Hamiltonian, and now we turn to the loss function. Since every final
observable of the dressed Hamiltonian is again a random Pauli, it has $2^{-N}$
probability of having non-zero expectation value. Therefore, the distribution
of Fourier terms by level $l_{M}(m)$ is simply

$l_{M}(m)=2^{-N}n_{M}(m).$ (33)

The expected number of all non-zero terms in the loss function is

$\displaystyle l_{M}=2^{-N}\left(\frac{3}{2}\right)^{M}\approx 2^{0.59M-N}\ .$
(34)

Now, we recall that each term at level $m$ contributes exactly $2^{-m}$ to the
$L^{2}$ norm. From the distribution of Fourier coefficients by level, we can
derive the distribution of norm by level $\nu_{M}(m)=2^{-m}l_{M}(m)$,
explicitly given by

$\displaystyle\nu_{M}(m)=2^{-N-M}\binom{M}{m}\ .$ (35)

Note that $\sum_{m=0}^{M}\nu_{M}(m)=2^{-N}$. At Fig. 3(a) we plot results of
numerical simulations for random circuits, which convincingly confirm our
estimates.
Accuracy of a truncated expansion.

We can now address the question of how many terms need to be included in the
loss function to give a good $L^{2}$ approximation. Since binomial
distribution (35) is symmetric around $m=\frac{M}{2}$, including Fourier terms
up to level $\frac{M}{2}$ will on average account for $50\%$ of the norm. The
total number of nodes at the included levels can be estimated as
$l_{M}(\frac{M}{2})\sim 2^{\frac{M}{2}-N}$. While in the large $M$ limit this
is an exponentially small fraction of all terms (34), the number of relevant
terms itself is still an exponential in $M$.
Complexity of the algorithm and simulation limits.

Exponential increase in the number of relevant terms with $M$ clearly limits
the depth of the circuits we can address. Importantly, the property that most
final observables have zero expectation values in turn limits the number of
qubits $N$ we can meaningfully simulate. While in principle, the number of
qubits is only limited by the simulation cost of Clifford circuits, to yield a
non-zero loss function the number of Pauli rotation gates $M$ needs to
increase with $N$.

Requiring the number of non-zero terms in the loss function (34) to stay of
order $\mathcal{O}(1)$ as we increase $N$ requires to scale the depth of the
circuit as $M\sim N/\log_{2}\frac{3}{2}$. Thanks to the subroutine filtering
by the expectation value, the algorithm only branches during processing of the
first $M-N$ gates, leading to the number of computational nodes around
$\sim\left(\frac{3}{2}\right)^{M-N}\simeq 2^{0.59(M-N)}$. Therefore, the
number of computational nodes generated per non-zero term in the loss function
grows with the number of qubit roughly as

$\sim 2^{0.41N}\approx 10^{N/8}$ (36)

As reported in Fig.3(b), this scaling is well confirmed numerically.

With a computational budget to process $10^{6}$ nodes, which takes minutes
with a basic implementation run on a laptop, the Fourier expansion of the loss
function for a $50$-qubit random circuit with $85$ Pauli rotation gates can be
computed exactly. With resources to process $10^{12}$ nodes, which should be
feasible with an efficient implementation on a computational cluster,
$100$-qubit circuits with about $170$ Pauli gates can be handled.

###  V.2 QAOA

degree$\backslash$level | 1 | 2 | 3 | 4
---|---|---|---|---
2 | $10^{0.7\pm 0.08}$ | $10^{1.5\pm 0.9}$ | $10^{3.4\pm 1.7}$ | $10^{4.7\pm 2.6}$
3 | $10^{1.1\pm 0.1}$ | $10^{3.6\pm 0.4}$ | $10^{8.6\pm 0.9}$ | $10^{17.5\pm 2.2}$
4 | $10^{1.4\pm 0.3}$ | $10^{5.6\pm 0.5}$ | $10^{16.4\pm 1.3}$ | $10^{36.6\pm 3.0}$
Table 1: Estimated number of computational nodes to exactly compute Fourier
expansion for QAOA circuits of varying degree and level. Statistics is
collected over different graphs and observables, see App. C.2 for details.

Variational circuits appearing in practice are far from the random Pauli model
described above. Instead, they typically only involve local gates and
observables. If two local Pauli strings are supported on different subsets of
qubits they necessarily commute, and hence the probability of two generic
local operators anti-commuting is much smaller than $\frac{1}{2}$. Although we
expect the random Pauli model to describe well the large depth asymptotic,
circuits of shorter depth may behave quite differently. Indeed, until the
observables at computational nodes become sufficiently scrambled, the
branching probability is much smaller than $\frac{1}{2}$, and the complexity
growth is much slower. This allows to compute Fourier expansions for circuits
with much higher depths than anticipated for the random model.

As a case study, we consider instances of the Quantum Approximate Optimization
Algorithm (QAOA) for the MaxCut problem on regular graphs [14], which is the
most studied approach to combinatorial optimization. Two key characteristics
of a particular QAOA circuit is the degree of a graph $d$ and the number of
layers $p$, see App. C.2 for details.

Due to locality, shallow instances of QAOA allow for efficient classical
computation of the loss function independently of the number of qubits. To
every observable one associates a reverse light cone, containing all qubits
that are connected to the observable by the entangling gates. For graphs of
bounded degree, the size of the reverse light cone stays constant in the large
$N$ limit and classical computation of the loss function never involves
simulating quantum circuits larger than that size. For QAOA on a graph of
degree $d$ with $p$ layers the size of the reverse light cone is bounded by
[14]

$\displaystyle N_{c}=2\frac{(d-1)^{p+1}-1}{d-2}\ .$ (37)

Note that $N_{c}$ scales exponentially with the level $p$, so that large $p$
regime is the most difficult to simulate.

Large $p$ also implies large number of parametric gates, which is the key
limiting factor for our algorithm. At the same time, while the locality of
gates is not necessary in our approach, it certainly helps. Moreover, the
benefits of locality are incorporated automatically. Indeed, the Pauli
generators supported outside the reverse light cone of the Hamiltonian commute
with all observables at the computational nodes and are trivially eliminated.
Using a basic Monte-Carlo sampling, we estimate the expected complexity of our
algorithm to compute the full Fourier expansion of the loss function for
several values of $d$ and $p$. Results are reported in Tab. 1, details of
numerical simulations are specified in App.C.2. Though estimates are very
crude, they provide useful anchors for expected complexities.

We stress here that the special structure of QAOA circuits and observables
makes it possible to handle much larger number of parameters compared to
random circuits. Indeed, for a $d$-regular graph with $N$ nodes there are
exactly $|E|=Nd/2$ edges, so the number of parametric gates in our simulations
is given by

$\displaystyle M=p(N_{c}+|E|)=pN_{c}\left(\frac{d}{2}+1\right)\ .$ (38)

For instance, for $p=3,d=3$ this gives $N_{c}=30,M=225$, and a random circuit
in this setup would require processing
$\sim\left(\frac{3}{2}\right)^{M-N_{c}}\simeq 2.18\times 10^{34}$
computational nodes, clearly an unmanageable amount. In contrast, for the
actual QAOA instances, we find about $10^{9}$ computational nodes to be
sufficient. The distribution of Fourier terms in QAOA appears to be somewhat
irregular, although with clustering characteristics similar to the random
circuits of equivalent complexity, see App. C.2 for details.

###  V.3 Hardware-efficient circuits

One frequently studied design of variational algorithms is the hardware-
efficient form [15], where the circuit is constructed to give the maximal
expressivity with limited depth, efficiently using native hardware gates.

We consider hardware-efficient circuits with brick wall architecture, where
each block is built of an entangling CZ  gate and four single-qubit Pauli
rotations, and observables of weight $2$, see App. C.3 for details. In this
setting, we estimate the complexity of our algorithm to exactly compute the
loss function and report the results in Fig. 4. The takeaway is similar to the
QAOA case, locality of the circuit and the observable strongly reduces the
number of Fourier terms compared to the random case, so the computational
budget needed to process random circuits with $M\sim 180$ should be sufficient
to fully characterize the Fourier series for circuits with $M\sim 600$
parametric gates in the current local setup.

Figure 4: Estimated algorithm complexity for computing the full Fourier series
of a two-local hardware efficient circuit with 50 qubits and Pauli Hamiltonian
of weight two, as a function of the number of parametric gates. Details are
specified in App. C.3.

##  VI Discussion and outlook

We looked at some qualitative and quantitative aspects of Fourier series
expansion of VQA loss functions. The main observation is that restricting to
the class of Clifford+Pauli circuits allows giving a much finer picture than
possible for generic VQA. We presented an efficient classical algorithm for
computing the Fourier expansion level by level, with the worst case complexity
bounded by ${\cal O}(N2^{M})$ for a single Pauli observable. We estimated the
complexity of the algorithm and characterized the distribution of Fourier
terms in several examples, including random non-local Clifford+Pauli circuits
as well as more conventional local circuits such as QAOA and HEA. We
anticipate our findings to provide useful tools for further study of the
interplay between the properties of VQA and the structure of their Fourier
series expansion.

One major issue facing VQA is trainability, where two crucial obstacles are
barren plateaus [31, 32] and local minimums [13, 33]. Although many ansatz
structures, initialization and optimization heuristics have been proposed to
mitigate these problems (we refer to [34] for a nice summary, and to [35, 36]
for the discussion of the over-parameterized setting), they persist in many
practical scenarios. Interestingly, the Fourier series contains a wealth of
global information about the loss landscape, and may hence give a useful
perspective on these problems. For instance, we note that the variance of the
loss function gradient, frequently used to diagnose barren plateaus, is
naturally related to the Fourier coefficients

$\displaystyle\operatorname{var}_{\boldsymbol{\phi}}\left(\nabla
F(\boldsymbol{\phi})\right)^{2}=-\int
d\boldsymbol{\phi}\,\,F(\boldsymbol{\phi})\Delta
F(\boldsymbol{\phi})=\sum_{m=0}^{M}m\,l(m)\ .$ (39)

(To arrive at this expression first integrate by parts and then use the fact
that for every trigonometric monomial $t_{m}(\boldsymbol{\phi})$ of degree $m$
it holds $\Delta t_{m}(\boldsymbol{\phi})=-mt_{m}(\boldsymbol{\phi})$.)

Therefore, already coarse-grained characteristics of the Fourier expansion,
such as the distribution of terms by level $l(m)$, may provide a very useful
input. Note that this distribution can be estimated by a Monte-Carlo sampling,
even when the full computation of the loss function may be out of reach. Also,
it would be very interesting to evaluate the role of lower vs higher order
Fourier modes in shaping the loss landscape. Superficially, since the higher
order modes are more oscillatory, yet typically contribute less to the $L^{2}$
norm of the loss function, they may be a justified suspect in creating the
majority of spurious local minima.

Also, it looks promising to explore the potential of our algorithm for
classical computation of quantum mean values. For circuits with local gates
and constant depth, classical algorithms exist that scale favorably with the
number of qubits [37]. While gate locality helps in practice, our algorithm is
based on the properties of stabilizer circuits it is not principally
constrained by it. Neither it is constrained by the degree of entanglement
[38, 39] or circuit’s graph treewidth [40], which can be a bottleneck for
simulators based on tensor networks. Apparently, the flavor of our approach is
most similar to simulation schemes for circuits dominated by Clifford gates.
The parameter ranges that can be handled look similar, involving 50-100 qubits
and dozens to hundreds of non-Clifford gates [41, 42], but detailed benchmarks
are necessary to make a reasonable comparison. We should also stress that our
approach is neither a simulation algorithm (in a weak or strong sense) nor an
algorithm exclusively computing mean values. Instead, it returns the full (or
truncated) Fourier series of a VQA, and hence has an interesting and novel
character.

Acknowledgments. We thank Vsevolod Yashin and V Vijendran for useful
discussions. We thank the Priority 2030 program at the National University of
Science and Technology “MISIS” under the project K1-2022-027.

##  Appendix A Structure of a generic Fourier expansion

###  A.1 Level expansion and number of terms

For the sake of clarity, here we give a more detailed description of the
general structure of Fourier expansion for the loss function introduced in
Sec. II. Let us first illustrate trigonometric expansion of the unitary (3)
for the case with $M=2$ angles

$\displaystyle
U(\phi_{1},\phi_{2})=U_{00}\cos\frac{\phi_{1}}{2}\cos\frac{\phi_{2}}{2}+U_{01}\cos\frac{\phi_{1}}{2}\sin\frac{\phi_{2}}{2}+U_{10}\sin\frac{\phi_{1}}{2}\cos\frac{\phi_{2}}{2}+U_{11}\sin\frac{\phi_{1}}{2}\sin\frac{\phi_{2}}{2}\
.$ (40)

Note that trigonometric monomials here have the same degree and period $4\pi$.
Substituting such expansions into the definition of the loss function (7)
leads to the Fourier series expansion of the loss function. Applying
identities
$\cos^{2}\frac{\phi}{2}=\frac{1+\cos\phi}{2},\sin^{2}\frac{\phi}{2}=\frac{1-\cos\phi}{2},\cos\frac{\phi}{2}\sin\frac{\phi}{2}=\frac{\sin\phi}{2}$
leads to an expression involving trigonometric monomials of a smaller period
$2\pi$ and degrees up to $M$. For instance,

$\displaystyle
F(\phi_{1},\phi_{2})=\Braket{0}{U^{\dagger}(\phi_{1},\phi_{2})HU(\phi_{1},\phi_{2})}{0}=\frac{1}{4}\Braket{0}{U_{00}^{\dagger}HU_{00}}{0}\left(1+\cos\phi_{1}+\cos\phi_{2}+\cos\phi_{1}\cos\phi_{2}\right)+\dots$
(41)

and we wrote explicitly only the contribution from the first term.

More generally, terms in the expansion (9) assume the following form

$\displaystyle F_{0}=\text{const},\quad
F_{1}(\boldsymbol{\phi})=\sum_{i=1}^{M}F_{i}(\phi_{i}),\quad
F_{2}(\boldsymbol{\phi})=\sum_{i,j=1}^{M}F_{ij}(\phi_{i},\phi_{j}),\quad
F_{3}(\boldsymbol{\phi})=\sum_{i,j,k=1}^{M}F_{ijk}(\phi_{i},\phi_{j},\phi_{k}),\quad\dots$
(42)

Here

$\displaystyle\begin{split}&F_{i}(\phi_{i})=A_{i}\cos\phi_{i}+B_{i}\sin\phi_{i}\\\
&F_{ij}(\phi_{i},\phi_{j})=A_{ij}\cos\phi_{i}\cos\phi_{j}+B_{ij}\cos\phi_{i}\sin\phi_{j}+C_{ij}\sin\phi_{i}\cos\phi_{j}+D_{ij}\sin\phi_{i}\sin\phi_{j}\\\
&\dots\end{split}$ (43)

At each level $m$ there are $\binom{M}{m}$ subsets of parameters, enumerating
possible indices of homogeneous polynomials $F_{i_{1}\dots i_{m}}$. To define
the polynomial for each parameter configuration requires specifying $2^{m}$
coefficients. This leads to the counting (10) for the total number of
coefficients in the Fourier series.

###  A.2 Coefficients from averages

Trivially $F_{0}=\Braket{F(\boldsymbol{\phi})}_{\boldsymbol{\phi}}$. Higher
order terms can be obtained similarly. For instance,
$F_{i}(\phi_{i})=\Braket{F(\boldsymbol{\phi})-F_{0}}_{\boldsymbol{\phi}\neq\phi_{i}}$,
i.e. averaging over all angles except $\phi_{i}$ only leaves first-order
monomials involving a given angle $\phi_{i}$. Higher order terms can be found
recursively. Note that in this prescription, to compute terms at level $m$
requires to first compute and subtract the contribution of all levels below
$m$.

##  Appendix B Optimizing Pauli order

A Pauli form of a Clifford+Pauli circuit is not unique, because adjacent Pauli
rotations with commuting generators can be swapped. In principle, the number
of computational nodes may be sensitive to this ordering. We briefly mention
two possible optimizations along these lines.

Delayed branching. The first strategy may be to delay branching as long as
possible, by swapping the anti-commuting Pauli to the left. For illustration,
consider a Pauli circuit $(P_{1}\dots P_{N}P^{\prime}_{1}\dots
P^{\prime}_{N}|O)$ where all Pauli generators mutually commute, and in
addition $[P_{i},O]=0$. Generators $P^{\prime}_{i}$ may not commute with $O$.
While processing operators $P_{i}^{\prime}$ up to $2^{N}$ nodes may be
generated. After that, assuming that $P_{i}$ span an $X$-basis, no new nodes
will be produced. Had we started with the circuit $(P^{\prime}_{1}\dots
P^{\prime}_{N}P_{1}\dots P_{N}|O)$ instead, which is equivalent by assumption,
generators $P_{i}$ would be eliminated right away and no branching would be
required (assuming $P^{\prime}_{i}$ also span an $X$-basis).

Early pruning. The freedom to swap commuting Pauli generators can also be used
to enforce early pruning. As an illustration, consider circuit
$(PP^{\prime}_{1}\dots P^{\prime}_{m}|O)$, where $P$ commutes with all
$P^{\prime}_{i}$ and, moreover, $P$ is an independent generator, in the sense
that $P_{X}$ does not lie in the span of $(P_{i}^{\prime})_{X}$. If we process
the circuit directly, up to $2^{m}$ computational nodes are generated and then
tested against compatibility with the last generator $P$. However, in the
equivalent setting $(P^{\prime}_{1}\dots P^{\prime}_{m}P|O)$ compatibility
with $P$ is tested right away, and on average half of observables $O$ will be
pruned leading on average to $2^{m-1}$ computation nodes. If there are $n$
generators of $P$-type, savings up to $2^{-n}$ can be expected.

Note that both optimizations rely on the presence of large streaks of
commuting Pauli generators. Therefore, for the non-local random Pauli circuits
very little is to be gained, while, for structured local circuits (e.g. QAOA)
the savings might be substantial.

##  Appendix C Details of numerical computations

An implementation of the algorithm and the data presented in the paper are
available at GitHub repository [43].

###  C.1 Random circuits

Statistics in Fig. 3(a) is collected from 20 random non-local Pauli circuits
with $N=30$ qubits and depth $M=25$. Note that the distribution of norm
$\nu(m)$ is not an independent characteristic, but is computed from the
distribution of nodes $l(m)$ according to $\nu(m)=2^{-m}l(m)$. In Fig. 5(a) we
also plot the distribution of Fourier terms for random circuits, where only
the observable is non-local, while the circuit consists of random Pauli
exponentials of weight $2$. As expected, the average values are the same as
for the non-local random circuits, but the fluctuations around the mean are
much higher.

((a))

((b))

((c))

Figure 5: Normalized distribution of non-zero Fourier coefficients in (a)
local circuits with random observables (b) QAOA circuits with $d=3,p=2$ and
(c) brick wall HEA circuits with $N=50$ and $M=304$. Statistics is collected
for random circuits structures (except for HEA, where it is fixed) and random
observables. Solid lines are theoretical curves for random circuits of the
equivalent complexity.

###  C.2 QAOA

@C=1.0em @R=0.2em @!R
& H 1 ZZ (γ_0) 3 2 R_X (β_0)
H ZZ (γ_2) 1 ZZ (γ_3) 2 R_X (β_1)
H ZZ (γ_1) ZZ (γ_4) 1 ZZ (γ_5) R_X (β_2)
H R_X (β_3)

((a))

@C=1.0em @R=0.2em @!R
& 1 R_X R_Z 1 R_X R_Z
R_X R_Z 1 R_X R_Z R_X R_Z
1 R_X R_Z R_X R_Z 1 R_X R_Z
R_X R_Z 1 R_X R_Z R_X R_Z
R_X R_Z

((b))

Figure 6: Circuit layouts for (a) Single-layer QAOA instance on a 3-regular
graph with 4 nodes (b) Hardware-efficient brick wall circuit with 5 qubits
(all rotation gates have different parameters, not indicated explicitly).

Given a graph $G$ with edges $E_{ij}$, a single layer of the QAOA circuit
consists of two-qubit Pauli rotation gates
$Z_{i}Z_{j}(\gamma_{ij})=\exp(-\imath Z_{i}Z_{j}\gamma_{ij}/2)$ (here $Z_{i}$
stands for a Pauli string with $Z$ on $i$-th position and identities on all
others) for each edge $E_{ij}$, followed by the sequence of single-qubit
$X_{i}(\beta_{i})$ gates placed on every qubit, see Fig. 6(a) for an example.
Note that in the standard formulation of QAOA,
$\gamma_{ij}=\gamma,\beta_{i}=\beta$ i.e. all $ZZ$ gates and all $X$ gates
have the same parameters within each layer. However, taking this correlation
into account does not simplify our analysis and we will not impose it. Single
layer repeated $p$ times gives an instance of QAOA with $p$ layers. The
Hamiltonian is given by the sum of all $Z$-Pauli generators
$H=\sum_{E_{ij}}Z_{i}Z_{j}$.

To estimate the complexity of computing the Fourier series for instances of
QAOA of degree $d$ with $p$ layers, it is sufficient to consider circuits of
size $N_{c}$ (37). Statistics presented in Tab. 1 was collected in the
following way. For each configuration $(d,p)$, we generate 20 random regular
graphs of degree $d$ and for each graph choose a single observable
$Z_{i}Z_{j}$ at random. Then for each circuit-observable pair we estimate the
complexity of the algorithm using the Monte-Carlo technique C.4 with $10^{4}$
samples. Fluctuations are significant both due to inaccuracies of the Monte-
Carlo sampling, but more importantly, due to inhomogeneous nature of data
gathered over different graphs and observables. To represent large
fluctuations more meaningfully, we compute average values and deviations in
the log scale, i.e. for the exponents of the expected number of computations
nodes. Hence the format used in Tab. 1. We stress that the numbers reported
are estimated number of computational nodes, i.e. the expected complexity of
the algorithm. The number of non-zero terms in the loss function is less by
orders of magnitude, but also harder to estimate with reasonable precision.

Besides the overall complexity, it might be of interest to look at the
particular distributions of the Fourier terms in the loss function. To this
end, in Fig. 5(b) we plot distributions on non-zero Fourier coefficients by
level. Statistics is gathered for 100 different instances of $d=3,p=2$ QAOA
circuits, with a single randomly chosen QAOA observable for each circuit.
Theoretical curve for random circuits with $M=30$ is plotted for comparison.
We note that while the distribution of the coefficients looks close enough to
the random case, the distribution of norm has significant differences and
large fluctuations. In particular, the norm appears to be concentrated at
lower Fourier levels than expected for random circuits (after adjusting for
the distribution of terms), which might make the truncated schemes quite
useful.

###  C.3 Hardware-efficient circuits

An example of a brick wall hardware-efficient circuit that we consider is
shown in Fig. 6(b). Statistics for Fig.4 is collected for $N=50$ qubit
circuits with up to $M=600$ rotation gates. For each $M$ we take 10 random
observables of weight $2$ and estimate the number of computational nodes with
$10^{4}$ Monte-Carlo samples. Similarly to the QAOA case, due to large
individual fluctuations, we compute means and deviations in the logarithmic
scale.

In Fig. 5(c) we also plot the distribution of Fourier coefficients for $N=50$
circuit with $M=308$ parametric gates (corresponding to 77 CZ gates), averaged
over 100 random Pauli observables of weight 2. Similarly to the QAOA case, the
most interesting part is the norm distribution, which shows large fluctuations
and also the tendency to cluster at lower levels (compared to a random circuit
with a similar distribution of coefficients).

###  C.4 Estimating complexity from Monte-Carlo sampling

It is useful to estimate the runtime/complexity of the algorithm without
performing the full computation. To do this, we can use a simple version of
Monte-Carlo sampling to probe the structure of the computational tree Fig. 1.
First, we will describe a version that allows to estimate the number of terms
in the expansion of the dressed Hamiltonian. To this end, we can use the basic
algorithm to traverse the computational tree, but instead of keeping both
branches at each split, we choose one at random and disregard the other. This
procedure produces a single complete branch from the computational tree. The
probability to generate any particular branch is $2^{-m}$ with $m$ being the
number of splits performed in the process. Let $s(m)$ be the number of samples
at with $m$ splits. We estimate the total number of terms $n$ as

$\displaystyle n\approx\frac{\sum_{m=0}^{M}2^{m}s(m)}{\sum_{m=0}^{M}s(m)}\ .$
(44)

Note that in this case the number of splits $m$ is the same as the Fourier
level of the resulting term.

To estimate the complexity of the actual algorithm, i.e. the number of
computation nodes generated during the computation, we need to take into
account the pruning based on the expectation values, see Sec. IV.2. The
sampling prescription above is modified in a simple way. We traverse the
computational tree and choose branches at random. Importantly, we only choose
from the branches compatible with the pruning restrictions. If one of the
branches in not admissible, the actual algorithm will not generate additional
nodes and the probability of the sample does not need to be updated. In this
case, the Fourier level of the resulting term does not correspond to the
probability of sampling. The prescription (44) remains unchanged.

## References

  * Cerezo  _et al._ [2021a] 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 2021 3:9 3, 625 (2021a), arXiv:2012.09265 .
  * Preskill [2018] J. Preskill, Quantum 2, 1 (2018), arXiv:1801.00862 .
  * Bharti  _et al._ [2021] K. Bharti, A. Cervera-Lierta, T. H. Kyaw, T. Haug, S. Alperin-Lea, A. Anand, M. Degroote, H. Heimonen, J. S. Kottmann, T. Menke, W.-K. Mok, S. Sim, L.-C. Kwek, and A. Aspuru-Guzik, Reviews of Modern Physics 94 (2021), 10.1103/RevModPhys.94.015004, arXiv:2101.08448v2 .
  * Biamonte  _et al._ [2017] J. Biamonte, P. Wittek, N. Pancotti, P. Rebentrost, N. Wiebe, and S. Lloyd, Nature 549, 195 (2017), arXiv:1611.09347 .
  * Schuld  _et al._ [2020] M. Schuld, R. Sweke, and J. J. Meyer, Physical Review A 103 (2020), 10.1103/PhysRevA.103.032430, arXiv:2008.08605v2 .
  * Gil Vidal and Theis [2020] F. J. Gil Vidal and D. O. Theis, Frontiers in Physics 8, 297 (2020), arXiv:1901.11434 .
  * AtchadeAdelomou and Larson [2023] P. AtchadeAdelomou and K. Larson,  , 1 (2023), arXiv:2302.00105 .
  * Casas and Cervera-Lierta [2023] B. Casas and A. Cervera-Lierta,  , 1 (2023), arXiv:2302.03389 .
  * Schreiber  _et al._ [2022] F. J. Schreiber, J. Eisert, and J. J. Meyer,  (2022), 10.48550/arxiv.2206.11740, arXiv:2206.11740 .
  * Landman  _et al._ [2023] J. Landman, S. Thabet, C. Dalyac, H. Mhiri, and E. Kashefi, “Classically Approximating Variational Quantum Machine Learning with Random Fourier Features,” (2023).
  * Fontana  _et al._ [2022a] E. Fontana, I. Rungger, R. Duncan, and C. C. Cˆırstoiu,  (2022a), 10.48550/arxiv.2208.05958, arXiv:2208.05958 .
  * Fontana  _et al._ [2022b] E. Fontana, I. Rungger, R. Duncan, and C. C. Cˆırstoiu,  (2022b), 10.48550/arxiv.2206.08811, arXiv:2206.08811 .
  * Bittel and Kliesch [2021] L. Bittel and M. Kliesch,  (2021), 10.1103/PhysRevLett.127.120502, arXiv:2101.07267 .
  * Farhi  _et al._ [2014] E. Farhi, J. Goldstone, and S. Gutmann,  (2014), 10.48550/arxiv.1411.4028, arXiv:1411.4028 .
  * 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), arXiv:1704.05018 .
  * Romero  _et al._ [2018] J. Romero, R. Babbush, J. R. McClean, C. Hempel, P. J. Love, and A. Aspuru-Guzik, Quantum Science and Technology 4, 014008 (2018), arXiv:1701.02691 .
  * Gottesman [1997] D. Gottesman,  2008 (1997), arXiv:9705052 [quant-ph] .
  * Aaronson and Gottesman [2004] S. Aaronson and D. Gottesman, Physical Review A - Atomic, Molecular, and Optical Physics 70 (2004), 10.1103/PhysRevA.70.052328, arXiv:0406196v5 [quant-ph] .
  * Cao  _et al._ [2019] Y. Cao, J. Romero, J. P. Olson, M. Degroote, P. D. Johnson, M. Kieferová, I. D. Kivlichan, T. Menke, B. Peropadre, N. P. Sawaya, S. Sim, L. Veis, and A. Aspuru-Guzik, Chemical Reviews 119, 10856 (2019), arXiv:1812.09976 .
  * Bauer  _et al._ [2020] B. Bauer, S. Bravyi, M. Motta, and G. K.-L. Chan, Chemical Reviews 120, 12685 (2020), arXiv:2001.03685v2 .
  * Sapova and Fedorov [2022] M. D. Sapova and A. K. Fedorov, Communications Physics 5, 199 (2022).
  * Mitarai  _et al._ [2020] K. Mitarai, Y. Suzuki, W. Mizukami, Y. O. Nakagawa, and K. Fujii, Physical Review Research 4, 033012 (2020), arXiv:2011.09927 .
  * Cheng  _et al._ [2022] M. H. Cheng, K. E. Khosla, C. N. Self, M. Lin, B. X. Li, A. C. Medina, and M. S. Kim,  , 1 (2022), arXiv:2207.01539 .
  * Ravi  _et al._ [2022] G. S. Ravi, P. Gokhale, Y. Ding, W. M. Kirby, K. N. Smith, J. M. Baker, P. J. Love, H. Hoffmann, K. R. Brown, and F. T. Chong,  , 1 (2022), arXiv:2202.12924 .
  * Ryabinkin  _et al._ [2020] I. G. Ryabinkin, R. A. Lang, S. N. Genin, and A. F. Izmaylov, Journal of Chemical Theory and Computation 16, 1055 (2020), arXiv:1906.11192 .
  * Brown  _et al._ [2022] J. Brown, M. P. Coons, E. Lloyd, A. Fleury, K. Bieniasz, V. Senicourt, and A. Zaribafiyan,  , 1 (2022), arXiv:2211.10501 .
  * Schleich  _et al._ [2023] P. Schleich, J. Boen, L. Cincio, A. Anand, J. S. Kottmann, S. Tretiak, P. A. Dub, and A. Aspuru-Guzik,  (2023), 10.48550/arxiv.2303.01221, arXiv:2303.01221 .
  * Hadfield [2018] S. Hadfield,  (2018), 10.48550/arxiv.1805.03265, arXiv:1805.03265 .
  * Dinur [2020] I. Dinur, Proceedings of the Annual ACM-SIAM Symposium on Discrete Algorithms , 2550 (2020), arXiv:2005.04800 .
  * Barbero  _et al._ [2022] S. Barbero, E. Bellini, C. Sanna, and J. Verbel, Discrete Applied Mathematics 309, 13 (2022).
  * McClean  _et al._ [2018] J. R. McClean, S. Boixo, V. N. Smelyanskiy, R. Babbush, and H. Neven, Nature Communications 9, 1 (2018), arXiv:1803.11173 .
  * Cerezo  _et al._ [2021b] M. Cerezo, A. Sone, T. Volkoff, L. Cincio, and P. J. Coles, Nature Communications 2021 12:1 12, 1 (2021b), arXiv:2001.00550 .
  * Anschuetz and Kiani [2022] E. R. Anschuetz and B. T. Kiani,  , 1 (2022), arXiv:2205.05786 .
  * Wang  _et al._ [2023] Y. Wang, B. Qi, C. Ferrie, and D. Dong,  (2023), arXiv:2302.06858 .
  * Larocca  _et al._ [2021] M. Larocca, N. Ju, D. García-Martín, P. J. Coles, and M. Cerezo,  (2021), arXiv:2109.11676 .
  * Kim  _et al._ [2020] J. Kim, J. Kim, and D. Rosa, Physical Review Research 3, 023203 (2020), arXiv:2010.00157 .
  * Bravyi  _et al._ [2019a] S. Bravyi, D. Gosset, and R. Movassagh, Nature Physics 2021 17:3 17, 337 (2019a), arXiv:1909.11485 .
  * Vidal [2003] G. Vidal, Physical Review Letters 91, 147902 (2003), arXiv:0301063 [quant-ph] .
  * Zhou  _et al._ [2020] Y. Zhou, E. M. Stoudenmire, and X. Waintal, Physical Review X 10, 041038 (2020), arXiv:2002.07730 .
  * Markov and Shi [2008] I. L. Markov and Y. Shi, SIAM Journal on Computing 38, 963 (2008), arXiv:0511069 [quant-ph] .
  * Bravyi and Gosset [2016] S. Bravyi and D. Gosset, Physical Review Letters 116 (2016), 10.1103/PhysRevLett.116.250501, arXiv:1601.07601v3 .
  * Bravyi  _et al._ [2019b] S. Bravyi, D. Browne, P. Calpin, E. Campbell, D. Gosset, and M. Howard, Quantum 3, 181 (2019b), arXiv:1808.00128v2 .
  * Nemkov  _et al._ [2023] N. Nemkov, E. Kiktenko, and A. Fedorov, “Fouriervqa,” https://github.com/idnm/FourierVQA (2023).