††thanks: Former affiliation was IBM Quantum, IBM Japan.

# Recursive Quantum Relaxation for Combinatorial Optimization Problems

Ruho Kondo  Toyota Central R&D Labs., Inc., 41-1, Yokomichi, Nagakute, Aichi
480-1192, Japan  Quantum Computing Center, Keio University, 3-14-1 Hiyoshi,
Kohoku-ku, Yokohama, Kanagawa 223-8522, Japan<EMAIL_ADDRESS>Yuki
Sato  Toyota Central R&D Labs., Inc., 41-1, Yokomichi, Nagakute, Aichi
480-1192, Japan  Quantum Computing Center, Keio University, 3-14-1 Hiyoshi,
Kohoku-ku, Yokohama, Kanagawa 223-8522, Japan  Rudy Raymond  Quantum Computing
Center, Keio University, 3-14-1 Hiyoshi, Kohoku-ku, Yokohama, Kanagawa
223-8522, Japan  Department of Computer Science, The University of Tokyo,
7-3-1, Hongo, Bunkyo-ku, Tokyo 113-0033, Japan  Naoki Yamamoto  Quantum
Computing Center, Keio University, 3-14-1 Hiyoshi, Kohoku-ku, Yokohama,
Kanagawa 223-8522, Japan  Department of Applied Physics and Physico-
Informatics, Keio University, Hiyoshi 3-14-1, Kohoku-ku, Yokohama 223-8522,
Japan

###### Abstract

Quantum optimization methods use a continuous degree-of-freedom of quantum
states to heuristically solve combinatorial problems, such as the MAX-CUT
problem, which can be attributed to various NP-hard combinatorial problems.
This paper shows that some existing quantum optimization methods can be
unified into a solver that finds the binary solution that is most likely
measured from the optimal quantum state. Combining this finding with the
concept of quantum random access codes (QRACs) for encoding bits into quantum
states on fewer qubits, we propose an efficient recursive quantum relaxation
method called recursive quantum random access optimization (RQRAO) for MAX-
CUT. Experiments on standard benchmark graphs with several hundred nodes in
the MAX-CUT problem, conducted in a fully classical manner using a tensor
network technique, show that RQRAO outperforms the Goemans–Williamson method
and is comparable to state-of-the-art classical solvers. The codes will be
made available soon.

##  1 Introduction

Figure 1:  Concept of the proposed method, recursive quantum random access
optimization (RQRAO). (a) One/two/three classical bits are embedded into one
qubit using a quantum random access code (QRAC) [1, 2]. A total of
two/four/eight variants of the information of the one/two/three classical bits
are represented as a vector pointing to the corner of the line/square/cube
inscribed in the Bloch sphere. These specific quantum states are called (one-
qubit) magic states. (b) Using a high compression QRAC, the size of the
Hilbert space corresponding to the given problem is reduced and each (multi-
qubit) magic state is no longer necessarily orthogonal to the others. Note
that only four of the $2^{6}$ magic states are illustrated. After optimizing
the objective function defined on the reduced Hilbert space, the bit string is
decoded by searching for the nearest magic state in the fidelity, which is
called rounding. To accomplish the rounding without strong approximation, a
recursive strategy similar to the RQAOA [3] is adopted.

The MAX-CUT problem is one of Karp’s 21 NP-complete problems [4] and is
frequently used as a benchmark of combinatorial optimization algorithms
because of its simplicity. Namely, given a graph that consists of nodes
connected by weighted edges, the optimization variant of the MAX-CUT problem
asks for the labeling of nodes into two classes such that the sum of the
weighted edges connecting nodes with different labels is maximized. Obtaining
good solutions for the MAX-CUT problem is in great demand, because various NP-
hard combinatorial problems, e.g., partitioning, covering, and coloring
problems, can be relaxed to the MAX-CUT problem by treating the constraints as
penalty terms [5, 6]. Although the MAX-CUT problem is a discrete optimization
problem, it can be translated to a continuous variant by treating the discrete
solution in a probabilistic manner. Classically, such a continuous variant of
the MAX-CUT problem is relaxed to a semidefinite program by approximating the
probability [7]. However, the quantum optimization does not take such an
approximation because the probabilistic characteristics are naturally
incorporated into the quantum states. In the quantum MAX-CUT , all solutions
are encoded in the eigenstates of the problem Hamiltonian, and then the
quantum state is optimized such that the expected value of the problem
Hamiltonian is maximized. From a theoretical viewpoint, it has been proved
that the eigenstate with the maximum eigenvalue can be obtained by the
infinite-level quantum approximate optimization algorithm (QAOA) [8, 9] with
the appropriate parameters. However, in a noisy quantum device, it is
intractable to carry out even a sufficiently large level QAOA, let alone an
infinite level QAOA. For this reason, the field of combinatorial problems with
quantum computing is dominated by the development of heuristics, that mainly
focus on the development of the ansatz variants and improving the trainability
of QAOA [10].

Rather than developing the QAOA ansatz, an alternative approach is to develop
a non-diagonal Hamiltonian instead of the diagonal Ising-type one used in the
QAOA. In [11], the MAX-CUT Hamiltonian is constructed using quantum random
access codes (QRACs) [1, 2] (see Sec. 5.3 for a review). Using QRACs, the
$n$-qubit quantum state can be used to encode $m$ binary variables for $m>n$.
The QRAC Hamiltonian is still a 2-local Hamiltonian, and therefore it is QMA-
hard to find the ground state [12], but the QRAC formulation has the advantage
that the number of qubits can be reduced, making it easier to optimize the
variational state.

Recently, [3] proposed the recursive QAOA (RQAOA), showing that the recursive
application of level-1 QAOA gives an empirically higher approximation ratio
than the naïve application of level-1 QAOA and the Goemans–Williamson (GW)
algorithm [7], where the latter is the best generic classical algorithm with a
theoretical guarantee for the MAX-CUT problem under the unique games
conjecture [13]. This result is interesting because in [3], it was proved that
the approximation ratio of the naïve level-1 QAOA cannot outperform the GW for
some specific graphs. The application of the RQAOA has been extended to the
MAX-$k$-CUT problem in [14], which also shows that the RQAOA outperforms the
best-known classical approximation algorithms.

Here, the following questions naturally arise: Is it possible to apply the
QRAC in the RQAOA framework? Furthermore, if it is possible, are there any
benefits of using the QRAC in the RQAOA framework, such as improved
approximation ratios? In this paper, we propose a general framework for the
MAX-CUT problem that incorporates the QRAC [11] into the RQAOA [3] framework,
namely recursive quantum random access optimization (RQRAO, see Fig. 1). Our
contributions are summarized as follows:

  1. 1.

We theoretically show that incorporating QRACs into the RQAOA framework is a
natural expansion of [11] and [3], i.e., they are all uniformly expressed as a
method that (i) maximizes the expectation of the problem Hamiltonian that all
candidate solutions are embedded in using QRACs, and (ii) searches for the
nearest embedding in the fidelity (Theorem 1).

  2. 2.

We give numerical experiments to show the benefit of incorporating QRACs into
the RQAOA framework. Demonstrations on an open graph dataset, Gset [15], were
carried out using the tensor network simulation conducted on a classical
computer. The results showed that the approximation ratio of the proposed
method outperforms not only GW but also the methods from which inspiration was
taken [11, 3], and is comparable to existing state-of-the-art classical
heuristics [16].

The remainder of the paper is organized as follows: Section 2 presents the
mathematical definition of the MAX-CUT problem and describes the RQAOA and
QRAC in brief. Section 3 describes the proposed algorithm and its theoretical
background. Section 4 presents the results of numerical experiments. Section 5
summarizes the related works. Section 6 concludes the paper.

##  2 Preliminaries

###  2.1 MAX-CUT Problem

As stated in Sec. 1, the goal of the MAX-CUT problem is to find binary node
labels such that the sum of the weights of edges with different node labels is
maximized. It is mathematically formulated as follows:

######  Definition 1 (MAX-CUT problem)

Given a graph $G=(V,E)$, where $V$ is a node set and $E$ is an edge set with
edge weights $w_{jk}\in\mathbb{R}$ for $(j,k)\in E$, the MAX-CUT problem is
defined as

$\max_{{\bf b}\in\\{0,1\\}^{|V|}}\texttt{CW}({\bf b})$ (1)

where

$\texttt{CW}({\bf b}):=\sum_{(j,k)\in E}w_{jk}\
\frac{1-(-1)^{b_{j}+b_{k}}}{2}$ (2)

is the cut weight corresponding to ${\bf b}$.

###  2.2 Recursive QAOA

Recursive QAOA (RQAOA) [3] is a variant of QAOA, which is an iterative method
that determines the parity of one edge per iteration, then deletes a node and
modifies a graph. Here, edge parity is defined as follows:

######  Definition 2 (Edge parity)

Given edge $(j,k)$, positive and negative edges are defined as $b_{j}=b_{k}$
and $b_{j}\neq b_{k}$, respectively, where $b_{i}\in\\{0,1\\}$ is an attribute
of node $i$. Define edge parity as being a positive or negative edge.

In the first step of RQAOA, a candidate edge whose parity is to be determined
is calculated using optimization. Let $|\psi({\boldsymbol{\theta}})\rangle$ be
a QAOA ansatz with learnable parameters $\boldsymbol{\theta}$ and

$H_{\rm Ising}:=\sum_{(j,k)\in E}w_{jk}\,\frac{I-Z_{j}Z_{k}}{2}\\\ $ (3)

be the Ising-type MAX-CUT Hamiltonian, where $Z_{j}$ is a Pauli $Z$ matrix
corresponding to the $j$th qubit. The optimized parameters
${\boldsymbol{\theta}}^{*}$ are obtained by maximizing
$\langle\psi({\boldsymbol{\theta}})|H_{\rm
Ising}|\psi({\boldsymbol{\theta}})\rangle$. After optimization, the energy of
edge $(j,k)\in E$ is calculated as

$\mathcal{E}_{jk}:=\langle\psi({\boldsymbol{\theta}}^{*})|Z_{j}Z_{k}|\psi({\boldsymbol{\theta}}^{*})\rangle.$
(4)

The candidate edge $(j^{*},k^{*})$ is defined as a maximal energy edge as

$(j^{*},k^{*})=\mathop{\rm argmax}_{(j,k)\in E}|\mathcal{E}_{jk}|.$ (5)

Then, the parity of edge $(j^{*},k^{*})$ is determined as

$\displaystyle\left\\{\begin{array}[]{ll}b_{j^{*}}=b_{k^{*}}&\mathcal{E}_{j^{*}k^{*}}>0\\\\[14.22636pt]
b_{j^{*}}\neq b_{k^{*}}&\mathcal{E}_{j^{*}k^{*}}<0\end{array}\right..$ (8)

Note that the probability that all edge energy is equal to zero
($\mathcal{E}_{j^{*}k^{*}}=0$) is negligibly small and can be ignored.

In the next step, the problem Hamiltonian is modified to satisfy Eq. (8) as a
constraint. This modification can be interpreted as modifying the graph such
that the number of its nodes is reduced by one, as shown in Fig. 2.

Figure 2: Graph modification procedure of RQAOA. (a) Initial graph. (b) The
red edge indicates ${(j^{*},k^{*})=\mathop{\rm
argmax}_{(j,k)}}|\mathcal{E}_{jk}|$, where
$|\mathcal{E}_{jk}|:=|\langle\psi|Z_{j}Z_{k}|\psi\rangle|$. One of the nodes
in edge $(j^{*},k^{*})$ is kept (the green node) while the other one is
deleted (the red node). (c) After node removal. (d) Edges are re-connected to
the retained node by multiplying ${\rm sgn}(\mathcal{E}_{j^{*}k^{*}})$ by its
weight.

Node $j^{*}$ is deleted from the graph and edge $(j^{*},l)$ $(l\neq k^{*})$ is
re-connected as $(k^{*},l)$. At the same time, edge weight $w_{k^{*}l}$ is
updated to $w_{k^{*}l}+{\rm sgn}(\mathcal{E}_{j^{*}k^{*}})\cdot w_{j^{*}l}$
$(l\neq k^{*})$. As a result, one node is deleted from the graph by a single
iteration. Subsequently, the variational state is re-initialized and optimized
using the Hamiltonian defined with respect to the modified graph.

The number of nodes is gradually decreased by repeating the edge parity
determination and graph size reduction. When the number of nodes becomes
smaller than a threshold value $M$, the parities of the remaining nodes are
determined by brute-force search.

###  2.3 Quantum Random Access Optimization

In the Quantum Random Access Optimization (QRAO) [11], binary node attributes
are encoded in fewer qubits than the number of nodes using quantum random
cccess code (QRAC) [1, 2]. For each node $j\in V$, binary node attribute
$b_{j}\in\\{0,1\\}$ is randomly mapped to $(q^{(j)},P^{(j)})$, where $q^{(j)}$
is a qubit index and $P^{(j)}$ is a Pauli matrix with the following
constraints:

$\left\\{\begin{array}[]{ll}\displaystyle{P^{(j)}\neq P^{(k)}}&{\rm for\ \ \
}q^{(j)}=q^{(k)}\\\\[14.22636pt] \displaystyle{q^{(j)}\neq
q^{(k)}}&{\rm\forall\ \ \ }(j,k)\in E\end{array}\right..$ (9)

In the quantum random access coding, there are variations that depend on the
number of variants of the Pauli matrices. If only $Z$ is available,
$P^{(j)}=Z$, it reduces to the Ising-type formulation, which we refer to
$(1,1)$-QRAC. For $(2,1)$-QRAC, two variants of 1-local Pauli matrices, e.g.,
$P^{(j)}\in\\{X,Z\\}$, are available. For $(3,1)$-QRAC,
$P^{(j)}\in\\{X,Y,Z\\}$ are available. Then, the relaxed MAX-CUT Hamiltonian,
which we refer to as the $(m,1)$-QRAC Hamiltonian $(m=1,2,3)$, is defined as

$H_{m}:=\sum_{(j,k)\in E}w_{jk}\,\frac{I-mP_{\langle j\rangle}P_{\langle
k\rangle}}{2},$ (10)

where $P_{\langle j\rangle}:=P^{(j)}_{q^{(j)}}$ is a 1-local Pauli matrix
corresponding to the $q^{(j)}$th qubit. If the given graph is sparse,
constraint Eq. (9)2 rarely restricts the assignment, and hence, the total
number of qubits can be reduced by up to $1/m$th the number of nodes.

Let $\mu^{[q]}_{m}\in\mathbb{C}^{2\times 2}$ be the $q$th single qubit pure
state in the form of a density matrix defined as

$\mu^{[q]}_{m}(\\{b^{[q]}_{\bullet}\\}):=\frac{1}{2}\Bigg{(}I+\frac{1}{\sqrt{m}}\sum_{P\in\mathfrak{P}}(-1)^{b^{[q]}_{P}}P\Bigg{)},$
(11)

where $\mathfrak{P}\subset\\{X,Y,Z\\}$, $|\mathfrak{P}|=m$, and
$b^{[q]}_{P}\in\\{0,1\\}$ is a binary node attribute mapped to $(q,P)$.

Figure 3: Representation of the $(3,1)$-magic state in the Bloch sphere. The
coordinates of the vertex of the cube inscribed in the sphere correspond to
the embedded classical bits $(b^{[q]}_{X},b^{[q]}_{Y},b^{[q]}_{Z})$.

Because $\mu^{[q]}_{m}(\\{b^{[q]}_{\bullet}\\})$ is a one-qubit pure state, it
can be drawn in the Bloch sphere (Fig. 3). For example, in the case of $m=3$,
the states $\mu^{[q]}_{3}(\\{b^{[q]}_{\bullet}\\})$ for
$(b^{[q]}_{X},b^{[q]}_{Y},b^{[q]}_{Z})\in\\{0,1\\}^{3}$ correspond to the
vectors pointing the vertices of the cube inscribed in the Bloch sphere. The
$n$-qubit product state of $\mu^{[q]}_{m}$ is defined as

$\mu_{m}({\bf
b}):=\mu_{m}^{[n-1]}(\\{b^{[n-1]}_{\bullet}\\})\otimes\cdots\otimes\mu_{m}^{[0]}(\\{b^{[0]}_{\bullet}\\}).$
(12)

According to [11, Proposition 1], the following property holds.

${\rm tr}(H_{m}\,\mu_{m}({\bf b}))=\texttt{CW}({\bf b}).$ (13)

Of course, the $n$-qubit quantum state is not required to be a product state
such as Eq. (12). If $\rho$ is assumed to be an arbitrary $n$-qubit state that
includes the aforementioned product state, ${\rm tr}(H_{m}\,\rho)$ could be
higher than ${\rm tr}(H_{m}\,\mu({\bf b}))$ when $m>1$, where $\rho$ may not
be directly related to an optimal binary solution. This is why $H_{m}$ $(m>1)$
is called the relaxed Hamiltonian. Note that it has not been clarified what
exactly ${\rm tr}(H_{m}\,\rho)$ represents when $\rho$ is an arbitrary state
in [11], but this will be clarified in Corollary 1 in Sec. 3.2.

After maximizing ${\rm tr}(H_{m}\rho)$ with respect to $\rho$, the classical
bit string ${\bf b}$ is decoded from the optimized $\rho$. This process is
called rounding. In [11], two kinds of rounding were proposed, magic state
rounding and Pauli rounding. Note that magic state rounding is a measurement
method rather than a decoding method, and hence magic state rounding is
henceforth referred to as random magic measurements. The random magic
measurement is a projective measurement whose basis is the magic state
selected uniformly at random. The detailed protocol is defined in Definition 3
of Supplementary Information A. Note that $(1,1)$-random magic measurement
with $\mathfrak{P}=\\{Z\\}$ is equivalent to measurements on a computational
basis. Pauli rounding is a decoding method based on the value of ${\rm
tr}(P_{\langle j\rangle}\rho)$. The value of ${\bf b}$ is determined as
$b_{j}=0$ if ${\rm tr}(P_{\langle j\rangle}\rho)>0$ and $b_{j}=1$ if ${\rm
tr}(P_{\langle j\rangle}\rho)<0$. If ${\rm tr}(P_{\langle j\rangle}\rho)=0$,
$b_{j}$ is determined as $0$ or $1$ uniformly at random. In Sec. 3.2, we
provide theorems to show the relationship between random magic measurements
and Pauli rounding and reveal what assumption is made in Pauli rounding.

##  3 Proposed Method

The proposed method, named RQRAO, is introduced in this section.

Section 3.1 describes the concrete algorithm of the proposed method. In brief,
the proposed method incorporates the $(m,1)$-QRAC Hamiltonian (Eq. (10)) into
the RQAOA framework (Sec. 2.2). The edge parity is recursively determined
based on the value of the edge energy, which is defined using the $(m,1)$-QRAC
Hamiltonian. Two heuristics are also incorporated into the proposed method.
One is the ensemble method of defining the edge energy, and the other
determines several edge parities at once using the maximum spanning tree. The
former improves the resulting cut weight whereas the latter reduces the
runtime of the proposed method.

In Sec. 3.2, we show that not only the proposed method but also RQAOA and QRAO

  1. 1.

try to maximize the expected cut weight using $(m,1)$-random magic
measurements with respect to the quantum state, and

  2. 2.

try to decode the bit string embedded in the nearest magic state in the view
of the fidelity against the optimized quantum state.

###  3.1 Algorithm

The pseudocode of RQRAO is presented in Algorithm 1. The algorithm describes
how the ensemble edge energy $\mathfrak{E}_{jk}$ of edge $(j,k)$ is calculated
and how the edge parity is calculated from $\mathfrak{E}_{jk}$. The details
are described below.

input : Graph $G=(V,E)$, number of ensembles $N$, scale factor $S$, brute-
force search threshold $M$

output : MAX-CUT solution ${\bf b}^{*}$

initialize : $\mathcal{P}\leftarrow\emptyset$,
$G^{\prime}=(V^{\prime},E^{\prime})\leftarrow G=(V,E)$

1 $w^{\prime}_{jk}\leftarrow w_{jk}+\xi,\ \ \ \xi\sim{\rm
Uniform}([-10^{-5},10^{-5}])$;

2 while  _$|V^{\prime}| >M$_ do

/* Compute the ensemble edge energy */

3 for  _$t=1,\ldots,N$_ do

4 Assign Pauli $P^{(t)}$ to each node $j\in V^{\prime}$ and make the
$(m,1)$-QRAC Hamiltonian $H_{m}^{(t)}(G^{\prime})$;

5 $\theta^{(t)*}\leftarrow$ Result of maximization of
$\langle\psi(\theta)|H_{m}^{(t)}(G^{\prime})|\psi(\theta)\rangle$;

6
$\displaystyle{\mathcal{E}^{(t)}_{jk}\leftarrow\langle\psi(\theta^{(t)*})|P^{(t)}_{\langle
j\rangle}P^{(t)}_{\langle k\rangle}~{}|\psi(\theta^{(t)*})\rangle}$ for all
$(j,k)\in E^{\prime}$;

7

8 end for

9 for  _$(j,k)\in E^{\prime}$_ do

10 $\mu_{jk}\leftarrow\texttt{Mean}(\\{\mathcal{E}_{jk}^{(t)}\\})$;

11
$\sigma_{jk}\leftarrow\texttt{StandardDeviation}(\\{\mathcal{E}_{jk}^{(t)}\\})$;

12 $\displaystyle{\mathfrak{E}_{jk}\leftarrow}\mu_{jk}-{\rm
sign}(\mu_{jk})\cdot{\rm min}(S\sigma_{jk},|\mu_{jk}|)$;

13

14 end for

/* Compute the edge parity and modify the graph */

15 $E_{\rm S}\leftarrow\\{(j,k,|\mathfrak{E}_{jk}|):(j,k)\in
E^{\prime},|\mathfrak{E}_{jk}|>0\\}$;

16 $E_{\rm T}\leftarrow\texttt{MaximumSpanningTree}(E_{\rm S})$;

17 for  _$(j,k)\in E_{\rm T}$_ do

18 $\mathcal{P}\leftarrow\mathcal{P}\cup\\{(j,k,{\rm
sign}(\mathfrak{E}_{jk}))\\}$;

19
$G^{\prime}=(V^{\prime},E^{\prime})\leftarrow\texttt{GetReducedGraph}(G^{\prime},(j,k))$

20 end for

21

22 end while

23${\bf b}^{*}\leftarrow\texttt{BruteForceSearch}(G,\mathcal{P})$;

Algorithm 1 Algorithm for recursive quantum random access optimization
(RQRAO).

##### Line 1: Add Noise to Edges

Add small perturbation noise to all edge weights to mitigate the isolated
nodes problem (see Supplementary Information B).

##### Line 4: Define the Relaxed Hamiltonian

Assign 1-local Pauli $P^{(t)}$ to each node included in the current graph
$G^{\prime}=(V^{\prime},E^{\prime})$, where superscript $(t)$ describes the
$t$th trial. Note that $t$ is the sample index of the ensemble, not to be
confused with the recursive iteration index.

The assignment is randomly carried out so as not to violate constraints Eq.
(9). Then, the $t$th relaxed Hamiltonian is defined as

$H_{m}^{(t)}:=\sum_{(j,k)\in E^{\prime}}\frac{I-mP^{(t)}_{\langle
j\rangle}P^{(t)}_{\langle k\rangle}}{2}w^{\prime}_{jk}.$ (14)

##### Line 5: Optimization

Define ansatz $|\psi({\boldsymbol{\theta}})\rangle$, where
${\boldsymbol{\theta}}\in\mathbb{R}^{K}$ are $K$ learnable parameters. The
objective function is defined as

$\mathcal{L}^{(t)}({\boldsymbol{\theta}}):=\langle{\psi}(\boldsymbol{\theta})|H_{m}^{(t)}|{\psi}(\boldsymbol{\theta})\rangle.$
(15)

The optimized quantum state $|\psi({\boldsymbol{\theta}}^{*})\rangle$ is
obtained by maximizing $\mathcal{L}^{(t)}({\boldsymbol{\theta}})$. Any
optimization method can be used here.

##### Line 6: Compute the Edge Energy

The $t$th edge energy is computed as

$\mathcal{E}^{(t)}_{jk}:=\langle\psi({\boldsymbol{\theta}}^{*})|P^{(t)}_{\langle
j\rangle}P^{(t)}_{\langle k\rangle}|\psi({\boldsymbol{\theta}}^{*})\rangle.$
(16)

The resulting edge energy strongly depends on both the assignment of 1-local
Pauli $P^{(t)}$ and the reached quantum state, which is in general not a
eigenstate corresponding to the maximum eigenvalue of $H_{m}^{(t)}$.

##### Lines 7–11: Compute the Ensemble Edge Energy

Using $N$ edge energies, $\\{\mathcal{E}^{(t)}_{jk}\\}^{N}_{t=1}$, the
ensemble edge energy is calculated as

$\mathfrak{E}_{jk}:=\mu_{jk}-{\rm sign}(\mu_{jk})\cdot{\rm
min}(S\sigma_{jk},|\mu_{jk}|)$ (17)

where $\mu_{jk}$ and $\sigma_{jk}$ are the mean and standard deviation of
$\\{\mathcal{E}^{(t)}_{jk}\\}^{N}_{t=1}$, which are defined as
$\mu_{jk}:=\frac{1}{N}\sum^{N}_{t=1}\mathcal{E}^{(t)}_{jk}$ and
$\sigma_{jk}:=\sqrt{\frac{1}{N}\sum^{N}_{t=1}(\mathcal{E}^{(t)}_{jk}-\mu_{jk})^{2}}$,
respectively. Equation (17) means that $\mathfrak{E}_{jk}$ has a nonzero value
if $\mu_{jk}\pm S\sigma_{jk}$ have the same sign as $\mu_{jk}$.

##### Lines 12–16: Compute the Maximum Spanning Tree

After computing the ensemble edge energy $\mathfrak{E}_{jk}$, the parities of
edges with nonzero $\mathfrak{E}_{jk}$ are determined according to the sign of
$\mathfrak{E}_{jk}$. However, sometimes inconsistency arises. For example,
this occurs when three negative edges form a cycle
$\\{(b_{0},b_{1}),(b_{1},b_{2}),(b_{2},b_{0})\\}$, one of which is
inconsistent because $(b_{0}\neq b_{1})\land(b_{1}\neq b_{2})\land(b_{2}\neq
b_{0})$ cannot hold. To avoid this inconsistency, we adopt the maximum
spanning tree to remove any cycles from the candidates.

Edges with nonzero $\mathfrak{E}_{jk}$ are collected and $|\mathfrak{E}_{jk}|$
are treated as edge weights. For each subgraph $E_{\rm S}$ made by this
collection that is disconnected from the others, the maximum spanning tree
$E_{\rm T}$ can be calculated in linear time in the number of edges by Prim’s
algorithm [17] or by the randomized KKT algorithm [18]. The obtained tree
$E_{\rm T}$ can be easily converted into a rooted tree by setting an arbitrary
node as the root.

##### Lines 17–21: Modify the Graph

Now, we have candidates $E_{\rm T}$ that can be used to determine the parity
and can be deleted from graph $G^{\prime}$ without inconsistency. The parity
for each edge in $E_{\rm T}$ is determined from the leaf to the root one by
one according to its sign. At the same time, the node on the leaf side is
deleted and the edge weights are updated. The graph modification procedure is
the same as the RQAOA, as depicted in Fig. 2.

##### Line 23: Brute-force Search

The recursive procedure is stopped when the number of nodes of $G^{\prime}$ is
less than or equal to $M$. Here, the number of nodes for which parity with
other nodes has not yet been determined is less than or equal to $M$. That is,
the number of candidate solutions is less than or equal to $2^{M}$. If $M$ is
sufficiently small, e.g., $M\leq 20$, a brute-force search among the
candidates is not time-consuming, and hence the best solution ${\bf b}^{*}$
can be quickly obtained. The RQAOA [3] also uses brute-force search.

##### Example

Figure 4:  Schematic illustration of the RQRAO procedure. (a) Initial graph
whose edge weights are $+1$ (red) or $-1$ (blue). (b) Three samples of edge
energy, $\mathcal{E}^{(1)}_{jk}$, $\mathcal{E}^{(2)}_{jk}$, and
$\mathcal{E}^{(3)}_{jk}$, are evaluated using optimized states whose
corresponding Hamiltonians are $H^{(1)}_{3}$, $H^{(2)}_{3}$, and
$H^{(3)}_{3}$, respectively. Node labels are assigned 1-local Paulis. (c)
Ensemble edge energy $\mathfrak{E}_{jk}$ using the three samples illustrated
in (b). Closed loop $0-4-2-8-0$ is inconsistent because
$(b_{0}=b_{4})\land(b_{4}\neq b_{2})\land(b_{2}=b_{8})\land(b_{8}=b_{0})$
cannot be satisfied. (d) Maximum spanning tree with $|\mathfrak{E}_{jk}|$ edge
weights. Red and blue indicate positive and negative edges, respectively. Edge
$(0,4)$ is ignored at this step. (e) Modified graph after all leaf nodes of
(d) have been removed and the edges in the tree have been reconnected to the
remaining nodes. RQRAO repeats steps (a) to (e) until the number of nodes
becomes smaller than $M$.

The whole RQRAO procedure is schematically shown in Fig. 4. The problem graph
$G$ is shown in Fig. 4, where the red and blue edges have weights $+1$ and
$-1$, respectively. The random Pauli assignments and optimization are carried
out to obtain the edge energy for each trial $t$, $\mathcal{E}^{(t)}_{jk}$
(Fig. 4. The ensemble edge energy $\mathfrak{E}_{jk}$ is then calculated using
$\\{\mathcal{E}^{(t)}_{jk}\\}^{N}_{t=1}$ (Fig. 4. The inconsistent edges are
removed by solving the maximum spanning tree problem, where
$|\mathfrak{E}_{jk}|$ are regarded as the edge weights. The rooted tree is
constructed by giving the direction from root to leaf, where the root node has
been randomly selected (Fig. 4. The edge parities in the rooted tree are
determined from leaf to root according to the sign of corresponding
$\mathfrak{E}_{jk}$. At the same time, the graph is modified to satisfy the
determined parity. Figure 4 represents the modified graph after all edge
parities in the rooted tree have been determined. The protocol from Fig. 4 to
4 is repeated until the number of nodes of the modified graph becomes smaller
than $M$.

###  3.2 Theoretical Analysis

The proposed method, RQRAO, combines the RQAOA [3] and QRAO [11]. These
existing methods work well experimentally, but there are some theoretical
ambiguities: (i) Why does the recursive approach improve the cut weight
performance in the RQAOA? (ii) What is the meaning of maximizing the
expectation value of the QRAC Hamiltonian, and what is the theoretical
background of Pauli rounding in QRAO? To answer these questions, we prove the
following theorem:

###### Theorem 1

The objective of the RQAOA, QRAO, and RQRAO can be unifiedly expressed as the
problem of solving

$\rho^{*}=\mathop{\rm argmax}_{\rho}\mathbb{E}_{{\bf
b}\sim\mathcal{P}_{m}({\bf b};\rho)}\Big{[}\texttt{CW}({\bf b})\Big{]},$ (18)
${\bf b}^{\rm R}=\mathop{\rm argmax}_{{\bf
b}\in\\{0,1\\}^{|V|}}\mathcal{P}_{m}({\bf b};\rho^{*}),$ (19)

where $\mathcal{P}_{m}({\bf b};\rho)$ is the probability of obtaining ${\bf
b}$ by $(m,1)$-random magic measurements and is proportional to the fidelity
between $\rho$ and the $(m,1)$-magic state embedding $\mu({\bf b})$.

From Theorem 1, it is clarified that both the RQAOA and QRAO, let alone RQRAO,
are models searching for the highest probability sample of the probability
distribution that maximizes the expectation value of the cut weight through
the quantum state as the hidden variable. It will be proved later that
maximizing the expectation value of the QRAC Hamiltonian, let alone that of
the Ising-type Hamiltonian, can be seen as carrying out Eq. (18). In addition,
it will also be proved later that both the recursive approach in the RQAOA and
Pauli rounding in the QRAO are just methods that approximate methods of Eq.
(19).

Of course, Theorem 1 does not say the MAX-CUT problem is efficiently solvable
because both maximizing the expected cut weight and obtaining the highest
probability sample is, in general, intractable when ${\bf b}$ lives in an
extremely high dimensional space, $\\{0,1\\}^{|V|}$ for $|V|\gg 1$. However,
we can say that it is natural to take the ${\bf b}$ with the highest
probability because, after solving Eq. (18), which maximizes the expected cut
weight, the probability for larger $\texttt{CW}({\bf b})$ is expected to be
large. It is therefore reasonable calculate Eqs. (18) and (19) to obtain a
better cut weight.

In Sec. 3.2.1, Eq. (18) is proved. Then, Eq. (19) is proved in Sec. 3.2.2. All
the proofs of the following theorems and corollaries are deferred to
Supplementary Information C.

####  3.2.1 Proof of Eq. (18)

The objectives of the RQAOA, QRAO, and RQRAO are to maximize ${\rm
tr}(H_{m}\rho)$, where $m=1$ for the RQAOA whereas $m\in\\{1,2,3\\}$ for QRAO
and RQRAO. Hence, to prove Eq. (18), we need to show that maximizing ${\rm
tr}(H_{m}\rho)$ is equivalent to maximizing $\mathbb{E}_{{\bf
b}\sim\mathcal{P}_{m}({\bf b};\rho)}\big{[}\texttt{CW}({\bf b})\big{]}$. To
show this, let us start with the following theorem and corollary, which are
proved for the first time in this paper.

###### Theorem 2

Let $P^{\otimes k}:=\prod^{k}_{i=1}P_{q^{(i)}}\in{\mathbb{C}}^{2^{n}\times
2^{n}}$ be the $k$-local Pauli observable, where $P_{q}$ represents $1$-local
Pauli $P\in\mathfrak{P}\subset\\{X,Y,Z\\}$ acting on the $q$th qubit, and
$q^{(i)}$ is the $i$th qubit index, where $q^{(i)}\neq q^{(j)}$ for $i\neq j$.
Suppose the case $|\mathfrak{P}|=m$. Let $\rho\in{\mathbb{C}}^{2^{n}\times
2^{n}}$ be the arbitrary $n$ qubit state, $\mu_{m}({\bf
b})\in{\mathbb{C}}^{2^{n}\times 2^{n}}$ be the $n$ qubit $(m,1)$-magic state
defined in Eq. (12), in which ${\bf b}\in\\{0,1\\}^{mn}$ is embedded, and
$f_{m}({\bf b},\rho):={\rm tr}(\mu_{m}({\bf b})\rho)$ be the fidelity between
$\mu_{m}({\bf b})$ and $\rho$. The following relationship holds:

${\rm tr}\big{(}P^{\otimes k}\rho\big{)}=m^{k/2}\,\mathbb{E}_{{\bf
b}\sim\mathcal{P}_{m}({\bf
b};\rho)}\Big{[}(-1)^{\sum^{k}_{i=1}b^{[q^{(i)}]}_{P}}\Big{]},$ (20)

where $b^{[q]}_{P}\in\\{0,1\\}$ is a binary mapped to the Pauli $P$ of the
$q$th qubit and

$\mathcal{P}_{m}({\bf b};\rho):=\frac{f_{m}({\bf b},\rho)}{\sum_{{\bf
b}^{\prime}\in\\{0,1\\}^{mn}}f_{m}({\bf b}^{\prime},\rho)}.$ (21)

Note that ${\bf b}$ following the probability $\mathcal{P}_{m}$ is the outcome
of the $(m,1)$-random magic measurements, whose details are described in
Supplementary Information A.

###### Corollary 1

Let $\rho\in{\mathbb{C}}^{2^{n}\times 2^{n}}$ be the arbitrary $n$ qubit state
and $H_{m}$ be the $(m,1)$-QRAC Hamiltonian defined as Eq. (10). The following
relationship holds:

$\displaystyle{{\rm tr}(H_{m}\,\rho)}$ $\displaystyle{=\sum_{(j,k)\in
E}\frac{1-m^{2}\,\mathbb{E}_{{\bf b}\sim\mathcal{P}_{m}({\bf
b};\rho)}\big{[}(-1)^{b_{j}+b_{k}}\big{]}}{2}w_{jk}}$
$=m^{2}\,\mathbb{E}_{{\bf b}\sim\mathcal{P}_{m}({\bf
b};\rho)}\big{[}\texttt{CW}({\bf b})\big{]}-\frac{m^{2}-1}{2}\sum_{(j,k)\in
E}w_{jk},$ (22)

where $\mathcal{P}_{m}$ is defined as Eq. (21).

According to Eq. (22), the expectation value of the $(m,1)$-QRAC Hamiltonian
is proportional to the expectation value of the cut weight, where the
expectation is taken by the bit strings obtained by $(m,1)$-random magic
measurements. Hence, maximizing the expectation value of the $(m,1)$-QRAC
Hamiltonian can be regarded as maximizing the expected cut weight. Note that
the lower bounds of $\mathbb{E}_{{\bf b}\sim\mathcal{P}_{m}({\bf
b};\rho)}\big{[}\texttt{CW}({\bf b})\big{]}$ for $w_{jk}=1$ when ${\rm
tr}(H_{m}\rho)\geq\texttt{CW}({\bf b}^{*})$ have been given in [11, Theorem 2
and Supplementary Information VII]. They are $\mathbb{E}_{{\bf
b}\sim\mathcal{P}_{m}({\bf b};\rho)}\big{[}\texttt{CW}({\bf
b})\big{]}\geq\frac{5}{9}\texttt{CW}({\bf b}^{*})$ for $m=3$ and
$\mathbb{E}_{{\bf b}\sim\mathcal{P}_{m}({\bf b};\rho)}\big{[}\texttt{CW}({\bf
b})\big{]}\geq\frac{5}{8}\texttt{CW}({\bf b}^{*})$ for $m=2$. These bounds are
reproduced by substituting ${\rm tr}(H_{m}\rho)\geq\texttt{CW}({\bf b}^{*})$
and $\sum_{(j,k)\in E}w_{jk}\geq\texttt{CW}({\bf b}^{*})$ into Eq. (22).

Note that not only the expectation value but also the variance of the cut
weight can be obtained using the $(m,1)$-QRAC Hamiltonian as follows.

###### Corollary 2

Let $\rho\in{\mathbb{C}}^{2^{n}\times 2^{n}}$ be the arbitrary $n$ qubit state
and $H_{m}$ be the $(m,1)$-QRAC Hamiltonian defined as Eq. (10). The variance
of the cut weight is

$\displaystyle{{\rm Var}_{{\bf b}\sim\mathcal{P}_{m}({\bf
b};\rho)}\big{[}\texttt{CW}({\bf b})\big{]}}$ $=\frac{1}{m^{4}}\Big{(}{\rm
tr}(H_{m}^{2}\rho)-{\rm tr}(H_{m}\rho)^{2}\Big{)},$ (23)

where $\mathcal{P}_{m}$ is defined as Eq. (21).

####  3.2.2 Proof of Eq. (19)

First, Eq. (19) is rewritten in the form of the conditional distribution to
make it easier to understand. A joint distribution $\mathcal{P}_{m}({\bf
b};\rho^{*})$ can be written as the product of the conditional distribution as

$\mathcal{P}_{m}({\bf
b};\rho^{*})=\prod^{|V|}_{a=1}\mathcal{P}_{m}(b_{j_{a}}|{\bf
b}_{<j_{a-1}};\rho^{*}),$ (24)

where ${\bf b}_{<j_{a-1}}:=\\{b_{j_{1}},\ldots,b_{j_{a-1}}\\}$. Using this
relationship, Eq. (19) can be written as

$\displaystyle{b^{\rm R}_{j_{a}}=\mathop{\rm
argmax}_{b_{j_{a}}\in\\{0,1\\}}\mathcal{P}_{m}(b_{j_{a}}|{\bf b}^{\rm
R}_{<j_{a}};\rho^{*})}$ ${\rm for}\ \ \ a=1,2,\ldots,|V|.$ (25)

Because the cut weight is invariant to the bit flip, we can fix one of the
$b^{\rm R}_{j}$ to $0$ or $1$. In this case, there exists a one-to-one
relationship between $\mathcal{P}_{m}({\bf b};\rho^{*})$ and

${\mathcal{P}}_{m}(\mathfrak{p};\rho^{*},T):=\prod^{|V|-1}_{b=1}\mathcal{P}_{m}(\mathfrak{p}_{e_{b}}|\mathfrak{p}_{<e_{b}};\rho^{*},T),$
(26)

where

$\mathfrak{p}_{e=(c,p)}:=\left\\{\begin{array}[]{ll}+&(b_{c}=b_{p})\\\\[5.69054pt]
-&(b_{c}\neq b_{p})\end{array}\right.$ (27)

is a random variable of an edge parity, where “$+$” and “$-$” indicate
positive and negative parities, respectively,
$T=\\{e_{j_{b}}=(c_{j_{b}},p_{j_{b}})\\}^{|V|-1}_{b=1}$ is an arbitrary
spanning tree of given graph $G$, and
$\mathfrak{p}_{<e_{b}}:=\\{\mathfrak{p}_{e_{1}},\ldots,\mathfrak{p}_{e_{b-1}}\\}$.
Here, the reason for considering the probability of the edge parity on the
tree is to exclude inconsistent parities such as a triangle with negative
parities. Excluding such an inconsistency does not lose the generality because
all possible configurations of solutions can be treated by considering only
the probability on the tree. From the view of edge parity, Eq. (19) can be
rewritten as

$\displaystyle{\mathfrak{p}^{\rm R}_{e_{b}}=\mathop{\rm
argmax}_{\mathfrak{p}_{e_{b}}\in\\{+,-\\}}\mathcal{P}_{m}(\mathfrak{p}_{e_{b}}|\mathfrak{p}^{\rm
R}_{<e_{b}};\rho^{*},T^{*})}$ ${\rm for}\ \ \ b=1,2,\ldots,|V|-1,$ (28)

where

$T^{*}:=\mathop{\rm
argmax}_{T}\max_{\mathfrak{p}\in\\{+,-\\}^{|V|-1}}{\mathcal{P}}_{m}(\mathfrak{p};\rho^{*},T).$
(29)

In the following, we show that QRAO uses Eq. (25), whereas both the RQAOA and
RQRAO use Eq. (28) to obtain the solution. All of them assume a moderate
approximation on a probability distribution. The following corollary is
frequently used for the proof.

###### Corollary 3

Define $\mathcal{E}_{j}:={\rm tr}(P_{\langle j\rangle}\,\rho)$ and
$\mathcal{E}_{jk}:={\rm tr}(P_{\langle j\rangle}P_{\langle k\rangle}\,\rho)$.
When the bit string ${\bf b}=b_{|V|-1}\ldots b_{1}b_{0}$ is obtained by the
$(m,1)$-random magic measurements, ${\bf b}\sim\mathcal{P}_{m}({\bf b};\rho)$,
the probabilities of $b_{j}=0$ and $b_{j}=1$ are

$\begin{array}[]{lll}\mathcal{P}(b_{j}=0)&=&\displaystyle{\frac{1}{2}+\frac{1}{2\sqrt{m}}\mathcal{E}_{j}}\\\\[14.22636pt]
\mathcal{P}(b_{j}=1)&=&\displaystyle{\frac{1}{2}-\frac{1}{2\sqrt{m}}\mathcal{E}_{j}}\end{array},$
(30)

respectively, whereas the probability of $b_{j}=b_{k}$ and $b_{j}\neq b_{k}$
are

$\begin{array}[]{lll}\mathcal{P}(b_{j}=b_{k})&=&\left\\{\begin{array}[]{ll}\displaystyle{\frac{1}{2}+\frac{1}{2m}\mathcal{E}_{j}\mathcal{E}_{k}}&(q^{(j)}=q^{(k)})\\\\[14.22636pt]
\displaystyle{\frac{1}{2}+\frac{1}{2m}\mathcal{E}_{jk}}&(q^{(j)}\neq
q^{(k)})\end{array}\right.\\\\[28.45274pt] \mathcal{P}(b_{j}\neq
b_{k})&=&\left\\{\begin{array}[]{ll}\displaystyle{\frac{1}{2}-\frac{1}{2m}\mathcal{E}_{j}\mathcal{E}_{k}}&(q^{(j)}=q^{(k)})\\\\[14.22636pt]
\displaystyle{\frac{1}{2}-\frac{1}{2m}\mathcal{E}_{jk}}&(q^{(j)}\neq
q^{(k)})\end{array}\right.\end{array},$ (31)

respectively, where $q^{(j)}$ is a qubit with $b_{j}$ embedded.

Equation (30) is consistent with the upper bound of the decoding probability,
$\frac{1}{2}+\frac{1}{2\sqrt{m}}$ [1], which is recovered when
$|\mathcal{E}_{j}|=1$.

##### Pauli rounding

Assume all $b_{j_{a}}$ for $a=1,2,\ldots,|V|$ are independent. Under this
assumption, Eq. (24) is approximated as

${\mathcal{P}}^{\rm Pauli}_{m}({\bf
b};\rho^{*}):=\prod^{|V|}_{a=1}\mathcal{P}_{m}(b_{j_{a}};\rho^{*}).$ (32)

Substituting Eq. (32) into Eq. (25) gives

$\begin{array}[]{lll}b^{\rm Pauli}_{j_{a}}&:=&\displaystyle{\mathop{\rm
argmax}_{b_{j_{a}}\in\\{0,1\\}}\mathcal{P}_{m}(b_{j_{a}};\rho^{*})}\\\\[14.22636pt]
&=&\displaystyle{\mathop{\rm
argmax}_{b_{j_{a}}\in\\{0,1\\}}\Bigg{(}\frac{1}{2}+\frac{(-1)^{b_{j_{a}}}}{2\sqrt{m}}\mathcal{E}_{j_{a}}\Bigg{)}}\\\\[14.22636pt]
&=&\left\\{\begin{array}[]{ll}0&(\mathcal{E}_{j_{a}}>0)\\\
1&(\mathcal{E}_{j_{a}}<0)\end{array}\right.,\end{array}$ (33)

where the second equality uses Corollary 3. Interestingly, ${\bf b}^{\rm
Pauli}$ is equivalent to the result of the Pauli rounding of QRAO [11]. This
comes out that the Pauli rounding seeks the variable with the highest
probability of ${\mathcal{P}}_{m}({\bf b};\rho)$ under the assumption of Eq.
(32).

##### Tree rounding

Assume all $\mathfrak{p}_{e_{b}}$ for $b=1,2,\ldots,|V|-1$ are independent.
Under this assumption, Eq. (26) is approximated as

${\mathcal{P}}^{\rm
tree}_{m}(\mathfrak{p};\rho^{*},T):=\displaystyle{\prod_{b=1}^{|V|-1}\mathcal{P}_{m}(\mathfrak{p}_{e_{b}};\rho^{*},T)}.\\\
$ (34)

Substituting Eq. (34) into Eq. (28) gives

$\begin{array}[]{lll}\mathfrak{p}^{\rm
tree}_{e_{b}=(c,p)}&:=&\displaystyle{\mathop{\rm
argmax}_{\mathfrak{p}_{e_{b}}\in\\{+,-\\}}\mathcal{P}_{m}(\mathfrak{p}_{e_{b}};\rho^{*},T^{*})}\\\\[14.22636pt]
&=&\displaystyle{\mathop{\rm
argmax}_{(b_{c},b_{p})\in\\{0,1\\}^{2}}\Bigg{(}\frac{1}{2}+\frac{(-1)^{b_{c}+b_{p}}}{2m}\mathcal{E}_{cp}\Bigg{)}}\\\\[14.22636pt]
&=&\left\\{\begin{array}[]{ll}+&(\mathcal{E}_{cp}>0)\\\
-&(\mathcal{E}_{cp}<0)\end{array}\right.,\end{array}$ (35)

where the second equality uses Corollary 3, and

$\begin{array}[]{lll}T^{*}&=&\displaystyle{\mathop{\rm
argmax}_{T}\max_{\mathfrak{p}}{\mathcal{P}}^{\rm
tree}_{m}(\mathfrak{p};\rho^{*},T)}\\\\[14.22636pt]
&=&\displaystyle{\mathop{\rm
argmax}_{T\ni\\{(c,p)\\}}\max_{(b_{c},b_{p})\in\\{0,1\\}^{2}}\Bigg{(}\frac{1}{2}+\frac{(-1)^{b_{c}+b_{p}}}{2m}\mathcal{E}_{cp}\Bigg{)}}\\\\[14.22636pt]
&=&\displaystyle{\mathop{\rm
argmax}_{T\ni\\{(c,p)\\}}\big{|}\mathcal{E}_{cp}\big{|}}\end{array}$ (36)

is a maximum spanning tree. Computing Eq. (35) can be carried out by a two-
step calculation. In the first step, the maximum spanning tree $T^{*}$ is
computed with the edge weights $|\mathcal{E}_{cp}|$ for edge $(c,p)$. In the
second step, the edge parity is determined by the sign of corresponding
$\mathcal{E}_{cp}$. We refer to this rounding method as tree rounding.

RQRAO with $N=1$, which has no recursive steps, is equivalent to using tree
rounding. Because both Pauli and tree rounding assume that all random
variables are independent, it is expected that QRAO and RQRAO ($N=1$) will
yield a similar cut weight.

##### Recursive rounding

Both Pauli and tree rounding impose the strong approximation on
$\mathcal{P}_{m}$ that all random variables are independent. This is because
conditional sampling is in general intractable to carry out. For this reason,
let us impose the condition $\mathfrak{p}_{<e_{b}}=\mathfrak{p}^{*}_{<e_{b}}$
on the optimization problem as the constraint instead of on $\mathcal{P}_{m}$.
Of course, these two approaches are not exactly equivalent, but imposing the
condition as the constraint of the optimization should be a good approximation
to the conditional sampling. The specific steps are given below.

First, one parity of $e_{1}$ is determined using the optimization result,
$\rho^{*}_{1}$ as

$\mathfrak{p}^{*}_{e_{1}}=\mathop{\rm
argmax}_{\mathfrak{p}_{e_{1}}\in\\{+,-\\}}\mathcal{P}_{m}(\mathfrak{p}_{e_{1}};\rho^{*}_{1}),$
(37)

where

$e_{1}=\mathop{\rm
argmax}_{e}\max_{\mathfrak{p}_{e}\in\\{+,-\\}}\mathcal{P}_{m}(\mathfrak{p}_{e};\rho^{*}_{1}).$
(38)

According to Corollary 3, this can be rewritten as

$\mathfrak{p}^{\rm
rec}_{e_{1}=(c_{1},p_{1})}=\left\\{\begin{array}[]{ll}+&(\mathcal{E}_{c_{1}p_{1}}>0)\\\
-&(\mathcal{E}_{c_{1}p_{1}}<0)\end{array}\right.,$ (39)

where

$\begin{array}[]{lll}(c_{1},p_{1})&=&\displaystyle{\mathop{\rm
argmax}_{(c,p)\in
E}\Bigg{(}\frac{1}{2}+\frac{1}{2m}\big{|}\mathcal{E}_{cp}\big{|}\Bigg{)}}\\\\[14.22636pt]
&=&\displaystyle{\mathop{\rm argmax}_{(c,p)\in
E}\big{|}\mathcal{E}_{cp}\big{|}}.\end{array}$ (40)

After that, the given graph $G$ is modified to satisfy
$\mathfrak{p}_{e_{1}}=\mathfrak{p}^{*}_{e_{1}}$, which refers to $G_{2}$.
Next, the second optimization is carried out with the $(m,1)$-QRAC Hamiltonian
based on $G_{2}$, and the optimized state $\rho^{*}_{2}$ is obtained. Because
the modified graph is forced to be
$\mathfrak{p}_{e_{1}}=\mathfrak{p}^{*}_{e_{1}}$, it can be seen that
$\rho^{*}_{2}$ is conditioned on
$\mathfrak{p}_{e_{1}}=\mathfrak{p}^{*}_{e_{1}}$. Hence, the determined parity
of the second step can be written as

$\mathfrak{p}^{*}_{e_{2}}=\mathop{\rm
argmax}_{\mathfrak{p}_{e_{2}}\in\\{+,-\\}}\mathcal{P}_{m}(\mathfrak{p}_{e_{2}};\rho^{*}_{2}(\mathfrak{p}_{e_{1}}=\mathfrak{p}^{*}_{e_{1}}))$
(41)

where

$e_{2}=\mathop{\rm
argmax}_{e}\max_{\mathfrak{p}_{e}\in\\{+,-\\}}\mathcal{P}_{m}(\mathfrak{p}_{e};\rho^{*}_{2}(\mathfrak{p}_{e_{1}}=\mathfrak{p}^{*}_{e_{1}})).$
(42)

Of course, this could also be rewritten using the edge energy
$\mathcal{E}_{cp}$, but we will leave that out as it is self-explanatory.
Repeating this procedure gives

$\displaystyle{\mathfrak{p}^{*}_{e_{b}}=\mathop{\rm
argmax}_{\mathfrak{p}_{e_{b}}}\mathcal{P}_{m}(\mathfrak{p}_{e_{b}};\rho^{*}_{b}(\mathfrak{p}_{<e_{b}}=\mathfrak{p}^{*}_{<e_{b}}))}$
$\ \ \ {\rm for}\ \ \ b=1,\ldots,|V|-1.$ (43)

Finally,
$\mathfrak{p}^{*}=\\{\mathfrak{p}^{*}_{1},\ldots,\mathfrak{p}^{*}_{|V|-1}\\}$
is obtained. The corresponding bit string ${\bf b}^{\rm rec}$ can be easily
recovered from $\mathfrak{p}^{*}$ by fixing one of the $b^{\rm rec}_{j}$ to
$0$ or $1$. The condition $\mathfrak{p}_{<e_{b}}=\mathfrak{p}^{*}_{<e_{b}}$ is
indirectly imposed through $\rho^{*}_{b}$ to $\mathcal{P}_{m}$, and hence this
rounding method imposes an approximation that is weaker than that of Pauli and
tree rounding. We refer to this rounding method as recursive rounding.

In recursive rounding, the edge that determines the parity is selected using
the absolute value of its edge energy (Eq. (42)), and its parity is determined
by the sign of the corresponding edge energy (Eq. (41)). This process is the
same as explained in Sec. 2.2 and hence, it can be said that the RQAOA adopts
recursive rounding. The main idea behind RQAOA is to approximate $\mathop{\rm
max}_{{\bf b}\in\\{0,1\\}^{|V|}}{\rm tr}(H_{1}|{\bf b}\rangle\langle{\bf b}|)$
[3], but in actuality, it approximates $\mathop{\rm max}_{{\bf
b}\in\\{0,1\\}^{|V|}}\mathcal{P}_{1}({\bf b};\rho^{*})=\mathop{\rm max}_{{\bf
b}\in\\{0,1\\}^{|V|}}{\rm tr}(\rho^{*}|{\bf b}\rangle\langle{\bf b}|)$.111
Since $\mu_{1}({\bf b})=|{\bf b}\rangle\langle{\bf b}|$, $f_{1}({\bf
b},\rho)={\rm tr}(\mu_{1}({\bf b})\rho)$, and $\sum_{{\bf
b}\in\\{0,1\\}^{|V|}}f_{1}({\bf b},\rho)=1$, $\mathcal{P}_{1}({\bf
b};\rho)=f_{1}({\bf b},\rho)/\sum_{{\bf b}^{\prime}}f_{1}({\bf
b}^{\prime},\rho)={\rm tr}(|{\bf b}\rangle\langle{\bf b}|\rho)$.  By contrast,
the rounding for RQRAO with $N>1$ can be considered to be an intermediate
approach between tree rounding and recursive rounding because the edge
parities that are being determined simultaneously in one step are assumed to
be independent.

##  4 Numerical Experiments

In this section, we show the results of the numerical experiments of the MAX-
CUT problem using RQRAO and other existing algorithms. The main claim is that
RQRAO achieves a better cut weight than existing algorithms and is comparable
to the state-of-the-art classical heuristic.

##### Datasets

We used three graph datasets.

The first dataset is used to show the validity of Theorem 1. That is, it shows
that the bit string with the highest probability yields a good cut weight
after maximizing the expectation of the $(m,1)$-QRAC Hamiltonian. Only one
graph instance is generated by a machine-independent graph generator rudy [19]
using -rnd_graph 14 50 0 -random 0 1 0 -times 2 -plus -1, where the number of
nodes is $14$, the graph density is $0.5$, and the edge weights are set to
$\pm 1$ uniformly at random. This instance has only one MAX-CUT solution,
${\bf b}=00001101101010$, except for its bit-flip solution, where the maximum
cut weight is 12. We call this graph instance Rnd14.

The second dataset is the benchmark graph dataset Gset [15]
(http://www.stanford.edu/~yyye/yyye/Gset/). Gset includes random, toroidal,
and almost planar graphs with $|V|=800$ to $20000$ indexed by G1 to G81, which
are generated by rudy. Only G1, G6, G14, and G18 are used for evaluation, all
of which have 800 nodes. G1 and G6 are random graphs, G11 is a toroidal graph,
and G14 and G18 are the combined graphs of two planar graphs. Graphs G1 and
G6, and graphs G14 and G18 have the same edge structure except for the
weights, where G1 and G14 have $+1$ edge weights whereas G6 and G18 have edge
weights of $\pm 1$ assigned at random. The edge weights of G11 are $\pm 1$
assigned at random. We call this graph dataset Gset800.

The third dataset consists of 3-regular graphs with $\pm 1$ edge weights
generated at random. The number of nodes varies from $10$ to $10^{4}$. Ten
graphs are generated for each number of nodes. These graphs are used to
evaluate the scalability with respect to the cut weight and the runtime of the
algorithms. We call this graph dataset Rnd3R.

##### Baselines

We compared our model with GW [7], CirCut [16], QRAO [11], and RQAOA [3]. GW
is a well-known classical approximation algorithm based on semidefinite
programming. CirCut is the Fortran program name of the rank-two relaxation
algorithm that is a classical heuristic algorithm and empirically outperforms
other classical heuristics on various graph instances [20]. The details of GW
and CirCut are described in Sec. 5.1. The RQAOA is a QAOA-based algorithm
whose details are explained in Sec. 2.2. QRAO is a quantum relaxation-based
optimization method whose details are explained in Sec. 2.3.

##### Implementation Details

The matrix product state (MPS) ansatz is used for
$|\psi({\boldsymbol{\theta}})\rangle$, where all the matrices in the MPS are
treated as trainable parameters (see Supplementary Information D for details).
The default hyperparameters for RQRAO used in the numerical experiments are
listed in Table 1.

We implemented our model in Python, and the MPS optimization is implemented
using PyTorch version 1.10.0 [21]. For the optimizer, we used L-BFGS [22] with
line search using strong Wolfe conditions. The termination tolerance for
changing the function value and parameters is set to $10^{-2}$. Implementation
details for GW, CirCut, the RQAOA, and QRAO are described in Supplementary
Information E.

All experiments were carried out on a single Intel(R) Core(TM) i9-9900X CPU @
3.50GHz with 32GB (8GBx4) DDR4-2666 Quad-Channel.

Table 1: Default hyperparameters for RQRAO. Hyperparameter | Symbol | Value
---|---|---
# embeddings / qubits | $m$ | $3$
# ensemble | $N$ | $20$
Scale factor | $S$ | $2$
Bond dimension | $\chi$ | $2$
Brute-force search threshold | $M$ | $10$

##### Details of the Experiments

Using the Rnd14 dataset, Theorem 1 was verified. First, $(m,1)$-QRAC
Hamiltonians were randomly generated 100 times each for $m=1,2,$ and $3$. For
each Hamiltonian, the eigenstate with the maximum eigenvalue and optimized
state with bond dimension $\chi=2$ were then calculated. Using each state, the
cut weight with the highest probability was calculated. We call this
experiment ExpVer.

Using the Gset800 dataset, the procedure was run 10 times on the same graph
instance and the best cut of the trials was recorded. RQRAO and all baseline
methods were used. We call this experiment ExpGset.

Using the Rnd3R dataset, the procedure was run once for each graph instance
and the cut weight and elapsed time were recorded. The cut weight is
normalized with respect to that of the GW algorithm and is referred to as the
relative cut weight. Using Rnd3R, two experiments were carried out. The first
one evaluated the scalability. RQRAO and all baseline methods were used. We
call this experiment ExpScale. The second one analyzed the sensitivity of the
RQRAO hyperparameters. Only one hyperparameter of RQRAO was varied from its
default value (Table 1). We call this experiment ExpHyp.

###  4.1 Results and Discussion

##### Verification of Theorem 1 (ExpVer)

(a) MaxEig

(b) $\chi=2$

Figure 5:  Solutions with the highest probability of 100 trials for the same
graph instance, Rnd14, where the state used was (a) the eigenstate
corresponding to the maximum eigenvalue and (b) the state obtained through
optimization with $\chi=2$.

Figure 5 summarizes the results of ExpVer. If one can obtain the eigenstate
with the maximum eigenvalue, $(1,1)$-QRAC is the best choice because the cut
weight with the highest probability for $(1,1)$-QRAC is always the MAX-CUT
solution (Fig. 5(a)). By contrast, in the case of $(2,1)$ and $(3,1)$-QRACs,
the solution with the highest probability is sometimes not the MAX-CUT
solution. This occurs when low cut weight solutions are embedded near the MAX-
CUT solution, and higher but not the maximum cut weight solutions are
aggregated somewhere else in Hilbert space. In such a case, the state with the
maximum eigenvalue is placed in its aggregation, which could be far from the
MAX-CUT solution. This never occurs in the case of $(1,1)$-QRAC, because all
embeddings are orthogonal to each other. There are drawbacks due to these
accidental cases. However, $(2,1)$ and $(3,1)$-QRACs still have relatively
higher cut weights. Therefore, it makes sense to choose a sample with the
highest probability as a solution to obtain a higher cut weight solution,
which is what the RQAOA, QRAO, and RQRAO attempt to implement according to
Theorem 1, if the state corresponding to the maximum eigenvalue can be
obtained. However, in practice, the state obtained through optimization is not
the eigenstate with the maximum eigenvalue. In this case, the cut weight with
the highest probability for both $(2,1)$\- and $(3,1)$-QRACs is much better
than that for $(1,1)$-QRAC (Fig. 5(b)). This suggests that $(2,1)$\- or
$(3,1)$-QRAC is a better choice for obtaining a higher cut weight solution
when the state is obtained through optimization. One reason for this
difference comes from the difficulty of optimization. Optimization using the
$(2,1)$\- and $(3,1)$-QRAC Hamiltonians would be easier than using the
$(1,1)$-QRAC because the numbers of qubits of $(2,1)$\- and $(3,1)$-QRAC
Hamiltonians are less than that of $(1,1)$-QRAC. The other reason is the
difference in the landscape of the objective functions. However, the details
have not yet been clarified and are left for future work.

Originally, the MAX-CUT problem in an Ising-type formulation was to find the
eigenstate with the maximum eigenvalue, which is known to be QMA-hard [12].
Remember that the variable with the highest probability is embedded in the
magic state with the highest fidelity to the current quantum state. This
suggests that the MAX-CUT problem has now been changed into a problem of
finding the magic state with the highest fidelity to the optimized quantum
state. Of course, finding the magic state with the largest fidelity is also
intractable, and the largest fidelity sample is not always the MAX-CUT
solution, as seen in Fig. 5(a). However, in the practical case, we can expect
that a better solution can be obtained using $(2,1)$\- and $(3,1)$-QRACs
rather than using $(1,1)$-QRAC, as shown in Fig. 5(b). In other words, the
MAX-CUT problem could be successfully relaxed by the QRAC Hamiltonian.

##### Gset (ExpGset)

Table 2: Comparison with existing methods on MAX-CUT problems from the Gset. All instances have $|V|=800$. All best-known values were reported in [23] with the breakout local search algorithm. GW: Goemans–Williamson algorithm [7], RQAOA: recursive quantum approximate optimization algorithm [3], QRAO: quantum random access optimization [11], RQRAO: recursive quantum random access optimization (proposed). Instance | Best-known | GW | CirCut | QRAO | RQAOA | ${\rm RQRAO}_{N=1}$ | ${\rm RQRAO}_{N=20}$
---|---|---|---|---|---|---|---
G1 | 11624 | 11467 | 11622 | 11453 | 11441 | 11466 | 11562
G6 | 2178 | 2013 | 2178 | 1980 | 1862 | 1980 | 2148
G11 | 564 | 536 | 556 | 540 | 562 | 544 | 564
G14 | 3064 | 2999 | 3052 | 2971 | 2965 | 2965 | 3043
G18 | 992 | 924 | 987 | 929 | 882 | 920 | 980

The cut weight on the Gset dataset for each algorithm is summarized in Table
2. RQRAO with $N=20$ exceeds both the classical approximation algorithm (GW)
and quantum algorithms (QRAO and the RQAOA), but is slightly worse than the
classical heuristic algorithm (CirCut) on all graphs except for G11.

As expected in Sec. 3.2.2, the cut weight of RQRAO with $N=1$ is very similar
to that of QRAO because both assume all random variables are independent. The
cut weight performance of RQRAO increases when the number of ensembles $N$ is
increased, which confirms the effectiveness of using an ensemble, which is
heuristically introduced in our method.

##### Scalability (ExpScale)

Figure 6: Comparison of scalability with existing methods. Ten 3-regular
graphs with the same number of nodes were randomly generated whose edge
weights were $\pm 1$ at random. The markers indicate the means. (a) Relative
cut weight with respect to the cut weight obtained by the GW algorithm; (b)
Algorithm run time. The gray dashed line is a fitted curve using $({\rm
Elapsed\ time\ [sec]})=a({\rm Number\ of\ nodes})^{b}$

.

Figure 6 summarizes the results on $3$-regular graphs. Although the relative
cuts of QRAO and the RQAOA are lower than 1, that of RQRAO exceeds 1 when the
number of nodes is larger than 200 and is larger than that of CirCut (Fig. 6).
In addition, although the absolute runtime of RQRAO is longer than that of
CirCut, the order of empirical time complexity of these methods is close (Fig.
6). To quantitatively compare the empirical time complexity of the methods,
the runtime curves were fitted using

$\displaystyle{({\rm Wall\ time\ [sec]})=a({\rm Number\ of\ nodes})^{b}}$

using only data with a large number of nodes. The fitted values of $b$ are
$1.34$ for CirCut, $2.78$ for the RQAOA, $1.73$ for QRAO, $1.44$ for RQRAO,
and $2.44$ for GW. That is, the empirical time complexity of the proposed
method is much lower than that of GW and the RQAOA. Moreover, it is comparable
to that of CirCut. Because most of the computational time is spent on
optimizing the variational states, the computational time of the proposed
method can be reduced by improving the optimization algorithm. For example,
the runtime cost of RQRAO can be reduced by parallel computation in the
ensemble part. Further reductions in the computational costs of RQRAO are
future work.

##### Sensitivity Analysis (ExpHyp)

Figure 7: Results of the ablation study for the number of embeddings per qubit
$m$. The problem graph instances are $3$-regular graphs with $\pm 1$ edge
weights. For each number of nodes, $10$ graph instances were randomly
generated.

Only the results for varying the number of embeddings per qubit $m$ are
described here. Other hyperparameter sensitivities are summarized in
Supplemental Information F. In Fig. 7, we can see the clear advantage of using
$m>1$ over $m=1$, which is that only the cut weight in the cases of $m>1$
overcomes that of GW. When $m=1$, the optimized state may be stacked into a
bad local minimum, as shown in Fig. 5(b), and this is the reason for the lower
cut weight. Almost no runtime dependency on $m$ is seen in Fig. 7, which means
there is no disadvantage to using $m>1$ as long as the given graph is
sufficiently sparse and the MPS ansatz is adopted.

##  5 Related Work

###  5.1 Classical Algorithm for the MAX-CUT problem

Before reviewing the classical algorithms for the MAX-CUT problem, we
reformulate Definition 1 to make the existing classical algorithms easier to
understand. Because $(-1)^{b_{j}+b_{k}}=v_{j}v_{k}$ where
$v_{j}:=(-1)^{b_{j}}$, the MAX-CUT problem can be reformulated as

$\max_{{\bf X}}\texttt{CW}_{\rm X}({\bf X}),$ (44)

where

$\texttt{CW}_{\rm X}({\bf X}):={\rm tr}({\bf w}{\bf F}({\bf X})),$ (45)

${\bf F}({\bf X}):=({\bf 1}-{\bf X})/2$, and ${\bf X}:={\bf v}{\bf
v}^{\mathsf{T}}$. By definition, ${\bf X}$ is restricted as ${\rm diag}({\bf
X})={\bf 1}$, ${\rm rank}({\bf X})=1$, and ${\bf X}\succeq 0$.

The best known theoretically guaranteed generic algorithm for the MAX-CUT
problem is the GW [7] as stated in Sec. 1. In the GW algorithm, the constraint
that ${\rm rank}{(\bf X})=1$ is removed. This transforms the MAX-CUT problem
to a semidefinite program (SDP), which can be solved in polynomial time for
the number of nodes. After solving the SDP, the solution ${\bf v}\in\\{\pm
1\\}^{|V|}$ needs to be extracted from the obtained ${\bf X}$. When ${\bf X}$
is written in the form ${\bf X}={\bf vv}^{\mathsf{T}}$, the solution is ${\bf
v}$. If not, as is generally the case here, the obtained ${\bf X}$ is
decomposed as ${\bf X}={\bf L}{\bf L}^{\mathsf{T}}$ using, for example,
Cholesky decomposition, where ${\bf L}=\begin{bmatrix}{\bf l}_{1}&\dots&{\bf
l}_{|V|}\end{bmatrix}^{\mathsf{T}}\in\mathbb{R}^{|V|\times|V|}$ is a lower
triangular matrix and ${\bf l}_{i}\in{\mathbb{R}}^{|V|}$ is a unit vector.
Then, the solution $\\{v_{i}\\}$ is determined by the random projection of
${\bf l}_{i}$ as $v_{i}=1$ if ${\bf l}_{i}^{\mathsf{T}}{\bf r}\geq 0$ and
$v_{i}=-1$ otherwise. Here, ${\bf r}\in{\mathbb{R}}^{|V|}$ is a random vector
uniformly distributed on the $|V|$-dimensional unit sphere. This can be viewed
as dividing the set $\\{{\bf l}_{i}\\}^{|V|}_{i=1}$ into two sets by a
hyperplane with ${\bf r}$ as the normal vector. When the hyperplane cutting is
performed uniformly at random, the lower bound of the expected cut weight is
$0.878\,\texttt{CW}({\bf b}^{*})$ [7]. This bound is explained as follows.
According to [7, Theorem 3.1], the MAX-CUT problem can be further reformulated
as

$\\{{\bf l}^{\rm opt}_{i}\\}:=\mathop{\rm argmax}_{\\{{\bf
l}_{i}\\}}{\mathbb{E}}_{{\bf b}\sim\mathcal{P}^{\rm MC}({\bf b};\\{{\bf
l}_{i}\\})}[\texttt{CW}({\bf b})]$ (46)

where

$\mathcal{P}^{\rm MC}(b_{j}\neq b_{k};\\{{\bf l}_{i}\\}):=\frac{{\rm
arccos}({\bf l}_{j}^{\mathsf{T}}{\bf l}_{k})}{\pi}.$ (47)

In the GW algorithm, the objective is rewritten as

$\\{{\bf l}^{*}_{i}\\}:=\mathop{\rm argmax}_{\\{{\bf
l}_{i}\\}}{\mathbb{E}}_{{\bf b}\sim\mathcal{P}^{\rm GW}({\bf b};\\{{\bf
l}_{i}\\})}[\texttt{CW}({\bf b})]$ (48)

where

$\displaystyle{\mathcal{P}^{\rm GW}(b_{j}\neq b_{k};\\{{\bf
l}_{i}\\}):=\frac{1-{\bf l}_{j}^{\mathsf{T}}{\bf l}_{k}}{2}}$
$\leq\frac{1}{0.878}\mathcal{P}^{\rm MC}(b_{j}\neq b_{k};\\{{\bf l}_{i}\\}).$
(49)

The last inequality of Eq. (49) comes from [7, Lemma 3.4, Lemma 3.5]. 222
$\min_{x\in[-1,1]}\frac{{\rm
arccos}(x)}{\pi}\big{(}\frac{1-x}{2}\big{)}^{-1}\approx 0.878$.  The right-
hand side of Eq. (48) is an SDP, and hence it can be solved in polynomial
time. The resulting ${\mathbb{E}}_{{\bf b}\sim\mathcal{P}^{\rm MC}({\bf
b};\\{{\bf l}^{*}_{i}\\})}[\texttt{CW}({\bf b})]$, which is the expected cut
weight of random hyperplane cutting, is guaranteed to be greater than
$0.878\,{\mathbb{E}}_{{\bf b}\sim\mathcal{P}^{\rm MC}({\bf b};\\{{\bf l}^{\rm
opt}_{i}\\})}[\texttt{CW}({\bf b})]=0.878\max_{{\bf b}}\texttt{CW}({\bf b})$
because

$\displaystyle{\mathcal{P}^{\rm MC}(b_{j}\neq b_{k};\\{{\bf l}^{*}_{i}\\})\geq
0.878\,\mathcal{P}^{\rm GW}(b_{j}\neq b_{k};\\{{\bf l}^{*}_{i}\\})}$
$\displaystyle{\geq 0.878\,\mathcal{P}^{\rm GW}(b_{j}\neq b_{k};\\{{\bf
l}^{\rm opt}_{i}\\})}$ $=0.878\,\mathcal{P}^{\rm MC}(b_{j}\neq b_{k};\\{{\bf
l}^{\rm opt}_{i}\\}).$ (50)

Equation (48) is very similar to Eq. (18), but the rounding process of the GW
is not the same as that of Eq. (19). Instead of carrying out ${\bf b}^{\rm
R}=\mathop{\rm argmax}_{{\bf b}}\mathcal{P}^{\rm GW}({\bf b};\\{{\bf
l}^{*}_{i}\\})$ or ${\bf b}^{\rm R}=\mathop{\rm argmax}_{{\bf
b}}\mathcal{P}^{\rm MC}({\bf b};\\{{\bf l}^{*}_{i}\\})$, a number of ${\bf b}$
obeying $\mathcal{P}^{\rm MC}({\bf b};\\{{\bf l}^{*}_{i}\\})$ are sampled by
random hyperplane cutting and then, the best $\texttt{CW}({\bf b})$ is adopted
as the candidate solution.

To solve the GW algorithm faster with less memory, various SDP solvers have
been developed [24]. The runtime of the state-of-the-art SDP solvers are
$\widetilde{O}(\sqrt{n}(mn^{2}+m^{\omega}+n^{\omega})\log(1/\epsilon)))$ [25]
and
$\widetilde{O}((\sqrt{n}(m^{2}+n^{4})+m^{\omega}+n^{2\omega})\log(1/\epsilon)$
[26], where $n$ is the problem size, $m$ the number of constraints, $\omega$
the matrix multiplication constant $\leq 2.372927$, and $\epsilon$ the
relative accuracy. In the case of the MAX-CUT problem, the runtime is
$\widetilde{O}(|V|^{3.5}\log(1/\epsilon))$ because $n=m=|V|$. The overall time
complexity for the GW algorithm depends on the runtime of SDP because the
naïve Cholesky decomposition takes $O(n^{3})=O(|V|^{3})$ to run, which is
faster than the state-of-the-art SDP solver, as long as the number of trials
of random hyperplane cutting is less than the time complexity of the SDP. Note
that the computational complexity of SDP with respect to the number of nodes
can be further reduced to $\widetilde{O}(|V|/\epsilon^{3.5})$ using a specific
algorithm for the MAX-CUT , by taking a disadvantage with respect to the
accuracy [27].

Although there is no approximation ratio guarantee nor computational
complexity bound, there are several classical heuristics that are empirically
faster and have higher approximation ratios than the GW algorithm. The rank-
two relaxation algorithm [16] relaxed ${\rm rank}({\bf X})=1$ to ${\bf X}={\bf
AA}^{\mathsf{T}}$, where ${\bf A}=\begin{bmatrix}{\bf a}_{1}&\dots&{\bf
a}_{|V|}\end{bmatrix}^{\mathsf{T}}$, ${\bf
a}_{j}=\begin{bmatrix}\cos\theta_{j}&\sin\theta_{j}\end{bmatrix}^{\mathsf{T}}$,
and $\theta_{j}\in[0,2\pi)$. Similar to the GW, the objective of the rank-two
relaxation can be reformulated as

$\\{\theta^{*}_{j}\\}:=\mathop{\rm
argmax}_{\\{\theta_{j}\\}}\mathbb{E}_{\mathcal{P}^{\rm R2}({\bf
b};\\{\theta_{j}\\})}[\texttt{CW}({\bf b})],$ (51)

where

$\mathcal{P}^{\rm R2}(b_{j}\neq
b_{k};\\{\theta_{j}\\})=\sin^{2}\Big{(}\frac{\theta_{j}-\theta_{k}}{2}\Big{)}.$
(52)

As with the GW algorithm, the rank-two relaxation does not take the sample
with the highest probability as a candidate. Moreover, because the hyperplane
is only two-dimensional, random hyperplane cutting is not used. Instead,
exhaustive search is used, expressed as

${\bf b}^{\rm R2}=\texttt{CW}({\bf b}^{*}(\alpha^{*})),$ (53)

where

$b^{*}_{j}(\alpha)=\left\\{\begin{array}[]{lll}0&{\rm
for}&\cos(\theta^{*}_{j}+\alpha)\geq 0\\\\[5.69054pt] 1&{\rm
for}&\cos(\theta^{*}_{j}+\alpha)<0\\\ \end{array}\right.$ (54)

and

$\alpha^{*}=\mathop{{\rm argmax}}_{\alpha}\texttt{CW}({\bf b}^{*}(\alpha)).$
(55)

To optimize ${\theta_{j}}$, a gradient-based optimizer is used. In addition,
$\\{1,2\\}$-local search and restarting optimization are incorporated. Here,
$k$-local search is an exhaustive method that searches the entire solution
within Hamming distance $k$ of the current solution. Breakout local search
[23] repeats the 1-local search and perturbation of the current solution using
a tabu list. Once a local search falls into a local optimum, perturbation is
added to the current solution, and then the local search is restarted. In the
perturbation phase, one of three kinds of perturbation is applied with a
specific probability. The tabu list is used to prohibit the perturbation of
the listed nodes so that the same local minimum will not be visited again.

###  5.2 Quantum Algorithms for the MAX-CUT problem

The quantum algorithms for the MAX-CUT problem are categorized into two
research streams. One is quantum approximate optimization algorithms (QAOAs)
[8, 9] and the other is quantum algorithms for SDP [28, 29, 30, 31, 32, 33,
34, 35, 36]. Here, quantum annealing is excluded.

QAOA variants are well summarized in a recent review [10] and are not re-
summarized here.

The quantum SDP algorithms are categorized into two types. One encodes the
solution to a density matrix and the other uses QRAM [37] to store the
solution.

For the density matrix encoding approach, there are two types of encoding
methods. One uses a Gibbs state to represent a solution [28, 29, 30, 31, 33],
which can be regarded as a solution that is encoded in the corresponding
Hamiltonian. The Hamiltonian used in the Gibbs state is updated by the matrix
multiplicative weight update method [38] or the matrix exponentiated gradient
update method [39]. The other density matrix encoding approach uses a
parametrized quantum circuit to represent a solution [34, 35, 36]. The
parameter is optimized by frameworks using the variational quantum algorithm
[40, 41].

The state-of-the-art quantum SDP algorithm [42] is based on the primal-dual
central path method, which uses quantum linear algebra and stores the solution
in QRAM. The runtime of this quantum SDP algorithm is $O((mn^{1.5}+n^{3}){\rm
poly}(\kappa,\log(mn/\epsilon))$ where $\kappa$ is the condition number of the
matrices appearing in the algorithm, and $\epsilon$ is the accuracy parameter.
In the case of the MAX-CUT problem, we have $\widetilde{O}(|V|^{3}{\rm
poly}(\kappa,\log(|V|^{2}/\epsilon))$, which is a slight improvement on its
classical counterpart, which is $\widetilde{O}(|V|^{3.5}\log(1/\epsilon))$.

Although these algorithms have the potential to solve an SDP faster than
classical ones, the approximation ratio is the same as that of the GW
algorithm, which is empirically worse than that of classical heuristics.

Very recently, [43] proposed the quantum-inspired algorithm based on the
observation that ADAPT-QAOA [44] can be approximately expressed as the
Clifford circuit, which can be efficiently simulatable using a classical
computer. Their model showed the comparable cut weight performance to the GW
with several thousand hyperplane cuttings.

###  5.3 QRACs

The concept that encodes $m$ bits into $n$ (qu)bits was originally proposed by
Wiesner as conjugate coding [45] and later rediscovered by Ambainis et al. [1]
as (quantum) random access codes ((Q)RACs). $(m,n,p)$ random access encoding
is defined as the function that maps $m$ classical bits into $n$ (qu)bits with
a decoding probability at least $p$ [1, Definition 1.1]. The minimum number of
$n$ with $p>1/2$ is bounded as $(1-\eta(p))m$ for both classical [1, Theorem
2.1] and quantum [46, Theorem 2.3] settings, where $\eta(p)$ is the binary
entropy function. In the case of $(2,1,p)$\- and $(3,1,p)$-QRACs, there exists
$p=\frac{1}{2}+\frac{1}{2\sqrt{2}}$ [1, Lemma 3.1] and
$p=\frac{1}{2}+\frac{1}{2\sqrt{3}}$ [1, attributed to Chuang], respectively,
which are consistent with our results (Eq. (30)) in the case of
$\mathcal{E}_{j}\in\\{\pm 1\\}$. It was proved that $(2^{n},n,>1/2)$ classical
random access encoding [47, Theorem 5.5] and $(2^{2n},n,>1/2)$ quantum random
access encoding [48, Theorem 2] do not exist. From these observations,
$(3,1)$-QRAC shows the maximum number of classical bits that can be embedded
in one qubit with a decoding probability that is larger than $1/2$. This is
why we only consider $m\in\\{1,2,3\\}$ in this study. In addition to
$(m,1,p)$-QRAC, $(m,2,p)$-QRAC is particularly of interest and well studied in
terms of the theoretical bound [48, 49, 50, 51] and concrete construction
scheme [52]. For arbitrary $(m,n,p)$-QRAC, the upper bound of decoding
probability $p$ is improved by considering shared randomness [53], shared
entanglement [54], and $d$-level system [55, 49]. The information that can be
decoded is not limited to the $m$ bits that are encoded, but also extends to
functions that use $k$-out-of-$m$ bits [56, 57]. Using this extension, QRACs
have also been applied in machine learning to efficiently encode discrete
features [58].

##  6 Conclusion

In this study, we proposed a recursive and quantum-relaxed algorithm to solve
the MAX-CUT problem. The quantum relaxation is achieved using the QRAC, which
is incorporated into the existing recursive MAX-CUT algorithm. To show the
validity of the proposed algorithm, new properties of the quantum-relaxed
Hamiltonian were revealed. According to these findings, we showed that the
MAX-CUT problem can be approximately reduced to the problem finding the
variable with the highest probability, which is determined by the quantum
state optimized through the quantum-relaxed Hamiltonian. We confirmed that the
proposed algorithm is consistent with this objective and includes the existing
methods as special cases. The numerical experiments on graphs with several
hundred nodes showed that the proposed algorithm outperforms the GW algorithm,
which is the best-known classical approximation algorithm, on graphs with more
than 300 nodes both in the approximation ratio and the empirical time
complexity.

## Acknowledgement

This work is supported by MEXT Quantum Leap Flagship Program Grant Number
JPMXS0118067285 and JPMXS0120319794, and JSPS KAKENHI Grant Number 20H05966.
We would like to thank Yohichi Suzuki and Michihiko Sugawara of Keio Quantum
Computing Center and Hiroshi C. Watanabe of Kyushu University for discussing
the theoretical and experimental justification of RQRAO. RR would like to
acknowledge Prof. Hiroshi Imai for information on Gset instances and MAX-CUT .

## References

  * Ambainis et al. [1999] Ambainis, A., Nayak, A., Ta-Shma, A. and Vazirani, U.  Dense quantum coding and a lower bound for 1-way quantum automata.  In _Proceedings of the thirty-first annual ACM symposium on Theory of computing_ , pages 376–383, 1999.  doi: 10.1145/301250.301347.
  * Ambainis et al. [2002] Ambainis, A., Nayak, A., Ta-Shma, A. and Vazirani, U.  Dense quantum coding and quantum finite automata.  _Journal of the ACM (JACM)_ , 49(4):496–511, 2002\.  doi: 10.1145/581771.581773.
  * Bravyi et al. [2020] Bravyi, S., Kliesch, A., Koenig, R. and Tang, E.  Obstacles to variational quantum optimization from symmetry protection.  _Physical review letters_ , 125(26):260505, 2020\.  doi: 10.1103/PhysRevLett.125.260505.
  * Karp [2010] Karp, R.M.  _Reducibility among combinatorial problems_.  Springer, 2010.  doi: 10.1007/978-1-4684-2001-2˙9.
  * Lucas [2014] Lucas, A.  Ising formulations of many NP problems.  _Frontiers in physics_ , 2:5, 2014.  doi: 10.3389/fphy.2014.00005.
  * Glover et al. [2022] Glover, F., Kochenberger, G., Hennig, R. and Du, Y.  Quantum bridge analytics I: a tutorial on formulating and using QUBO models.  _Annals of Operations Research_ , 314(1):141–183, 2022.  doi: 10.1007/s10479-022-04634-2.
  * Goemans and Williamson [1995] Goemans, M.X. and Williamson, D.P.  Improved approximation algorithms for maximum cut and satisfiability problems using semidefinite programming.  _Journal of the ACM (JACM)_ , 42(6):1115–1145, 1995.  doi: 10.1145/227683.227684.
  * Farhi et al. [2014] Farhi, E., Goldstone, J. and Gutmann, S.  A quantum approximate optimization algorithm.  _arXiv preprint arXiv:1411.4028_ , 2014.  doi: 10.48550/arXiv.1411.4028.
  * Wang et al. [2018] Wang, Z., Hadfield, S., Jiang, Z. and Rieffel, E.G.  Quantum approximate optimization algorithm for MaxCut: A fermionic view.  _Physical Review A_ , 97(2):022304, 2018.  doi: 10.1103/PhysRevA.97.022304.
  * Blekos et al. [2023] Blekos, K. et al.  A Review on Quantum Approximate Optimization Algorithm and its Variants.  _arXiv preprint arXiv:2306.09198_ , 2023.  doi: 10.48550/arXiv.2306.09198.
  * Fuller et al. [2021] Fuller, B. et al.  Approximate solutions of combinatorial problems via quantum relaxations.  _arXiv preprint arXiv:2111.03167_ , 2021.  doi: 10.48550/arXiv.2111.03167.
  * Kempe et al. [2006] Kempe, J., Kitaev, A. and Regev, O.  The complexity of the local Hamiltonian problem.  _Siam journal on computing_ , 35(5):1070–1097, 2006.  doi: 10.1007/978-3-540-30538-5˙31.
  * Khot [2002] Khot, S.  On the power of unique 2-prover 1-round games.  In _Proceedings of the thiry-fourth annual ACM symposium on Theory of computing_ , pages 767–775, 2002.  doi: 10.1145/509907.510017.
  * Bravyi et al. [2022] Bravyi, S., Kliesch, A., Koenig, R. and Tang, E.  Hybrid quantum-classical algorithms for approximate graph coloring.  _Quantum_ , 6:678, 2022.  doi: 10.22331/q-2022-03-30-678.
  * Helmberg and Rendl [2000] Helmberg, C. and Rendl, F.  A spectral bundle method for semidefinite programming.  _SIAM Journal on Optimization_ , 10(3):673–696, 2000.  doi: 10.1137/S1052623497328987.
  * Burer et al. [2002] Burer, S., Monteiro, R.D. and Zhang, Y.  Rank-two relaxation heuristics for max-cut and other binary quadratic programs.  _SIAM Journal on Optimization_ , 12(2):503–521, 2002.  doi: 10.1137/S1052623400382467.
  * Prim [1957] Prim, R.C.  Shortest Connection Networks And Some Generalizations.  _Bell System Technical Journal_ , 36(6):1389–1401, November 1957.  doi: 10.1002/j.1538-7305.1957.tb01515.x.
  * Karger et al. [1995] Karger, D.R., Klein, P.N. and Tarjan, R.E.  A randomized linear-time algorithm to find minimum spanning trees.  _Journal of the ACM (JACM)_ , 42(2):321–328, 1995\.  doi: 10.1145/201019.201022.
  * Rinaldi [1998] Rinaldi, G.  Rudy, 1998.  https://www-user.tu-chemnitz.de/~helmberg/rudy.tar.gz.
  * Dunning et al. [2018] Dunning, I., Gupta, S. and Silberholz, J.  What works best when? A systematic evaluation of heuristics for Max-Cut and QUBO.  _INFORMS Journal on Computing_ , 30(3):608–624, 2018.  doi: 10.1287/ijoc.2017.0798.
  * Paszke et al. [2019] Paszke, A. et al.  PyTorch: An imperative style, high-performance deep learning library.  _Advances in neural information processing systems_ , 32, 2019.  doi: 10.48550/arXiv.1912.01703.
  * Liu and Nocedal [1989] Liu, D.C. and Nocedal, J.  On the limited memory BFGS method for large scale optimization.  _Mathematical programming_ , 45(1-3):503–528, 1989.  doi: 10.1007/BF01589116.
  * Benlic and Hao [2013] Benlic, U. and Hao, J.K.  Breakout local search for the max-cutproblem.  _Engineering Applications of Artificial Intelligence_ , 26(3):1162–1173, 2013.  doi: 10.1016/j.engappai.2012.09.001.
  * Majumdar et al. [2020] Majumdar, A., Hall, G. and Ahmadi, A.A.  Recent scalability improvements for semidefinite programming with applications in machine learning, control, and robotics.  _Annual Review of Control, Robotics, and Autonomous Systems_ , 3:331–360, 2020.  doi: 10.1146/annurev-control-091819-074326.
  * Jiang et al. [2020] Jiang, H., Kathuria, T., Lee, Y.T., Padmanabhan, S. and Song, Z.  A faster interior point method for semidefinite programming.  In _2020 IEEE 61st annual symposium on foundations of computer science (FOCS)_ , pages 910–918. IEEE, 2020.  doi: 10.1109/FOCS46700.2020.00089.
  * Huang et al. [2022a] Huang, B., Jiang, S., Song, Z., Tao, R. and Zhang, R.  Solving SDP faster: A robust IPM framework and efficient implementation.  In _2022 IEEE 63rd Annual Symposium on Foundations of Computer Science (FOCS)_ , pages 233–244. IEEE, 2022a.  doi: 10.1109/FOCS54457.2022.00029.
  * Lee and Padmanabhan [2020] Lee, Y.T. and Padmanabhan, S.  An $\widetilde{\mathcal{O}}(m/\varepsilon^{3.5})$-cost algorithm for semidefinite programs with diagonal constraints.  In Abernethy, J. and Agarwal, S., editors, _Proceedings of Thirty Third Conference on Learning Theory_ , volume 125 of _Proceedings of Machine Learning Research_ , pages 3069–3119. PMLR, 09–12 Jul 2020.  URL https://proceedings.mlr.press/v125/lee20c.html.
  * Brandao and Svore [2017] Brandao, F.G. and Svore, K.M.  Quantum speed-ups for solving semidefinite programs.  In _2017 IEEE 58th Annual Symposium on Foundations of Computer Science (FOCS)_ , pages 415–426. IEEE, 2017.  doi: 10.1109/FOCS.2017.45.
  * Van Apeldoorn et al. [2020] Van Apeldoorn, J., Gilyén, A., Gribling, S. and de Wolf, R.  Quantum SDP-solvers: Better upper and lower bounds.  _Quantum_ , 4:230, 2020.  doi: 10.22331/q-2020-02-14-230.
  * Brandão et al. [2019] Brandão, F.G.S.L., Kalev, A., Li, T., Lin, C.Y.Y., Svore, K.M. and Wu, X.  Quantum SDP Solvers: Large Speed-Ups, Optimality, and Applications to Quantum Learning.  In Baier, C., Chatzigiannakis, I., Flocchini, P. and Leonardi, S., editors, _46th International Colloquium on Automata, Languages, and Programming (ICALP 2019)_ , volume 132 of _Leibniz International Proceedings in Informatics (LIPIcs)_ , pages 27:1–27:14, Dagstuhl, Germany, 2019\. Schloss Dagstuhl–Leibniz-Zentrum fuer Informatik.  ISBN 978-3-95977-109-2.  doi: 10.4230/LIPIcs.ICALP.2019.27.  URL http://drops.dagstuhl.de/opus/volltexte/2019/10603.
  * van Apeldoorn and Gilyén [2019] van Apeldoorn, J. and Gilyén, A.  Improvements in Quantum SDP-Solving with Applications.  In Baier, C., Chatzigiannakis, I., Flocchini, P. and Leonardi, S., editors, _46th International Colloquium on Automata, Languages, and Programming (ICALP 2019)_ , volume 132 of _Leibniz International Proceedings in Informatics (LIPIcs)_ , pages 99:1–99:15, Dagstuhl, Germany, 2019\. Schloss Dagstuhl–Leibniz-Zentrum fuer Informatik.  ISBN 978-3-95977-109-2.  doi: 10.4230/LIPIcs.ICALP.2019.99.  URL http://drops.dagstuhl.de/opus/volltexte/2019/10675.
  * Kerenidis and Prakash [2020] Kerenidis, I. and Prakash, A.  A quantum interior point method for LPs and SDPs.  _ACM Transactions on Quantum Computing_ , 1(1):1–32, 2020.  doi: 10.1145/3406306.
  * Brandao et al. [2022] Brandao, F.G.L., Kueng, R. and França, D.S.  Faster quantum and classical SDP approximations for quadratic binary optimization.  _Quantum_ , 6:625, 2022.  doi: 10.22331/q-2022-01-20-625.
  * Bharti et al. [2022] Bharti, K., Haug, T., Vedral, V. and Kwek, L.C.  Noisy intermediate-scale quantum algorithm for semidefinite programming.  _Physical Review A_ , 105(5):052445, 2022.  doi: 10.1103/PhysRevA.105.052445.
  * Patel et al. [2021] Patel, D., Coles, P.J. and Wilde, M.M.  Variational quantum algorithms for semidefinite programming.  _arXiv preprint arXiv:2112.08859_ , 2021.  doi: 10.48550/arXiv.2112.08859.
  * Patti et al. [2023] Patti, T.L., Kossaifi, J., Anandkumar, A. and Yelin, S.F.  Quantum Goemans-Williamson Algorithm with the Hadamard Test and Approximate Amplitude Constraints.  _Quantum_ , 7:1057, 2023.  doi: https://doi.org/10.22331/q-2023-07-12-1057.
  * Giovannetti et al. [2008] Giovannetti, V., Lloyd, S. and Maccone, L.  Quantum random access memory.  _Physical review letters_ , 100(16):160501, 2008\.  doi: 10.1103/PhysRevLett.100.160501.
  * Arora and Kale [2007] Arora, S. and Kale, S.  A combinatorial, primal-dual approach to semidefinite programs.  In _Proceedings of the thirty-ninth annual ACM symposium on Theory of computing_ , pages 227–236, 2007.  doi: 10.1145/2837020.
  * Tsuda et al. [2005] Tsuda, K., Rätsch, G. and Warmuth, M.K.  Matrix exponentiated gradient updates for on-line learning and Bregman projection.  _Journal of Machine Learning Research_ , 6(Jun):995–1018, 2005.  Retrieved from https://dl.acm.org/doi/10.5555/1046920.1088706.
  * Peruzzo et al. [2014] Peruzzo, A. et al.  A variational eigenvalue solver on a photonic quantum processor.  _Nature communications_ , 5(1):1–7, 2014.  doi: 10.1038/ncomms5213.
  * Cerezo et al. [2021] Cerezo, M. et al.  Variational quantum algorithms.  _Nature Reviews Physics_ , 3(9):625–644, 2021\.  doi: 10.1038/s42254-021-00348-9.
  * Huang et al. [2022b] Huang, B., Jiang, S., Song, Z., Tao, R. and Zhang, R.  A faster quantum algorithm for semidefinite programming via robust IPM framework.  _arXiv preprint arXiv:2207.11154_ , 2022b.  doi: 10.48550/arXiv.2207.11154.
  * Muñoz-Arias et al. [2023] Muñoz-Arias, M.H., Kourtis, S. and Blais, A.  Low-depth clifford circuits approximately solve maxcut.  _arXiv preprint arXiv:2310.15022_ , 2023.  doi: 10.48550/arXiv.2310.15022.
  * Zhu et al. [2022] Zhu, L. et al.  Adaptive quantum approximate optimization algorithm for solving combinatorial problems on a quantum computer.  _Physical Review Research_ , 4(3):033029, 2022\.  doi: 10.1103/PhysRevResearch.4.033029.
  * Wiesner [1983] Wiesner, S.  Conjugate coding.  _ACM Sigact News_ , 15(1):78–88, 1983.  doi: 10.1145/1008908.1008920.
  * Nayak [1999] Nayak, A.  Optimal lower bounds for quantum automata and random access codes.  In _40th Annual Symposium on Foundations of Computer Science (Cat. No. 99CB37039)_ , pages 369–376. IEEE, 1999.  doi: 10.1109/SFFCS.1999.814608.
  * Iwama et al. [2007] Iwama, K., Nishimura, H., Raymond, R. and Yamashita, S.  Unbounded-error one-way classical and quantum communication complexity.  In _Automata, Languages and Programming: 34th International Colloquium, ICALP 2007, Wrocław, Poland, July 9-13, 2007. Proceedings 34_ , pages 110–121. Springer, 2007.  doi: 10.1007/978-3-540-73420-8˙12.
  * Hayashi et al. [2006] Hayashi, M., Iwama, K., Nishimura, H., Raymond, R. and Yamashita, S.  (4, 1)-quantum random access coding does not exist—one qubit is not enough to recover one of four bits.  _New Journal of Physics_ , 8(8):129, 2006.  doi: 10.1088/1367-2630/8/8/129.
  * Liabøtrø [2017] Liabøtrø, O.  Improved classical and quantum random access codes.  _Physical Review A_ , 95(5):052315, 2017.  doi: 10.1103/PhysRevA.95.052315.
  * Imamichi and Raymond [2018] Imamichi, T. and Raymond, R.  Constructions of quantum random access codes.  In _Asian Quantum Information Symposium (AQIS)_ , volume 66, 2018\.  Retrieved from https://research.ibm.com/publications/constructions-of-quantum-random-access-codes.
  * Mančinska and Storgaard [2022] Mančinska, L. and Storgaard, S.A.  The geometry of Bloch space in the context of quantum random access codes.  _Quantum Information Processing_ , 21(4):143, 2022\.  doi: 10.1007/s11128-022-03470-4.
  * Teramoto et al. [2023] Teramoto, K., Raymond, R., Wakakuwa, E. and Imai, H.  Quantum-Relaxation Based Optimization Algorithms: Theoretical Extensions.  _arXiv preprint arXiv:2302.09481_ , 2023.  doi: 10.48550/arXiv.2302.09481.
  * Ambainis et al. [2008] Ambainis, A., Leung, D., Mancinska, L. and Ozols, M.  Quantum random access codes with shared randomness.  _arXiv preprint arXiv:0810.2937_ , 2008.  doi: 10.48550/arXiv.0810.2937.
  * Pawłowski and Żukowski [2010] Pawłowski, M. and Żukowski, M.  Entanglement-assisted random access codes.  _Physical Review A_ , 81(4):042326, 2010.  doi: 10.1103/PhysRevA.81.042326.
  * Tavakoli et al. [2015] Tavakoli, A., Hameedi, A., Marques, B. and Bourennane, M.  Quantum random access codes using single d-level systems.  _Physical review letters_ , 114(17):170502, 2015\.  doi: 10.1103/PhysRevLett.114.170502.
  * Ben-Aroya et al. [2008] Ben-Aroya, A., Regev, O. and De Wolf, R.  A hypercontractive inequality for matrix-valued functions with applications to quantum computing and ldcs.  In _2008 49th Annual IEEE Symposium on Foundations of Computer Science_ , pages 477–486. IEEE, 2008.  doi: 10.1109/FOCS.2008.45.
  * Doriguello and Montanaro [2021] Doriguello, J.F. and Montanaro, A.  Quantum random access codes for boolean functions.  _Quantum_ , 5:402, 2021.  doi: 10.22331/q-2021-03-07-402.
  * Yano et al. [2021] Yano, H., Suzuki, Y., Itoh, K.M., Raymond, R. and Yamamoto, N.  Efficient discrete feature encoding for variational quantum classifier.  _IEEE Transactions on Quantum Engineering_ , 2:1–14, 2021\.  doi: 10.1109/QCE49297.2020.00012.
  * Huang et al. [2020] Huang, H.Y., Kueng, R. and Preskill, J.  Predicting many properties of a quantum system from very few measurements.  _Nature Physics_ , 16(10):1050–1057, 2020.  doi: 10.1038/s41567-020-0932-7.
  * Schollwöck [2011] Schollwöck, U.  The density-matrix renormalization group in the age of matrix product states.  _Annals of physics_ , 326(1):96–192, 2011.  doi: 10.1016/j.aop.2010.09.012.
  * Toh et al. [1999] Toh, K.C., Todd, M.J. and Tütüncü, R.H.  SDPT3—a MATLAB software package for semidefinite programming, version 1.3.  _Optimization methods and software_ , 11(1-4):545–581, 1999.  doi: 10.1080/10556789908805762.
  * Tütüncü et al. [2003] Tütüncü, R.H., Toh, K.C. and Todd, M.J.  Solving semidefinite-quadratic-linear programs using SDPT3.  _Mathematical programming_ , 95:189–217, 2003.  doi: 10.1007/s10107-002-0347-5.

##  Appendix A Classical Shadow with Random Magic Measurements

Let us define $(m,1)$-random magic measurements, which is referred to as magic
state rounding in [11]. In this section, only $(3,1)$-random magic
measurements are introduced, as $(2,1)$\- and $(1,1)$-random magic
measurements are easily derived similarly.

######  Definition 3 (Random magic measurements [11])

Let $\mu^{[q]\pm}_{i}$ be the single qubit $(3,1)$-magic state corresponding
to the $q$th qubit defined as follows:

$\displaystyle{\mu_{1}^{[q]+}:=\mu_{3}^{[q]}(0,0,0),\ \ \
\mu_{1}^{[q]-}:=\mu_{3}^{[q]}(1,1,1),}$
$\displaystyle{\mu_{2}^{[q]+}:=\mu_{3}^{[q]}(0,1,1),\ \ \
\mu_{2}^{[q]-}:=\mu_{3}^{[q]}(1,0,0),}$
$\displaystyle{\mu_{3}^{[q]+}:=\mu_{3}^{[q]}(1,0,1),\ \ \
\mu_{3}^{[q]-}:=\mu_{3}^{[q]}(0,1,0),}$
$\displaystyle{\mu_{4}^{[q]+}:=\mu_{3}^{[q]}(1,1,0),\ \ \
\mu_{4}^{[q]-}:=\mu_{3}^{[q]}(0,0,1),}$

where $\mu_{3}^{[q]}(b^{[q]}_{X},b^{[q]}_{Y},b^{[q]}_{Z})$ is defined as Eq.
(11). Random magic measurement is the protocol that obtains ${\bf
b}\in\\{0,1\\}^{3n}$ as follows:

  1. 1.

Select $i\in\\{1,2,3,4\\}$ uniformly at random for each $q$.

  2. 2.

Measure the state $\rho$ in the basis $\\{\mu^{[q]\pm}_{i}\\}$ for all qubits
$q\in[n]$ with the selected $i$.

  3. 3.

Convert the measurement outcomes to the corresponding classical bits as listed
in Table 3.

Table 3: Measurement outcome of the random magic measurement and corresponding classical bits. Here, $i_{q}$ is the selected $i\in\\{1,2,3,4\\}$ for qubit $q$. $i_{q}$ | Outcome | $\big{(}b^{[q]}_{X},b^{[q]}_{Y},b^{[q]}_{Z}\big{)}$
---|---|---
1 | 0 | $(0,0,0)$
1 | 1 | $(1,1,1)$
2 | 0 | $(0,1,1)$
2 | 1 | $(1,0,0)$
3 | 0 | $(1,0,1)$
3 | 1 | $(0,1,0)$
4 | 0 | $(1,1,0)$
4 | 1 | $(0,0,1)$

The bit-strings ${\bf b}$ obtained by a random magic measurement follows the
probability distribution $\mathcal{P}_{\rm f}({\bf b})$ (see the proof of
Theorem 2).

Because the expectation of the $k$-local Pauli observable in Eq. (20) is
written as the expectation of the random variable

$\hat{o}:=3^{k/2}(-1)^{\sum^{k}_{i=1}b^{[q^{(i)}]}_{P}},$ (S.1)

the classical shadow [59] based on random magic measurements can be
considered. That is, the expectation value of the $k$-local Pauli observables
can be estimated by averaging $\hat{o}$ computed from the corrected bits
$\\{b^{[q]}_{P}\\}$ by the random magic measurements. The upper bound of the
variance of random variable $\hat{o}$ is readily derived as

${\rm
Var}[\hat{o}]=\mathbb{E}[\hat{o}^{2}]-(\mathbb{E}[\hat{o}])^{2}\leq\mathbb{E}[\hat{o}^{2}]=3^{k},$
(S.2)

which is equivalent to that of random Pauli measurements [59]. The classical
snapshot based on random magic measurement is

$\hat{\rho}=\bigotimes^{n-1}_{q=0}\Big{(}3U^{[q]\dagger}|\hat{b}^{[q]}\rangle\langle\hat{b}^{[q]}|U^{[q]}-I\Big{)},\
\ \ U^{[q]}\sim\mathcal{U}_{\rm ms},$ (S.3)

where $\hat{b}^{[q]}$ is the measurement outcome of the $q$th qubit. Equation
(S.3) is the same as that of random Pauli measurements except for the unitary
ensemble $\mathcal{U}_{\rm ms}$. In the case of random magic measurements, the
unitary ensemble is $\mathcal{U}_{\rm ms}=\\{U_{i}\\}^{4}_{i=1}$, where
$\mu_{i}^{+}:=U_{i}|0\rangle\langle 0|U_{i}^{\dagger}$ and
$\mu_{i}^{-}:=U_{i}|1\rangle\langle 1|U_{i}^{\dagger}$. Note that $U_{1}$ can
be expressed as $U_{1}^{\dagger}=e^{{\rm i}tX}e^{{\rm i}sZ}$ where $s=\pi/8$
and $t={\rm arccos}(\sqrt{(1+1/\sqrt{3})/2})$, and
$U_{2}^{\dagger}=U^{\dagger}_{1}X$, $U^{\dagger}_{3}=U^{\dagger}_{1}Y$,
$U^{\dagger}_{4}=U^{\dagger}_{1}Z$ [11].

Figure S.1: Accuracy comparison of a classical shadow obtained using random
Pauli measurements and that obtained using $(3,1)$-random magic measurements.
(a) Absolute error of the expectation value of a $2$-local Pauli observable in
a ten qubit system; (b) Fidelity between the true state and the one obtained
using state tomography in a five qubit system.

To demonstrate the performance of random magic measurements, two numerical
experiments were carried out.

##### Expectation Estimation

A state $\rho$ is sampled from ten-qubit Haar’s random states and is used
throughout the experiments. A hundred 2-local Pauli observables are randomly
generated and their expectations are evaluated using random Pauli measurements
and random magic measurements, respectively. Although the expectation of the
$k$-local Pauli observable is within $[-1,1]$, the range of the expectation
estimated by random magic measurements is $[-3^{k/2},3^{k/2}]$. To mitigate
the anomaly estimation, the expectation value estimated by the random magic
measurements is truncated so that the value is within $[-1,1]$. In Fig. S.1,
the accuracy of the expectation evaluation of 2-local Pauli observable is
shown. The dependence of the accuracy of the expectation evaluation on the
number of shots is almost the same for both the random Pauli measurements and
random magic measurements, which is as expected because the variance upper
bounds of these methods are the same.

##### State Tomography

A state $\rho$ is sampled from five-qubit Haar’s random states and used
throughout the experiments. $N$ snapshots, $\\{\hat{\rho}_{i}\\}^{N}_{i=1}$,
are collected using random Pauli measurements and the random magic
measurements. The candidate matrix is then defined as
$\hat{\rho}=\sum^{N}_{i=1}\hat{\rho}_{i}/N$. Here, in general, $\hat{\rho}$
does not satisfy the property of the density matrix, i.e., ${\rm
tr}(\hat{\rho})=1$ and $\hat{\rho}\succeq 0$. To avoid improper estimation of
the fidelity, the nearest state to $\hat{\rho}$ that satisfied the property of
the density matrix is obtained by the optimization, and is defined as the
result of the state tomography. The optimization process is as follows. The
$n$-qubit nearest state is parametrized as
$|\hat{\psi}({\boldsymbol{\theta}})\rangle_{k}=f_{k}/\sqrt{\sum^{2^{n}}_{k^{\prime}=1}f^{*}_{k^{\prime}}f_{k^{\prime}}}$
where $f_{k}=\theta^{\rm Re}_{k}+{\rm i}\theta^{\rm Im}_{k}$ and
${\boldsymbol{\theta}}=\\{\theta^{\rm Re}_{k},\theta^{\rm
Im}_{k}\\}^{2^{n}}_{k=1}$. Here,
$|\hat{\psi}({\boldsymbol{\theta}})\rangle_{k}$ denotes the $k$th component of
$|\hat{\psi}({\boldsymbol{\theta}})\rangle$ in the computational basis. The
objective function is defined as
$\mathcal{L}({\boldsymbol{\theta}}):=\|\hat{\rho}-|\hat{\psi}({\boldsymbol{\theta}})\rangle\langle\hat{\psi}({\boldsymbol{\theta}})|\|_{\rm
F}$ where $\|\bullet\|_{\rm F}$ is the Frobenius norm. Parameter
${\boldsymbol{\theta}}$ is optimized by L-BFGS to minimize
$\mathcal{L}({\boldsymbol{\theta}})$. Then, the tomography result is defined
as
$\overline{\rho}:=|\hat{\psi}({\boldsymbol{\theta}}^{*})\rangle\langle\hat{\psi}({\boldsymbol{\theta}}^{*})|$,
where ${\boldsymbol{\theta}}^{*}$ is the optimized parameter. According to
Fig. S.1, the dependence of the state tomography accuracy on the number of
shots is almost the same for both the random Pauli measurements and random
magic measurements.

##  Appendix B Node Isolation Problem

Figure S.2:  Schematic illustration of the node isolation process. (a) The red
and green nodes represent the node to be deleted and the node to be kept,
respectively. The red edge indicates the candidate edge. The ensemble edge
energy of the candidate edge $\mathcal{E}$ is negative. (b) After removing the
candidate edge. (c) Multiplying the sign of the ensemble edge energy of the
candidate edge and reconnecting it to the green node. (d) After removing the
zero-weight edge.

When the edge weights of problem graph $G$ are integers, the number of
isolated nodes, which are defined as nodes with undetermined parity with
respect to other nodes, will increase as the iteration progresses. An example
is illustrated in Fig. S.2. In Fig. S.2, the node to be deleted and the node
to be kept are shown in red and green, respectively. The red edge indicates
the candidate whose ensemble edge energy $\mathcal{E}$ is negative. After
removing the candidate edge (Fig. S.2), the edge that was connected to the red
node at one end is reconnected to the green node (Fig. S.2). At this time, the
edge weight is multiplied by the sign of $\mathcal{E}$. As a result, the
weight of the edge the bottom of the illustration is updated to zero (Fig.
S.2), which means that the black node becomes isolated. If there are $N_{\rm
iso}$ isolated nodes after the parity determination procedure, a brute-force
search must be carried out on the $M+N_{\rm iso}$ nodes, which is intractable
when there are several tens of isolated nodes. One easy way to mitigate the
node isolation problem is to apply a small perturbation to all edge weights at
random. When perturbed edge weights are used, the probability that a
reconnected edge has zero weight becomes negligibly small.

##  Appendix C Proofs

  * Proof of Theorem 2.

Consider a quantum channel $\mathcal{E}$ that measures an $n$ qubit state
$\rho$ with the measurement basis $\mu_{m}({\bf b})$, which is an $n$ qubit
product state of magic states selected uniformly at random:

$\mathcal{E}(\rho):=\frac{1}{2^{(m-1)n}}\sum_{{\bf b}\in\\{0,1\\}^{mn}}{\rm
tr}(\mu_{m}({\bf b})\rho)\mu_{m}({\bf b}),$ (S.4)

where $2^{(m-1)n}$ is the total number of $(m,1)$-magic state bases living in
$n$-qubit quantum state space. According to [11], the following holds:

${\rm tr}\big{(}P^{\otimes k}\mathcal{E}(\rho)\big{)}=\frac{1}{m^{k}}{\rm
tr}\big{(}P^{\otimes k}\rho\big{)}$ (S.5)

In addition, it is easily confirmed that

${\rm tr}(P^{\otimes k}\mu_{m}({\bf
b}))=\frac{(-1)^{\sum^{k}_{i=1}b^{[q^{(i)}]}_{P}}}{m^{k/2}}.$ (S.6)

Substituting Eqs. (S.4) and (S.6) into Eq. (S.5) yields

$\displaystyle{{\rm tr}(P^{\otimes k}\rho)=m^{k}{\rm tr}(P^{\otimes
k}\mathcal{E}(\rho))}$ $\displaystyle{=m^{k}{\rm tr}\Bigg{(}P^{\otimes
k}\frac{1}{2^{(m-1)n}}\sum_{{\bf b}\in\\{0,1\\}^{mn}}{\rm tr}(\mu_{m}({\bf
b})\rho)\mu_{m}({\bf b})\Bigg{)}}$
$\displaystyle{=\frac{m^{k}}{2^{(m-1)n}}\sum_{{\bf b}\in\\{0,1\\}^{mn}}{\rm
tr}(P^{\otimes k}\mu_{m}({\bf b}))f_{m}({\bf b},\rho)}$
$=\frac{m^{k/2}}{2^{(m-1)n}}\sum_{{\bf
b}\in\\{0,1\\}^{mn}}(-1)^{\sum^{k}_{i=1}b^{[q^{(i)}]}_{P}}f_{m}({\bf
b},\rho),$ (S.7)

where

$f_{m}({\bf b},\rho):={\rm tr}(\mu_{m}({\bf b})\rho)$ (S.8)

is a fidelity of $\mu_{m}({\bf b})$ and $\rho$. Meanwhile, it is easily
confirmed that

$1={\rm tr}(\mathcal{E}(\rho))=\frac{1}{2^{(m-1)n}}\sum_{{\bf
b}^{\prime}\in\\{0,1\\}^{mn}}f({\bf b}^{\prime},\rho).$ (S.9)

Then, the following holds:

$\displaystyle{{\rm tr}(P^{\otimes k}\rho)=m^{k/2}\,\frac{\sum_{{\bf
b}\in\\{0,1\\}^{mn}}(-1)^{\sum^{k}_{i=1}b^{[q^{(i)}]}_{P}}f_{m}({\bf
b},\rho)}{\sum_{{\bf b}^{\prime}\in\\{0,1\\}^{mn}}f_{m}({\bf
b}^{\prime},\rho)}}$ $=m^{k/2}\,\mathbb{E}_{{\bf b}\sim\mathcal{P}m({\bf
b};\rho)}\Big{[}(-1)^{\sum^{k}_{i=1}b^{[q^{(i)}]}_{P}}\Big{]},$ (S.10)

where

$\mathcal{P}_{m}({\bf b};\rho):=\frac{f_{m}({\bf b},\rho)}{\sum_{{\bf
b}^{\prime}\in\\{0,1\\}^{mn}}f_{m}({\bf b}^{\prime},\rho)}$ (S.11)

$\Box$

  * Proof of Corollary 1.

Consider that a binary variable of node $j$, $b_{j}$, is mapped to
$(q^{(j)},P^{(j)})$, where $q^{(j)}$ is a qubit index,
$P^{(j)}\in\mathfrak{P}\subset\\{X,Y,Z\\}$ and
$|\mathfrak{P}|=m\in\\{1,2,3\\}$ (see Sec. 2.3). Equation (S.7) in the case of
$k=2$ reduces

${\rm
tr}\big{(}P^{(j)}_{q^{(j)}}P^{(k)}_{q^{(k)}}\rho\big{)}=m\,\mathbb{E}_{{\bf
b}\sim\mathcal{P}_{m}({\bf b};\rho)}\big{[}(-1)^{b_{j}+b_{k}}\big{]}.$ (S.12)

The left-hand side of Eq. (S.12) is written as ${\rm tr}\big{(}P_{\langle
j\rangle}P_{\langle k\rangle}\rho\big{)}$ in the main body. Substituting the
$(m,1)$-QRAC Hamiltonian (Eq. (10)) into Eq. (S.12) leads to

${\rm tr}(H_{m}\,\rho)=\sum_{(j,k)\in E}\frac{{\rm tr}(\rho)-m{\rm
tr}(P_{\langle j\rangle}P_{\langle k\rangle}\,\rho)}{2}w_{jk}=\sum_{(j,k)\in
E}\frac{1-m^{2}\mathbb{E}_{{\bf b}\sim\mathcal{P}_{m}({\bf
b};\rho)}\big{[}(-1)^{b_{j}+b_{k}}\big{]}}{2}w_{jk}.$ (S.13)

According to the definition of $\texttt{CW}({\bf b})$ (Eq. (2)),

$\sum_{(j,k)\in E}\frac{(-1)^{b_{j}+b_{k}}}{2}w_{jk}=\sum_{(j,k)\in
E}\frac{1}{2}w_{jk}-\texttt{CW}({\bf b}).$ (S.14)

Substituting Eq. (S.14) into (S.13) gives

${\rm tr}(H_{\rm QRAC}\,\rho)=m^{2}\mathbb{E}_{{\bf b}\sim\mathcal{P}_{m}({\bf
b};\rho)}\big{[}\texttt{CW}({\bf b})]-\frac{m^{2}-1}{2}\sum_{(j,k)\in
E}w_{jk}.$ (S.15)

$\Box$

  * Proof of Corollary 2.

According to Eq. (10),

$\begin{array}[]{ll}\displaystyle{H^{2}_{m}}&\displaystyle{=\Big{(}\sum_{(j,k)\in
E}\frac{I-mP_{\langle j\rangle}P_{\langle
k\rangle}}{2}w_{jk}\Big{)}\cdot\Big{(}\sum_{(j^{\prime},k^{\prime})\in
E}\frac{I-mP_{\langle j^{\prime}\rangle}P_{\langle
k^{\prime}\rangle}}{2}w_{j^{\prime}k^{\prime}}\Big{)}}\\\\[14.22636pt]
&\displaystyle{=\frac{m^{2}}{4}\sum_{(j,k)\in
E}\sum_{(j^{\prime},k^{\prime})\in E}P_{\langle j\rangle}P_{\langle
k\rangle}P_{\langle j^{\prime}\rangle}P_{\langle
k^{\prime}\rangle}w_{jk}w_{j^{\prime}k^{\prime}}-\frac{mW}{2}\sum_{(j,k)\in
E}P_{\langle j\rangle}P_{\langle
k\rangle}w_{jk}}+\frac{1}{4}W^{2}\end{array},$

where $W:=\sum_{(j,k)\in E}w_{jk}$. The expectation value of $H^{2}_{m}$ for
an arbitrary state $\rho$ is calculated as

$\begin{array}[]{ll}\displaystyle{{\rm
tr}(H^{2}_{m}\rho)}&\displaystyle{=\frac{m^{2}}{4}\sum_{(j,k)\in
E}\sum_{(j^{\prime},k^{\prime})\in E}{\rm tr}(P_{\langle j\rangle}P_{\langle
k\rangle}P_{\langle j^{\prime}\rangle}P_{\langle
k^{\prime}\rangle}w_{jk}w_{j^{\prime}k^{\prime}}\rho)-\frac{mW}{2}\sum_{(j,k)\in
E}{\rm tr}(P_{\langle j\rangle}P_{\langle
k\rangle}w_{jk}}\rho)+\frac{1}{4}W^{2}\\\\[14.22636pt]
&\displaystyle{=\frac{m^{2}}{4}\sum_{(j,k)\in
E}\sum_{(j^{\prime},k^{\prime})\in E}{\rm tr}(P_{\langle j\rangle}P_{\langle
k\rangle}P_{\langle j^{\prime}\rangle}P_{\langle
k^{\prime}\rangle}w_{jk}w_{j^{\prime}k^{\prime}}\rho)+W{\rm
tr}(H_{m}\rho)-\frac{1}{4}W^{2}}.\end{array}$

Meanwhile, from Eq. (22),

$\begin{array}[]{ll}\displaystyle{(\texttt{CW}({\bf
b}))^{2}}&\displaystyle{=\Big{(}\sum_{(j,k)\in
E}\frac{1-m^{2}(-1)^{b_{j}+b_{k}}}{2m^{2}}w_{jk}+\frac{m^{2}-1}{2m^{2}}W\Big{)}\Big{(}\sum_{(j^{\prime},k^{\prime})\in
E}\frac{1-m^{2}(-1)^{b_{j^{\prime}}+b_{k^{\prime}}}}{2m^{2}}w_{j^{\prime}k^{\prime}}+\frac{m^{2}-1}{2m^{2}}W\Big{)}}\\\\[14.22636pt]
&\displaystyle{=\frac{1}{4}\sum_{(j,k)\in E}\sum_{(j^{\prime},k^{\prime})\in
E}(-1)^{b_{j}+b_{k}+b_{j^{\prime}}+b_{k^{\prime}}}w_{jk}w_{j^{\prime}k^{\prime}}-\frac{1}{2}W\sum_{(j,k)\in
E}(-1)^{b_{j}+b_{k}}w_{jk}+\frac{1}{4}W^{2}}.\end{array}$

Using Theorem 2 and substituting the expression of ${\rm tr}(H^{2}_{m}\rho)$
derived above yields

$\begin{array}[]{ll}\displaystyle{\mathbb{E}_{{\bf b}\sim\mathcal{P}_{m}({\bf
b};\rho)}\big{[}(\texttt{CW}({\bf
b}))^{2}\big{]}}&\displaystyle{=\frac{1}{m^{4}}{\rm
tr}(H^{2}_{m}\rho)+\frac{m^{2}-1}{m^{4}}W{\rm
tr}(W\rho)+\frac{(m^{2}-1)^{2}}{4m^{2}}W^{2}}.\end{array}$

Combining with Eq. (22), the variance of the cut weight can be written as

$\begin{array}[]{ll}\displaystyle{{\rm Var}_{{\bf b}\sim\mathcal{P}_{m}({\bf
b};\rho)}\big{[}\texttt{CW}({\bf b})\big{]}}&\displaystyle{=\mathbb{E}_{{\bf
b}\sim\mathcal{P}_{m}({\bf b};\rho)}\big{[}\texttt{CW}({\bf
b})^{2}\big{]}-\mathbb{E}_{{\bf b}\sim\mathcal{P}_{m}({\bf
b};\rho)}\big{[}\texttt{CW}({\bf b})\big{]}^{2}}\\\\[14.22636pt]
&\displaystyle{=\frac{1}{m^{4}}{\rm
tr}(H_{m}^{2}\rho)+\frac{m^{2}-1}{m^{4}}W{\rm
tr}(H_{m}\rho)+\frac{(m^{2}-1)^{2}}{4m^{4}}W^{2}}\\\\[5.69054pt]
&\displaystyle{{\hskip 28.45274pt}-\Big{(}\frac{1}{m^{2}}{\rm
tr}(H_{m}\rho)+\frac{m^{2}-1}{2m^{2}}W\Big{)}^{2}}\\\\[14.22636pt]
&\displaystyle{=\frac{1}{m^{4}}\Big{(}{\rm tr}(H_{m}^{2}\rho)-{\rm
tr}(H_{m}\rho)^{2}\Big{)}}.\end{array}$

$\Box$

  * Proof of Corollary 3.

According to Eq. (20),

$\begin{array}[]{ll}\displaystyle{\mathcal{E}_{j}}&\displaystyle{:={\rm
tr}(P_{\langle j\rangle}\rho)}\\\\[5.69054pt]
&\displaystyle{=\sqrt{m}\,\mathbb{E}_{{\bf b}\sim\mathcal{P}_{\rm f}({\bf
b};\rho)}[(-1)^{b_{j}}]}\\\\[5.69054pt]
&\displaystyle{=\sqrt{m}\,\big{[}\mathcal{P}(b_{j}=0)\cdot(+1)+\mathcal{P}(b_{j}=1)\cdot(-1)\big{]}}\\\\[5.69054pt]
&\displaystyle{=\sqrt{m}\,\big{[}\mathcal{P}(b_{j}=0)-(1-\mathcal{P}(b_{j}=0))\big{]}}\\\\[5.69054pt]
&\displaystyle{=2\sqrt{m}\mathcal{P}(b_{j}=0)-\sqrt{m}}.\end{array}$

Then,

$\displaystyle{\mathcal{P}(b_{j}=0)=\frac{1}{2}+\frac{1}{2\sqrt{m}}\mathcal{E}_{j}}$
$\displaystyle{\mathcal{P}(b_{j}=1)=1-\mathcal{P}(b_{j}=0)=\frac{1}{2}-\frac{1}{2\sqrt{m}}\mathcal{E}_{j}.}$

Similarly, when $q^{(j)}\neq q^{(k)}$,

$\begin{array}[]{ll}\displaystyle{\mathcal{E}_{jk}}&\displaystyle{:={\rm
tr}(P_{\langle j\rangle}P_{\langle k\rangle}\rho)}\\\\[5.69054pt]
&\displaystyle{=m\,\mathbb{E}_{{\bf b}\sim\mathcal{P}_{\rm f}({\bf
b};\rho)}[(-1)^{b_{j}+b_{k}}]}\\\\[5.69054pt]
&\displaystyle{=m\,\big{[}\mathcal{P}(b_{j}=b_{k})\cdot(+1)+\mathcal{P}(b_{j}\neq
b_{k})\cdot(-1)\big{]}}\\\\[5.69054pt]
&\displaystyle{=m\,\big{[}\mathcal{P}(b_{j}=b_{k})-(1-\mathcal{P}(b_{j}=b_{k}))\big{]}}\\\\[5.69054pt]
&\displaystyle{=2m\mathcal{P}(b_{j}=b_{k})-m}.\end{array}$

Then,

$\displaystyle{\mathcal{P}(b_{j}=b_{k})=\frac{1}{2}+\frac{1}{2m}\mathcal{E}_{jk}}$
$\displaystyle{\mathcal{P}(b_{j}\neq
b_{k})=1-\mathcal{P}(b_{j}=b_{k})=\frac{1}{2}-\frac{1}{2m}\mathcal{E}_{jk}.}$

However, when $q^{(j)}=q^{(k)}$

$\begin{array}[]{ll}\displaystyle{\mathcal{P}(b_{j}=b_{k})}&\displaystyle{=\mathcal{P}(b_{j}=0)\cdot\mathcal{P}(b_{k}=0)+\mathcal{P}(b_{j}=1)\cdot\mathcal{P}(b_{k}=1)}\\\\[5.69054pt]
&\displaystyle{=\frac{1}{2}+\frac{1}{2m}\mathcal{E}_{j}\mathcal{E}_{k}}\end{array}$

and

$\begin{array}[]{ll}\displaystyle{\mathcal{P}(b_{j}\neq
b_{k})}&\displaystyle{=1-\mathcal{P}(b_{j}=b_{k})}\\\\[5.69054pt]
&\displaystyle{=\frac{1}{2}-\frac{1}{2m}\mathcal{E}_{j}\mathcal{E}_{k}}.\end{array}$

$\Box$

##  Appendix D Matrix Product State

Let $\Psi\in\mathbb{C}^{2^{n}}$ be $n$-qubit quantum state. Generally, $\Psi$
can be represented as the matrix multiplication of $n$ matrices as

$\Psi=A^{[n-1]}A^{[n-2]}\cdots A^{[1]}A^{[0]},$ (S.16)

where

$A^{[i]}\in\mathbb{C}^{\min(2^{i+1},2^{n-i-1},\chi)\times
2\times\min(2^{i},2^{n-i-1},\chi)}$ (S.17)

and

$(A^{[i]}A^{[j]})_{k(ab)m}=\sum^{\min(2^{i},2^{n-i-1},\chi)}_{l=1}(A^{[i]})_{kal}(A^{[j]})_{lbm}.$
(S.18)

The subscript in the parentheses, $(ab)$, summarizes two indices in one index,
i.e., $(a,b)=(1,1),(1,2),\ldots,(2,1),(2,2),\ldots$. The expression of Eq.
(S.16) is called the MPS [60]. The dimension between matrices is called the
bond dimension and can be restricted to $\chi$ at the expense of
representation power. The number of components of an $n$-qubit state is
$2^{n}$, whereas the total number of matrix components of MPS is $\sim
2n\chi^{2}$. For example, a $4$-qubit state
$|\psi\rangle\in\mathbb{C}^{2^{4}}$ with bond dimension $\chi=2$ can be
represented as

$\displaystyle{\psi_{(abcd)}=\sum_{\alpha,\beta,\gamma,\delta,\epsilon}\begin{bmatrix}\begin{bmatrix}\bullet\\\
\bullet\end{bmatrix}_{a}&\begin{bmatrix}\bullet\\\
\bullet\end{bmatrix}_{a}\end{bmatrix}_{\alpha\beta}\begin{bmatrix}\begin{bmatrix}\bullet\\\
\bullet\end{bmatrix}_{b}&\begin{bmatrix}\bullet\\\
\bullet\end{bmatrix}_{b}&\begin{bmatrix}\bullet\\\
\bullet\end{bmatrix}_{b}&\begin{bmatrix}\bullet\\\
\bullet\end{bmatrix}_{b}\\\\[14.22636pt] \begin{bmatrix}\bullet\\\
\bullet\end{bmatrix}_{b}&\begin{bmatrix}\bullet\\\
\bullet\end{bmatrix}_{b}&\begin{bmatrix}\bullet\\\
\bullet\end{bmatrix}_{b}&\begin{bmatrix}\bullet\\\
\bullet\end{bmatrix}_{b}\end{bmatrix}_{\beta\gamma}\begin{bmatrix}\begin{bmatrix}\bullet\\\
\bullet\end{bmatrix}_{c}&\begin{bmatrix}\bullet\\\
\bullet\end{bmatrix}_{c}\\\\[14.22636pt] \begin{bmatrix}\bullet\\\
\bullet\end{bmatrix}_{c}&\begin{bmatrix}\bullet\\\
\bullet\end{bmatrix}_{c}\\\\[14.22636pt] \begin{bmatrix}\bullet\\\
\bullet\end{bmatrix}_{c}&\begin{bmatrix}\bullet\\\
\bullet\end{bmatrix}_{c}\\\\[14.22636pt] \begin{bmatrix}\bullet\\\
\bullet\end{bmatrix}_{c}&\begin{bmatrix}\bullet\\\
\bullet\end{bmatrix}_{c}\end{bmatrix}_{\gamma\delta}\begin{bmatrix}\begin{bmatrix}\bullet\\\
\bullet\end{bmatrix}_{d}\\\\[14.22636pt] \begin{bmatrix}\bullet\\\
\bullet\end{bmatrix}_{d}\end{bmatrix}_{\delta\epsilon}}$
$\approx\sum_{\alpha,\beta,\gamma,\delta,\epsilon}\begin{bmatrix}\begin{bmatrix}\bullet\\\
\bullet\end{bmatrix}_{a}&\begin{bmatrix}\bullet\\\
\bullet\end{bmatrix}_{a}\end{bmatrix}_{\alpha\beta}\begin{bmatrix}\begin{bmatrix}\bullet\\\
\bullet\end{bmatrix}_{b}&\begin{bmatrix}\bullet\\\
\bullet\end{bmatrix}_{b}\\\\[14.22636pt] \begin{bmatrix}\bullet\\\
\bullet\end{bmatrix}_{b}&\begin{bmatrix}\bullet\\\
\bullet\end{bmatrix}_{b}\end{bmatrix}_{\beta\gamma}\begin{bmatrix}\begin{bmatrix}\bullet\\\
\bullet\end{bmatrix}_{c}&\begin{bmatrix}\bullet\\\
\bullet\end{bmatrix}_{c}\\\\[14.22636pt] \begin{bmatrix}\bullet\\\
\bullet\end{bmatrix}_{c}&\begin{bmatrix}\bullet\\\
\bullet\end{bmatrix}_{c}\end{bmatrix}_{\gamma\delta}\begin{bmatrix}\begin{bmatrix}\bullet\\\
\bullet\end{bmatrix}_{d}\\\\[14.22636pt] \begin{bmatrix}\bullet\\\
\bullet\end{bmatrix}_{d}\end{bmatrix}_{\delta\epsilon},$ (S.19)

where $\bullet$ indicates the matrix element and the subscripts represent the
index symbols of the inner submatrix elements.

Let $P$ be the identity or Pauli matrix and

$H=\sum^{K}_{k=1}C_{k}\bigotimes^{n-1}_{l=0}P_{kl}$ (S.20)

be a $n$-qubit Hamiltonian that consists of a weighted sum of $K$ Pauli
products. Similar to the MPS, Eq. (S.20) can be represented as matrix
multiplication of $n$ matrices as

$H=B^{[n-1]}B^{[n-2]}\cdots B^{[1]}B^{[0]},$ (S.21)

where

$B^{[i]}\in\left\\{\begin{array}[]{ll}\displaystyle{\mathbb{C}^{K\times
2\times 2\times 1}}&i=0\\\\[14.22636pt] \displaystyle{\mathbb{C}^{1\times
2\times 2\times K}}&i=n-1\\\\[14.22636pt] \displaystyle{\mathbb{C}^{K\times
2\times 2\times K}}&{\rm else}\end{array}\right.,$ (S.22)
$\left\\{\begin{array}[]{ll}\displaystyle{(B^{[i]})_{kaa^{\prime}1}=C_{k}(P_{ki})_{aa^{\prime}}}&i=0\\\\[14.22636pt]
\displaystyle{(B^{[i]})_{1aa^{\prime}k}=(P_{ki})_{aa^{\prime}}}&i=n-1\\\\[14.22636pt]
\displaystyle{(B^{[i]})_{kaa^{\prime}k}=(P_{ki})_{aa^{\prime}}}&{\rm
else}\end{array}\right.$ (S.23)

and

$(B^{[i]}B^{[j]})_{l(ab)(a^{\prime}b^{\prime})m}=\sum^{K}_{k=1}(B^{[i]})_{laa^{\prime}k}(B^{[j]})_{kbb^{\prime}m}.$
(S.24)

For example,

$\displaystyle{H_{(abcd)(a^{\prime}b^{\prime}c^{\prime}d^{\prime})}=(C_{1}IXYI+C_{2}IZII+C_{3}IYIX)_{(abcd)(a^{\prime}b^{\prime}c^{\prime}d^{\prime})}}$
$=\sum_{\alpha,\beta,\gamma,\delta,\epsilon}\begin{bmatrix}[I]_{aa^{\prime}}&[I]_{aa^{\prime}}&[I]_{aa^{\prime}}\end{bmatrix}_{\alpha\beta}\begin{bmatrix}[X]_{bb^{\prime}}\\\
&[Z]_{bb^{\prime}}\\\
&&[Y]_{bb^{\prime}}\end{bmatrix}_{\beta\gamma}\begin{bmatrix}[Y]_{cc^{\prime}}\\\
&[I]_{cc^{\prime}}\\\
&&[I]_{cc^{\prime}}\end{bmatrix}_{\gamma\delta}\begin{bmatrix}[C_{1}I]_{dd^{\prime}}\\\
[C_{2}I]_{dd^{\prime}}\\\
[C_{3}X]_{dd^{\prime}}\end{bmatrix}_{\delta\epsilon}$ (S.25)

The expression of Eq. (S.21) is called the matrix product operator (MPO) [60].

The MPO is a set of sparse matrices and thus can be efficiently stored in
classical memory in COOrdinate (COO), Compressed Sparse Row (CSR), or
Compressed Sparse Column (CSC) formats.

Using the MPS and MPO, the expectation value $\langle\Psi|H|\Psi\rangle$ can
be efficiently calculated by the matrix multiplication of $\\{A^{[i]}\\}$ and
$\\{B^{[i]}\\}$.

All edge energy can be simultaneously evaluated by modifying the MPO of the
problem Hamiltonian.

$\left\\{\begin{array}[]{ll}\displaystyle{(B^{[i]})_{kaa^{\prime}1}=(P_{ki})_{aa^{\prime}}}&i=0\\\\[14.22636pt]
\displaystyle{(B^{[i]})_{kaa^{\prime}k}=(P_{ki})_{aa^{\prime}}}&{\rm
else}\end{array}\right.$ (S.26)

For example,

$\displaystyle{\begin{bmatrix}C_{1}IXYI_{(abcd)(a^{\prime}b^{\prime}c^{\prime}d^{\prime})}\\\
C_{2}IZII_{(abcd)(a^{\prime}b^{\prime}c^{\prime}d^{\prime})}\\\
C_{3}IYIX_{(abcd)(a^{\prime}b^{\prime}c^{\prime}d^{\prime})}\end{bmatrix}}$
$=\sum_{\alpha,\beta,\gamma,\delta,\epsilon}\begin{bmatrix}[I]_{aa^{\prime}}\\\
&[I]_{aa^{\prime}}\\\
&&[I]_{aa^{\prime}}\end{bmatrix}_{\alpha\beta}\begin{bmatrix}[X]_{bb^{\prime}}\\\
&[Z]_{bb^{\prime}}\\\
&&[Y]_{bb^{\prime}}\end{bmatrix}_{\beta\gamma}\begin{bmatrix}[Y]_{cc}\\\
&[I]_{cc^{\prime}}\\\
&&[I]_{cc^{\prime}}\end{bmatrix}_{\gamma\delta}\begin{bmatrix}[C_{1}I]_{dd^{\prime}}\\\
[C_{2}I]_{dd^{\prime}}\\\
[C_{3}X]_{dd^{\prime}}.\end{bmatrix}_{\delta\epsilon}$ (S.27)

##  Appendix E Implementation Details for GW, CirCut, the RQAOA, and QRAO

For GW, an official MATLAB implementation of SDPT3 [61, 62] is used to solve
the SDP. The gap tolerance is set to $10^{-9}$ for high-precision calculation.
For hyperplane cutting, which is a post-processing step in GW performed after
computing the SDP solution, in-house Python code is used. The number of
hyperplane cuttings is set to $10,000$.

For CirCut [16], the official Fortran implementation (https://www.cmor-
faculty.rice.edu/~zhang/circut/index.html) is used. The hyperparameters for
CirCut are not changed from the default ones.

For RQAOA, in-house Fortran code is used. The analytic form of
$\langle\psi_{\rm QAOA}({\boldsymbol{\theta}})|Z_{j}Z_{k}|\psi_{\rm
QAOA}({\boldsymbol{\theta}})\rangle$ [3] is used to evaluate the expectation
value of the problem Hamiltonian and edge energy, where $|\psi_{\rm
QAOA}({\boldsymbol{\theta}})\rangle=e^{i\beta B}e^{i\gamma
C}|+\rangle^{\otimes n}$ is a level-1 QAOA ansatz, $B=\sum^{n-1}_{i=0}X_{i}$,
$C=\sum_{(j,k)\in E}w_{jk}Z_{j}Z_{k}$, and
${\boldsymbol{\theta}}=\\{\beta,\gamma\\}$. Only the optimal $\gamma$ is
determined using a grid search; the optimal $\beta$ is calculated for each
$\gamma$ using the same method used in [14]. We derived the concrete
representation of the optimal $\beta$ as

$\beta^{*}=-\frac{1}{4}{\rm
arctan}\left(\frac{\frac{1}{2}(R_{11}-R_{22}-R_{33}+R_{44})+R_{14}-R_{23}}{-I_{12}-I_{13}+I_{24}+I_{34}}\right)+\frac{\pi}{8},$
(S.28)

where

$\begin{bmatrix}R_{11}&R_{12}+iI_{12}&R_{13}+iI_{13}&R_{14}+iI_{14}\\\
R_{12}-iI_{12}&R_{22}&R_{23}+iI_{23}&R_{24}+iI_{24}\\\
R_{13}-iI_{13}&R_{23}-iI_{23}&R_{33}&R_{34}+iI_{34}\\\
R_{14}-iI_{14}&R_{24}-iI_{24}&R_{34}-iI_{34}&R_{44}\end{bmatrix}:=\sum_{(u,v)\in
E}w_{uv}\rho_{uv}$ (S.29)

and

$\rho_{uv}:={\rm tr}_{uv}\Big{(}e^{i\gamma C}|+\rangle^{\otimes
n}\langle+|^{\otimes n}e^{-i\gamma C}\Big{)}.$ (S.30)

Here, ${\rm tr}_{uv}(\sigma)$ denotes the reduced state obtained by tracing
out, excepting the $u$th and $v$th qubits from $\sigma$. For an algorithm to
calculate $\rho_{uv}$, see [14]. The parameter spaces are set to
$\beta\in[0,\frac{\pi}{2}]$ and $\gamma\in[0,\pi]$. Once the optimal $\gamma$
is found in $50$ equidistant grids of $[0,\pi]$, the more precise solution is
searched for in the fine-grained grids around the course-grained solution with
a mesh size is $1/50$. Fine-graining is repeated twice, resulting in a mesh
size of $1/250$ of the parameter space.

For QRAO, in-house Python code is used. The MPS ansatz with bond dimension
$\chi=2$ is used for $|\psi({\boldsymbol{\theta}})\rangle$ and Pauli rounding
is used for decoding. The same optimizer as used for RQRAO, L-BFGS with the
strong Wolfe condition, is used to optimize
$|\psi({\boldsymbol{\theta}})\rangle$.

##  Appendix F Additional Sensitivity Analysis Results (ExpHyp)

Figure S.3: Results of the ablation study. Only one hyperparameter is varied
for each trial. The problem graph instances are $3$-regular graphs with $\pm
1$ edge weights. For each number of nodes, $10$ graph instances are randomly
generated. The varying hyperparameters are the (a) number of ensembles $N$;
(b) scale factor $S$; and (c) bond dimension $\chi$.

Here, the sensitivities of all RQRAO hyperparameters except for the number of
embeddings per qubit $m$ on the performance of the cut weight and runtime are
reported. The results are summarized in Fig. S.3. In brief, a higher cut
weight is obtained through a longer runtime.

Except for the number of embeddings per qubit $m$, the sensitive
hyperparameters with respect to cut weight performance are the number of
ensembles $N$ and scale factor $S$. In Fig. S.3, the cut weight increases as
the number of ensembles $N$ increases. In addition, as the number of ensembles
increases, the number of edges whose parities are simultaneously determined in
one optimization step tends to decrease. This is because the edges, which are
a candidate for determining the parity when the number of ensembles is small,
drop out as the accuracy of the parity confidence improves, which is achieved
by increasing the number of ensembles. This trend (i.e., the number of
candidates decreases with as the accuracy of the parity confidence increases),
can be seen in the larger scale factor $S$ in Fig. S.3. A larger $S$ tends to
decrease the number of candidates because the number of edges with nonzero
$\mathcal{E}_{jk}$ decreases.

The effect of the bond dimension $\chi$ on the cut weight performance needs to
be carefully considered. Recall that a large entanglement can be treated using
a large bond dimension. The $\chi=1$ result is worse both in terms of relative
cut weight and runtime than the $\chi>1$ results. In the case of $\chi=1$,
$|\psi({\boldsymbol{\theta}})\rangle$ is restricted to a direct product state
that tends to be trapped in a lower energy state than in the $\chi>1$ case
because the search space is much smaller. As a result, each edge energy in an
ensemble, $\mathcal{E}^{(t)}_{jk}$, is no longer biasedly distributed in
either positive or negative sides because the lower energy state includes
smaller information of a large cut weight. By contrast, when $\chi>1$, no
significant difference in the relative cut weight can be seen. This does not
mean, however, that a large entanglement will not affect the cut weight
performance because even the value of $\chi=10$, which was the maximum value
considered in the experiment, is small compared with the maximum bond
dimension that the system can take, which is $\chi_{\rm max}=2^{\lfloor
n/2\rfloor}$. For example, consider the case of $|V|=300$. The number of
qubits is assumed to be $300/3=100$, which can be achieved by a $(3,1)$-QRAC
formulation if the graph is sufficiently sparse. In this case, $\chi_{\rm
max}=2^{\lfloor 99/2\rfloor}\approx 10^{14}$, which is much larger than the
experiments considered in this study. To take into account a large bond
dimension within a reasonable runtime, computation using a real quantum device
would be required.