\section{Introduction}

The Quantum Approximate Optimization Algorithm (\textsc{qaoa}\xspace) \cite{farhi2014quantum} is a variational quantum algorithm \cite{moll2018quantum} that is able to find approximate solutions to combinatorial-optimization problems \cite{crooks2018performance}, and can also used to heuristically find optimal solutions \cite{zhou2018quantum}.

The basic structure of the \textsc{qaoa}\xspace, and some precursor quantum-optimization algorithms \cite{hogg2000quantum,trugenberger2002quantum}, is the application of an alternating series of unitary operations. These unitary operations correspond to applications of the cost function (cost Hamiltonian) that one is attempting to optimize for, and of a \emph{mixer} or \emph{driver} Hamiltonian. The application of these unitaries is intended to evolve the state of the quantum computer such that it has substantial overlap with computational-basis states representing low-energy configurations of the variables being optimized over. In the \textsc{qaoa}\xspace, each unitary is parameterized by a single real number (often referred to in the \textsc{qaoa}\xspace literature as an \emph{angle}), and these parameters are chosen to optimize the quality of the solutions output by the algorithm. A crucial question is how the choice of the \textsc{qaoa}\xspace parameters affects the solutions the \textsc{qaoa}\xspace is able to find, and hence how to best pick the parameters to be used when running the \textsc{qaoa}\xspace.

In this paper we will discuss the \textsc{qaoa}\xspace in the setting of solving (classical) Ising optimization problems, i.e., where the optimization problem's cost function takes the form $\sum_{(i,j)\in E} J_{ij} s_i s_j + \sum_{i\in V} h_i s_i$, and $s_i \in \{-1,+1\}$.
We note that \textsc{Max}-\textsc{Cut}\xspace problems are a subclass of Ising-optimization problems, for which the external-field terms $h_i$ are zero and the couplings $J_{ij}$ are either $0$ or $-1$. Our results are however applicable to the solution of any NP-hard problem by the \textsc{qaoa}\xspace, via an appropriate mapping \cite{lucas2014ising}.

The \textsc{qaoa}\xspace has recently been experimentally demonstrated with up to 23 qubits on Ising-optimization problems with no external-field terms \cite{google_qaoa}.
The experimental results in \cite{google_qaoa} include optimization landscapes for the case where the \textsc{qaoa}\xspace has only a single \emph{layer} (a single application of the cost-function Hamiltonian, followed by a single application of the mixer Hamiltonian).
These landscapes are plots of the expected value of the cost-function Hamiltonian as a function of the two \textsc{qaoa}\xspace parameters.
In our work we present an analytical formula that is able to compute the energy landscape of the single-layer \textsc{qaoa}\xspace for any Ising-problem instance.
We compare our analytical results with the numerical results from \cite{google_qaoa} (which include both numerical results from classical simulations of the \textsc{qaoa}\xspace, and measured data from running the \textsc{qaoa}\xspace on a quantum processor), and find excellent agreement---showing that our analytical formula is able to predict the energy landscape.
Our analytical formula gives us the ability to analytically optimize the \textsc{qaoa}\xspace parameters---this enables one to compute the optimal \textsc{qaoa}\xspace parameters classically, and then only run the \textsc{qaoa}\xspace on a quantum processor with the optimal set of angles \cite{streif2020training}.
This yields a large speedup versus the current standard mode of operating the single-layer \textsc{qaoa}\xspace, where the parameters are optimized numerically via an iterative process or via performing a brute-force search of the landscape.

In addition to replicating the landscapes for the three problem instances given in \cite{google_qaoa}, we also explore the landscapes our analytical formula obtains for a range of other Ising problems, some of which we discover have qualitatively different features than the landscapes published in the \textsc{qaoa}\xspace literature so far.

One of the key benefits of having an analytical formula is that we can study the performance of the \textsc{qaoa}\xspace on problems using far more qubits than we can classically simulate or that we have experimental hardware containing.
We use our analytical formula to predict the performance of the \textsc{qaoa}\xspace on the G-set \textsc{Max}-\textsc{Cut}\xspace benchmark instances \cite{G_set}, which includes problems having up to $\num{7000}$ vertices (qubits), as well as on Ising instances with up to $\num{100000}$ vertices.
These results give us an indication of how well the \textsc{qaoa}\xspace may perform on problem instances with sizes far beyond the $\sim 50$ spins that is the current limit for direct testing of the \textsc{qaoa}\xspace.

There has been previous work \cite{wang2018quantum} on studying the \textsc{qaoa}\xspace parameters and energy-expectation landscapes analytically in the single-layer setting. \cite{wang2018quantum} derives a formula for the energy expectation as a function of the \textsc{qaoa}\xspace parameters for the specific case of \textsc{Max}-\textsc{Cut}\xspace on unweighted graphs, and \cite{Bravyi2019ObstaclesTS} generalized this formula for the energy expectation for the \textsc{qaoa}\xspace to the case of weighted graphs. Our work builds on this: we include in our analysis the case of non-zero external-field terms in the Ising Hamiltonian. This allows us to assess the energy expectation for the single-layer \textsc{qaoa}\xspace solving any Ising problem.

\section{Analytical Formulae for Expected Energies}

\begin{figure}[t]
\centering
\includegraphics[width=0.5\textwidth]{subgraph2_pic.pdf}
\caption{\label{fig:subgraph}
A graph with $5$ vertices $V = \{1,\dots,5\}$ and $5$ edges $E=\{(12), (13),\dots, (24)\}$.
The product $\prod_{(1k)\in E, k\neq 2}$ involves the vertices connected to $1$ that are not $2$, i.e.\ the vertices $3$ and $5$.
The product $\prod_{(1k)\in E, (2k)\notin E}$ involves the vertices that are connected to $1$ but that are not connected to $2$, which is the single vertex $5$. (Note that it is implicit that $k=2$ should be excluded.)
The product $\prod_{(1k)\in E, (2k)\in E}$ involves the vertices that are connected to both $1$ and $2$, which is the vertex $3$.
}
\end{figure}

Given an undirected graph $G = (V, E)$ defined by vertices $V=\{1,\dots,n\}$ and edges $E\subset V\times V$, as well an external field $h_i$ associated with each vertex, and a coupling strength (edge weight) $J_{ij}$ associated with each edge, we aim to find the configuration $s\in\{-1,+1\}^n$ that minimizes the (Ising) cost function
\begin{align}
C & =  \sum_{i\in V} h_i s_i + \sum_{(i,j)\in E} J_{ij} s_i s_j.
\end{align}
This cost function is equivalent to the cost Hamiltonian
\begin{align}
H_C &= \sum_{i\in V} C_i + \sum_{(i,j)\in E} C_{ij}
\end{align}
where
\begin{align}
\left\{
\begin{aligned}
C_{i} & = h_i \sigma^z_i \\
C_{ij} & = J_{ij} \sigma^z_i \sigma^z_j
\end{aligned}
\right.
.
\end{align}
The unitary transformation $U$ for a \textsc{qaoa}\xspace circuit with depth $p=1$ is defined as
\begin{align}
U & = U_B^\beta U_C^\gamma
\end{align}
where
\begin{align}
U_C^\gamma & = e^{-\mathrm{i}\gamma H_C}= e^{-\mathrm{i}\gamma (\sum_{ij} J_{ij} \sigma^z_i \sigma^z_j + \sum h_i \sigma^z_i)} \\
& = \prod_{(i,j)\in E} (\cos(J_{ij} \gamma)\sigma_i^0\sigma_j^0 - \mathrm{i} \sin( J_{ij}\gamma) \sigma^z_i \sigma^z_j) \prod_{i\in V} (\cos(h_i \gamma)\sigma_i^0 - \mathrm{i} \sin(h_i\gamma) \sigma^z_i) \\
U_B^\beta & = e^{-\mathrm{i}\beta H_B}= e^{-\mathrm{i}\beta \sum_i \sigma^x_i} \\
& = \prod_{i\in V} ( \cos(\beta)\sigma_i^0 - \mathrm{i} \sin(\beta) \sigma^x_i )  .
\end{align}
Applying $U = U_B^\beta U_C^\gamma$ to an initial, uniform superposition of all bit strings, $\ket{\psi_0}$, we obtain the final \textsc{qaoa}\xspace state
\begin{align}
\ket{\beta,\gamma}& =U \ket{\psi_0}.
\end{align}
The expectation value of $H_C$ for this final state is
\begin{align}
F(\beta, \gamma) & \coloneqq \sum_{i\in V} \langle C_{i} \rangle  + \sum_{(i,j)\in E} \langle C_{ij} \rangle
\end{align}
where
\begin{align}
\left\{
\begin{aligned}
\langle C_{i} \rangle
& = \bra{\beta,\gamma}C_i\ket{\beta,\gamma}
= \bra{\psi_0}U^\dagger C_i U \ket{\psi_0}
= h_i \Tr [\rho_0 U^{\dagger} \sigma^z_{i} U]\\
\langle C_{ij} \rangle
& = \bra{\beta,\gamma}C_{ij}\ket{\beta,\gamma}
= \bra{\psi_0}U^\dagger C_{ij} U \ket{\psi_0}
= J_{ij}\Tr [\rho_0 U^{\dagger} \sigma^z_i\sigma^z_j U]\\
\end{aligned}
\right.
\label{eq:trace_ci}
\end{align}
and $\rho_0 =\ket{\psi_0} \bra{\psi_0} $ is the initial density matrix.

\subsection{Formula for General Case with Variable Couplings}
In Appendix~\ref{app:traces} we derive explicit expressions for the two traces of Equation~\ref{eq:trace_ci}, giving us the results
\begin{align}
\left\{\begin{aligned}
\label{eq:main_res}
\langle C_i \rangle & =
h_i \sin(2\beta)\sin(2\gamma h_i)\prod_{(ik)\in E}\cos(2\gamma J_{ik}) \\
\begin{split}
\langle C_{ij} \rangle & =
\frac{J_{ij} \sin(4\beta)}{2} \sin(2 \gamma J_{ij})
\Big[ \cos(2 \gamma h_i) \prod_{\substack{(ik) \in E\\ k \neq j}} \cos(2 \gamma J_{ik} ) + \cos(2 \gamma h_j) \prod_{\substack{(jk) \in E \\ k \neq i}} \cos(2 \gamma J_{jk})\Big]  \\
& \quad -\frac{J_{ij}}{2} (\sin(2\beta))^2  \prod_{\substack{(ik)\in E\\ (jk)\notin E}} \cos(2 \gamma J_{ik})  \prod_{\substack{(jk)\in E\\ (ik)\notin E}}\cos(2\gamma J_{jk})  \\
& \qquad \times \Big[\cos(2 \gamma (h_i{+}h_j)) \prod_{\substack{(ik)\in E\\(jk)\in E}} \cos(2 \gamma (J_{ik}{+}J_{jk})) -
\cos(2 \gamma (h_i{-}h_j)) \prod_{\substack{(ik)\in E\\(jk)\in E}} \cos(2 \gamma (J_{ik}{-}J_{jk}))\Big].
\end{split}
\end{aligned}
\right.
\end{align}
Equation~\ref{eq:main_res} is the contribution of a single vertex and connected subgraphs in the graph. Connected subgraphs are formed by two vertices $(i,j)$ connected by an edge and the nearest neighbours of the $i$ and $j$ vertices, as shown in Fig.~\ref{fig:subgraph}. From this general formula we will derive the expressions corresponding to simpler problem instances and graphs, including previously derived results.

Without loss of generality we can assume the graph to the complete graph $K_n$ as we can use $J_{ij}=0$ for absent edges.
Doing so simplifies the above expressions  to
\begin{align}
\left\{\begin{aligned}
\label{eq:K_res}
\langle C_i \rangle_K & =
h_i \sin(2\beta)\sin(2\gamma h_i)\prod_{k\neq i}\cos(2\gamma J_{ik}) \\
\begin{split}
\langle C_{ij} \rangle_K & =
\frac{J_{ij} \sin(4\beta)}{2} \sin(2 \gamma J_{ij})
\Big[ \cos(2 \gamma h_i) \prod_{k \neq i,j} \cos(2 \gamma J_{ik} ) + \cos(2 \gamma h_j) \prod_{k \neq i,j} \cos(2 \gamma J_{jk})\Big]  \\
& \quad -\frac{J_{ij}}{2} (\sin(2\beta))^2 \Big[\cos(2 \gamma (h_i{+}h_j)) \prod_{k \neq i,j} \cos(2 \gamma (J_{ik}{+}J_{jk})) -
\cos(2 \gamma (h_i{-}h_j)) \prod_{k \neq i,j} \cos(2 \gamma (J_{ik}{-}J_{jk}))\Big].
\end{split}
\end{aligned}
\right.
\end{align}
If, additionally, one assumes the absence of local fields, i.e., $h_i=0$ and $h_j=0$, this further simplifies to
\begin{align}
\begin{split}
\langle C_{ij} \rangle_{K,h=0} & =
\frac{J_{ij} \sin(4\beta)}{2}\sin(2 \gamma J_{ij})
\Big[  \prod_{k \neq i,j} \cos(2 \gamma J_{ik} ) + \prod_{k \neq i,j} \cos(2 \gamma J_{jk})\Big]  \\
& \quad -\frac{J_{ij}}{2} (\sin(2\beta))^2 \Big[ \prod_{k \neq i,j} \cos(2 \gamma (J_{ik}{+}J_{jk})) -
 \prod_{k \neq i,j} \cos(2 \gamma (J_{ik}{-}J_{jk}))\Big].
\end{split}
\end{align}
This last expression has been derived previously and appears in Lemma~C.1 of \cite{Bravyi2019ObstaclesTS}.

Returning to the general Equation~\ref{eq:main_res}, we can also consider the case of triangle-free graphs, i.e.,  graphs where $i$ and $j$ share no common vertices $k$ with both $(i,k)\in E$ and $(j,k)\in E$, for which the expectation reduces to
\begin{align}
\label{eq:main_res2_notrian]}
\langle C_{ij} \rangle_{\mathrm{tf}} & =
\frac{J_{ij} \sin(4\beta)}{2} \sin(2 \gamma J_{ij})  \Big[ \cos(2 \gamma h_i) \prod_{(ik) \in E} \cos(2 \gamma J_{ik} ) + \cos(2 \gamma h_j) \prod_{(jl) \in E} \cos(2 \gamma  J_{jl})\Big]  \\
& \quad +J_{ij} (\sin(2\beta))^2 \sin(2 \gamma h_i) \sin(2 \gamma h_j)
\prod_{\substack{(ik)\in E\\ (jk)\notin E}} \cos(2 \gamma J_{ik})  \prod_{\substack{(jk)\in E\\ (ik)\notin E}}\cos(2\gamma J_{jk}).
\end{align}

\subsection{Formulae for Ising Problems with Constant Couplings}
Next, we will derive the expressions of the expectation value for some special-case Hamiltonians where the couplings are constant ($J_{ij}{=}\pm 1$). For the NP-hard ``P5'' Ising problem of \cite{Barahona_1982} for a graph $(V,E)$ the goal is to minimize the Hamiltonian
\begin{align}
    H_{\mathrm{P5}} &= \sum_{(i,j)\in E} \sigma^z_i \sigma^z_j + \sum_{i\in V} \sigma^z_i
\end{align}
With $d_i$ and $d_j$ the degrees of $i$ and $j$, and $f_{ij}$ the number of shared vertices between $i$ and $j$, it is easy to see that for the expectation we now have
\begin{align}
\label{eq:Cij_P5}
\left\{
\begin{aligned}
\langle C_{i} \rangle_{\mathrm{P5}}
& = \sin(2\beta) \sin(2\gamma)(\cos(2\gamma))^{d_i} \\
\begin{split}
\langle C_{ij} \rangle_{\mathrm{P5}}
& =
\frac{\sin(4\beta)}{2} \sin(2 \gamma )
\Big[ (\cos(2 \gamma ))^{d_i} + (\cos(2 \gamma))^{d_j}\Big]\\
& \quad + \frac{1}{2} (\sin(2\beta))^2  (\cos(2 \gamma ))^{d_i+d_j-2f_{ij}-2}  (1-(\cos(4 \gamma ))^{f_{ij}+1}).
\end{split}
\end{aligned}\right.
\end{align}

For the \textsc{Max}-\textsc{Cut}\xspace problem, with no local fields and couplings $J_{ij}=-1$, the expectation reads
\begin{align}
\left\{
\begin{aligned}
H_{\mathrm{MC}}& =-\sum_{(i,j)\in E} \sigma^z_i \sigma^z_j  \\
\langle C_{i} \rangle_{\mathrm{MC}} & = 0 \\
\langle C_{ij} \rangle_{\mathrm{MC}} &
= \frac{\sin(4\beta)}{2} \sin(2\gamma) [(\cos(2\gamma))^{d_i-1}+(\cos(2\gamma))^{d_j-1}] \\
& \quad -\frac{1}{2}(\sin(2\beta))^2 (\cos(2 \gamma))^{d_i+d_j-2f_{ij}-2} (1-(\cos(4\gamma))^{f_{ij}}) .
\end{aligned}
\right.
\end{align}
This is the result obtained by Wang et al.\ \cite{PhysRevA.97.022304}.

The expression in the absence of spin-spin couplings reads,
\begin{align}
\label{eq:main_result_nocoupling1}
H_{\mathrm{NC}} & = \sum_{i\in V}h_i  \sigma^z_i  \\
\langle C_{i} \rangle_{\mathrm{NC}} & =
h_i \sin(2\beta)  \sin(2\gamma h_i) \label{eq:no_coupling} \\
\label{eq:main_result_nocoupling2}
\langle C_{ij} \rangle_{\mathrm{NC}} & = 0.
\end{align}
The solution to this optimization problem is trivial, however, the determination of the angles for the optimal expectation value is not.

We will use the analytical formulae derived in this section to gain insight into the behavior of expectation value landscapes for \textsc{qaoa}\xspace $p=1$. First, we will study \textsc{Max}-\textsc{Cut}\xspace instances, as the paper progresses we will gradually increase the complexity of the instances we study until we reach general Ising problems.

\subsection{Formulae for Standard Sherrington-Kirkpatrick Model without Local Field}
\label{sec:SK_maxcut_expectation}
To analyze the Sherrington-Kirkpatrick model we assume  no local field, i.e.\ $h_i=0$, and that the all-to-all coupling values $J_{ij}$ are distributed according to an independent and identical distribution with mean $0$.
Taking the average over this distribution and using the facts that all $J_{ij}$ distributions are independent and that $\EE_J[J_{ij}]=0$ we are interested in the expectation of expectations
\begin{align}
\EE_J[\langle C_{ij} \rangle] & =
\sin(4\beta)\cdot \EE_J[J  \sin(2 \gamma J)]
\cdot(\EE_J[\cos(2 \gamma J)])^{n-2}.
\end{align}
Often one is interested in the $n\rightarrow \infty$ behaviour of this model, which is easiest to discuss when we scale the $J_{ij}$ random variables by $1/\sqrt{n-1}$ such that we have
\begin{align}
J_{ij} & = \frac{J'}{\sqrt{n-1}},
\end{align}
where $J'$ is a random variable with mean $0$ and
variance $\sigma^2 \in O(1)$ independent of $n$.
Doing so will result in an $O(1)$ energy per spin,  i.e.\ $\langle C/n\rangle = (n-1)/2\cdot \langle C_{ij}\rangle$ will be of order $1$.
We use $\sqrt{n-1}$ instead of the more common $\sqrt{n}$ to make the finite $n$ expressions work out nicer.

For the normal distribution $J'\sim N(0,\sigma^2)$ we have $J\sim N(0,\sigma^2/(n-1))$, and hence
\begin{align}
\EE_{J\sim N(0,\sigma^2/(n-1))}[J\sin(2\gamma J)] & = 2\gamma\sigma^2\cdot (n-1)^{-1}\cdot e^{-2\gamma^2\sigma^2/(n-1)}\\
\EE_{J \sim N(0,\sigma^2/(n-1))}[\cos(2\gamma J)] & = e^{-2\gamma^2\sigma^2/(n-1)}.
\end{align}
The expectation of the expectations of the energy per spin therefore equals
\begin{align} \label{eq:expSK}
\EE_{J'\sim N(0,\sigma^2)}[\langle C/n \rangle] & =
\sin(4\beta)\cdot\gamma \sigma^2 e^{-2\gamma^2\sigma^2}.
\end{align}
Note that this is an exact expression for finite $n$ and it allows us to compute the optimal angles $\beta$, $\gamma$ that minimize the expectation as a function of $n$ and $\sigma$:
\begin{align} \label{eq:expSKopt}
\beta_\mathrm{min} = \frac{-\pi}{8}, \quad
\gamma_\mathrm{min} =  \frac{1}{2\sigma}, \quad
\text{with}\quad
\EE_{J'\sim N(0,\sigma^2)}[\langle C/n \rangle]  = - \frac{\sigma}{2\sqrt{e}} \approx - 0.303 \sigma.
\end{align}

All of the above results agree with those in \cite[Section~4]{farhi2019quantum}.
Moreover, \cite{farhi2019quantum} shows that in the $n\rightarrow \infty$ limit this result is independent of the details of the distribution $J'$.
As an example, consider the bimodal distribution $B$ where $\Pr[J'{=}{\pm}\sigma] = 1/2$ for some fixed $\sigma$, such that $\EE[J']=0$ and $\EE[J'^2] = \sigma^2$.
As a result,  we have $\Pr[J{=}{\pm}\sigma/\sqrt{n-1}] = 1/2$ with $\EE_{J}[J \sin(2\gamma J)] = (\sigma/\sqrt{n-1})\sin(2\gamma \sigma/\sqrt{n-1})$ and $\EE_{J}[\cos(2\gamma J)] = \cos(2\gamma \sigma/\sqrt{n-1})$, and hence
\begin{align}
\EE_{J'\sim B}[\langle C/n\rangle] & =
\frac{1}{2}\sin(4\beta) \cdot
\sigma\sqrt{n-1}\cdot \sin(2\gamma \sigma/\sqrt{n-1})\cdot(\cos(2\gamma \sigma/\sqrt{n-1}))^{n-2}.
\end{align}
In the $n\rightarrow \infty$ limit for large systems we thus get again
\begin{align}
\lim_{n\rightarrow \infty}\EE_{J'\sim B}[\langle C/n \rangle] & =
\sin(4\beta)\cdot \gamma \sigma^2 e^{-2\gamma^2 \sigma^2},
\end{align}
which is identical to Equation~\ref{eq:expSK}.

It is important to realize, though, that there are distributions over $J'$ with mean $0$ and variance $\sigma^2$ that give different results for the optimal angles and expectations.
Consider the setting where we want to describe a spin glass where each vertex has an expected number of $d$ edges for which it has $J{=}{\pm}1$, while for the other $n-1-d$ edges it has $J{=}0$.
To capture this setting  we define the random variable $J$ according to a trimodal distribution $T$ over coupling values $J{\in} \{-1,0,+1\}$ with $\Pr[J{=}{\pm}1] = d/(2n-2)$ and $\Pr[J{=}0]=1-d/(n-1)$ such that $\EE_{J\sim T}[J]=0$ and $\EE_{J\sim T}[J^2]= d/(n-1)$.
For this distribution we have the expectations $\EE_{J\sim T}[J\sin(2\gamma J)] = d/(n-1) \sin(2\gamma)$, and $\EE_{J\sim T}[\cos(2\gamma J)] = 1+d(\cos(2\gamma)-1)/(n-1)$ and hence
\begin{align}
\EE_{J\sim T}[\langle C/n\rangle] & =
\frac{1}{2}\sin(4\beta) \cdot  d \cdot \sin(2\gamma) \cdot (1-d/(n-1)+d/(n-1)\cos(2\gamma))^{n-2}.
\end{align}
In the $n\rightarrow \infty$ limit this gives
\begin{align}\label{eq:SK-d-n-inf}
\EE_{J\sim T}[\langle C/n \rangle] & =
\frac{1}{2}\sin(4\beta)\cdot d \cdot \sin(2\gamma)\cdot e^{-2d(\sin(\gamma))^2}.
\end{align}

If we set $d=1$ the Equation \ref{eq:SK-d-n-inf} is maximized by $\beta = -\pi/8$, $\gamma = 0.452278$, such that $\EE_{J\sim T}[\langle C/n\rangle] = -0.268281$.
Note now that this derived optimum is different from the results of Equation~\ref{eq:expSKopt}.
This difference is explained by the fact that the related distribution of $J' = J\sqrt{n-1}$ (with mean $0$ and variance $1$) is no longer independent of $n$, as it is defined by $\Pr[J'{=}{\pm}\sqrt{n-1}] = 1/(2n-2)$ and $\Pr[J'{=}0] = (n-2)/(n-1)$.
Because of this $n$-dependency it no longer holds that the $n\rightarrow \infty$ properties are described by  Equations~\ref{eq:expSK} and \ref{eq:expSKopt}.

More specifically,  Equation~\ref{eq:expSK} will hold as long as we have the expectations in the $n\rightarrow \infty$ limit:
\begin{align}\label{eq:Jmoments}
 \EE_J[nJ^2] \rightarrow \sigma^2\text{, and~}
 \EE_J[nJ^r] \rightarrow 0, \text{~for $r \in \{3,4,\dots\}$.}
\end{align}
That way we have
\begin{align} \label{eq:ninfinitylimit}
\lim_{n\rightarrow \infty}\EE_J[\langle C/n\rangle] & =
\lim_{n\rightarrow \infty}\frac{n-1}{2}\cdot\sin(4\beta)\cdot \EE_J[J\sin(2\gamma J)]\cdot(\EE_J[\cos(2\gamma J)])^{n-2} \\
& = \sin(4\beta)\cdot \gamma\sigma^2 \cdot e^{-2\gamma^2 \sigma^2}.
\end{align}
The above trimodal distribution does not fit this mold as it has $\EE_J[nJ^3] \rightarrow d$.

\subsection{Formulae for the Sherrington-Kirkpatrick Model with Constant Local Field}
Consider now the SK model with a constant local field $h$ and with $J_{ij}$ couplings that again are i.i.d.\ with mean zero. The local energies can now be expressed as
\begin{align}
\left\{\begin{aligned}
\EE_J[\langle C_i \rangle] & =
h \sin(2\beta)\sin(2\gamma h)\cdot \EE_J[\cos(2\gamma J)]^{n-1} \\
\EE_J[\langle C_{ij} \rangle] & =
\sin(4\beta)\cos(2 \gamma h)\cdot \EE_J[J  \sin(2 \gamma J)] \cdot \EE_J[\cos(2 \gamma J)]^{n-2}.
\end{aligned}
\right.
\end{align}
The per spin energy thus equals
\begin{align}
\EE_J[\langle C/n\rangle] & =
h \sin(2\beta)\sin(2\gamma h)\cdot \EE_J[\cos(2\gamma J)]^{n-1} + \frac{n-1}{2}\sin(4\beta)\cos(2 \gamma h)\cdot \EE_J[J  \sin(2 \gamma J)] \cdot \EE_J[\cos(2 \gamma J)]^{n-2}.
\end{align}

Assuming again that the $J$ values are normally distributed according to $J=J'/\sqrt{n-1}$ with $J'\sim N(0,\sigma^2)$ we get for the energy per spin
\begin{align}
\EE_{J'\sim N(0,\sigma^2)}[\langle C/n \rangle]
& = (\gamma \sigma^2 \cos(2 h \gamma) \sin(4 \beta) + h \sin(2 \beta) \sin(2 h \gamma))e^{-2 \gamma^2\sigma^2}.
\end{align}

As an example, if we take $\sigma=1$ and $h=1$ we have the optimal angles $\beta \approx -0.54854$ and $\gamma \approx 0.3962$, such that $\EE[\langle C/n\rangle] \approx -0.62791$.
Compare this with previous result for $(\sigma,h) = (1,0)$ for which the optimal angles are $\beta = \frac{-\pi}{8} \approx -0.3927$ and $\gamma = 0.5$, with
$\EE[\langle C/n\rangle]= - \frac{1}{2\sqrt{e}} \approx -0.3033$.

\subsection{Formulae for the Sherrington-Kirkpatrick Model with Normal Local Field}
Next, we will study the SK model with local fields distributed according to an independent and identical distribution with mean 0. The local energies can be expressed as
\begin{align}
\left\{\begin{aligned}
\EE_{J,h}[\langle C_i \rangle] & =
 \sin(2\beta)\cdot \EE_h[h \sin(2\gamma h)]\cdot \EE_J[\cos(2\gamma J)]^{n-1} \\
\EE_{J,h}[\langle C_{ij} \rangle] & =
\sin(4\beta)\cdot \EE_h[\cos(2 \gamma h)]\cdot \EE_J[J  \sin(2 \gamma J)] \cdot \EE_J[\cos(2 \gamma J)]^{n-2}.
\end{aligned}
\right.
\end{align}
The per spin energy thus equals
\begin{align}
\EE_{J,h}[\langle C/n\rangle] & =
\sin(2\beta)\cdot \EE_h[h \sin(2\gamma h)]\cdot \EE_J[\cos(2\gamma J)]^{n-1}\\& \quad + \frac{n-1}{2}\sin(4\beta)\cdot\EE_h[\cos(2 \gamma h)]\cdot \EE_J[J  \sin(2 \gamma J)] \cdot \EE_J[\cos(2 \gamma J)]^{n-2}.
\end{align}
Assume now that the  $h_i$ random variables will have the same normal distribution $h\sim N(0,\sigma^2)$ as $J'$:
\begin{align}
 \EE_{h\sim N(0,\sigma^2)}[h\sin(2\gamma h)]  &= 2\gamma\sigma^2\cdot e^{-2\gamma^2\sigma^2}\\
\EE_{h \sim N(0,\sigma^2)}[\cos(2\gamma h)]  &= e^{-2\gamma^2\sigma^2}.
\end{align}
The expectation of the expectations of the energy per spin equals
\begin{align}
\EE_{J'\sim N(0,\sigma^2), h \sim N(0,\sigma^2)}[\langle C/n \rangle] & =
(2\sin(2\beta) + \sin(4\beta))\cdot \gamma \sigma^2 e^{-4\gamma^2\sigma^2}.
\end{align}
We then compute the optimal angles $\beta$, $\gamma$,
\begin{align}
\beta_\mathrm{min} = \frac{-\pi}{6}, \quad
\gamma_\mathrm{min} =  \frac{1}{2 \sqrt{2} \sigma}, \quad
\text{with}\quad
\EE_{J'\sim N(0,\sigma^2), h \sim N(0,\sigma^2)}[\langle C/n \rangle]  = - \frac{3 }{4} \sqrt{\frac{3}{2 e}} \sigma \approx -0.557 \sigma.
\end{align}

\subsection{Formulae for Random Spin Glasses on Regular Graphs}
\label{sec:reg_graphs_expectation}
When considering random $J$ couplings between vertices of a constant degree $d$ graph $G=(V,E)$, the expectation of expectations for edges $(ij)\in E$ will be
\begin{align}
\EE_J[\langle C_{ij} \rangle] & =
\sin(4\beta)\cdot \EE_J[J  \sin(2 \gamma J)]
\cdot(\EE_J[\cos(2 \gamma J)])^{d-1},
\end{align}
where $d$ is the degree of the graph.
For non-edges $(ij)\notin E$ we obviously have $\EE_J[\langle C_{ij} \rangle] =  0$.
If we assume a normal distribution of $J$ couplings over the edges and for each vertex $i$ the expectation $\EE_J[\sum_j J_{ij}^2] = \sigma^2$, we should take $J\sim N(0,\sigma^2/d)$.
With this distribution we get
\begin{align}
\EE_{J\sim N(0,\sigma^2/d)}[J\sin(2\gamma J)] & = 2\gamma\sigma^2\cdot d^{-1}\cdot e^{-2\gamma^2\sigma^2/d}\\
\EE_{J \sim N(0,\sigma^2/d)}[\cos(2\gamma J)] & = e^{-2\gamma^2\sigma^2/d}
\end{align}
and hence an expectation of the expectations of the energy per spin of
\begin{align}
\EE_{J\sim N(0,\sigma^2/d)}[\langle C/n \rangle] & =
\sin(4\beta)\cdot \frac{n-1}{d}\cdot \gamma\sigma^2\cdot e^{-2\gamma^2\sigma^2}.
\end{align}
The corresponding optimal values are thus
\begin{align}
\beta_\mathrm{min} = \frac{-\pi}{8}, \quad
\gamma_\mathrm{min} =  \frac{1}{2\sigma}, \quad
\text{and}\quad
\EE_{J'\sim N(0,\sigma^2)}[\langle C/n \rangle]  = - \frac{n-1}{d}\cdot \frac{\sigma}{2\sqrt{e}},
\end{align}
which for $d=n-1$ returns the standard SK result.

Next, we add local fields. We obtain for degree $d$ graphs with normal distribution $h\sim N(0,\sigma^2)$,
\begin{align}
\left\{\begin{aligned}
\EE_{J,h}[\langle C_i \rangle] & =
 \sin(2\beta)\cdot \EE_h[h \sin(2\gamma h)]\cdot \EE_J[\cos(2\gamma J)]^{d} \\
\EE_{J,h}[\langle C_{ij} \rangle] & =
\sin(4\beta)\cdot \EE_h[\cos(2 \gamma h)]\cdot \EE_J[J  \sin(2 \gamma J)] \cdot \EE_J[\cos(2 \gamma J)]^{d-1}.
\end{aligned}
\right.
\end{align}

The expectation of the expectations of the energy per spin equals
\begin{align}
\EE_{J'\sim N(0,\sigma^2), h \sim N(0,\sigma^2)}[\langle C/n \rangle] & =
(2\sin(2\beta) + \frac{n-1}{d} \sin(4\beta))\cdot \gamma \sigma^2 e^{-4\gamma^2\sigma^2}.
\end{align}

Here the optimal angle $\gamma_\mathrm{min}$ remains constant as we modify $d$, while $\beta_\mathrm{min}$ changes.

\section{Numerical Results and discussion}

\begin{figure}[t]
\centering
\includegraphics[width=0.6\textwidth]{plot1.pdf}
\caption{\label{fig:google_plots} Negative expectation-value landscapes, obtained using our analytical formulae, Eq.~\ref{eq:main_res}, for the same problem instances studied experimentally in \cite{google_qaoa}. The graphs for each instance are shown above their respective landscape.
(A) Sycamore grid problem with $23$ qubits.
(B) $3$-regular graph with $14$ qubits.
(C) Fully connected SK-model with $11$ qubits.}
\end{figure}

\begin{figure}[t]
\centering
\includegraphics[width=0.6\textwidth]{plot2.pdf}
\caption{\label{fig:google_plots_morequbits} Negative expectation-value landscapes for two $100$-qubit instances.
(A) 3-regular graph with $100$ qubits.
(B) Fully connected SK-model with $100$ qubits.}
\end{figure}

\begin{table}[t]
 \centering
 \begin{tabular}{||c S[table-format=4.0] S[table-format=5.0] S[table-format=5.0] S[table-format=4.0] S[table-format=4.3] S[table-format=1.3] S[table-format=5.3] S[table-format=1.3] c||}
 \hline
 Name & 	{$|V|$}    & {$|E|$}     & {Sum of weights} & {Best-known Ising} & {\textsc{qaoa}\xspace expectation} & {Ratio Ising} & {$1.18\sqrt{|E| \times |V|}$}    & {Ratio exp.} & Weights \\ \hline\hline
 G1   & 800  & 19176 & 19176          & 4072             & 1482.034         & 0.364       & 4621.745  & 0.321 & $\{0,1\}$     \\ \hline
 G2   & 800  & 19176 & 19176          & 4064             & 1480.775          & 0.364       & 4621.745  & 0.32 & $\{0,1\}$     \\ \hline
 G3   & 800  & 19176 & 19176          & 4068             & 1478.236         & 0.363       & 4621.745  & 0.32  & $\{0,1\}$    \\ \hline
 G4   & 800  & 19176 & 19176          & 4116             & 1481.303         & 0.36       & 4621.745  & 0.321 & $\{0,1\}$      \\ \hline
 G5   & 800  & 19176 & 19176          & 4086             & 1479.496         & 0.362       & 4621.745  & 0.32 & $\{0,1\}$     \\ \hline
 G6   & 800  & 19176 & 154            & 4202             & 1679.595         & 0.4       & 4621.745  & 0.363 & $\{0,\pm 1\}$     \\ \hline
 G7   & 800  & 19176 & -150           & 4162             & 1674.41         & 0.402       & 4621.745  & 0.362 & $\{0,\pm 1\}$      \\ \hline
 G8   & 800  & 19176 & -170           & 4180             & 1676.846         & 0.401       & 4621.745  & 0.363 & $\{0,\pm 1\}$      \\ \hline
 G9   & 800  & 19176 & -64            & 4172             & 1676.17        & 0.402       & 4621.745  & 0.363  & $\{0,\pm 1\}$     \\ \hline
 G10  & 800  & 19176 & -160           & 4160             & 1675.396        & 0.403       & 4621.745  & 0.363  & $\{0,\pm 1\}$     \\ \hline
 G11  & 800  & 1600  & 34             & 1094             & 519.61           & 0.475       & 1335.018  & 0.389 & $\{0,\pm 1\}$      \\ \hline
 G12  & 800  & 1600  & -4             & 1116             & 519.61          & 0.466       & 1335.018  & 0.389 & $\{0,\pm 1\}$      \\ \hline
 G13  & 800  & 1600  & 34             & 1130             & 519.61           & 0.46        & 1335.018  & 0.389 & $\{0,\pm 1\}$      \\ \hline
 G14  & 800  & 4694  & 4694           & 1434             & 577.546          & 0.403      & 2286.644  & 0.253 & $\{0,1\}$     \\ \hline
 G15  & 800  & 4661  & 4661           & 1439             & 575.252         & 0.4        & 2278.592  & 0.252 & $\{0,1\}$     \\ \hline
 G16  & 800  & 4672  & 4672           & 1432             & 574.294          & 0.401       & 2281.279  & 0.252 & $\{0,1\}$     \\ \hline
 G17  & 800  & 4667  & 4667           & 1427             & 573.562          & 0.402      & 2280.058  & 0.252 & $\{0,1\}$     \\ \hline
 G18  & 800  & 4694  & 64             & 1920             & 714.097          & 0.372      & 2286.644  & 0.312 & $\{0,\pm 1\}$     \\ \hline
 G19  & 800  & 4661  & -113           & 1925             & 710.511          & 0.369       & 2278.592  & 0.312 & $\{0,\pm 1\}$     \\ \hline
 G20  & 800  & 4672  & -46            & 1928             & 709.332          & 0.368       & 2281.279  & 0.311 & $\{0,\pm 1\}$     \\ \hline
 G21  & 800  & 4667  & -67            & 1929             & 710.908          & 0.369       & 2280.058  & 0.312 & $\{0,\pm 1\}$     \\ \hline
 G27  & 2000 & 19990 & -42            & 6724             & 2693.905         & 0.401         & 7461.109  & 0.361 & $\{0,\pm 1\}$     \\ \hline
 G28  & 2000 & 19990 & -104           & 6700             & 2695.314         & 0.402       & 7461.109  & 0.361 & $\{0,\pm 1\}$     \\ \hline
 G29  & 2000 & 19990 & 80             & 6730             & 2693.685         & 0.4         & 7461.109  & 0.361 & $\{0,\pm 1\}$     \\ \hline
 G30  & 2000 & 19990 & 100            & 6726             & 2693.264          & 0.4       & 7461.109  & 0.361 & $\{0,\pm 1\}$     \\ \hline
 G31  & 2000 & 19990 & -80            & 6700             & 2691.868         & 0.402       & 7461.109  & 0.361 & $\{0,\pm 1\}$     \\ \hline
 G59  & 5000 & 29570 & 96             & 12076            & 4452.877          & 0.369       & 14348.043 & 0.31 & $\{0,\pm 1\}$      \\ \hline
 G61  & 7000 & 17148 & 362            & 11230            & 4590.429         & 0.409       & 12928.191 & 0.355 & $\{0,\pm 1\}$     \\ \hline
 G64  & 7000 & 41459 & 527            & 16975            & 6235.328         & 0.367       & 20102.054 & 0.31 & $\{0,\pm 1\}$      \\ \hline
 \end{tabular}
\caption{Computational results of \textsc{qaoa}\xspace $p=1$ expectation values in Ising energies for a subset of the G-set instances.
Column \emph{Sum of weights} shows the sum of all $J_{ij}$ values.
Column \emph{Best-known Ising} shows the best classical results for the Ising energy, column \emph{\textsc{qaoa}\xspace expectation} gives the best \textsc{qaoa}\xspace expectation value.
The best-known classical results were obtained from \cite{toshiba_benchmarking}.
The columns \emph{Ratio Ising} and \emph{Ratio exp.} show the corresponding ratio of the \textsc{qaoa}\xspace value versus the best classical result and versus the analytical formula value, respectively.}
\label{table:G-set_Ising}
\end{table}

\begin{table}[t]
\begin{tabular}{||c S[table-format=4.0] S[table-format=5.0] S[table-format=5.0] S[table-format=5.0] S[table-format=5.3] S[table-format=0.3] S[table-format=5.0] c||}
\hline
Name & {$|V|$} & {$|E|$}  & {Sum of weights} & {Best-known cut} & {\textsc{qaoa}\xspace expectation cut} & {Ratio cut} &  {Time (s)} & {Weights} \\ \hline\hline
G1   & 800  & 19176 & 19176          & 11624          & 10329.017                  & 0.889 & 28322 & $\{0,1\}$     \\ \hline
G2   & 800  & 19176 & 19176          & 11620          & 10328.387                   & 0.889 & 28562 & $\{0,1\}$     \\ \hline
G3   & 800  & 19176 & 19176          & 11622          & 10327.118                  & 0.889 & 28332 & $\{0,1\}$     \\ \hline
G4   & 800  & 19176 & 19176          & 11646          & 10328.651                  & 0.887 & 28430 & $\{0,1\}$     \\ \hline
G5   & 800  & 19176 & 19176          & 11631          & 10327.748                  & 0.888 & 29039 & $\{0,1\}$     \\ \hline
G6   & 800  & 19176 & 154            & 2178           & 916.797                    & 0.42 & 31953 & $\{0,\pm 1\}$      \\ \hline
G7   & 800  & 19176 & -150           & 2006           & 762.205                     & 0.38 & 32223 & $\{0,\pm 1\}$     \\ \hline
G8   & 800  & 19176 & -170           & 2005           & 753.423                    & 0.376 & 32685 & $\{0,\pm 1\}$     \\ \hline
G9   & 800  & 19176 & -64            & 2054           & 806.085                    & 0.392 & 33056 & $\{0,\pm 1\}$    \\ \hline
G10  & 800  & 19176 & -160           & 2000           & 757.698                    & 0.379 & 33681 & $\{0,\pm 1\}$    \\ \hline
G11  & 800  & 1600  & 34             & 564            & 276.805                    & 0.491 & 1193 & $\{0,\pm 1\}$    \\ \hline
G12  & 800  & 1600  & -4             & 556            & 257.805                    & 0.464 & 1188 & $\{0,\pm 1\}$     \\ \hline
G13  & 800  & 1600  & 34             & 582            & 276.805                    & 0.476 & 1194 & $\{0,\pm 1\}$     \\ \hline
G14  & 800  & 4694  & 4694           & 3064           & 2635.773                   & 0.86 & 4808 & $\{0,1\}$     \\ \hline
G15  & 800  & 4661  & 4661           & 3050           & 2618.126                   & 0.858 & 4982 & $\{0,1\}$     \\ \hline
G16  & 800  & 4672  & 4672           & 3052           & 2623.147                   & 0.859 & 5076 & $\{0,1\}$     \\ \hline
G17  & 800  & 4667  & 4667           & 3047           & 2620.281                   & 0.86 & 5056 & $\{0,1\}$      \\ \hline
G18  & 800  & 4694  & 64             & 992            & 389.049                    & 0.392 & 5125 & $\{0,\pm 1\}$     \\ \hline
G19  & 800  & 4661  & -113           & 906            & 298.755                    & 0.33 & 5187 & $\{0,\pm 1\}$     \\ \hline
G20  & 800  & 4672  & -46            & 941            & 331.666                    & 0.352 & 5186 & $\{0,\pm 1\}$     \\ \hline
G21  & 800  & 4667  & -67            & 931            & 321.954                    & 0.346 & 5237 & $\{0,\pm 1\}$     \\ \hline
G27  & 2000 & 19990 & -42            & 3341           & 1325.952                   & 0.397 & 22103 & $\{0,\pm 1\}$     \\ \hline
G28  & 2000 & 19990 & -104           & 3298           & 1295.657                   & 0.393 & 21916 & $\{0,\pm 1\}$     \\ \hline
G29  & 2000 & 19990 & 80             & 3405           & 1386.082                   & 0.407 & 21741 & $\{0,\pm 1\}$     \\ \hline
G30  & 2000 & 19990 & 100            & 3413           & 1396.632                    & 0.409 & 22292 & $\{0,\pm 1\}$     \\ \hline
G31  & 2000 & 19990 & -80            & 3310           & 1305.934                   & 0.395 & 22460 & $\{0,\pm 1\}$     \\ \hline
G59  & 5000 & 29570 & 96             & 6086           & 2274.439                   & 0.374 & 39462 & $\{0,\pm 1\}$     \\ \hline
G61  & 7000 & 17148 & 362            & 5796           & 2476.214                   & 0.427 & 14794 & $\{0,\pm 1\}$     \\ \hline
G64  & 7000 & 41459 & 527            & 8751           & 3381.164                   & 0.386 & 57647 & $\{0,\pm 1\}$     \\ \hline
\end{tabular}
\caption{Computational results of \textsc{qaoa}\xspace $p=1$ expectation values in cut values for a subset of the G-set instances. Column \emph{Sum of weights} shows the sum of all $J_{ij}$ values. Column \emph{Best-known cut} shows the best classical results for the cut value, column \emph{\textsc{qaoa}\xspace expectation cut} gives the best \textsc{qaoa}\xspace expectation value. The best-known classical results were obtained from \cite{toshiba_benchmarking}. The column \emph{Ratio cut} shows the corresponding ratio of the \textsc{qaoa}\xspace value versus the best classical result for the cut value. \emph{Time (s)} column shows the time required to obtain the full landscape and optimal angles, in seconds, on a \emph{ac5n.large} AWS instance.}
\label{table:G-set_cut}
\end{table}

\begin{table}[t]
\begin{tabular}{||c S[table-format=6.0] S[table-format=6.0] S[table-format=3.0] S[table-format=5.3] S[table-format=6.3] S[table-format=1.3] S[table-format=6.0] c||}
\hline\hline
Name & {$|V|$} & {$|E|$} & {Sum of weights} & {\textsc{qaoa}\xspace expectation} & {$1.18\sqrt{|V| \times (|E|+|V|)}$} & {Ratio exp.} & {Time (s)} & {Weights} \\ \hline \hline
I10  & 800  & 19176 & -144           & 1706.929          & 4717.167          & 0.362 & 31376 & $\{0,\pm 1\}$       \\ \hline
I20  & 800  & 4672  & -128           & 787.059          & 2468.88           & 0.319 & 5269 & $\{0,\pm 1\}$        \\ \hline
I30  & 2000 & 19990 & 86             & 2820.98         & 7825.455          & 0.36 & 22648 & $\{0,\pm 1\}$       \\ \hline
I59  & 5000 & 29570 & 156            & 4940.085         & 15513.747         & 0.318 & 40740 & $\{0,\pm 1\}$       \\ \hline
I64  & 7000 & 41459 & 661            & 6914.15         & 21732.929         & 0.318 & 59894 & $\{0,\pm 1\}$       \\ \hline
I99 & 70000 & 140000 & -62 & 54809.093 & 143067.397 & 0.383 & 125363 & $\{0,\pm 1\}$ \\ \hline
I100 & 100000 & 150000 & 956 & 72847.211 & 186574.382 & 0.39 & 137774 & $\{0,\pm 1\}$ \\ \hline
\end{tabular}
\caption{Computational results of \textsc{qaoa}\xspace $p=1$ expectation values in Ising energies for a subset of the G-set instances with external fields, the I-set. Plus two more custom generated instances, I99, and I100. Column \emph{Sum of weights} shows the sum of all $h_i$ and $J_{ij}$ values,  column \emph{\textsc{qaoa}\xspace expectation} gives the best \textsc{qaoa}\xspace expectation value. The column \emph{Ratio exp.} corresponds to the ratio of the \textsc{qaoa}\xspace energy expectation value versus the optimal energy predicted by the  formula $1.18\sqrt{|V| \times (|E|+|V|)}$. \emph{Time (s)} column shows the time required to obtain the full landscape and optimal angles, in seconds, on an \emph{ac5n.large} AWS instance.}
\label{table:G-set_ising_field}
\end{table}

\begin{figure}[t]
\centering
\includegraphics[width=0.4\textwidth]{plot3.pdf}
\caption{\label{fig:G_set} Expectation-value landscape for the G64 instance with $\num{7000}$ vertices and $\num{41459}$ edges, obtained using the analytical formula Eq.~\ref{eq:main_res}.}
\end{figure}

\begin{figure}[t]
\centering
\includegraphics[width=0.6\textwidth]{plot4.pdf}
\caption{\label{fig:google_plots_withfields} Expectation-value landscapes obtained for problem instances with $h_i, J_{ij} \in \{-1,1\}$, using the analytical formula Eq.~\ref{eq:main_res}.
Problem graphs are shown above, with blue nodes representing those with $h_i=-1$, and red for those with $h_i=1$.
(A) Sycamore grid problem with $23$ qubits.
(B) $3$-regular graph with $22$ qubits.
(C) Fully connected SK-model with $17$ qubits.}
\end{figure}

\begin{figure}[t]
\centering
\includegraphics[width=0.6\textwidth]{plot5.pdf}
\caption{\label{fig:google_plots_withfields2} Expectation-value landscapes obtained for problem instances with $h_i, J_{ij} \in \{-2,2\}$, using the analytical formula Eq.~\ref{eq:main_res}.
Problem graphs are shown above, with blue nodes representing those with $h_i=-2$, and red for those with $h_i=2$.
(A) Sycamore grid problem with $23$ qubits.
(B) $3$-regular graph with $22$ qubits.
(C) Fully connected SK-model with $17$ qubits.}
\end{figure}

\begin{figure}[t]
\centering
\includegraphics[width=0.6\textwidth]{plot6.pdf}
\caption{\label{fig:google_plots_hgeJ} Expectation-value landscapes obtained for problem instances with $h_i \in \{-2,2\}$ and $ J_{ij} \in \{-1,1\}$, using the analytical formula Eq.~\ref{eq:main_res}.
Problem graphs are shown above, with blue nodes representing those with $h_i=-2$, and red for those with $h_i=2$.
(A) Sycamore grid problem with $23$ qubits.
(B) $3$-regular graph with $22$ qubits.
(C) Fully connected SK-model with $17$ qubits.}
\end{figure}

\begin{figure}[t]
\centering
\includegraphics[width=0.6\textwidth]{plot7.pdf}
\caption{\label{fig:google_plots_hgeJ10} Expectation-value landscapes obtained for problem instances with $h_i \in \{-10,10\}$ and $ J_{ij} \in \{-1,1\}$, using the analytical formula Eq.~\ref{eq:main_res}.
Problem graphs are shown above, with blue nodes representing those with $h_i=-10$, and red for those with $h_i=10$.
(A) Sycamore grid problem with $23$ qubits.
(B) $3$-regular graph with $22$ qubits.
(C) Fully connected SK-model with $17$ qubits.}
\end{figure}

\begin{figure}[t]
\centering
\includegraphics[width=0.6\textwidth]{plot8.pdf}
\caption{\label{fig:google_plots_hleJ} Expectation-value landscapes obtained for problem instances with $h_i \in \{-1,1\}$ and $ J_{ij} \in \{-2,2\}$, using the analytical formula Eq.~\ref{eq:main_res}.
Problem graphs are shown above, with blue nodes representing those with $h_i=-1$, and red for those with $h_i=1$.
(A) Sycamore grid problem with $23$ qubits.
(B) $3$-regular graph with $22$ qubits.
(C) Fully connected SK-model with $17$ qubits.}
\end{figure}

\begin{figure}[t]
\centering
\includegraphics[width=0.6\textwidth]{plot9.pdf}
\caption{\label{fig:google_plots_hleJ10} Expectation-value landscapes obtained for problem instances with $h_i \in \{-1,1\}$ and $ J_{ij} \in \{-10,10\}$, using the analytical formula Eq.~\ref{eq:main_res}.
Problem graphs are shown above, with blue nodes representing those with $h_i=-1$, and red for those with $h_i=1$.
(A) Sycamore grid problem with $23$ qubits.
(B) $3$-regular graph with $22$ qubits.
(C) Fully connected SK-model with $17$ qubits.}
\end{figure}

\begin{figure}[t]
\centering
\includegraphics[width=0.4\textwidth]{plot10.pdf}
\caption{\label{fig:google_plot_hJwild10} Expectation-value landscape, obtained using the analytical formula Eq.~\ref{eq:main_res}, for a 3-regular graph problem instance that has $16$ qubits, and $h_i \in \{-30, -29,\dots,29,30\}$ and $J_{ij} \in \{-30,30\}$.}
\end{figure}

\subsection{Numerical Results on Google \textsc{Max}-\textsc{Cut}\xspace and Sherrington-Kirkpatrick (SK) Instances}

Eq.~\ref{eq:main_res} allows us to replicate the expectation-value landscapes of the recent experimental work of Ref.~\cite{google_qaoa}, which provides numerical support for the validity of our analytical results.
In this section, we present our investigation of the problems studied in Ref.~\cite{google_qaoa}, as well as some similar but larger problem instances.

The problems studied in \cite{google_qaoa} were weighted \textsc{Max}-\textsc{Cut}\xspace instances with weights $-1$ and $+1$ assigned at random; we obtained the three specific instances used in \cite{google_qaoa} from the \emph{ReCirq} repository~\cite{ReCirq}.
We numerically evaluated Eq.~\ref{eq:main_res} for these three problem instances. This gave us the energy-expectation landscapes for the single-layer \textsc{qaoa}\xspace for these instances, i.e., the numerical evaluation of the function $F(\beta,\gamma)$ for varying values of both $\beta$ and $\gamma$. We have plotted these energy-expectation landscapes in Fig.~\ref{fig:google_plots}. We note that these landscapes, which we obtained using the analytical formula Eq.~\ref{eq:main_res}, are in excellent agreement with the brute-force-classical and quantum-hardware evaluations of the landscapes shown in \cite{google_qaoa}. Observe, also, that we plot $-F(\beta,\gamma)$ due to the normalization of the previous evaluations.

The analytical formula Eq.~\ref{eq:main_res} allows us to explore the landscapes of problem instances far larger than those possible to study on current experimental hardware.
We have computed the landscapes for two $100$-spin (i.e., $100$-qubit) Ising problem instances: one instance has a $3$-regular graph structure, and the other is a Sherrington-Kirkpatrick instance. The graphs and their corresponding landscapes are shown in Fig.~\ref{fig:google_plots_morequbits}.

In relation with the previous section, note how the models shown in Fig.~\ref{fig:google_plots}.C and  Fig.~\ref{fig:google_plots_morequbits}.B correspond to the cases studied in Section~\ref{sec:SK_maxcut_expectation}. We predicted, as also shown in Ref.~\cite{farhi2019quantum}, that $\sqrt{n} \gamma_\textrm{min} = 0.5$ for this bimodal distribution, and how that quantity will remain constant for any SK instances. For $n=11$ we obtain $\sqrt{11} \gamma^{(11)}_\textrm{min} \approx \sqrt{11} \cdot 0.158 \approx 0.5$, and for $n=100$, $\sqrt{100} \gamma^{(100)}_\textrm{min} \approx 10 \cdot 0.0527 \approx 0.5$, confirming our theoretical prediction.

We will now comment briefly on the features that one can observe in the landscapes shown in Fig.~\ref{fig:google_plots} and Fig.~\ref{fig:google_plots_morequbits}. In Fig.~\ref{fig:google_plots}, the peaks are clearly narrower and closer to the edge for the $11$-spin SK instance than for the two ($23$-spin and $14$-spin) sparsely-connected instances. Furthermore, for the SK instance, additional peaks appear at $\gamma=\pi/4$. These peaks will vary in height depending on the distribution of the $\pm 1$ values in the graph.

In Fig.~\ref{fig:google_plots_morequbits}, we see that the landscape for the sparsely-connected ($3$-regular) instance with $100$ spins is the same as the landscape for the $3$-regular instance with $14$ spins---i.e., as the number of spins has increased, the landscape has not changed. In contrast, the landscape for the $100$-spin SK instance is qualitatively different than the landscape for the $11$-spin SK instance: the peaks in the landscape for the $100$-spin instance are narrower and closer to the edges of the $\gamma$ range, and the peaks at $\gamma=\pi/4$ for the $11$-spin instance are not present for the $100$-spin instance. The peaks at $\gamma=\pi/4$ also become less pronounced as we increase the size of the SK graph. From the results derived in theory Sections~\ref{sec:SK_maxcut_expectation} and \ref{sec:reg_graphs_expectation}, we showed above that $\gamma_\textrm{min}$ will remain constant for $3$-regular graphs, while decreasing with the number of spins for SK instances. The landscapes that we observe are a direct effect of this phenomena.

An SK instance with $n$ spins is an $n$-regular graph. We can understand the narrowing and shifting of the peaks between the $n=11$ and $n=100$ SK instances we studied as an example of a more general phenomenon for regular graphs: the larger the degree of a regular graph, the narrower the peaks in the landscape will be.
This is explored in Sec.~\ref{sec:analytical_landscapes}.
This phenomenon has an important experimental consequence: the peak sharpness places a constraint on how accurately the angle $\gamma$ needs to be implemented in hardware, and so the necessary precision for specifying $\gamma$ will increase with the {\it degree} of a regular graph.

While the expectation-value landscapes presented in Fig.~\ref{fig:google_plots} already show some qualitative differences from instance to instance, we will show later in this paper that landscapes can vary far more dramatically depending on the parameters of the cost Hamiltonian.
Furthermore, we will study how the optimal $\gamma$ value, $\gamma_\textrm{min}$, that minimizes the expectation value, changes with those parameters.
So far we have seen evidence that increased degree of regular graphs can cause $\gamma_\textrm{min}$ to be closer to $0$ or $\pi/2$, but we will find that the value of $\gamma_\textrm{min}$ can have a rather complex dependence on the cost-Hamiltonian parameters.

One of the potential benefits of using an analytical formula to determine the optimal angles instead of performing a brute-force search is speed, since the former requires on classical-computing calculations whereas the latter uses quantum simulations or runs on quantum hardware. This potential benefit is seen already with the instances studied in this section. To illustrate this, we simulated the expectation-value landscape of fully connected graphs with $23$ qubits using \emph{Cirq}~\cite{Cirq}; these calculations took approximately $1$ minute to complete.
In comparison, using the analytical formula we could compute the landscape for $100$-qubit problems in approximately $4$ seconds.
Both calculations were performed on a \emph{c5n.large} AWS instance.

\subsection{Numerical Results on G-Set \textsc{Max}-\textsc{Cut}\xspace Instances and Ising Instances with External Fields}

The G-set \cite{G_set} collection of benchmark instances for \textsc{Max}-\textsc{Cut}\xspace is widely used in the comparison of classical \textsc{Max}-\textsc{Cut}\xspace heuristics \cite{Benlic_2013,leleu2019destabilization,toshiba_benchmarking}. We endeavored to predict how well the single-layer \textsc{qaoa}\xspace would perform on these instances, which affords a comparison between the single-layer \textsc{qaoa}\xspace running on a (future) quantum computer and the best-known classical methods for heuristically solving \textsc{Max}-\textsc{Cut}\xspace problems.

We calculated the optimal \textsc{qaoa}\xspace expectation values for a subset of the G-set instances having $\num{800}$, $\num{2000}$, $\num{5000}$, or $\num{7000}$ vertices and up to $\num{41459}$ edges.
The results are divided into two tables; Table~\ref{table:G-set_Ising}, where results are expressed in terms of Ising energies, and Table~\ref{table:G-set_cut}, where we show them as cut values. The instances have weights that are drawn from two possible sets, depending on the instance: $\{0,1\}$ and $\{-1,0,+1\}$. For Table~\ref{table:G-set_Ising}, we derive in the appendix an analytical formula that allows us to estimate the optimal energy of a G-set instance, $E_\textrm{opt} \approx 1.18 \sqrt{|E|\times |V|}$. The ratio between the theoretical maximum and the \textsc{qaoa}\xspace expectation, $r_\textrm{exp.}=F_\textrm{max}/E_\textrm{opt}$, ranges between $0.252-0.389$. This ratio is similar to the formula of the approximation ratio, $r=F_\textrm{max}/C_\textrm{max}$, that is frequently used in the literature to determine the efficiency of the \textsc{qaoa}\xspace~\cite{farhi2014quantum,google_qaoa, wang2018quantum}.
We just substitute the real maximum of the cost function by our theoretical approximation.
The ratio for the best-known classical solution to the \textsc{qaoa}\xspace expectation, $r_\textrm{Ising}=F_\textrm{max}/E^{cl}_\textrm{max}$, ranges between $0.36-0.475$. Note how the theoretical maximum is close to the classical best-known solution for most instances except G14-G17. For this set the \emph{Ratio Ising} and \emph{Ratio cut} are also particularly high. This leads us to believe that the best-known classical solution is not as good for these particular instances. This same logic applies to other instances with high values of \emph{Ratio Ising} e.g.\ G11. Additionally, the \textsc{qaoa}\xspace expectation value is expected to be higher for instances with small amount of edges and common neighbours between vertices. These ratios can be used as a method to quickly determine the \textsc{qaoa}\xspace expectation of a given instance; one can use the analytical formula to determine the theoretical maximum, and use the ratio range to determine the \textsc{qaoa}\xspace expectation. The expectation value landscapes are very similar for the G-set instances, with similar optimal angles as predicted in \cite{brandao2018fixed}. Fig.~\ref{fig:G_set} shows the landscape of instance G64.

We also calculated the best \textsc{qaoa}\xspace expectation values for $7$ Ising instances in Table~\ref{table:G-set_ising_field}.
These instances were generated by adding weight $\{-1,+1\}$ external fields at random to instances of the G-set.
We denote these instances as I$x$ where $x$ is the G instance label on which it is based.
Additionally, two more instances, I99 and I100, were custom generated by assigning $\{-1,+1\}$ weights at random.
The number of vertex and edges goes up to $\num{100000}$ and $\num{150000}$, respectively, in I100.
An approximation for the optimal Ising energy of these instances was obtained using the analytical formula $E_\textrm{opt} \approx 1.18 \sqrt{|V| \times (|E|+|V|)}$, that is derived in the appendix.
We obtained a ratio $r_\textrm{exp.}$ that ranges between $~0.318-0.39$ when comparing the theoretical optima and the \textsc{qaoa}\xspace expectation.
The values of the $r_\textrm{exp.}$ ratios are very similar to those obtained for \textsc{Max}-\textsc{Cut}\xspace instances.
The \emph{Time} column displays in seconds the time required to obtain the full expectation value landscape and optimal angles.

\subsection{Numerical and Analytical Results on Ising Problem Expectation Value Landscapes}
\label{sec:analytical_landscapes}

In this section we will show the expectation-value landscapes of simple Ising problem instances with different weights, and derive analytic expressions that describe some of the observed phenomena.

First, we add external fields and study the case in which $h_i$ and $J_{ij}$ take only the values $\{-h,h\}$ assigned at random, where $h$ can take any (real) value, and the spin connectivity is that of a $d$-regular graph.
For this particular case we can approximate $\gamma_\textrm{min}$ as
\begin{align}
\gamma_\textrm{min} & \approx \frac{1}{2h}\arctan{\frac{1}{\sqrt{d}}},\frac{1}{2h}( \pi - \arctan{\frac{1}{\sqrt{d}}}) ~.
\end{align}
We previously observed how the regularity of the graph affects the $\gamma_\textrm{min}$ value, and now we see that increasing the $h$ value will also bring $\gamma_\textrm{min}$ closer to the edges ($\gamma = 0$ and $\gamma = \pi/2$). The $h$ value also affects the $\gamma$-periodicity of the landscape: the period of the oscillations in the landscape as a function of $\gamma$ is proportional to $1/h$. For example, peaks that for $h=1$ were located at $\gamma=\pi/2$, for $h=2$ will be located at $\gamma=\pi/4$. The expectation value of $\gamma_\textrm{min}$ for the single-node contribution in Eq.~\ref{eq:main_res} reads,
\begin{align}
\langle C_{i} \rangle(\beta_\textrm{min}, \gamma_\textrm{min}) & =\frac{\sqrt{d}}{d+1}(1+\frac{1}{d})^{-(d-1)/2}.
\end{align}
This value decreases monotonically with $d$. Here we operate under the assumption that the single-node term will be the main contribution to the expectation value.
We observe the predicted phenomena in Fig.~\ref{fig:google_plots_withfields} and \ref{fig:google_plots_withfields2} where we plot the expectation landscapes for the cases $h=1$ and $h=2$, respectively. The first plot is, as expected, similar to what we encountered in the absence of external fields. Due to the presence of local minima and maxima with different expectation values, and additional features, the $\beta$ values of the plots now range between $-\pi/2$ and $\pi/2$. For $h=2$ we observe that the periodicity in $\gamma$ is divided by half, and there is a narrowing of the minimum and maximum peaks in $\gamma$. We had already seen (e.g., in Figures \ref{fig:google_plots} and \ref{fig:google_plots_morequbits}) that increasing number of edges per node causes the peaks in the landscape to become narrower in $\gamma$, and we now see that the peaks also become narrower when $h$ is increased.

The difference in expectation values for different minima is going to arise from the summation terms in Eq.~\ref{eq:main_res}. For the general case these terms have the form
\begin{equation}
[\cos(2 \gamma h)^o \cos(2 \gamma J)^p+\cos(2\gamma h)^q \cos(2 \gamma J)^r].
\end{equation}
If by our choice of $o,p,q,r,h,J$ parameters we can make one term positive and another negative for the same value of $\gamma$, this will cause different expectation values for different minima.

Consider the case of problem instances with $h_i \in \{-rJ,rJ\}$ and $J_{ij} \in \{-J,J \}$ assigned at random, where the graph is again $d$-regular. In this case $J_{ij} \le h_i$ for all $i$ and $j$ (since $r$ is an integer). The following are two different approximations for $\gamma_\textrm{min}$ one can make depending on the value of $r$. For $r \le 5$,
\begin{align}
\gamma_\textrm{min} & \approx \frac{1}{2J}\arctan{(\frac{1}{\sqrt{d}}(1-0.15(r-1)))}.
\end{align}
For $r > 5$,
\begin{align}
\gamma_\textrm{min} & \approx \frac{1}{2r^{0.6}J}\arctan{(\frac{1}{\sqrt{d}})}.
\end{align}
These equations show that the relative increase of $h_i$ with regards to $J_{ij}$ is going to make $\gamma_\textrm{min}$ closer to the edges. This phenomenon is observed in the $r=2, 10$ cases shown in Fig.~\ref{fig:google_plots_hgeJ} and \ref{fig:google_plots_hgeJ10}, respectively. The $r=2$ case shows the peaks taking different shapes to what we previously observed. And for $r=10$, the generation of additional peaks at different locations in the landscape, and the narrowing of those. By evaluating the expectation value in $\gamma_\textrm{min}$ we observe that it decreases monotonically with $r$. Thus, an increase in $r$ will not only reduce $\gamma_\textrm{min}$ but it will generate additional peaks, the narrowing of those, a reduced optimal expectation value, and generate new shapes.

Alternatively, we can study the case $J_{ij} \ge h_i$, where we set $h_i \in \{ -h,h \}$ and $J_{ij} \in \{ -rh,rh \}$ assigned at random in a $d$ regular graph. For this case we can approximate $\gamma_\textrm{min}$,
\begin{align}
\gamma_\textrm{min} & \approx \frac{1}{2rh}\arctan{\frac{1}{\sqrt{d}}}
\end{align}
This equation shows how an increase in $r$ will again lead to $\gamma_\textrm{min}$ being closer to the edge. The expectation value of $\gamma_\textrm{min}$ for the single node contribution Eq.~\ref{eq:main_res} reads,
\begin{align}
\langle C_{i} \rangle(\beta_\textrm{min},\gamma_\textrm{min}) & =(1+\frac{1}{d})^{-\frac{d}{2}} \sin(\frac{\arctan(1/\sqrt{d)})}{r})~.
\end{align}
This expression decreases monotonically with $r$, as well as $d$, leading to a worse expectation value with increasing $r$.

Fig.~\ref{fig:google_plots_hleJ} and \ref{fig:google_plots_hleJ10} show the $r=2, 10$ cases respectively. The $r=2$ case shows additional peaks, narrowing of those, and new shapes. While for $r=10$ we have a complex landscape with a big number of narrow peaks with different shapes. Thus, an increase of $r$ for these kind of problem instances will lead to a $\gamma_\textrm{min}$ closer to the edge, a reduced optimal expectation value, additional peaks, narrowing of these, and new shapes.

All these phenomena are very similar to those observed in the $J_{ij} \le h_i$ case. However, unlike previously, the number of additional peaks is larger, the narrowing is much more noticeable, and the new shapes are easily observable. Thus, the increase of the relative $J_{ij}$ values has a bigger impact in generating wild expectation value landscapes. Notice also that the expectation value of some of these additional peaks is very similar to that of the optima. For the instance displayed in Fig.~\ref{fig:google_plots_hleJ10} the peaks are already very narrow and difficult to distinguish. This can cause difficulties when trying to set the optimal angles in an experimental setup due to precision issues. Our method can be useful in this case to identify the most suitable set of optimal angles that are experimentally more feasible.

We have shown that varying the parameters of the cost Hamiltonian the expectation value landscapes can look very different to the ones presented at the beginning of this paper; one example of this is shown in  Fig.~\ref{fig:google_plots_hleJ10}. This means that for certain Ising instances the computation of optimal angles is not straightforward and further study is required. As another example, consider the class of Ising instances where $h_i$ can take any integer value at random between $h$ and $-h$, i.e., $h_i  \in \{-h,\dots,h\}$, and $J_{ij} \in \{ -h, h \}$. The landscape of one such instance is shown in Fig.~\ref{fig:google_plot_hJwild10} for $h=30$. Here we no longer have any of the recognisable patterns that we had in the first landscapes, and the optimal angles are not at the edges.

More general and practically relevant problems may have complicated landscapes and require their own study. As we have shown, \textsc{qaoa}\xspace landscapes can vary substantially depending on the class of problem instance, and finding the optimal angles using classical methods requires some care to be taken.

\section{Conclusions}

We introduced an analytical formula for the expectation value of single-layer ($p=1$) \textsc{qaoa}\xspace circuits solving general Ising Hamiltonians, which allows for an efficient classical calculation of the optimal angles $\beta_\textrm{min}$, $\gamma_\textrm{min}$.
We validated our analytical formula by numerically evaluating the \textsc{qaoa}\xspace expectation values of several Ising problems studied in a recent experimental paper \cite{google_qaoa}.
Furthermore, we predicted the energy landscapes for single-layer \textsc{qaoa}\xspace with $100$ qubits, and benchmarked the expectation values for the G-set \textsc{Max}-\textsc{Cut}\xspace problems and custom Ising problems with hundreds to thousands of vertices/spins/qubits.
For the G-set benchmark instances the ratio between the theoretical estimated maximum Ising energy and the \textsc{qaoa}\xspace expectation for Ising energy was in the range $~0.252-0.389$; for the Ising problems with external fields that we tested, the ratio was in the range $~0.318-0.39$.
In order to compute these ratios---since we don't know the true optimal energy for each of such large Ising instances---we derived an analytical formula that allows one to estimate the optimal energy of an Ising or \textsc{Max}-\textsc{Cut}\xspace instance with weights $\{-1,0,+1\}$.
The largest instance that we estimated the accuracy of the single-layer \textsc{qaoa}\xspace for had $\num{100000}$ vertices and $\num{150000}$ edges, and found that it achieved a ratio of $\approx 0.39$.

For the G-set \textsc{Max}-\textsc{Cut}\xspace instances, we also computed the ratio between the expected cut value from the single-layer \textsc{qaoa}\xspace and the best-known cut value (obtained using classical heuristic methods) for each instance.
For the G-set instances that are {\it unweighted} (i.e., $J_{ij} \in \{0,1\}$), the ratios are all $\geq 0.85$---meaning that we can expect the \textsc{qaoa}\xspace to find \textsc{Max}-\textsc{Cut}\xspace solutions for these instances that have at least $85\%$ of the maximum (known) cut value.
The accuracy of the solutions one would obtain by running the single-layer \textsc{qaoa}\xspace on such large instances in practice (under the assumption of noiseless quantum hardware) should be even higher, because each run of the \textsc{qaoa}\xspace will yield a single sample from a distribution that has a {\it mean} solution value with $\sim 85\%$ accuracy, and so multiple runs of the \textsc{qaoa}\xspace are likely to yield some samples that have solution values better than that of the mean value. Two related open questions of great interest are {\it what the shape of the output-sample distribution from the \textsc{qaoa}\xspace is} and {\it what the probabilities of measuring solutions with energies much better than the mean energy are} \cite{pagano2020quantum}.

We also note that the multi-layer \textsc{qaoa}\xspace cannot perform worse than the single-layer \textsc{qaoa}\xspace, since one can always choose to set the angles for all layers after the first one to be equal to $0$, and hence have the layers after the first layer perform the identity operation., so our results provide a lower bound for the performance of the multi-layer \textsc{qaoa}\xspace.

While the numerical evaluation of our analytical formula has polynomial scaling, the wall-clock times for our current classical code implementation are typically $\gg 1000$~seconds for Ising instances comprising $|V| \gg 100$ spins (exact times are given in Tables \ref{table:G-set_cut} and \ref{table:G-set_ising_field}). An area for future work is to increase the speed of the classical code for computing optimal angles.

Finally, we studied the effect that the values of the external-field terms $h_i$ and spin-spin-coupling terms $J_{ij}$ have on the expectation-value landscape and the optimal $\gamma$ parameter, $\gamma_\textrm{min}$. We showed that an increase of the external fields or spin-spin couplings appear, in general, to lead to the generation of additional peaks in the landscape, the narrowing of these peaks, new shapes, a reduction in optimal expectation value, and $\gamma_\textrm{min}$ being closer to the edges. A relative increase in $J_{ij}$ has a bigger impact in the generation of complex expectation landscapes. We found that a broad class of Ising Hamiltonians have highly complex landscapes that are very dissimilar to the ones that we are aware of in the literature. Our analytical methods for finding optimal angles may be particularly useful in these cases, especially since angle searches performed in situ using quantum hardware may easily ``miss'' very narrow features in the landscapes.

{\it Note added in preparation:} Related work studying the single-layer \textsc{qaoa}\xspace numerically appeared on the arXiv on 29 November 2020: \cite{bravyi2020hybrid}, which investigates the performance of the \textsc{qaoa}\xspace on a particular subclass of MAX-CUT instances related to graph coloring.

\section{Acknowledgements}

The authors thank V.~Putz, J.\,I.~Adame, R.~Parrish and E.~Farhi for useful discussions, and
J.\,I.~Adame for a thorough reading of a draft of this paper. PLM acknowledges membership of the CIFAR Quantum Information Science Program as an Azrieli Global Scholar.

\bibliographystyle{ieeetr}


\section{Introduction}
Training deep neural networks is well known to be highly compute-intensive, with the multiplications typically being the most demanding in terms of the required hardware resources.
For this reason, various reduced-precision methods have been proposed for neural network training\citep{Dett, Gins, Mici, Cour}.
The authors of \citep{Abad} and \citep{Dean} introduced BFLOAT16 format (1+8+7 format).
Despite these proposals, the multipliers still take the largest part of the circuit

AdderNet \citep{Chen} reduced the number of multiplication by replacing the convolutions with $\ell_1$ norms, which only uses addition and absolute value, but it showed a mild ($\sim 1\%$) accuracy degradation. It is unclear how to generalize this method to layers other than convolution and fully-connected layers, which may limit applicable neural network architectures. It also required a special training procedure to make it stable.

The logarithmic number representation (LNR) is shown to give better accuracy than linear number representation with the same bit-width for deep neural networks \citep{miyashita2016convolutional, lee2017lognet}. These methods are mainly used for inference and have not been used for training of large networks and datasets such as ImageNet
(Our proposal may be considered to be a natural generalization of LNR, which linearly interpolates it, and makes the training of large networks possible.)

To avoid hardware-intense multiplication, \citep{kim2018efficient} proposed using Mitchell's approximate multiplication \citep{mitchell1962computer} to convolutional neural networks and had shown only a slight ($\sim 0.2\%$) degradation in accuracy for VGG16, but the authors did not suggest using it for training.

In this article, we propose replacing the multiplication with integer addition (detail in Figure~\ref{list1}), and call it ‘addition-as-int’.
We point out that adding two IEEE 754 floating-point representations with an integer addition instruction actually results in Mitchell's approximate multiplication.
We show that this operation may be used as a drop-in replacement of floating-point multiplication. We found only a small ($\sim 0.5\%$) or no accuracy degradation in training ResNet \citep{He} on the ImageNet dataset.

In addition to using our addition-as-int as an approximation to multiplication, we propose using it as a different but exact operation, which is expected to be more stable.

\begin{figure}[t]
\begin{center}
\begin{tabular}{l}
\hline \\
{\tt float int\_mul(const float a, const float b)\{} \\
{\tt\hspace{2em} int c = *(int*)\&a + *(int*)\&b - 0x3f800000;} \\
{\tt\hspace{2em} return *(float*)\&c;} \\
{\tt\}} \\
\hline
\end{tabular}
\end{center}
\caption{\label{list1}
Our pseudo-multiplication in C++. This function gives a close ($\sim 12.5\%$) approximation of the FP32 multiplication.
In detail, the sign bit is also computed correctly with this simple code if no underflow happens. We should supply rounding to zero if the underflow of the exponent is possible.
Exponent bias is 127 in IEEE754 representation. Then the exponent bias in 32bit integer representation is $127\times 2^{23} = {\tt 0x3F800000}$, since the LSB of the exponent is at the 23rd bit.}
\end{figure}

\section{Addition-as-int as float multiplication}
How can adding IEEE754 floating-point representations as mere integers (addition-as-int) give a close approximation to multiplication?
Intuitively, the reason is that it gives exact multiplication when either of the multiplicands is powers of 2. Otherwise, it will give linear interpolation between those exact values because the result is linear in the mantissa. As the result, the error is only 12.5\% in the worst case. ($1.5 \times 1.5$ gives $2$ in this operation, while $(1.5)^2$ is $2.25$.)

We may more easily understand our addition-as-int in terms of logarithms and exponentials. Let us consider two functions defined by
$e(x) =	1 + x$ for $0 =< x < 1$, otherwise
$e(x+1) = 2 e(x)$,
and
$l(x) = x - 1$ for $1 =< x < 2$, otherwise
$l(2x) = l(x) + 1$,
which give exponential of base 2 for integer $x$ and logarithm of base 2 when x is a power of 2 respectively, and give linearly interpolated values in between. The function $e(x)$ may easily be implemented reading an FP32 (IEEE754) representation into an integer register, subtracting the exponent bias, and considering the fixed decimal point to be at the left of the 22nd bit. The function $l(x)$ may easily be implemented by doing the same in the opposite direction. Then the addition-as-int as a substitute of multiplication of $a$ and $b$ may be obtained as $    f(a, b) = e(l(a) + l(b))$.

Backpropagation should be calculated using
$	\partial f(a, b) / \partial a = e’(l(a) + l(b)) l’(a),$
where the derivative $e’$ and $l’$ may easily be calculated by zeroing the mantissa. One drawback in implementing this operation is that when a forward path is $y = ab$ then backward gradient $da$ requires not only $b$ and $dy$ (like $da = b\ dy$) but also $a$, which results in an increase of bandwidth about $a$.
The ``a-operations", which will be introduced later, does not have this drawback.

In the light of logarithm and exponential, a linear layer of neural networks
$	y_i = \sum_j w_{ij} x_j$
may be considered as
$	y_i = \sum_j \exp( \log w_{ij} + \log x_j ),$
which is a log activation layer before addition before an exp activation layer.  Therefore replacing exp and log with our $e$ and $l$ function is just like changing the activation functions a little, and we may expect that the change of the result should be as little as replacing the ReLu activation function with another activation function.

\section{Experimental settings}
We used Chainer \citep{Toku} and CuPy \citep{Okut} to implement the layers and Nvidia PASCAL P100 GPUs to run the experiments.
The ResNet-50 was trained on ImageNet (ILSVRC 2015) dataset for 100-epoch as our benchmark. The standard augmentation methods, i.e. flip, crop, resize, were applied as \citep{He}. The training hyperparameters are lr=0.05*32/256, batchsize=32, the learning rate was decreased at 30, 60, and 90 epoch by the factor of 10, weight decay rate was $10^{-4}$, and momentum SGD with $\mu=0.9$, which are standard. All the FP32 baseline and other training are done under the same random seed.

In the MNIST training, a MLP of linear(784, 1000)-ReLu-linear(1000,1000)-ReLu-linear(1000,10) structure with an ordinary softmax-cross entropy layer was used without augmentation with Adam for 20 epochs with batchsize=100.

In implementing softmax cross-entropy, we have to be careful about the shift-invariance of the logsumexp function ${\rm lse}(\vec x)$ in it.
We have relationship ${\rm lse}(\vec x + a) = {\rm lse}(\vec x) + a$, and then we usually subtract the maximum element of $\vec x$ from $x_i$ before the calculation and add it afterward. In our approximate exponential function, this invariance holds only for integer $a$, and the implementation should respect that. We observed performance degradation without this consideration.

\begin{table}[t]
\caption{Top-1 accuracy (\%) on ImageNet of ResNet-50 with e-operations, a-operations.
The meaning of the abbreviation like ``cE" is explained at the beginning of the result section.}
\label{table1}
\begin{center}
\begin{tabular}{llll}
\toprule
\multicolumn{2}{c}{\bf e-operations}  &\multicolumn{2}{c}{\bf a-operations} \\
\cmidrule(r){1-2} \cmidrule(r){3-4}
baseline(FP32)  & 75.6 & & \\
cE  & 75.1 & ca.fE  & 75.6\\
cE.fE  & 75.0 & cE.fa  & fail\\
cE.fE.bE  & 73.6 & ca.fE.ba  & 75.5\\
cE.fE.bE.eE  & 73.9 & ca.fE.ba.ea  & 75.6\\
\bottomrule
\end{tabular}
\end{center}
\end{table}

\section{Results}
We present results for two variants of our addition-as-in (Figure~\ref{list1}). The first is the exact operation (``e-operations") scheme  and the second is the approximate operation (``a-operations") scheme.
Also, the results of using BFLOAT16 together will be presented.

We denote the network configuration as follows: `cE' denotes using e-operations for all the convolution layers, `ca' denotes using a-operations for all the convolution layers, and ordinary multiplication is used if none is specified. Other than `c' for convolution layers, `f' denotes a fully-connected layer, `b' for batch normalization layers, and `e' for a softmax cross-entropy layer. Thus `cE.ba' denotes that the network has convolution layers with e-operations, batch normalization layers with a-operations, and other layers with ordinary definition, for example.

\subsection{Exact operations}

Here we use our addition-as-int instead of multiplication and consider the exact differentiation of this function. We call such operations by ``e-operations", and each of the operations will be called ``e-mult", ``e-exp", ``e-div", and so on.
This e-mult is not approximating multiplication but replacing multiplication with a similar but different operation.

For the MNIST dataset, applying e-mult to the simple MLP gave the accuracy of 0.980 and 0.982 for two runs of different initialization, where the ordinary multiplication gave 0.981 and 0.981, which may be considered to be the same.

Then we began ImageNet training with replacing the multiplications with our e-mult in all the convolution layers. We saw only a little ($\sim 0.5\%$) degradation in classification accuracy (cE in Table~\ref{table1}). Further replacing the multiplications in the fully-connected layer did not change the test accuracy (cE.fE in Table~\ref{table1}) from the cE case.

\subsection{Approximate operations}

On the other hand, we may also use ${\rm int\_mul}(x, y) + \gamma$ as an approximation to multiplication, where the correction constant $\gamma$ will be explained shortly. The gradient of multiplication is another multiplication, and it can again be approximated with ${\rm int\_mul}(x, y)+\gamma$.
We call this use of operations as “a-operations”, and each of the operations will be called as “a-mult”, “a-exp”, “a-div”, and so on.

The constant $\gamma$ was determined as follows. If we want to use $l(x)$ as an approximation to $\log_2(x)$, we need to use $l(x)+\gamma$  because $l(x)$ is consistently below $\log_2(x)$. Likewise, the same applies to $e(x+\gamma)$ as an approximation to $2^x$.  We assume the distribution of the incoming numbers are flat and the average difference of $\log_2(x)$ and $l(x)$ should be canceled, so we have
$	\gamma = \int_1^2 (\log x - (x-1)) dx = 3/2 - 1/\log 2.$
We found this assumption to be sufficient for training neural networks

We first replaced multiplications of the convolution layers with a-mults but kept e-mult for the fully-connected layer of ResNet-50. We observed equal accuracy (ca.fE in Table~\ref{table1}) in reference to the baseline FP32 model, which is better than e-operations. Replacing with a-mult for the fully-connected layer decreased accuracy more than 30\% at 10th epoch and no further training was performed. From this, we considered using a-operation in fully-connected layers is not good.

\subsection{Replacing all the other layers with our operations}

In addition to convolution and fully-connected layers, we replaced operations in batch normalization with e-operations.
We saw mild accuracy degradation ($\sim 1.4\%$) without other trouble (cE.fE.bE in Table~\ref{table1}). We further replaced the softmax cross-entropy layer using e-log, e-exp and observed no further drop of accuracy (cE.fE.bE.eE in Table~\ref{table1}).
The experimental result for a-operations gave the same accuracy as FP32. We replaced the batch normalization using a-mult to find no difference (ca.fE.ba in Table~\ref{table1}) from FP32. Additionally, we replaced the softmax cross-entropy layer using a-operations to find equal accuracy (ca.fE in Table~\ref{table1}) as FP32.
Therefore, if one is considering new hardware that lacks multipliers, a-operation seems better for batch normalization.

\begin{figure}[t]
\begin{center}
\includegraphics[width=0.5\linewidth]{bfloat.pdf}
\end{center}
\caption{\label{bfloat-lc}
Learning curves using BFLOAT16 for ImageNet.
We see no deviation from the FP32 result.
``ca.fE" means convolution layers with a-mult and a fully connected layer with e-mult.
}
\end{figure}

\subsection{Compatibility with a 16-bit float format}

Is it possible for our operation to cope with reduced-bit representation such as BFLOAT16?
If it is successful we would also obtain the merits of both ours and BFLOAT16's. Methodologically, we did almost the same as \citep{Kala}, i.e. pseudo multiplications were done for and resulted in BFLOAT16, then they are accumulated in FP32, and weights were kept in FP32. Our method was only different in that weight gradients were rounded after accumulation before adding to the master weights

As the result of ImageNet training, our BFLOAT16 baseline implementation did show the same accuracy (75.7\%) as the FP32 baseline, and we saw equal accuracy (75.7\%) for our approximate operations (ca.fE scheme) compared from the FP32 baseline (Figure~\ref{bfloat-lc}).
Here we showed that our method works with BFLOAT16 without any accuracy degradation.

\section{Conclusion and discussion}
In this article, we showed that replacing floating-point multiplications with integer additions gives competitive or sometimes equal performances as the FP32 baseline in training ResNet-50. We showed that all the multiplications in all the layers of the neural network can be replaced with integer additions. The experimental results also suggest that we may freely combine layers with our addition-as-int and ordinary layers. In addition, we showed that our method can work with a shorter floating-point format, i.e. BFLOAT16.

The future hardware that has fused instruction for our operation will run faster, though our integer-add-float-add operation ran slower for the current GPU, which is highly optimized for the fused-multiply-add.

Our results raise the possibility of further research on optimal activation functions. The multiplication may be considered as a combination of logarithmic activation function and exponential activation function. We have suggested only one way of changing these activation functions.
It is then interesting to consider future research directions that focus on finding the most hardware efficient activation functions that also attain the highest accuracy.

\begin{ack}
The author thanks Brian Vogel and Kentaro Minami for reading the manuscript and their suggestions, and Tanvir Ahmed for comments.
\end{ack}

\section*{Broader impact}
The impact of this study is expected to make deep neural-network training more efficient.
No one would get any disadvantage out of this study, but every person training neural networks will potentially benefit from improved efficiency, and every person using neural network inference can also benefit.
No failure modes nor leverage of bias specific to this method is expected to exist since it works equivalently to the existing methods as far as the authors know.

\bibliographystyle{plain}