tcb@breakable

# Connecting the Hamiltonian structure to the QAOA energy and Fourier
landscape structure

Michał Stęchły<EMAIL_ADDRESS>Zapata Computing Canada Inc., 25
Adelaide St E, Toronto, ON, M5C3A1  Lanruo Gao  Zapata Computing Canada Inc.,
25 Adelaide St E, Toronto, ON, M5C3A1  Boniface Yogendran  Zapata Computing
Canada Inc., 25 Adelaide St E, Toronto, ON, M5C3A1  Enrico Fontana  University
of Strathclyde, 16 Richmond St, Glasgow G1 1XQ, UK  National Physical
Laboratory, Hampton Rd, Teddington TW11 0LW, UK  Manuel S. Rudolph  Zapata
Computing Canada Inc., 25 Adelaide St E, Toronto, ON, M5C3A1  Institute of
Physics, Ecole Polytechnique Fédérale de Lausanne (EPFL), CH-1015 Lausanne,
Switzerland

###### Abstract

In this paper, we aim to expand the understanding of the relationship between
the composition of the Hamiltonian in the Quantum Approximate Optimization
Algorithm (QAOA) and the corresponding cost landscape characteristics. QAOA is
a prominent example of a Variational Quantum Algorithm (VQA), which is most
commonly used for combinatorial optimization. The success of QAOA heavily
relies on parameter optimization, which is a great challenge, especially on
scarce noisy quantum hardware. Thus understanding the cost function landscape
can aid in designing better optimization heuristics and therefore potentially
provide eventual value. We consider the case of 1-layer QAOA for Hamiltonians
with up to 5-local terms and up to 20 qubits. In addition to visualizing the
cost landscapes, we calculate their Fourier transform to study the
relationship with the structure of the Hamiltonians from a complementary
perspective. Furthermore, we introduce metrics to quantify the roughness of
the landscape, which provide valuable insights into the nature of high-
dimensional parametrized landscapes. While these techniques allow us to
elucidate the role of Hamiltonian structure, order of the terms and their
coefficients on the roughness of the optimization landscape, we also find that
predicting the intricate landscapes of VQAs from first principles is very
challenging and unlikely to be feasible in general.

##  I Introduction

Variational Quantum Algorithms (VQAs) are a family of hybrid quantum-classical
algorithms designed for near-term intermediate scale quantum (NISQ) [1]
devices. The goal of a VQA is to prepare a quantum state that optimizes a
given cost function (commonly also called objective function or loss
function). A cost function can be, for example, the expectation value of a
hermitian operator such as an Hamiltonian [2], or the distance between
generated and target distributions of measurement outcomes [3]. Unfortunately,
it has become evident that the highly non-convex cost landscapes are hard to
reliably navigate– especially with noise incurred from finite measurement
statistics. The particular shape of each cost landscape and its interaction
with the chosen optimization algorithm are key factors that determine the
reliability and efficiency of a VQA [4]. Numerous studies (such as Refs. [5,
6, 7]) that go beyond random parameter initialization and optimization naive
gradient descent were published in the past years to cleverly explore the
parametrized cost landscapes, but a universal solution has, perhaps
unsurprisingly, not yet been found. We argue that it is essential to
understand and appreciate each particular family of cost landscapes defined by
problem instances of interest in order to select an informed parameter
initialization strategy and an appropriate optimization algorithm.

The Quantum Approximate Optimization Algorithm (QAOA) [8] is a type of VQA,
which is typically used to solve combinatorial optimization problems. This is
done by approximately finding the ground state of a diagonal Ising Hamiltonian
which encodes the problem [9]. In QAOA, the cost function is usually the
expectation value of the problem Hamiltonian on a parametrized state, i.e.,
its energy. The circuit ansatz for the quantum state is highly structured and
problem dependent– a property that QAOA does not share with many types of
VQAs. As such, there are significantly fewer design choices that go into
utilizing QAOA compared to more generic VQA implementations. One of the
remaining being the choice of the optimizer for the particular problem
instance.

While working on analyzing the Variational Quantum Factoring algorithm (VQF)
[10], we became aware of the following pair of Ising Hamiltonians, which
challenged our understanding of QAOA landscapes:

$H_{1}=-2.75Z_{0}-3.25Z_{1}+3.75Z_{0}Z_{1}$ (1)

$\displaystyle H_{2}=$
$\displaystyle\frac{1}{4}\sum_{i=0}^{5}Z_{i}+\frac{3}{4}\sum_{i=\\{0,2,4\\}}Z_{i}Z_{i+1}+\frac{1}{8}\sum_{i=\\{0,2\\}}Z_{i}Z_{i+2}Z_{4}Z_{5}$
(2) $\displaystyle+$
$\displaystyle\frac{1}{8}\sum_{i=2}^{5}(-1)^{i+1}Z_{0}Z_{i}+\frac{1}{8}\sum_{i=2}^{5}(-1)^{i}Z_{1}Z_{i}+\frac{1}{8}\sum_{i=2}^{3}(-1)^{i}Z_{i}Z_{4}+\frac{1}{8}\sum_{i=2}^{3}(-1)^{i+1}Z_{i}Z_{5}$
$\displaystyle-$
$\displaystyle\frac{1}{8}\sum_{i=2}^{5}Z_{0}Z_{1}Z_{i}-\frac{1}{8}\sum_{i=0}^{3}Z_{i}Z_{4}Z_{5}-\frac{1}{8}\sum_{i=\\{0,1,4,5\\}}Z_{2}Z_{3}Z_{i},$

where $Z_{i}$ denote Pauli $Z$ operators acting on qubit $i$. These
Hamiltonians were the initial motivation for this work. One may naively assume
that the Hamiltonian with three 1- and 2-body terms has a significantly
simpler QAOA cost landscape than one with 36 up to 4-body terms. However,
after visually analyzing the cost landscapes of these Hamiltonians, we
discovered that this assumption is in fact incorrect (see Fig. 1). This led to
the desire to pin down and understand which properties of the Hamiltonian
influence the shape of its cost landscape and to what extent.

Figure 1: Cost landscape scans of the Hamiltonians $H_{1}$ (a) and $H_{2}$ (b)
with $p=1$ layer QAOA. The scans visually demonstrate that even though $H_{1}$
has an apparently simpler structure than $H_{2}$, its cost landscapes exhibit
many more local minima which are likely to lead to difficult optimization.
$H_{1}$ and $H_{2}$ are defined in Eqs. (1) and (2), respectively.

This work aims to study how the structure of a Hamiltonian influences the
resulting QAOA cost landscape, and by this contributes to our general
understanding. Our approach is based on visualization of the cost landscape,
but additionally, utilization of their Fourier transform, i.e., the Fourier
landscape, as a key tool. In order to more quantitatively evaluate the
measured landscapes, we introduce roughness metrics on the cost and Fourier
landscapes which can be correlated with how “simple” or “complicated” both
landscapes appear. We demonstrate the predictability of certain features of
the landscapes by systematically manipulating the QAOA Hamiltonians. While it
is unlikely that these or similar roughness metrics can be used in general to
predict which optimizer will be most successful in minimizing the cost, they
may be valuable to study a particular family of QAOA instances of interest.

Section II recapitulates the QAOA method, how it can be studied from the
perspective of the Fourier spectrum, and introduces the roughness metrics used
throughout this work. In Sec. III we present our findings on how the
particular form of a Hamiltonian affects the QAOA cost and Fourier landscapes.
Finally, we conclude in Sec. IV and point out possible future research
directions.

##  II Background

###  II.1 QAOA and its cost function

QAOA is a type of VQA that aims to find the ground state of an Ising
Hamiltonian. The QAOA ansatz, i.e., the architecture of the circuit, is
defined by the digitized time evolution under the problem Hamiltonian (cost
Hamiltonian) as well as a second mixer Hamiltonian. It consists of a total of
$p$ cost Hamiltonian layers (parametrized by parameters
$\gamma\in\mathbb{R}^{p}$) and mixer layers (parametrized by parameters
$\beta\in\mathbb{R}^{p}$), respectively, where $p$ is a user-controlled
hyperparameter which controls the expressive power of the circuit. As such,
the circuit ansatz is highly problem-specific and only sparsely parametrized
by $2p$ parameters.

In this work, we restrict ourselves to studying QAOA with $p=1$, as well as
the mixer Hamiltonian $\sum_{i=0}^{n}X_{i}$, which follows the commonly-used
original QAOA ansatz [8]. Here $n$ is the number of qubits and $X_{i}$ are
Pauli X operators. In this form, and with a 2-body problem Hamiltonian, the
cost function value $\mathcal{C}(\beta,\gamma)$, i.e., the energy, has a
closed-form expression which was defined in Ref. [11]. Such a Hamiltonian can
be interpreted as an undirected graph with its vertices $V={1,...,n}$ and
edges $E\subset V\times V$, where each vertex has an associated bias $h_{i}$
and each edge has coupling strength $J_{ij}$. The Hamiltonian consequently
consists of the sum of 1- and 2-body terms ($C_{i}$ and $C_{ij}$,
respectively)

$H=\sum_{i\in V}C_{i}+\sum_{(i,j)\in
E}C_{ij},\;\;C_{i}:=h_{i}Z_{i},\;C_{ij}:=J_{ij}Z_{i}Z_{j}\;,$ (3)

where the notation $(i,j)\in E$ denotes neighboring vertices on the graph that
are connected by an edge.

The cost function can be written as

$\mathcal{C}(\beta,\gamma):=\sum_{i\in V}\left\langle
C_{i}\right\rangle+\sum_{(i,j)\in E}\left\langle C_{ij}\right\rangle.$ (4)

The closed-form formulas for the individual cost terms have been derived in
[11] and are then given by:

$\displaystyle\left\langle C_{i}\right\rangle$
$\displaystyle=h_{i}\sin(2\beta)\sin\left(2\gamma h_{i}\right)\prod_{(i,k)\in
E}\cos\left(2\gamma J_{ik}\right)$ (5) $\displaystyle\left\langle
C_{ij}\right\rangle$
$\displaystyle=\frac{J_{ij}\sin(4\beta)}{2}\sin\left(2\gamma
J_{ij}\right)\biggl{[}\cos\left(2\gamma
h_{i}\right)\prod_{\begin{subarray}{c}(i,k)\in E\\\ k\neq
j\end{subarray}}\cos\left(2\gamma J_{ik}\right)+\cos\left(2\gamma
h_{j}\right)\prod_{\begin{subarray}{c}(j,k)\in E\\\ k\neq
i\end{subarray}}\cos\left(2\gamma J_{jk}\right)\biggr{]}$ (6)
$\displaystyle-\frac{J_{ij}}{2}(\sin(2\beta))^{2}\prod_{\begin{subarray}{c}(i,k)\in
E\\\ (j,k)\notin E\end{subarray}}\cos\left(2\gamma
J_{ik}\right)\prod_{\begin{subarray}{c}(j,k)\in E\\\ (i,k)\notin
E\end{subarray}}\cos\left(2\gamma J_{jk}\right)$
$\displaystyle\times\biggl{[}\cos\left(2\gamma\left(h_{i}+h_{j}\right)\right)\prod_{\begin{subarray}{c}(i,k)\in
E\\\ (j,k)\in
E\end{subarray}}\cos\left(2\gamma\left(J_{ik}+J_{jk}\right)\right)-\cos\left(2\gamma\left(h_{i}-h_{j}\right)\right)\prod_{\begin{subarray}{c}(i,k)\in
E\\\ (j,k)\in
E\end{subarray}}\cos\left(2\gamma\left(J_{ik}-J_{jk}\right)\right)\biggr{]}$

By inspecting Eqs. (5) and (6), one can notice that all terms are
trigonometric functions whose frequencies depend on the parameters $h_{i}$ and
$J_{ij}$ of the problem Hamiltonian. Thus one can immediately derive
properties of the landscape such as its maximum frequency and parameter
symmetries like parity. Unfortunately, it applies only to Hamiltonians with 1-
and 2-body terms, as closed-form formulas for higher-order terms are not
known. Furthermore, as we elaborate in the Sec. II.2, the trigonometric form
of QAOA cost functions justifies the use of Fourier analysis to characterize
their landscape.

###  II.2 Fourier analysis of the landscape

The Fourier transform is the quintessential classical signal processing tool.
It can be used to isolate the sinusoidal components of any given function. In
the particular case of VQAs, the cost function is evaluated on a parameterized
quantum circuit (PQC): $\mathcal{C}(\theta)=\langle
0|U^{{\dagger}}(\theta)HU(\theta)|0\rangle$, where $H$ is the cost
Hamiltonian, and $\theta=(\theta_{1},...,\theta_{N})$ the parameters that
control the variable gates in the quantum circuit. Typically, these gates are
generated by Hermitian operators. It has been demonstrated multiple times
across the literature that the cost functions of quantum algorithms can be
described as Fourier series (see e.g., Refs. [12, 13, 14, 15, 16, 17, 18, 19,
20]). In order words, we can write

$\mathcal{C}(\theta)=\sum_{\omega\in\Omega}c_{\omega}e^{i\omega\cdot\theta}$
(7)

where the coefficients $c_{\omega}$ above can be recovered by performing a
$N$-dimensional Fourier transform

$c_{\omega}=\dfrac{1}{(\sqrt{2\pi})^{N}}\underset{[0,2\pi]^{N}}{\int}e^{-i\omega\cdot\theta}\mathcal{C}(\theta)\,d^{N}\theta\,,$
(8)

where $\omega\cdot\theta$ is an inner product. The frequencies
$\omega\in\Omega\subset\mathbb{R}^{N}$ depend only on the Hamiltonians that
generate the parameterized unitary operations. Suppose we have a single gate
$e^{-i\theta M}$ generated by a Hermitian operator $M$, $N=1$. Then the
corresponding frequencies of the cost function obey

$\omega\in\\{\lambda_{i}(M)-\lambda_{j}(M)\\}_{i,j}\,,$ (9)

meaning all possible differences between pairs of eigenvalues $\lambda_{i}(M)$
of $M$ [12]. Conversely, the values of the coefficients
$c_{\omega}\in\mathbb{C}$ depend on the circuit, the input state and the cost
Hamiltonian $H$ [13, 12, 19], and in general are hard to estimate classically.
In the case of $p=1$ QAOA, $N=2$, the parameters are $\theta=(\beta,\gamma)$
and $\omega=(f_{\beta},f_{\gamma})$ are the corresponding frequencies. Since
the unitary parameterized by $\gamma$ is generated by the problem Hamiltonian
$H$, the frequencies $f_{\gamma}$ will depend on the differences between the
eigenvalues of $H$.

As outlined in Refs. [19, 20], the Fourier coefficients contain information
about general features of the cost landscape. Most relevant for our
discussion, in Ref. [20] it was shown that the number of nonzero Fourier
coefficients, also known as the Fourier sparsity, is an important metric for
optimization. For instance, it was shown that narrow gorges (an obstacle to
optimization) imply many nonzero coefficients. Roughly speaking, this occurs
because a narrow gorge is a concentration of the cost function in parameter
space, and by the uncertainty principle of Fourier analysis this implies a
frequency spectrum that is spread out. Conversely, the same paper showed that
a very sparse spectrum means that the cost function be efficiently
reconstructed classically via compressed sensing, thus making its optimization
much simpler.

In summary, Fourier analysis provides us with a powerful and suitable tool for
characterizing the cost landscape of QAOA instances.

###  II.3 Roughness metrics

Humans have an intuitive understanding of what is meant by a “rough” cost
function landscape. Upon comparing the two cost landscapes in Fig. 1, one may
conclude that the landscape in panel (a) is “rougher” than the one in panel
(b). This conclusion, however, is subjective and not quantifiable.
Furthermore, simple inspection becomes problematic for cost landscapes in
higher dimensions (i.e., with more parameters), which are significantly harder
to inspect and prone to the existence of artifacts due to the choice of
visualization technique.

One potential precise definition of roughness could be the density of local
optima (minima or maxima depending on the algorithm) and saddle points in a
given landscape. If the landscape presents numerous suboptimal local optima,
not only will it appear rough to a human observer, but optimization algorithms
will likely struggle to find good solutions. Both gradient-based and gradient-
free optimizers have the tendency to remain trapped in local optima instead of
reaching global optima, but the exact behavior depends on the particular
optimization algorithm of choice. Therefore, we seek objective metrics that
are able to quantify roughness, which can be fed back into the decision which
optimizer to use. To this end, various metrics were explored. Here we present
two selected roughness metrics that we have found to be useful, which
correlate well and yet are complimentary in certain aspects. The advantages
and disadvantages of the other metrics we considered are presented in Appendix
C.

####  II.3.1 Metric 1: Total variation

Total variation (TV) has previously been used to estimate roughness in neural
network loss landscapes [21]. Its name refers to the fact that it captures how
much a function varies across a finite domain. It can be calculated by
integrating the absolute value of the gradient:

$TV(f)=\int_{a}^{b}\lvert f^{\prime}(x)\rvert\,dx\ $ (10)

If a function oscillates more over a fixed domain, then it will have a steeper
average slope. TV thus represents a good measure for the oscillatory behavior
of a function, which makes it a good candidate for a roughness index. In our
implementation we normalize the expression in Eq. (10) maximum range of
function values

$[f]=\max\limits_{a\leq x\leq b}f(x)-\min\limits_{a\leq x\leq b}f(x)\,.$ (11)

This gives a normalized TV

$TV_{\text{norm}}(f)=\frac{1}{[f]}\int_{a}^{b}\lvert
f^{\prime}(x)\rvert\,dx\,,$ (12)

which is invariant to scaling of the cost function ($f\rightarrow\alpha f$),
$\alpha\in\mathbb{R}$. Because one does in practice not have access to a
closed-form expression for the function gradient, we use a first-order
numerical integration approximation of the TV which utilizes a first-order
numerical gradient approximation on a grid of parameters. The TV then takes
the simple form

$TV_{\text{discrete}}(f)=\frac{1}{[f]}\sum_{i=1}^{m}\left|f\left(a+(i+1)r\right)-f\left(a+ir\right)\right|\,,$
(13)

where $r=\frac{b-a}{m}$ and $m$ is the number of steps in the discrete
integral.

To generalize to multivariate functions, we take 1D slices of the landscape by
picking random vectors originating from a common point and sample the
landscape along those directions. The normalized TV is then calculated for
each section with the mean value representing the final roughness index. We
measure along 200 random directions with $m=200$ samples along each of these
directions. This value provided a good compromise between computation time and
the reliability of the results (for more information, see Appendix C). For
each section, the size of the sampling domain $[a,b]$ is set such that it
covers exactly one period of the cost function.

####  II.3.2 Metric 2: Fourier density

As discussed in Sec. II.2, the frequency modes revealed by the Fourier
transform offers a natural perspective to study quantum landscapes. In
particular, in Ref. [20] the Fourier sparsity defined as the number of nonzero
Fourier coefficients was shown to be an important metric for the purpose of
optimization. Mathematically, if $\textbf{c}_{\omega}$ is the vector of
Fourier coefficients, the Fourier sparsity is defined as its 0-norm
$\|\textbf{c}_{\omega}\|_{0}$.

In this work, we adapt Fourier sparsity to become a metric of landscape
roughness. We will refer to the metric by the name Fourier density (FD) to
avoid the potential confusion in interpreting phrases such as “large Fourier
sparsity” (which should describe that Fourier coefficients are not sparse).
However, this metric as defined above through the 0-norm poses a challenge for
numerical implementations such as ours, as the finite precision of floating
point calculations may lead to an artificially large value. One potential
solution is to set a threshold, below which a coefficient is considered to be
zero and is not counted. This, however, introduces an arbitrary free
parameter.

In Ref. [22], the numerical sparsity

$s(\textbf{c}_{\omega}):=\|\textbf{c}_{\omega}\|_{1}^{2}/\|\textbf{c}_{\omega}\|_{2}^{2}$
(14)

was proposed as a solution to the numerical instability of the 0-norm, where
$\|\textbf{c}_{\omega}\|_{1}^{2}$ and $\|\textbf{c}_{\omega}\|_{2}^{2}$ are
the 1- and 2-norm, respectively. $s(\textbf{c}_{\omega})$ is in fact a sharp
lower bound on $\|\textbf{c}_{\omega}\|_{0}$ with equality if and only if all
non-zero coefficients are of equal magnitude. Because of its lack of arbitrary
threshold and its numerical stability, we utilize the expression in Eq. (14)
as the second roughness metric and we refer to it as Fourier Density. Other
variants of this metric that we considered are discussed in Appendix C.

##  III Results

For visualization purposes, we use the orqviz package [23], which is a Python
library for visualizing the cost function landscapes in parametrized
algorithms. For details on the visualization techniques see Appendix A.

###  III.1 Toy problem

We begin with the study of a family of 2-qubit Hamiltonians,

$H(a,b,c)=aZ_{0}+bZ_{1}+cZ_{0}Z_{1}\,.$ (15)

The limited number of degrees of freedom allows for a closed-form analytical
expression of the cost function, which in turn can be compared to the
empirically measured cost and Fourier landscapes. The cost function for this
Hamiltonian can be derived from Eqs. (5) and (6) and reads

$\displaystyle\mathcal{C}(\beta,\gamma)$ $\displaystyle=\langle
C_{0}\rangle+\langle C_{1}\rangle+\langle C_{01}\rangle$ (16)
$\displaystyle=a\sin(2\beta)\sin(2a\gamma)\cos(2c\gamma)+b\sin(2\beta)\sin(2b\gamma)\cos(2c\gamma)$
$\displaystyle\hskip
14.22636pt+\frac{c}{2}\left[\sin(4\beta)\sin(2c\gamma)(\cos(2a\gamma)+\cos(2b\gamma))-\sin^{2}(2\beta)(\cos(2(a+b)\gamma)-\cos(2(a-b)\gamma))\right].$

Upon first examination, we see that each coefficient $a$, $b$ and $c$ affects
both the amplitudes and the frequencies of all partaking terms. The existence
of the two-qubit term in the Hamiltonian gives rise to cosine terms with
frequencies that are proportional to either the sum or the difference of the
coefficients $a$ and $b$. This already hints at the slew of frequencies that
contribute to larger problems with more coefficients. This form of the cost
function can be directly linked to the Fourier spectrum of $\mathcal{C}$, as
$\beta$ and $\gamma$ are present only as arguments of the trigonometric
functions. As described in Sec. II.2, the Fourier transform allows us to
define $f_{\beta}$ and $f_{\gamma}$ for every cost term. The formulation in
Eq. (16), however, allows us to observe that all the $\sin$ terms with $\beta$
in the argument have frequencies that do not depend on the values of the
coefficients. This implies that $f_{\beta}$ will be independent of the
particular choice of Hamiltonian coefficients.

To further simplify the problem, we first assign $a=b=1$ and calculate
expected frequencies for various values of $c$. Detailed calculations can be
found in Appendix B, the expected outcomes are presented in Table 1, and the
empirically observed cost and Fourier landscapes are depicted in Fig. 2.
Because there are terms that depend on the value of $c$ and others that do
not, as $c$ is increased some frequencies shift up, whereas others do not
change. This can also be observed in the cost landscape scans, where some
underlying features of the landscape stay consistent, while higher frequency
contributions are superimposed and differ in each case. Visually, one would
say that the cost landscapes with higher-frequency contributions are rougher–
a trend that is confirmed by the roughness metrics.

Note that in Table 1 we actually list the absolute value of the frequencies.
In reality, the Fourier transform, as defined in Eq. 8, will output both
positive and negative-valued frequencies. This is because each trigonometric
term gives rise to an equal number of positive and negative frequency terms
with the same absolute frequency values, e.g.
$\cos(4\beta)=\frac{1}{2}(e^{i4\beta}+e^{-i4\beta})$. Note that, while the
frequencies $\omega$ are the same in absolute value, the coefficients
$c_{\omega}$ at negative frequencies may have a completely different value
from those at positive ones. As a simple calculation shows, the cost function
being real implies that $c_{-\omega}=c^{*}_{\omega}$, but for the case of more
than one parameter (where $\omega$ is a vector), this still leaves some
freedom in the negative frequency coefficients. Indeed, take our case with 2
parameters $\theta=(\beta,\gamma)$ and frequencies
$\omega=(f_{\beta},f_{\gamma})$. The fact that the cost function is real-
valued implies that
$c_{(f_{\beta},f_{\gamma})}=c^{*}_{(-f_{\beta},-f_{\gamma})}$. In other words,
the Fourier landscape is point-symmetric around the origin. Hence we depict
the negative frequencies in the Fourier spectrum for the $f_{\beta}$
direction, only cutting the spectrum in half once to account for the symmetry
$c_{-\omega}=c^{*}_{\omega}$.

c | Frequencies from linear terms | Frequencies from quadratic term
---|---|---
| (2c$\pm$2, 2) | (2c$\pm$2, 4) (4, 4) (4, 0) (0, 4)
1 | (4, 2) | (4, 4) (4, 0) (0, 4)
5 | (8, 2) (12, 2) | (8, 4) (12, 4) (4, 4) (4, 0) (0, 4)
10 | (18, 2) (22, 2) | (18, 4) (22, 4) (4, 4) (4, 0) (0, 4)
20 | (38, 2) (42, 2) | (38, 4) (42, 4) (4, 4) (4, 0) (0, 4)
Table 1: Frequencies for Hamiltonian of form $H(c)=Z_{0}+Z_{1}+cZ_{0}Z_{1}$
(see Eq. (15)), written as $(|f_{\gamma}|,|f_{\beta}|)$. When multiple terms
with the same frequencies would be present in a particular row, their
amplitudes will be added together. This summation results in a reduced number
of points in the Fourier plot. The constant term $(0,0)$ is omitted. Appendix
B contains detailed steps on how these values have been obtained and why the
term $(0,2)$ for $c=1$ is missing.  Figure 2: Cost and Fourier landscapes
scans for the Hamiltonians ${H(c)=Z_{0}+Z_{1}+cZ_{0}Z_{1}}$, with
$c=\\{1,5,10,20\\}$. The analytically calculated frequencies are shown in
Table 1 and match the frequencies observed in the numerical experiments up to
finite-size artifacts.

We used the Hamiltonian in Eq. 15 also to study the effect of non-integer
coefficients. As explained in Appendix A.2, if the greatest common divisor
(GCD) of the coefficients is small, the period of the cost function can become
very large. For instance, if $a=0.171$, $b=0.340$, $c=0$, the GCD is 0.001 and
the period of the resulting cost function would be
$\frac{\pi}{0.001}=1000\pi$. If we round the coefficients to $a=0.17$,
$b=0.34$ and $c=0$ the period would be equal to $\frac{\pi}{0.17}$. However,
we can also note that this Hamiltonian is equivalent to one with $a=1$, $b=2$,
and $c=0$, up to a constant re-scaling factor in both frequency and amplitude,
which would further simplify the analysis. Since the GCD has been drastically
increased from 0.001 to 0.17, and consequently the period is reduced by orders
of magnitude, one might expect that the cost landscape is significantly
different. Interestingly, we show that this is not the case.

We generated multiple Hamiltonians with random coefficients, where the
coefficients were drawn from a uniform distribution ranging from -10 to 10.
For each of these Hamiltonians, a corresponding Hamiltonian was constructed
with coefficients rounded to the nearest integers. The cost landscapes of
these Hamiltonians were analyzed visually and through the use of roughness
metrics. We observe a high degree of similarity among the majority of the cost
landscapes generated between the non-integer and their corresponding rounded
integer-coefficient Hamiltonian. Naturally, the differences in the cost
landscapes are most noticeable when the rounding errors are high. A
representative example can be seen in Fig. 3. On the flip side, the Fourier
spectrum is significantly different between the pairs of Hamiltonian. By
nature of the discrete Fourier transform of the cost landscape, the
frequencies ($f_{\beta},f_{\gamma}$) are integer-valued. Therefore, if the
Hamiltonian eigenvalues are integer-valued, also their sum and difference
result in integers. This results in well-localized and distinct frequencies in
the Fourier landscape. In contrast, the original un-rounded Hamiltonian with
non-integer coefficients exhibits non-integer frequencies which need to be
distributed over adjacent integer frequencies. This is a phenomenon which is
captured well by the Fourier density metric, as it is highly sensitive to the
presence of additional frequencies in the spectrum. However, one can argue
that these are simply artifacts of the discrete Fourier transform. It follows
that the Fourier density should not be used to compare the roughness between
landscapes where one is integer-valued and the other is non-integer-valued.
Consequently, we constrain our study to Hamiltonians that do not fall into
this case.

Figure 3: Cost and Fourier landscapes for two similar Hamiltonians with
integer coefficients (a) and non-integer coefficients (b). The Hamiltonians
are ${H_{a}=1Z_{1}+4Z_{2}+6Z_{1}Z_{2}}$ and
${H_{b}=1.23Z_{1}+4.47Z_{2}+6.05Z_{1}Z_{2}}$, respectively. While the cost
landscapes appear to be very similar, the Fourier spectra differ more
strongly. This is a result of the fact that in discrete Fourier
transformation, non-integer frequencies do not result in well-localized peaks,
but rather in distributions of frequencies.

###  III.2 Higher k-body terms

Here we analyze the effect of the presence of $k$-body terms in the
Hamiltonian on the cost and Fourier landscapes in QAOA. Hamiltonians with
$k\geq 3$ body terms are rarely studied throughout the literature. One
apparent explanation is that such Hamiltonians do not natively map to problems
represented as graphs, which are typical problem instances for the QAOA
algorithm. However, we believe it is vital to understand the effect that such
terms have on the optimization landscape, as they are present in certain non-
graph-related problem instances (e.g., the VQF algorithm [10]).

To construct Hamiltonians for our experiments, we used as a base a 6-qubit
Hamiltonian consisting of all $1$-body Pauli $Z$ terms:

$\mathcal{H}_{1}=\sum_{i=1}^{6}Z_{i}\ .$ (17)

To control for the number of terms, we add one $k^{*}=2,3,4,5$ body term to
the Hamiltonian,

$\displaystyle\mathcal{H}_{2}=\mathcal{H}_{1}+Z_{1}Z_{2},\quad\quad\quad\;\mathcal{H}_{3}=\mathcal{H}_{1}+Z_{1}Z_{2}Z_{3}$
(18)
$\displaystyle\mathcal{H}_{4}=\mathcal{H}_{1}+Z_{1}Z_{2}Z_{3}Z_{4},\quad\mathcal{H}_{5}=\mathcal{H}_{1}+Z_{1}Z_{2}Z_{3}Z_{4}Z_{5}\
,$

and study the resulting landscapes. The results can be viewed in Fig. 4. The
cost landscapes clearly exhibit contributions of higher frequencies that are
generated by the higher-order terms, and which lead to the emergence of
additional local minima. The TV metric mirrors the observed trend and
increases steadily with $k^{*}$, while the Fourier Density decreases in the
last case with $k^{*}=5$. This highlights that the latter metric does not
necessarily correlate with the density of local minima in the cost landscape.

As it can be read off in the Fourier landscapes, higher order terms increase
the number of terms in the Fourier spectrum, as well as their possible
$f_{\beta}$ frequencies. The highest possible frequency in this dimension is
equal to twice the order of the highest term, i.e. $\pm 2\times k^{*}$. This
outcome is consistent with our expectations based on the analysis of equations
Eqs. (5) and (6). However, as we can see in the Fourier landscapes, replacing
the particular term also affects the relative weights of other Fourier
coefficients.

Figure 4: Cost and Fourier landscapes for the Hamiltonians $\mathcal{H}_{2}$
(a), $\mathcal{H}_{3}$ (b), $\mathcal{H}_{4}$ (c) and $\mathcal{H}_{5}$ (d) in
Eq. (18). We notice that as higher-order terms are introduced, Total Variation
tends to increase. The same clear trend is not reflected by the Fourier
Density. We can also see the highest possible frequency $f_{\beta}$ increase
as $\pm 2\times k^{*}$.

###  III.3 One large coefficient

In Sec. III.1, we show that the magnitude of Hamiltonian coefficients can
drastically impact the roughness of the cost landscape– an effect that we have
observed repeatedly throughout our studies. On the other hand, we show in Sec.
III.2 that introducing a single high-order term also tends to increase
roughness. A natural question to ask is which of the two factors has the
stronger impact? To study this, we use the Hamiltonian

$\mathcal{H}_{6}=\sum_{i=1}^{6}Z_{i}+\sum_{i<j}^{6}Z_{i}Z_{j}+\sum_{i<j<k}^{6}Z_{i}Z_{j}Z_{k}$
(19)

as the base and create three new Hamiltonians from it. For each we pick the
first term of the order of $k^{*}=1$ (the $Z_{1}$ term), $k^{*}=2$ (the
$Z_{1}Z_{2}$ term) or $k^{*}=3$ (the $Z_{1}Z_{2}Z_{3}$ term), and increase its
coefficient to 25. We label them as $\mathcal{H}_{7}$, $\mathcal{H}_{8}$ and
$\mathcal{H}_{9}$, respectively. Their cost and Fourier landscapes are shown
in Fig. 5.

A first observation is that the cost landscapes in fact become significantly
rougher when a single large coefficient is introduced. Furthermore, we observe
the addition of new dominant frequency at $(f_{\gamma},f_{\beta})=(\pm 2\times
25,2\times k^{*})$. However, while TV and FD correlate well with each other,
in this case we see that roughness metrics do not consistently increase with
$k^{*}$. Overall this clearly suggests that the effect of a single (or a few)
outlier coefficients can completely warp the cost landscape, but the interplay
with the order of the term is unclear. The coefficient manipulates
$f_{\gamma}$ frequencies and the order of the term influences $f_{\beta}$.
Their interference is what finally appears to determine the exact shape of the
cost and Fourier landscapes.

Figure 5: Cost and Fourier landscapes for the Hamiltonians $\mathcal{H}_{6}$
(a), $\mathcal{H}_{7}$ (b), $\mathcal{H}_{8}$ (c) and $\mathcal{H}_{9}$ (d)
which are a variation of the Hamiltonian in Eq. 19. For the latter three
Hamiltonians, one $k$-body term, i.e., $k^{*}=1,2,3\,,$ respectively, had its
coefficient increased to $25$.

###  III.4 Structure of the Hamiltonian

In order to better understand the influence of Hamiltonian structure on the
cost landscape, we analyze cases where we incrementally change Hamiltonians by
adding one term at a time. More specifically, we interpolate between two types
of Hamiltonians and study where the landscapes change most in terms of
roughness. A priori it is not evident whether the largest changes would happen
when the structure of the Hamiltonian is broken, i.e. when the first different
terms are added, or during the middle part of the interpolation when there is
arguably the least structure.

The Hamiltonians studied here form the sequence

$\displaystyle\mathcal{H}_{a}=\sum_{i=1}^{6}Z_{i}$ (20)
$\displaystyle\mathcal{H}_{a,1}=\mathcal{H}_{a}+Z_{0}Z_{1}Z_{2}$
$\displaystyle\mathcal{H}_{a,2}=\mathcal{H}_{a}+Z_{0}Z_{1}Z_{2}+Z_{0}Z_{1}Z_{3}$
$\displaystyle\quad\quad\quad\vdots$
$\displaystyle\mathcal{H}_{a,20}=\mathcal{H}_{a}+Z_{0}Z_{1}Z_{2}+Z_{0}Z_{1}Z_{3}+\cdots+Z_{3}Z_{4}Z_{5}=\mathcal{H}_{b}\
,$

and thus represent the interpolation between $\mathcal{H}_{a}$, the
Hamiltonian containing all $k=1$ terms, and $\mathcal{H}_{b}$, the Hamiltonian
containing all $k=1$ and $k=3$ terms.

Panel a) and b) of Fig. 6 show the TV and FD values along the interpolation
path where new terms are added one by one. These results imply that adding a
single term can have both minor or significant effects on the landscapes
depending on the current Hamiltonian. In particular, changes in roughness are
largest and most erratic in the intermediate region between $\mathcal{H}_{a}$
and $\mathcal{H}_{b}$ where the landscapes also tend to be rougher overall.
The fact that roughness is higher for the Hamiltonian $\mathcal{H}_{b}$, which
contains all $k=3$ terms, is consistent with our observations in Sec. III.2
that Hamiltonians containing higher-order terms tend to exhibit rougher
landscapes. In panels c) - f) we present representative cost and Fourier
landscape examples for the smallest and largest changes in roughness metrics.

Our findings suggest that predominant structures in the Hamiltonian (or when
applicable, its graph representation) stabilize the cost landscape. In this
regime, individual terms do not appear to contribute significantly, and it
requires a significant shift in structure for the landscape to qualitatively
change. We note that, as the interpolation can be understood in both
directions, adding or removing individual terms has equivalent effects.

Figure 6: Recorded roughness metrics (a, b) while interpolating between the
Hamiltonians $\mathcal{H}_{a}$ and $\mathcal{H}_{b}$ by adding one term at a
time (see Eq. (20)). We additionally show two pairs of cost Hamiltonians with
small and big differences in TV. Panels c) and e) show $\mathcal{H}_{a,1}$ and
$\mathcal{H}_{a,2}$, whereas panels d) and f) show $\mathcal{H}_{a,15}$ and
$\mathcal{H}_{a,16}$. The results demonstrate that adding or removing
individual terms does not have a drastic influence on the landscapes when the
dominant structure is mostly preserved. In contrast, large changes in the
landscapes occur when there is not one particular dominant structure.

###  III.5 Parameter concentration and barren plateaus

The phenomenon in QAOA that is most known for applying Hamiltonians with
particular structures is that of parameter concentration. It was shown in Ref.
[24] that, for MaxCut problems on random 3-regular graphs of the same size and
across graph sizes, local minima and maxima are located at similar parameter
values. Among other works, this finding was underpinned in Ref. [23], where it
was shown that in those cases the entire cost landscapes are very similar. It
was also shown that the phenomenon is not just limited to the MaxCut problem
or 3-regular graphs, and in fact persists for other types of problems and
graphs, even when the coefficients are selected from a random distribution.

In this work, we analyze the phenomenon of parameter concentration from the
perspective of the cost and Fourier landscapes. For this purpose, we study
MaxCut problems on 3-regular graphs, which are graphs where each node is
connected to exactly three other nodes. For a graph $G$ that is specified by
the sets of $n$ nodes $V(G)$ and edges $E(G)$, we label a node as $v_{i}\in
V(G)$, $i=1,2,\dots n$ and the edge connecting nodes $v_{i}$ and $v_{j}$ is
labeled as $(i,j)$. Each edge can have an individual weight $w_{ij}$
associated with it. The objective of the MaxCut problem is to find a partition
of $V(G)$ into two sets $A$ and $B$, such that the sum of weights on the edges
between $A$ and $B$ is maximized. We can encode whether a given node is in set
$A$ or $B$ using a binary string $z\in\\{-1,1\\}^{\otimes n}$, where
$z_{i}=-1$ if node $i$ is in $A$ and $z_{i}=1$ if it is in $B$. The MaxCut
cost function can then be written as the sum of weighted edges crossing $A$
and $B$,

$\mathcal{C}=\sum_{(i,j)\in E(G)}\frac{1}{2}w_{ij}(1-z_{i}z_{j})\,.$ (21)

This cost function can naturally be expressed as the expectation value of a
quantum state on the Hamiltonian

$\mathcal{H}_{\text{MaxCut}}=\sum_{(i,j)\in
E(G)}\frac{1}{2}w_{ij}(\mathbb{1}-Z_{i}Z_{j})\,,$ (22)

where $\mathbb{1}$ is the identity operator and $Z_{i}Z_{j}$ are Pauli $Z$
operators acting on qubits $i$ and $j$. For our simulations, we removed the
constant offset terms from the Hamiltonian that stem from the weighted
identity operators.

The features of the cost landscapes and their dependence on the values of the
Hamiltonian coefficients are extensively discussed in Ref. [11]. In our study,
we draw random 3-regular graphs of sizes $n=8,12,16$ and $20$, and with random
uniform edge weights in the domain $w_{ij}\in[-10,10]$. Representative results
are shown in Fig. 7. In all cases, both the cost and Fourier landscape
concentrate. It may be surprising that randomly connected graphs with a
3-regular structure and random weights concentrate in this manner, but our
observations are robust across many repetitions, coefficient ranges, other
types of random graphs (e.g., Erdős–Rényi graphs), and in agreement with the
findings in Ref. [23].

In the Fourier picture, the coefficients with $f_{\beta}=0$ tend to get
smaller as the size of the system increases, while the ones with $f_{\beta}=2$
form an approximately bell-shaped distribution with progressively widening
spread. The cancellation of most of the frequencies also leads to a decrease
in the value of the TV as the size of the problem increases, while values of
Fourier Density steadily increase.

One interpretation of the cost landscape scans and the decreasing TV is the
presence of barren plateaus [25]. In other words, the variation in the cost
function is concentrated in a shrinking region of parameter space, known as
narrow gorge [26]. One notices the inverse relation between the spread of the
Fourier coefficients and the size of the gorge. This is a manifestation of the
uncertainty principle in Fourier analysis as first mentioned in Ref. [20]. It
stands to reason that the non-local connectivity of the ansatz due to the
graph structure is sufficient to cause barren plateaus in an otherwise shallow
and expressivity-restricted circuit. We have found this to hold true even for
unweighted graphs, although random coefficients appear to amplify and
accelerate the phenomenon.

We can see that values of TV decrease as we increase the number of qubits,
while values of Fourier Density increase. This is a good example of why having
several complementary metrics can be useful. In this particular case, barren
plateaus seem to be characterized by very low Total Variation values and very
large Fourier Density.

Figure 7: Cost and Fourier landscapes of Hamiltonians for the MaxCut problem
on random 3-regular graphs. The system sizes are $8$ (a), $12$ (b), $16$ (c),
and $20$ (d) qubits. As the size of the problem increases, the cost and
Fourier landscapes concentrate. We can see the onset of barren plateaus in the
cost landscapes and the formation of an approximately bell-shaped distribution
of the frequencies in the Fourier landscape.

###  III.6 Revisiting the initial problem

At the end of this work, we wanted to apply the insights gained from the
previous cases to the examples presented in the introduction. Upon inspecting
the Hamiltonians, we can make some educated guesses about the cost function
landscapes and their Fourier spectrum. For both Hamiltonians, the resulting
cost function has the period of $2\pi$ in $f_{\gamma}$, which simplifies the
comparison.

The first factor we analyze is whether high coefficients in the Hamiltonian
lead to higher frequencies. $H_{1}$ has larger coefficients than $H_{2}$, and
given that both have the same periodicity, we would expect to see points with
larger frequencies $f_{\gamma}$ in the Fourier spectrum, which in turn cause
fast oscillations in the $\gamma$ direction. This prediction is indeed
correct, as we can see in Fig. 8.

The second factor is the order of the terms. Given that $H_{2}$ has 4-body
terms, we expect to see $f_{\beta}$ up to $\pm 8$, while for $H_{1}$ they will
be up to $\pm 4$. Consequently, cost function oscillations in the $\beta$
direction should be faster for $H_{2}$. However, while the maximal visible
frequencies in Fig. 8 agree with the prediction, the dominant frequencies in
terms of their amplitudes are on the two-body frequencies, i.e. $f_{\beta}=\pm
4$. This can be explained by looking at the relatively large coefficients of
those terms in the Hamiltonian. It again appears that the magnitude of
coefficients has a strong impact, because they not only directly control the
$f_{\gamma}$ frequencies, but they indirectly affect which $f_{\beta}$
frequencies shape the cost landscape the most.

The effect of the structure is much more difficult to quantify in this case.
$H_{1}$ consists of only 3 terms (all possible for two qubits), which makes it
is hard to talk about a particular structure. On the other hand, when
analyzing the structure of $H_{2}$ we can see that it has certain
regularities. All terms have coefficients that are multiples of
$\pm\frac{1}{8}$, all the 3-body terms are constructed in a similar way, and
more. This may hint that certain frequencies will cancel, which would lead to
a less rough landscape. However, this is speculative, as our work does not
provide tools to assess any given Hamiltonian structure in-depth.

Sadly, we see that our work does not provide all the necessary tools to
predict with certainty which of any two Hamiltonians produces a rougher
landscape. It does, however, provide tools that allow us to quantify the
roughness and to better understand which features of those Hamiltonians
influence the shape of cost and Fourier landscapes. As we have shown in the
previous section, our visualization techniques and roughness metrics are very
helpful in understanding each particular contribution to the QAOA landscapes,
but the sum of all components gives rise to a complexity that we do not know
whether it can be predicted without a known reference landscape.

In the introduction, we suggest that the landscape produced by $H_{1}$ is more
difficult to optimize compared to $H_{2}$, and that roughness is correlated
with the performance of the optimizer. While quantifying this effect
extensively is beyond the scope of this work, we test our hypothesis on this
pair of Hamiltonians. We performed 100 optimizations using the L-BFGS-B
optimizer [27, 28]. The initial parameters were drawn from a uniform
distribution from [-$0.9\pi$, $0.9\pi$]. As demonstrated in Fig. 8, the
optimization algorithm exhibits a higher likelihood of becoming trapped in
local minima for $H_{1}$, which is in line with the higher yielded values for
both roughness metrics. Additionally, we witness a higher density of low-lying
local minima, which can be both a curse and a blessing. If the absolute global
minimum needs to be found, these local minima present a big challenge because
gradients will most likely point the optimizer towards one of them. On the
flip side, if a good solution suffices, which of the local minima one
converges to may matter less than in a case similar to $H_{2}$, where any sub-
optimal local minimum is significantly higher in energy.

Figure 8: Cost and Fourier landscapes for the Hamiltonians $H_{1}$ (a) and
$H_{2}$ (b) as defined in Eqs. 1 and 2, as well as histograms on final
optimization outcomes. In accordance with the measured landscapes and recorded
roughness metrics, not only is the optimizer able to find the optimal solution
more than twice as often for $H_{2}$, but the density of sub-optimal local
minima is significantly higher for $H_{1}$.

##  IV Conclusion

This work was performed with the intention of understanding the relationship
between a Hamiltonian and its QAOA cost landscape. Towards this goal, we
applied visualization methods and utilized the Fourier transform of cost
landscape scans. We also introduced roughness metrics, which quantify apparent
changes in the landscapes and have been shown to correlate with optimization
performance in classical applications [21, 29]. These tools enabled us to
examine each component in a Hamiltonian of interest and how it affects the
cost function landscape.

Among the collection of metrics proposed by us, we conclude that Total
Variation is the best option, as its outcome closely resembles a subjective
interpretation of the cost landscape roughness. On the other hand, Fourier
density exhibited hypersensitivity to some features of the cost function, such
as the presence of fractional coefficients, but also to emergent properties in
barren landscapes, to name a positive example. How correlated these metrics
are with practical optimization performance remains unknown and is subject to
future work.

The visualization of the Fourier spectrum proved to be a useful tool in our
study. It enabled a more thorough understanding of specific phenomena, such as
the influence of high-order terms in the Hamiltonian, or the effect of their
coefficients. It also simplified the process of comparing cost landscapes. We
recommend it as a valuable tool to conduct analysis on cost function
landscapes. This capability has since been incorporated into the orqviz
visualization package to make it available for use by other researchers in the
field.

Finally, which properties of the Hamiltonian most affect the roughness? We
have identified three properties of the Hamiltonian that influence the shape
of the cost landscapes significantly. First are the relative magnitudes of the
coefficients of the Hamiltonian terms. In general, larger relative magnitudes
of the coefficients lead to rougher cost landscapes. Specifically, the
presence of a single term with a substantially larger coefficient can
significantly increase the roughness due to over-emphasizing a few frequencies
in the landscape. The second influencing factor is the existence of higher-
order terms in the Hamiltonian, i.e., operators which act on many qubits.
Their presence increases the maximum present frequencies $f_{\beta}$ in the
Fourier spectrum. Depending on the interplay with the coefficients, this can
lead to a landscape with an increased number of saddle points and local
minima. However, this increase in complexity can be counteracted by the third
factor– the structure of the Hamiltonian.

Hamiltonians exhibiting certain structures seem to have lower roughness scores
than those in which a structure has been disrupted by introducing or removing
terms. Additionally, Hamiltonians with structures (which can potentially be
traced back to the structure of a graph problem) appear robust to small
perturbation. This phenomenon is clearly exemplified by the phenomenon of
parameter concentration. This interestingly opens up the opportunity for
studying entire families of practical QAOA Hamiltonians if they share a
fundamental underlying structure. Given that the family of Hamiltonians
studied in Sec. III.5 are random 3-regular graphs, it is unclear how prevalent
barren plateaus will be in more realistic cases of structured problems.

Factors that we have not found to play a significant role in increasing the
roughness of the landscape are the number of qubits or the number of terms in
the Hamiltonian. The latter is oftentimes linked to the former, but an
increase in the number of qubits can, in fact, be shown to not only result in
concentrating landscapes but also barren plateaus. They diminish the magnitude
of smaller features and subsequently decrease the roughness of the cost
function landscapes, irrespective of the precise number of qubits or
individual terms.

It is important to note that there are limitations to the results presented in
this study. First of all, we performed our study for QAOA with one layer. This
limitation may affect the applicability of our findings to more practical
scenarios. Visualizing and analyzing multidimensional cost functions presents
a greater degree of difficulty. Fortunately, the orqviz visualization library
and the roughness metrics introduced in this work are designed for high-
dimensional cost functions.

Furthermore, we did not thoroughly study the relationship between roughness
metrics and the performance of various optimizers, nor did we present a study
on predicting the roughness of cost landscape for a given Hamiltonian.
Finally, it appears valuable to perform a similar analysis on other
variational algorithms, in which the Hamiltonian does not dictate the
structure of the circuit. In such cases, the choice of the circuit ansatz
plays a crucial role the shape of the optimization landscape. Notably, the
approaches used in this work appear to be good candidates to study phenomena
related to barren plateaus.

We believe that the development of tools to aid in “opening up the black box”
that variational quantum algorithms currently are will facilitate the creation
of more efficient algorithms in the future. There is significant value in the
ability to examine a given problem from multiple perspectives, and the
utilization of spectral analysis and metrics offers an alternate viewpoint to
the most widely-used methods. This approach allowed for the identification of
new insights, which to the best of our knowledge, have not yet been reported
in the literature.

###### Acknowledgements.

The authors would like to thank Pierre-Luc Dallaire-Demers for suggesting to
use Fourier transformation to analyze the landscapes of QAOA cost function to
investigate the VQF Hamiltonians and to Vladimir Vargas-Calderón for helpful
comments on the manuscript. The authors would also like to acknowledge the
ORQUESTRA® platform by Zapata Computing Inc. that was used for collecting the
data presented in this work.

## References

  * [1] John Preskill.  Quantum Computing in the NISQ era and Beyond.  Quantum, 2:79, 2018.
  * [2] 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):4213, 2014.
  * [3] Marcello Benedetti, Delfina Garcia-Pintos, Oscar Perdomo, Vicente Leyton-Ortega, Yunseong Nam, and Alejandro Perdomo-Ortiz.  A Generative Modeling Approach for Benchmarking and Training Shallow Quantum Circuits.  npj Quantum Information, 5(1):45, 2019.
  * [4] 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, 2020.
  * [5] Sukin Sim, Jonathan Romero, Jérôme F Gonthier, and Alexander A Kunitsa.  Adaptive Pruning-based Optimization of Parameterized Quantum Circuits.  Quantum Science and Technology, 6(2):025019, 2021.
  * [6] Max Wilson, Sam Stromswold, Filip Wudarski, Stuart Hadfield, Norm M. Tubman, and Eleanor Rieffel.  Optimizing Quantum Heuristics With Meta-learning.  Quantum Machine Intelligence, 2019.
  * [7] Jonas M. Kübler, Andrew Arrasmith, Lukasz Cincio, and Patrick J. Coles.  An Adaptive Optimizer for Measurement-Frugal Variational Algorithms.  Quantum, 4:263, 2020.
  * [8] Edward Farhi, Jeffrey Goldstone, and Sam Gutmann.  A Quantum Approximate Optimization Algorithm.  ArXiv, arXiv:1411.4028, 2014.
  * [9] Andrew Lucas.  Ising Formulations of Many NP Problems.  Frontiers in Physics, 2, 2014.
  * [10] Alán Aspuru-Guzik Eric Anschuetz, Jonathan Olson and Yudong Cao.  Variational Quantum Factoring.  ArXiv, arXiv/1808.08927, 2019.
  * [11] Asier Ozaeta, Wim van Dam, and Peter L McMahon.  Expectation Values From The Single-Layer Quantum Approximate Optimization Algorithm on Ising Problems.  Quantum Science and Technology, 7(4):045036, 2022.
  * [12] Maria Schuld, Ryan Sweke, and Johannes Jakob Meyer.  Effect of Data Encoding on The Expressive Power of Variational Quantum-Machine-Learning Models.  Physical Review A, 103(3):032430, 2021.
  * [13] Javier Gil Vidal and Dirk Oliver Theis.  Calculus on Parameterized Quantum Circuits.  ArXiv, arXiv/1812.06323, 2018.
  * [14] Francisco Javier Gil Vidal and Dirk Oliver Theis.  Input Redundancy For Parameterized Quantum Circuits.  Frontiers in Physics, 8:297, 2020.
  * [15] Ken M Nakanishi, Keisuke Fujii, and Synge Todo.  Sequential Minimal Optimization For Quantum-Classical Hybrid Algorithms.  Physical Review Research, 2(4):043158, 2020.
  * [16] Robert M Parrish, Joseph T Iosue, Asier Ozaeta, and Peter L McMahon.  A Jacobi Diagonalization and Anderson Acceleration Algorithm for Variational Quantum Algorithm Parameter Optimization.  ArXiv, arXiv/1904.03206, 2019.
  * [17] Mateusz Ostaszewski, Edward Grant, and Marcello Benedetti.  Structure Optimization For Parameterized Quantum Circuits.  Quantum, 5:391, 2021.
  * [18] David Wierichs, Josh Izaac, Cody Wang, and Cedric Yen-Yu Lin.  General Parameter-Shift Rules For Quantum Gradients.  Quantum, 6:677, 2022.
  * [19] Enrico Fontana, Ivan Rungger, Ross Duncan, and Cristina Cîrstoiu.  Spectral Analysis For Noise Diagnostics and Filter-Based Digital Error Mitigation.  ArXiv, arXiv/2206.08811, 2022.
  * [20] Enrico Fontana, Ivan Rungger, Ross Duncan, and Cristina Cîrstoiu.  Efficient Recovery of Variational quantum Algorithms Landscapes Using Classical Signal Processing.  ArXiv, arXiv/2208.05958, 2022.
  * [21] Keke Wu, Rui Du, Jingrun Chen, and Xiang Zhou.  Understanding Loss Landscapes of Neural Network Models in Solving Partial Differential Equations.  ArXiv, arXiv/2103.11069, 2021.
  * [22] Miles E. Lopes.  Estimating Unknown Sparsity in Compressed Sensing.  ArXiv, arXiv/1204.4227, 2012.
  * [23] Manuel S. Rudolph, Sukin Sim, Asad Raza, Michał Stęchły, Jarrod R. McClean, Eric R. Anschuetz, Luis Serrano, and Alejandro Perdomo-Ortiz.  Orqviz: Visualizing High-Dimensional Landscapes in Variational Quantum Algorithms.  arXiv, arXiv/2111.04695, 2021.
  * [24] Fernando G. S. L. Brandao, Michael Broughton, Edward Farhi, Sam Gutmann, and Hartmut Neven.  For Fixed Control Parameters The Quantum Approximate Optimization Algorithm’s Objective Function Value Concentrates for Typical Instances.  ArXiv, arXiv/1812.04170, 2018.
  * [25] 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):4812, 2018.
  * [26] Marco 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):1–12, 2021.
  * [27] Pauli Virtanen, Ralf Gommers, Travis E. Oliphant, Matt Haberland, Tyler Reddy, David Cournapeau, Evgeni Burovski, Pearu Peterson, Warren Weckesser, Jonathan Bright, Stéfan J. van der Walt, Matthew Brett, Joshua Wilson, K. Jarrod Millman, Nikolay Mayorov, Andrew R. J. Nelson, Eric Jones, Robert Kern, Eric Larson, C J Carey, İlhan Polat, Yu Feng, Eric W. Moore, Jake VanderPlas, Denis Laxalde, Josef Perktold, Robert Cimrman, Ian Henriksen, E. A. Quintero, Charles R. Harris, Anne M. Archibald, Antônio H. Ribeiro, Fabian Pedregosa, Paul van Mulbregt, and SciPy 1.0 Contributors.  SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python.  Nature Methods, 17:261–272, 2020.
  * [28] Richard H. Byrd, Peihuang Lu, Jorge Nocedal, and Ciyou Zhu.  A limited memory algorithm for bound constrained optimization.  SIAM Journal on Scientific Computing, 16(5):1190–1208, 1995.
  * [29] Matteo Aldeghi, David E. Graff, Nathan Frey, Joseph A. Morrone, Edward O. Pyzer-Knapp, Kirk E. Jordan, and Connor W. Coley.  Roughness of molecular property landscapes and its impact on modellability.  Journal of Chemical Information and Modeling, 62(19):4660–4671, 2022.  PMID: 36112568.
  * [30] Rafat N Siddiqi.  The Order of Fourier Coefficients of Function of Higher Variation.  Proceedings of the Japan Academy, 48(8):569–572, 1972.

##  Appendix A Visualization methods

For visualization purposes, we used the orqviz package [23], which is a Python
library for visualizing the cost function landscapes of parametrized
algorithms. While orqviz offers more sophisticated visualization techniques,
we utilize it only to easily perform 2D scans of the QAOA cost landscapes,
which can then be Fourier-transformed into their corresponding Fourier
landscape. The latter is a functionality which we added to orqviz during this
project.

###  A.1 Cost Landscape

To visualize the cost landscape, we perform 2D grid scans of the cost function
with a $201\times 201$ resolution, centered around the parameter values
$(\beta,\gamma)=(0,0)$ (except for Fig. 5, where we used $401\times 401$. The
basis vectors for the 2D grid can in principle be chosen to be any two
(preferably orthogonal) parameter vectors. However, since we restrict
ourselves to $p=1$ QAOA with two free parameters, the vectors are chosen in
the $\gamma$ and $\beta$ directions. We decided to use $\pi$ as the extent for
both axes for most of the figures presented here, which simplifies the
analysis and comparison.

###  A.2 Fourier Landscape

To visualize the Fourier space, we perform a discrete fast Fourier transform
(FFT) on a 2D scan of the cost function landscape. To illustrate how to
interpret the output of the Fourier transform, we will use the following
simple function: $\mathcal{C}(\beta,\gamma)=\cos(x\gamma)\cos(y\beta)$. When
we apply the Fourier transform using eq. 8, it produces a point with
coordinates $\omega=(f_{\gamma},f_{\beta})=(x,y)$ (as well as points at the
corresponding negative frequencies). Therefore the frequencies associated with
$\gamma$ and $\beta$ are being equal to $f_{\gamma}=x$ and $f_{\beta}=y$
accordingly (see Sec. II.2).

In our experiments, we plot only the absolute magnitude of the discrete
Fourier transform’s complex output. The magnitude of a coefficient determines
the amplitude of its corresponding frequency component and its influence on
the landscape. On the other hand, the phase dictates the location of crests
and troughs of frequencies relative to one another. For instance, $\sin(x)$
and $\cos(x)$, which are just shifted versions of each other, both have
coefficients of magnitude $1/2$ but with different phases. The analysis of
various plots has confirmed that the phases do not provide any significant
insight, so we decided not to include this in the data presented here.

As described in Sec. III.1, we decided to plot only half of the Fourier
spectrum, by cutting out the two quadrants corresponding to $f_{\gamma}<0$. To
summarise, the reason is that half of the spectrum is redundant due to the
cost function being real, which imposes a twofold symmetry in the coefficients
$c_{(f_{\beta},f_{\gamma})}=c^{*}_{(-f_{\beta},-f\gamma)}$. Finally, to ensure
that plots are easy to compare, we remove the constant offset from the
visualization, setting $c_{(0,0)}=0$.

####  A.2.1 Choosing period and maximum frequency

Before producing a Fourier plot, one must undertake some necessary
considerations. Firstly, it is important to ensure that the 2D scan grid
resolution is high enough when the intention is to perform a FT of the
landscape. Insufficient scan resolution causes spurious frequencies to appear
in the Fourier spectrum, a phenomenon known as “aliasing”. This behaviour is
very noticeable in lower resolution scans (linear resolution $<100$) and will
additionally inflate the value of the Fourier density roughness metric. The
nonzero artifacts tends to extend outwards from true frequencies like streaks
as shown in Fig. 9. Secondly, one must make sure that at least one period of
the function is scanned, otherwise, the key assumption of periodicity does not
hold, and taking a discrete Fourier transform does not return a faithful
representation of the spectrum.

Both the minimum resolution and the extent of the scans should depend directly
on the Hamiltonian that generates the rotations, i.e., the cost Hamiltonian.
Since the frequencies in the cost function are the pairwise differences of the
Hamiltonian’s eigenvalues, the minimum recommended scan resolution to resolve
landscape features is given by the highest frequency, i.e., the largest
difference of eigenvalues. Vice versa, the period is obtained from the
greatest common denominator (GCD) of the frequencies. Writing this as
$GCD(\\{\omega\\})$, the period is $\frac{2\pi}{GCD(\\{\omega\\})}$. This is
because the period is $2\pi$ times the smallest number that multiplies all
frequencies to integers, which is $\frac{1}{GCD(\\{\omega\\})}$. Since the
frequencies are the differences of eigenvalues, the GCD of the eigenvalues
will be a divisor of the GCD of the frequencies. For the case of diagonal
Hamiltonians such as the QAOA problem Hamiltonian, the eigenvalues are
themselves sums and differences of the Hamiltonian coefficients. Take an
Hamiltonian defined as

$H=\sum_{k=1}^{m}c_{k}H_{k}$ (23)

where the $H_{i}$ are diagonal Pauli strings (Pauli strings composed of $I$
and $Z$) with eigenvalues $\pm 1$. Then one can see that the eigenvalues of
$H$ obey $\lambda_{i}(H)\in\\{\pm c_{1}\pm c_{2}\cdots\pm c_{m}\\}$, meaning
that the frequency spectrum when this Hamiltonian is used as generator will be

$\Omega_{\gamma}=\\{\lambda_{i}(H)-\lambda_{j}(H)\\}_{i,j}\subset\\{\pm\\{2c_{1},0\\}\pm\\{2c_{2},0\\}\dots\pm\\{2c_{m},0\\}\\}.$
(24)

Therefore for QAOA, a suitable period in the $\gamma$ direction can be found
from the GCD of the coefficients of $H$, without the need to know the
eigenvalues. Since the $c_{k}$’s always appear in the spectrum in multiples of
2, the period is at most $\frac{\pi}{GCD(\\{c_{k}\\})}$. Regarding the maximum
theoretical frequency, which informs the sampling resolution, we can use Eq.
24 to infer that this is

$\max|f_{\gamma}|\leq 2\sum_{k}|c_{k}|.$ (25)

Conversely, for the $\beta$ direction the generator is always
$\sum_{i=1}^{n}X_{i}$, which gives a theoretical frequency spectrum

$\Omega_{\beta}\subset\\{0,\pm 2,\pm 4,\cdots\pm 2n\\}.$ (26)

The period is therefore simply $\pi$. The maximum frequency in theory is $2n$,
however for $p=1$ QAOA one finds (by backpropagating the measurement
Hamiltonian in the Heisenberg picture) that in practice

$\max|f_{\beta}|\leq 2w,$ (27)

where $w$ is the maximum weight of the Pauli strings composing the problem
operator.

Figure 9: Impact of cost landscape scan resolution (top row) and the
respective Fourier spectrum (bottom row). The scan resolutions are $20\times
20$ (a), $50\times 50$ (b) and $100\times 100$ (c). Increasing the resolution
of the cost landscape scans results in a reduction of artifacts in the Fourier
spectrum, revealing the true frequencies.

##  Appendix B Theoretical frequencies of selected Hamiltonians

In this appendix, we will show how we arrived at the theoretical frequencies
for Hamiltonians of form $Z_{0}+Z_{1}+cZ_{0}Z_{1}$ presented in Table 1. As
explained in the text, for brevity we write the absolute values of the
frequencies, as for each positive frequency the corresponding negative
frequencies are also present. We will start with equations for linear terms
($Z_{0}$, $Z_{1}$), derived from Eq. (16).

$\displaystyle\langle C_{0}\rangle=\langle C_{1}\rangle$
$\displaystyle=\sin(2\beta)\sin(2\gamma)\cos(2c\gamma)$ (28)
$\displaystyle=\frac{1}{2}\sin(2\beta)[\sin((2c+2)\gamma)+\sin((2c-2)\gamma)]$

In Eq. (28), the single $\beta$-containing factor $\sin(2\beta)$ means that
$f_{\beta}=2$. To find $f_{\gamma}$, we can use the identity
$\sin(x)\cos(y)=\frac{1}{2}[\sin(x+y)+\sin(x-y)]$. This tells us that
$\sin(2\gamma)\cos(2c\gamma)$ produces both $f_{\gamma}=2c+2$ and
$f_{\gamma}=2c-2$. Combined, the (absolute value) frequencies produced by
linear terms are $(2c\pm 2,2)$.

Exception to this rule: When $c=1$, $\langle
C_{0}\rangle=\frac{1}{2}\sin(2\beta)[\sin(4\gamma)+\sin(0\cdot\gamma)]=\frac{1}{2}\sin(2\beta)\sin(4\gamma)$
as given by Eq. (28). Contrary to the previous prediction, in this case the
only frequency left is (4, 2). The component (0,2), predicted to appear in the
term $\sin(2\beta)\sin(0\gamma)$, does not show up because
$\sin(0\cdot\gamma)=0$ and therefore the entire term is equal to 0.

$\displaystyle\langle C_{01}\rangle$
$\displaystyle=c\sin(4\beta)\sin(2c\gamma)\cos(2\gamma)-\frac{c}{2}\sin(2\beta)^{2}[\cos(4\gamma)-1]$
$\displaystyle=\frac{c\sin(4\beta)}{2}[\sin((2c+2)\gamma)+\sin((2c-2)\gamma)]-\frac{c}{2}[-\frac{\cos(4\beta)}{2}+\frac{1}{2}][\cos(4\gamma)-1]$
$\displaystyle=\frac{c}{2}[\sin(4\beta)\sin((2c+2)\gamma)+\sin(4\beta)\sin((2c-2)\gamma)]\;$
$\displaystyle\hskip
14.22636pt+\frac{c}{4}[\cos(4\beta)\cos(4\gamma)+\cos(4\gamma)-\cos(4\beta)+1]\,.$
(29)

In Eq. (29), each term represents one distinct frequency component. In the
respective order of terms, the frequencies $\omega=(f_{\gamma},f_{\beta})$
contributing to $\langle C_{01}\rangle$ are $(2c+2,4)$, $(2c-2,4)$, (4, 4),
(4, 0), (0, 4), and (0, 0). Note that terms missing a $\cos(\gamma)$ or
$\cos(\beta)$ component can be rewritten to include terms with a frequency of
$0$, for instance $\cos(4\gamma)=\cos(0\beta)\cos(4\gamma)$ and
$1=\cos(0\beta)\cos(0\gamma)$. The trigonometric identities
$\sin^{2}(x)=\frac{1-\cos(2x)}{2}$ and
$\sin(x)\cos(y)=\frac{1}{2}[\sin(x+y)+\sin(x-y)]$ were used in the
simplification of the above equations.

##  Appendix C Metrics

We experimented with several different roughness metrics before deciding on
the two metrics we ultimately used for our research (which are explained in
Section II.3). In this appendix, we discuss these other metrics and explore
their strengths and weaknesses.

###  C.1 Total variation family

We first attempted calculating TV using an algorithm defined in Section 3 of
Ref. [21], which we explained in Section II.3.1). In that work, TV is scaled
by a factor of $\frac{1}{b-a}$ such that its value is invariant to scaling of
the function’s domain. In the case of a periodic function, if one chooses the
sampling domain such that it covers exactly one period of the cost function,
this normalization can be omitted. To estimate the TV, we sample 200 random
directions, with 200 samples taken along each of these directions, which
usually gives an approximate std$/$mean ratio $\leq 10\%$ between runs. We
found this ratio a good compromise between the computation time and required
accuracy.

In Ref. [21], the final roughness index is defined as:

$I:=\frac{\sigma}{\mu}$ (30)

where $\mu$ and $\sigma$ are the mean and standard deviation of normalized TV
in different directions. This roughness index measures how much TV _varies in
different directions_ while seeming unaffected by the raw value of TV.
Experimentally, we found that that landscapes with higher $\mu$ also tend to
have higher values of $\frac{\sigma}{\mu}$. Since calculating $\mu$ is simpler
and easier to interpret than $\frac{\sigma}{\mu}$, we decided to use only it.
However, this choice might depend on the type of the landscapes analyzed,
since for some landscapes studied in Ref. [21], while the value of $\mu$ was
constant, the value of $\sigma$ was changing significantly, so it might be
valuable to monitor both values in future experiments.

In addition, we experimented with another version of estimating TV by
numerically approximating the $m$-dimensional integral of the cost function
(where $m$ is the number of parameters) through full discrete integration.
Especially on low dimensions, this provides similar results to using our
randomized estimation of $\mu$ and does not require setting metaparameters
such as the number of random directions. However, the sampling requirement
scales exponentially with the number of dimensions, while estimating $\mu$
doesn’t require more computation time as the number of dimensions increases,
as long as the number of sampled directions is fixed. This did not matter in
the context of this work, but we still decided to use $\mu$, as a potentially
more future-proof metric.

###  C.2 Fourier metric family

As discussed in Ref. II.3.2, large Fourier coefficients on higher frequencies
mean that a function is more oscillatory and hence should have a higher
roughness index. As observed in Ref. [30], the amplitudes that large
frequencies can have are bounded by total variation, which provides further
evidence that roughness metrics in the Fourier space can directly quantify
roughness of the cost landscape, but may contain complementary information. It
holds that

$\lvert\omega\rvert c_{\omega}\leq\frac{2}{\pi}TV(\mathcal{C}),$ (31)

where $\omega$ represents the frequency, $c_{\omega}$ is the coefficient at
that frequency, and $\mathcal{C}$ is the cost function. Inspired by the
relation in Eq. (31), we explore the following quantity as a potential metric
for the roughness of a multidimensional function:

$R_{f}^{\text{max}}:=\max_{\omega}[\lvert c_{\omega}\rvert\cdot\|\omega\|]\,,$
(32)

where we are now considering a $n$-dimensional Fourier transform and thus the
coefficients are labeled by a frequency vector $\omega\in\mathbb{Z}^{N}$. Up
to constant factors, this can be seen to upper bound the TV. We refer to this
quantity as “Fourier maximum”.

When the metric is defined like this, however, it is sensitive to small
changes in the value of the high-frequency coefficients and is blind to the
low-frequency modes. We, therefore, experimented with the index obtained by
taking the mean over the Fourier coefficients instead of the maximum,

$R_{f}^{\text{mean}}:=\frac{1}{\lvert\hat{\omega}\rvert}\sum_{\omega}\lvert
c_{\omega}\rvert\cdot\|\omega\|\,,$ (33)

where $\hat{\omega}$ is the maximal present frequency. We refer to this metric
as “Fourier mean”. We observe that Fourier density (described in Sec. II.3.2)
and Fourier mean correlate very well with each other. However, we decided to
use Fourier Density due to its numerical stability and ease of computation
(see II.3.2). It also does not come with the issues that Fourier maximum has
and therefore we decided to use it as the most reliable metric from this
family.