combinatorial optimization.

Finally, we remark that Grothendieck-type problems and inequalities have a
considerable historical connection to quantum theory. Tsirelson
[tsirelson1987quantum] employed Clifford algebras to reformulate the
commutative Grothendieck inequality into a statement about classical XOR games
with entanglement. Regev and Vidick [regev2015quantum] later introduced the
notion of quantum XOR games, which they studied through the generalization of
such ideas to noncommutative Grothendieck inequalities. The mathematical work
of Haagerup and Itoh [haagerup1995grothendieck] studied Grothendieck-type
inequalities as the norms of operators on $C^{*}$-algebras; their analysis
makes prominent use of canonical anticommutation relation algebras over
fermionic Fock spaces. Quadratic programming with orthogonality constraints
has also been applied for classical approximation algorithms for quantum many-
body problems, for instance by Bravyi _et al._ [bravyi2019approximation].
Recasting noncommutative Grothendieck problems into a quantum Hamiltonian
problem may therefore provide new insights into these connections.

The rest of this chapter is organized as follows: Section Problem statement
provides a formal description of the optimization problem that we study in
this chapter and reviews known complexity results of related problems. In
Section Applications of the LNCG problem we describe two well-known
applications of the problem: the group synchronization problem and the
generalized orthogonal Procrustes problem. Section Summary of results provides
a summary of our quantum relaxation which embeds the optimization problem into
a Hamiltonian, and two accompanying rounding protocols. In Section Quantum
formalism for optimization over orthogonal matrices we derive an embedding of
orthogonal matrices into quantum states via the Pin and Spin groups. We
elaborate on the connection to fermionic theories and provide a quantum
perspective on the convex hull of the orthogonal groups. From this embedding,
Section Quantum relaxation for the quadratic problem then establishes the
quantum Hamiltonian relaxation of the quadratic optimization problem. Section
Rounding algorithms describes both classical and quantum rounding protocols
for relaxations of the problem. Notably, for the classical SDP we derive an
approximation ratio for $\mathrm{SO}(n)$ rounding. Finally, in Section
Numerical experiments we demonstrate numerical experiments on random instances
of the group synchronization problem for $\mathrm{SO}(3)$ on three-regular
graphs and report the performance of various classical and quantum rounding
protocols. For our simulations of the quantum relaxation, we consider two
classes of quantum states: maximal eigenstates of the Hamiltonian and quasi-
adiabatically evolved states. We close in Section Discussion and future work
with a discussion on future lines of research.

### Problem statement

In this paper we consider the class of little noncommutative Grothendieck
(LNCG) problems over the orthogonal group, as studied previously by Bandeira
_et al._ [bandeira2016approximating].454545The authors also consider the
complex-valued problem over the unitary group, which is outside the scope of
this dissertation. Let $(V,E)$ be an undirected graph with $m=|V|$ vertices
and edge set $E$. For integer $n\geq 1$, let $C\in\mathbb{R}^{mn\times mn}$ be
a symmetric matrix, which for notation we partition into $n\times n$ blocks as

$C=\begin{bmatrix}C_{11}&\cdots&C_{1m}\\\ \vdots&\ddots&\vdots\\\
C_{m1}&\cdots&C_{mm}\end{bmatrix}.$ (780)

The quadratic program we wish to solve is of the form

$\max_{R_{1},\ldots,R_{m}\in G}\sum_{(u,v)\in E}\langle
C_{uv},R_{u}R_{v}^{\mathsf{T}}\rangle,$ (781)

where $G$ is either the orthogonal group

$\mathrm{O}(n)\coloneqq\\{R\in\mathbb{R}^{n\times n}\mid
R^{\mathsf{T}}R=\mathbb{I}_{n}\\},$ (782)

or the special orthogonal group

$\mathrm{SO}(n)\coloneqq\\{R\in\mathrm{O}(n)\mid\det R=1\\}$ (783)

on $\mathbb{R}^{n}$. Here, $\langle
A,B\rangle=\operatorname{tr}(A^{\mathsf{T}}B)$ denotes the Frobenius inner
product on the space of real matrices and $\mathbb{I}_{n}$ is the $n\times n$
identity matrix. Note that when $G=\mathrm{O}(1)=\\{\pm 1\\}$, Problem (781)
reduces to combinatorial optimization of the form

$\max_{x_{1},\ldots,x_{m}\in\\{\pm 1\\}}\sum_{(u,v)\in E}C_{uv}x_{u}x_{v},$
(784)

where now $C\in\mathbb{R}^{m\times m}$. This is sometimes referred to as the
commutative instance of the little Grothendieck problem. Problem (781) can
therefore be viewed as a natural generalization of quadratic binary
optimization to the noncommutative matrix setting.

We now comment on the known hardness results of these optimization problems.
The commutative problem (784) is already $\mathsf{NP}$-hard in general, as can
be seen by the fact that the Max-Cut problem can be expressed in this form. In
particular, Khot _et al._ [khot2007optimal] proved that, assuming the Unique
Games conjecture, it is $\mathsf{NP}$-hard to approximate the optimal Max-Cut
solution to better than a fraction of
$(2/\pi)\min_{\theta\in[0,\pi]}\frac{\theta}{1-\cos\theta}\approx 0.878$. This
value coincides with the approximation ratio achieved by the celebrated
Goemans–Williamson (GW) algorithm for rounding the semidefinite relaxation of
the problem [goemans1995improved]. More generally, consider the fully
connected graph and let $C\succeq 0$ be arbitrary. Nesterov
[nesterov1998semidefinite] showed that GW rounding guarantees an approximation
ratio of $2/\pi\approx 0.636$ in this setting, which Alon and Naor
[alon2004approximating] showed matches the integrality gap of the semidefinite
program. Khot and Naor [khot2009approximate] later demonstrated that this
approximation ratio is also Unique-Games-hard to exceed, and finally Briët _et
al._ [briet2017tight] strengthened this result to be unconditionally
$\mathsf{NP}$-hard.

For the noncommutative problem (781) that we are interested in, less is known
about its hardness of approximability. However, it is a subclass of more
general optimization problems for which some results are known. The most
general instance is the “big” noncommutative Grothendieck problem, for which
Naor _et al._ [naor2014efficient] provided a rounding procedure of its
semidefinite relaxation. Their algorithm achieves an approximation ratio of at
least $1/2\sqrt{2}\approx 0.353$ in the real-valued setting, and $1/2$ in the
complex-valued setting (wherein optimization is over the unitary group instead
of the orthogonal group). This $1/2$ result was later shown to be tight by
Briët _et al._ [briet2017tight] for both the real- and complex-valued
settings; in fact, they show that this is the $\mathsf{NP}$-hardness threshold
of a special case of the problem, called the _little_ noncommutative
Grothendieck problem.464646See Section 6 of Briët _et al._ [briet2017tight]
for the precise relation between the big and little NCG. However, the
threshold for Problem (781), which is an special case of LNCG, is not known.
Algorithmically, Bandeira _et al._ [bandeira2016approximating] demonstrated
constant approximation ratios for Problem (781) when $C\succeq 0$ and
$G=\mathrm{O}(n)$ or $\mathrm{U}(n)$ via an $(n\times n)$-dimensional
generalization of GW rounding, along with matching integrality gaps. These
approximation ratios exceed $1/2$, indicating that this subclass is
quantitatively less difficult than the general instance of the LNCG problem.
Although the optimization of rotation matrices is of central importance to
many applications, we are unaware of any general approximation ratio
guarantees for the $G=\mathrm{SO}(n)$ setting.

### Applications of the LNCG problem

Before describing our quantum relaxation, here we motivate the practical
interest in Problem (781) by briefly discussing some applications. Throughout,
let $G=\mathrm{O}(n)$ or $\mathrm{SO}(n)$ and $(V,E)$ be a graph as before.

#### Group synchronization

The group synchronization problem over orthogonal transformations has
applications in a variety of disciplines, including structural biology,
robotics, and wireless networking. For example, in structural biology the
problem appears as part of the cryogenic electron microscopy (cryo-EM)
technique. There, one uses electron microscopy on cryogenically frozen samples
of a molecular structure to obtain a collection of noisy images of the
structure. The images are noisy due to an inherently low signal-to-noise
ratio, and furthermore they feature the structure in different, unknown
orientations (represented by rotation matrices). One approach to solving the
group synchronization problem yields best-fit estimates for these orientations
via least-squares minimization [boumal2013robust], from which one can produce
a model of the desired 3D structure.474747Note that other loss functions are
also considered in the literature, which may not necessarily have a
reformulation as Problem (781). See Ref. [singer2018mathematics] for a further
overview, and Ref. [ozyecsil2017survey] for a survey of other applications of
group synchronization.

The formal problem description is as follows. To each vertex
$v\in[m]\coloneqq\\{1,\ldots,m\\}$ we assign an unknown but fixed element
$g_{v}\in G$. An interaction between each pair of vertices connected by an
edge $(u,v)\in E$ is modeled as $g_{uv}=g_{u}g_{v}^{\mathsf{T}}$. However,
measurements of the interactions are typically corrupted by some form of
noise. For instance, one may consider an additive noise model of the form
$C_{uv}=g_{uv}+\sigma W_{uv}$, where $\sigma\geq 0$ characterizes the strength
of the noise and each $W_{uv}\in\mathbb{R}^{n\times n}$ has independently,
normally distributed entries. We would like to recover each $g_{v}$ given only
access to the matrices $C_{uv}$. Therefore, as a proxy to the recovery problem
one may cast the solution as the least-squares minimizer

$\min_{\mathbf{R}\in G^{m}}\sum_{(u,v)\in
E}\|C_{uv}-R_{u}R_{v}^{\mathsf{T}}\|_{F}^{2},$ (785)

where $\|A\|_{F}=\sqrt{\langle A,A\rangle}$ is the Frobenius norm and we
employ the notation $\mathbf{R}\equiv(R_{1},\ldots,R_{m})$. It is
straightforward to see that the minimzer of this problem is equivalent to the
maximizer of

$\max_{\mathbf{R}\in G^{m}}\sum_{(u,v)\in E}\langle
C_{uv},R_{u}R_{v}^{\mathsf{T}}\rangle,$ (786)

which is precisely in the form of Problem (781).

#### Generalized orthogonal Procrustes problem

Procrustes analysis has applications in fields such as shape and image
recognition, as well as sensory analysis and market research on
$n$-dimensional data. In this problem, one has a collection of point clouds,
each representing for instance the important features of an image. One wishes
to determine how similar these images are to each other collectively. This is
achieved by simultaneously fitting each pair of point clouds to each other,
allowing for arbitrary orthogonal transformations on each cloud to best align
the individual points. We refer the reader to Ref. [gower2004procrustes] for a
comprehensive review.

Consider $m$ sets of $K$ points in $\mathbb{R}^{n}$,
$S_{v}=\\{x_{v,1},\ldots,x_{v,K}\\}\subset\mathbb{R}^{n}$ for each $v\in[m]$.
We wish to find an orthogonal transformation $R_{v}\in G$ for each $S_{v}$
that best aligns all sets of points simultaneously. That is, for each
$k\in[K]$ and $u,v\in[m]$ we wish to minimize the Euclidean distance
$\|R_{u}^{\mathsf{T}}x_{u,k}-R_{v}^{\mathsf{T}}x_{v,k}\|_{2}$. Taking least-
squares minimization as our objective, we seek to solve

$\min_{\mathbf{R}\in
G^{m}}\sum_{u,v\in[m]}\sum_{k\in[K]}\|R_{u}^{\mathsf{T}}x_{u,k}-R_{v}^{\mathsf{T}}x_{v,k}\|_{2}^{2}.$
(787)

From the relation between the vector 2-norm and matrix Frobenius norm, Eq.
(787) can be formulated as

$\max_{\mathbf{R}\in G^{m}}\sum_{u,v\in[m]}\langle
C_{uv},R_{u}R_{v}^{\mathsf{T}}\rangle,$ (788)

where each $C_{uv}\in\mathbb{R}^{n\times n}$ is defined as

$C_{uv}=\sum_{k\in[K]}x_{u,k}x_{v,k}^{\mathsf{T}}.$ (789)

### Summary of results

We now provide a high-level overview of the main contributions of this work.
We provide summary cartoon in Figure 23, depicting the quantum embedding of
the problem and the quantum rounding protocols. Let $(V,E)$ be a graph where
we label the vertices by $V=[m]$, and denote the objective function of Problem
(781) by

$f(\mathbf{R})\coloneqq\sum_{(u,v)\in E}\langle
C_{uv},R_{u}R_{v}^{\mathsf{T}}\rangle.$ (790)

Figure 23: A cartoon description of the quantum and classical encodings of an
LNCG problem, followed by classical and quantum rounding. (Left) The
description of the problem that we consider, which is described by a graph
$([m],E)$ and $n\times n$ matrices $C_{uv}$ for each edge $(u,v)$. We wish to
assign elements of $\mathrm{O}(n)$ or $\mathrm{SO}(n)$ to each vertex such
that the quadratic form of Eq. (790) is maximized. (Center top) The
description of the standard classical relaxation of the LNCG problem as an
$mn\times mn$ PSD matrix $M\succeq 0$, which is optimized using a semidefinite
program. (Right top) The classical rounding procedure, which returns a
collection of orthogonal matrices from $M$. (Center bottom) A description of
our quantum formulation of the LNCG problem as a two-body interacting
Hamiltonian. On each vertex we place a $d$-dimensional Hilbert space, and the
Hamiltonian corresponds to interaction terms $H_{uv}$ on the edges
($\mathcal{L}(\mathcal{H})$ is the set of linear operators on a Hilbert space
$\mathcal{H}$). The classical solution of the LNCG problem lies in a subset of
the full Hilbert space containing separable Gaussian states. (Right bottom)
Our proposed quantum rounding protocols. One protocol requires knowledge of
the two-body reduced density matrices across edges, while the other uses the
one-body reduced density matrices on each vertex.

#### Quantum Hamiltonian relaxation

First, consider the setting in which
$\mathbf{R}=(R_{1},\ldots,R_{m})\in\mathrm{O}(n)^{m}$. We embed this problem
into a Hamiltonian by placing $n$ qubits on each vertex $v\in[m]$, resulting
in a total Hilbert space $\mathcal{H}_{2^{n}}^{\otimes m}$ of $mn$ qubits.
Define the $n$-qubit Pauli operators

$P_{ij}\coloneqq\begin{cases}-X_{i}Z_{i+1}\cdots Z_{j-1}X_{j}&i<j,\\\
Z_{i}&i=j,\\\ -Y_{j}Z_{j+1}\cdots Z_{i-1}Y_{i}&i>j,\end{cases}$ (791)

where $Z_{i}\coloneqq\mathbb{I}_{2}^{\otimes(i-1)}\otimes
Z\otimes\mathbb{I}_{2}^{\otimes(n-i)}$ (similarly for $X_{i}$, $Y_{i}$). The
Hamiltonian

$H\coloneqq\sum_{(u,v)\in
E}\sum_{i,j\in[n]}[C_{uv}]_{ij}\sum_{k\in[n]}P_{ik}^{(u)}\otimes P_{jk}^{(v)}$
(792)

defines our quantum relaxation of the objective $f$ over $\mathrm{O}(n)^{m}$.
The notation $A^{(v)}$ denotes the operator $A$ acting only on the Hilbert
space of vertex $v$, and we overload this notation to indicate either the
$n$-qubit operator or $mn$-qubit operator acting trivially on the remaining
vertices. When the context is clear we typically omit writing the trivial
support.

For optimization over $\mathrm{SO}(n)^{m}$, we consider instead the
$(n-1)$-qubit Pauli operators

$\displaystyle\widetilde{P}_{ij}$
$\displaystyle\coloneqq\Pi_{0}P_{ij}\Pi_{0}^{\mathsf{T}},$ (793)
$\displaystyle\Pi_{0}$
$\displaystyle=\frac{1}{\sqrt{2}}\mathopen{}\left(\langle+|\otimes\mathbb{I}_{2}^{\otimes(n-1)}+\langle-|\otimes
Z^{\otimes(n-1)}\right)\mathclose{},$ (794)

where $\Pi_{0}:\mathcal{H}_{2^{n}}\to\mathcal{H}_{2^{n-1}}$ represents the
projection onto the even-parity subspace of $\mathcal{H}_{2^{n}}$. The
construction of the relaxed Hamiltonian for $\mathrm{SO}(n)$ is then analogous
to Eq. (792):

$\widetilde{H}\coloneqq\sum_{u,v\in
E}\sum_{i,j\in[n]}[C_{uv}]_{ij}\sum_{k\in[n]}\widetilde{P}_{ik}^{(u)}\otimes\widetilde{P}_{jk}^{(v)},$
(795)

where now the relaxed quantum problem is defined over $m(n-1)$ qubits.

These Hamiltonians serve as relaxations to Problem (781) in the following
sense. First, we show that for every $R\in\mathrm{O}(n)$, there is an
$n$-qubit state $|\phi(R)\rangle$ which is the maximum eigenstate of

$F(R)=\sum_{i,j\in[n]}R_{ij}P_{ij}.$ (796)

In particular, $F(R)$ is a free-fermion Hamiltonian, so $|\phi(R)\rangle$ is a
fermionic Gaussian state. If $R\in\mathrm{SO}(n)$, then furthermore
$|\phi(R)\rangle$ is an even-parity state, i.e., $\langle\phi(R)|Z^{\otimes
n}|\phi(R)\rangle=1$, so it is only supported on a subspace of dimension
$2^{n-1}$ (the image of $\Pi_{0}$). This correspondence establishes a
reformulation of the classical optimization problem as a constrained
Hamiltonian problem:

$\max_{\mathbf{R}\in
G^{m}}f(\mathbf{R})=\max_{\begin{subarray}{c}|\psi\rangle=\bigotimes_{v\in[m]}|\phi(R_{v})\rangle\\\
R_{v}\in G\end{subarray}}\langle\psi|H|\psi\rangle.$ (797)

Dropping these constraints on $|\psi\rangle$ implies the inequalities

$\displaystyle\max_{\mathbf{R}\in\mathrm{O}(n)^{m}}f(\mathbf{R})$
$\displaystyle\leq\max_{\rho\in\mathcal{D}(\mathcal{H}_{2^{n}}^{\otimes
m})}\operatorname{tr}(H\rho),$ (798)
$\displaystyle\max_{\mathbf{R}\in\mathrm{SO}(n)^{m}}f(\mathbf{R})$
$\displaystyle\leq\max_{\rho\in\mathcal{D}(\mathcal{H}_{2^{n-1}}^{\otimes
m})}\operatorname{tr}(\widetilde{H}\rho),$ (799)

where $\mathcal{D}(\mathcal{H})$ denotes the set of density operators on a
Hilbert space $\mathcal{H}$. This establishes the quantum Hamiltonian
relaxation.

#### Quantum rounding

Data: Quantum state $\rho\in\mathcal{D}(\mathcal{H}_{d}^{\otimes m})$ over a
graph of $m$ vertices, each with local Hilbert space of dimension $d=2^{n}$ if
$G=\mathrm{O}(n)$, or $d=2^{n-1}$ if $G=\mathrm{SO}(n)$

Result: Orthogonal matrices on each vertex, $R_{1},\ldots,R_{m}\in G$

$\mathcal{M}\leftarrow\mathbb{I}_{mn}$ ;

for  _$u\neq v\in[m]$_ do

for  _$(i,j)\in[n]^{2}$_ do

if  _$G=\mathrm{O}(n)$_ then

$[\mathcal{M}_{uv}]_{ij}\leftarrow\frac{1}{n}\operatorname{tr}(\Gamma_{ij}^{(u,v)}\rho)$
;

else if  _$G=\mathrm{SO}(n)$_ then

$[\mathcal{M}_{uv}]_{ij}\leftarrow\frac{1}{n}\operatorname{tr}(\widetilde{\Gamma}_{ij}^{(u,v)}\rho)$
;

end if

end for

end for

for  _$v\in[m]$_ do

$R_{v}\leftarrow\operatorname*{arg\,min}_{Y\in G}\|Y-\mathcal{M}_{1v}\|_{F}$ ;

end for

Algorithm 1 $\operatorname{conv}G$-based rounding of edge marginals.

Data: Quantum state $\rho\in\mathcal{D}(\mathcal{H}_{d}^{\otimes m})$ over a
graph of $m$ vertices, each with local Hilbert space of dimension $d=2^{n}$ if
$G=\mathrm{O}(n)$, or $d=2^{n-1}$ if $G=\mathrm{SO}(n)$

Result: Orthogonal matrices on each vertex, $R_{1},\ldots,R_{m}\in G$

for  _$v\in[m]$_ do

${Q}_{v}\leftarrow 0\in\mathbb{R}^{n\times n}$ ;

for  _$(i,j)\in[n]^{2}$_ do

if  _$G=\mathrm{O}(n)$_ then

$[{Q}_{v}]_{ij}\leftarrow\operatorname{tr}(P_{ij}^{(v)}\rho)$ ;

else if  _$G=\mathrm{SO}(n)$_ then

$[{Q}_{v}]_{ij}\leftarrow\operatorname{tr}(\widetilde{P}_{ij}^{(v)}\rho)$ ;

end if

end for

end for

for  _$v\in[m]$_ do

$R_{v}\leftarrow\operatorname*{arg\,min}_{Y\in G}\|Y-{Q}_{v}\|_{F}$ ;

end for

Algorithm 2 Rounding vertex marginals.

In order to recover orthogonal matrices from a relaxed quantum solution
$\rho$, we propose two rounding procedures, summarized in Algorithms 1 and 2.
These rounding procedures operate on local (i.e., single- or two-vertex
observables) expectation values of $\rho$ stored in classical memory, which
can be efficiently estimated, e.g., by partial state tomography.

Algorithm 1 is inspired by constructing a quantum analogue of the PSD variable
appearing in semidefinite relaxations to Problem (781). Consider the $mn\times
mn$ matrix of expectation values

$\mathcal{M}\coloneqq\begin{bmatrix}\mathbb{I}_{n}&T_{12}&\cdots&T_{1m}\\\
T_{21}&\mathbb{I}_{n}&\cdots&T_{2m}\\\ \vdots&\vdots&\ddots&\vdots\\\
T_{m1}&T_{m2}&\cdots&\mathbb{I}_{n}\end{bmatrix},$ (800)

where the off-diagonal blocks are defined as

$\displaystyle T_{uv}$
$\displaystyle\coloneqq\frac{1}{n}\begin{bmatrix}\operatorname{tr}(\Gamma_{11}^{(u,v)}\rho)&\cdots&\operatorname{tr}(\Gamma_{1n}^{(u,v)}\rho)\\\
\vdots&\ddots&\vdots\\\
\operatorname{tr}(\Gamma_{n1}^{(u,v)}\rho)&\cdots&\operatorname{tr}(\Gamma_{nn}^{(u,v)}\rho)\end{bmatrix}=T_{vu}^{\mathsf{T}},$
(801) $\displaystyle\Gamma_{ij}^{(u,v)}$
$\displaystyle\coloneqq\sum_{k\in[n]}P_{ik}^{(u)}\otimes P_{jk}^{(v)}$ (802)

when $G=\mathrm{O}(n)$, and we replace the operators $P_{ij}$ with
$\widetilde{P}_{ij}$ when $G=\mathrm{SO}(n)$. We show that $\mathcal{M}$
satisfies the following properties for all states $\rho$:

$\displaystyle\mathcal{M}\succeq 0,$ (803)
$\displaystyle\mathcal{M}_{uv}\in\operatorname{conv}G\quad\forall u,v\in[m],$
(804)

where $\operatorname{conv}G$ is the convex hull of $G$. Thus when
$G=\mathrm{SO}(n)$, $\mathcal{M}$ obeys the same constraints as the
$\operatorname{conv}\mathrm{SO}(n)$-based semidefinite relaxation proposed by
Saunderson _et al._ [saunderson2014semidefinite]. However, whereas the
classical representation of the $\operatorname{conv}\mathrm{SO}(n)$
constraints requires at least matrices of size $2^{n-1}\times 2^{n-1}$ for
each edge, our quantum state automatically satisfies these constraints (using
only $n-1$ qubits per vertex).

Algorithm 2 uses the single-vertex information
$\operatorname{tr}(P_{ij}^{(v)}\rho)$ of $\rho$, as opposed to the two-vertex
information $\operatorname{tr}(\Gamma_{ij}^{(u,v)}\rho)$. We consider this
rounding procedure due to the fact that, if $\rho$ is a pure Gaussian state
satisfying the constraint of Eq. (797), then the matrix of expectation values

$Q_{v}\coloneqq\begin{bmatrix}\operatorname{tr}(P_{11}^{(v)}\rho)&\cdots&\operatorname{tr}(P_{1n}^{(v)}\rho)\\\
\vdots&\ddots&\vdots\\\
\operatorname{tr}(P_{n1}^{(v)}\rho)&\cdots&\operatorname{tr}(P_{nn}^{(v)}\rho)\end{bmatrix}$
(805)

lies in $\mathrm{O}(n)$. On the other hand, for arbitrary density matrices we
have the relaxation $Q_{v}\in\operatorname{conv}\mathrm{O}(n)$, and again when
we replace $P_{ij}^{(v)}$ with $\widetilde{P}_{ij}^{(v)}$ then
$Q_{v}\in\operatorname{conv}\mathrm{SO}(n)$.

Both rounding procedures use the standard projection of the matrices
$X\in\operatorname{conv}G$ (e.g., the matrices $T_{uv}$ or $Q_{v}$ measured
from the quantum state) to some $R\in G$ by finding the nearest (special)
orthogonal matrix according to Frobenius-norm distance:

$R=\operatorname*{arg\,min}_{Y\in G}\|X-Y\|_{F}.$ (806)

This can be solved efficiently as a classical postprocessing step, essentially
by computing the singular value decomposition of $X=U\Sigma V^{\mathsf{T}}$.
When $G=\mathrm{O}(n)$, the solution is $R=UV^{\mathsf{T}}$. When
$G=\mathrm{SO}(n)$, we instead use the so-called special singular value
decomposition of $X=U\widetilde{\Sigma}\widetilde{V}^{\mathsf{T}}$, where
$\widetilde{\Sigma}=\Sigma J$ and $\widetilde{V}=VJ$, with $J$ being the
diagonal matrix

$J=\begin{bmatrix}\mathbb{I}_{n-1}&0\\\ 0&\det(UV^{\mathsf{T}})\end{bmatrix},$
(807)

assuming that the singular values $\sigma_{i}(X)$ are in descending order,
$\sigma_{1}(X)\geq\cdots\geq\sigma_{n}(X)$. Then the solution to Eq. (806) is
$R=U\widetilde{V}^{\mathsf{T}}\in\mathrm{SO}(n)$.

### Quantum formalism for optimization over orthogonal matrices

Our key insight into encoding orthogonal matrices into quantum states comes
from the construction of the orthogonal group from a Clifford algebra. We
review this mathematical construction in Appendix Clifford algebras and the
orthogonal group and only discuss the main aspects here. The Clifford algebra
$\mathrm{Cl}(n)$ is a $2^{n}$-dimensional real vector space equipped with an
inner product and multiplication operation satisfying the anticommutation
relation

$e_{i}e_{j}+e_{j}e_{i}=-2\delta_{ij}\mathbb{I},$ (808)

where $e_{1},\ldots,e_{n}$ is an orthonormal basis for $\mathbb{R}^{n}$ and
$\mathbb{I}$ is the multiplicative identity of the algebra. The orthogonal
group is then realized through a quadratic map
$Q:\mathrm{Cl}(n)\to\mathbb{R}^{n\times n}$ and the identification of a
subgroup $\mathrm{Pin}(n)\subset\mathrm{Cl}(n)$ such that
$Q(\mathrm{Pin}(n))=\mathrm{O}(n)$. Notably, the elements of $\mathrm{Pin}(n)$
have unit norm (with respect to the inner product on $\mathrm{Cl}(n)$). The
special orthogonal group, meanwhile, is constructed by considering only the
even-parity elements of $\mathrm{Cl}(n)$, denoted by $\mathrm{Cl}^{0}(n)$. The
group $\mathrm{Spin}(n)=\mathrm{Pin}(n)\cap\mathrm{Cl}^{0}(n)$ then yields
$Q(\mathrm{Spin}(n))=\mathrm{SO}(n)$.

Because the Clifford algebra $\mathrm{Cl}(n)$ is a $2^{n}$-dimensional vector
space, we observe that it can be identified with a Hilbert space of $n$
qubits.484848In fact, $n$ rebits suffice since $\mathrm{Cl}(n)$ is a real
vector space, but to keep the presentation straightforward we will not make
such a distinction. In this section we explore this connection in detail,
showing how to represent orthogonal matrices as quantum states and how the
mapping $Q$ acts as a linear functional on those states.

#### Qubit representation of the Clifford algebra

First we describe the canonical isomorphism between $\mathrm{Cl}(n)$ and
$\mathcal{H}_{2^{n}}\coloneqq(\mathbb{R}^{2})^{\otimes n}$ as Hilbert spaces.
We denote the standard basis of $\mathrm{Cl}(n)$ by $\\{e_{I}\coloneqq
e_{i_{1}}\cdots e_{i_{k}}\mid I=\\{i_{1},\ldots,i_{k}\\}\subseteq[n]\\}$. By
convention we assume that the elements of $I$ are ordered as
$i_{1}<\cdots<i_{k}$. Each basis element $e_{I}$ maps onto to a computational
basis state $|b\rangle$, where $b=b_{1}\cdots b_{n}\in\\{0,1\\}^{n}$, via the
correspondence

$e_{I}\equiv\bigotimes_{i\in[n]}|b_{i}\rangle,\quad\text{where
}b_{i}=\begin{cases}1&\text{if }i\in I,\\\ 0&\text{otherwise}.\end{cases}$
(809)

The inner products on both spaces coincide since this associates one
orthonormal basis to another. This correspondence also naturally equates the
grade $|I|$ of the Clifford algebra with the Hamming weight $|b|$ of the
qubits. The notion of parity, $|I|\>\mathrm{mod}\>2=|b|\>\mathrm{mod}\>2$, is
therefore preserved, so $\mathrm{Cl}^{0}(n)$ corresponds to the subspace of
$\mathcal{H}_{2^{n}}$ with even Hamming weight.

To represent the multiplication of algebra elements in this Hilbert space, we
use the fact that left- and right-multiplication are linear automorphisms on
$\mathrm{Cl}(n)$, which are denoted by

$\lambda_{x}(y)=xy,\quad\rho_{x}(y)=yx.$ (810)

The action of the algebra can therefore be represented on
$\mathcal{H}_{2^{n}}$ as linear operators. We shall use the matrix
representation provided in Ref. [saunderson2015semidefinite], as it precisely
coincides with the $n$-qubit computational basis described above. Because of
linearity, it suffices to specify left- and right-multiplication by the
generators $e_{i}$, which are the operators

$\displaystyle\lambda_{i}$ $\displaystyle\equiv
Z^{\otimes(i-1)}\otimes(-\mathrm{i}Y)\otimes\mathbb{I}_{2}^{\otimes(n-i)},$
(811) $\displaystyle\rho_{i}$
$\displaystyle\equiv\mathbb{I}_{2}^{\otimes(i-1)}\otimes(-\mathrm{i}Y)\otimes
Z^{\otimes(n-i)}.$ (812)

It will also be useful to write down the parity automorphism
$\alpha(e_{I})=(-1)^{|I|}e_{I}$ under this matrix representation. As the
notion of parity is equivalent between $\mathrm{Cl}(n)$ and
$\mathcal{H}_{2^{n}}$, $\alpha$ is simply the $n$-qubit parity operator,

$\alpha\equiv Z^{\otimes n}.$ (813)

It will also be useful to represent the subspace $\mathrm{Cl}^{0}(n)$
explicitly as an $(n-1)$-qubit Hilbert space. This is achieved by the
projection from $\mathrm{Cl}(n)$ to $\mathrm{Cl}^{0}(n)$, expressed in Ref.
[saunderson2015semidefinite] as the $2^{n-1}\times 2^{n}$ matrix

$\Pi_{0}\coloneqq\frac{1}{\sqrt{2}}\mathopen{}\left(\langle+|\otimes\mathbb{I}_{2}^{\otimes(n-1)}+\langle-|\otimes
Z^{\otimes(n-1)}\right)\mathclose{}.$ (814)

It is straightforward to check that $\Pi_{0}|b\rangle=0$ if
$|b|\>\mathrm{mod}\>2=1$, and that its image is a $2^{n-1}$-dimensional
Hilbert space.

#### The quadratic mapping as quantum expectation values

The quadratic map $Q:\mathrm{Cl}(n)\to\mathbb{R}^{n\times n}$ is defined as

$Q(x)(v)\coloneqq\pi_{\mathbb{R}^{n}}(\alpha(x)v\overline{x})\quad\forall
x\in\mathrm{Cl}(n),v\in\mathbb{R}^{n},$ (815)

where $\pi_{\mathbb{R}^{n}}$ is the projector from $\mathrm{Cl}(n)$ to
$\mathbb{R}^{n}$,

$\pi_{\mathbb{R}^{n}}(x)\coloneqq\sum_{i\in[n]}\langle e_{i},x\rangle
e_{i}\quad\forall x\in\mathrm{Cl}(n),$ (816)

and the conjugation operation $x\mapsto\overline{x}$ is defined as the linear
extension of $\overline{e_{I}}=(-1)^{|I|}e_{i_{k}}\cdots e_{i_{1}}$. This map
associates Clifford algebra elements with orthogonal matrices via the
relations $Q(\mathrm{Pin}(n))=\mathrm{O}(n)$ and
$Q(\mathrm{Spin}(n))=\mathrm{SO}(n)$ (see Appendix Clifford algebras and the
orthogonal group for a review of the construction). In the standard basis of
$\mathbb{R}^{n}$, the linear map $Q(x):\mathbb{R}^{n}\to\mathbb{R}^{n}$ has
the matrix elements

$\begin{split}[Q(x)]_{ij}&=\langle e_{i},Q(x)(e_{j})\rangle\\\ &=\langle
e_{i},\alpha(x)e_{j}\overline{x}\rangle.\end{split}$ (817)

Using the linear maps $\lambda_{i},\rho_{j}$ of left- and right-multiplication
by $e_{i},e_{j}$, as well as the conjugation identity $\langle
x,y\overline{z}\rangle=\langle xz,y\rangle$ in the Clifford algebra, these
matrix elements of $Q(x)$ can be rearranged as

$\begin{split}[Q(x)]_{ij}&=\langle e_{i},\alpha(x)e_{j}\overline{x}\rangle\\\
&=\langle e_{i}x,\alpha(x)e_{j}\rangle\\\
&=\langle\lambda_{i}(x),\rho_{j}(\alpha(x))\rangle\\\ &=\langle
x,\lambda_{i}^{\dagger}(\rho_{j}(\alpha(x)))\rangle.\end{split}$ (818)

We now transfer this expression to the quantum representation developed above.
First, define the following $n$-qubit Pauli operators as the composition of
the linear maps appearing in Eq. (818):

$P_{ij}\coloneqq\lambda_{i}^{\dagger}\rho_{j}\alpha=\begin{cases}-\mathbb{I}_{2}^{\otimes(i-1)}\otimes
X\otimes Z^{\otimes(j-i-1)}\otimes
X\otimes\mathbb{I}_{2}^{\otimes(n-j)}&i<j,\\\
\mathbb{I}_{2}^{\otimes(i-1)}\otimes
Z\otimes\mathbb{I}_{2}^{\otimes(n-i)}&i=j,\\\
-\mathbb{I}_{2}^{\otimes(j-1)}\otimes Y\otimes Z^{\otimes(i-j-1)}\otimes
Y\otimes\mathbb{I}_{2}^{\otimes(n-i)}&i>j,\end{cases}$ (819)

where the expressions in terms of Pauli matrices follow from Eqs. (811) to
(813). Then we may rewrite Eq. (818) as

$[Q(x)]_{ij}=\langle x|P_{ij}|x\rangle,$ (820)

where $|x\rangle\in\mathcal{H}_{2^{n}}$ is the quantum state identified with
$x\in\mathrm{Cl}(n)$. Hence, the matrix elements of
$Q(x)\in\mathbb{R}^{n\times n}$ possess the interpretation as expectation
values of a collection of $n^{2}$ Pauli observables
$\\{P_{ij}\\}_{i,j\in[n]}$. Furthermore, recall that $Q(x)\in\mathrm{O}(n)$ if
and only if $x\in\mathrm{Pin}(n)$, and $Q(x)\in\mathrm{SO}(n)$ if and only if
$x\in\mathrm{Spin}(n)$. Because
$\mathrm{Spin}(n)=\mathrm{Pin}(n)\cap\mathrm{Cl}^{0}(n)$, one can work in the
even-parity sector directly by projecting the operators as

$\widetilde{P}_{ij}\coloneqq\Pi_{0}P_{ij}\Pi_{0}^{\mathsf{T}}.$ (821)

These are $(n-1)$-qubit Pauli operators, and we provide explicit expressions
in Appendix Details for working in the even-parity subspace. When necessary,
we may specify another map
$\widetilde{Q}:\mathrm{Cl}^{0}(n)\to\mathbb{R}^{n\times n}$,

$[\widetilde{Q}(x)]_{ij}\coloneqq\langle x|\widetilde{P}_{ij}|x\rangle,$ (822)

for which $\widetilde{Q}(\mathrm{Spin}(n))=\mathrm{SO}(n)$.

In general, these double covers are only a subset of the unit sphere in
$\mathcal{H}_{d}$ ($d=2^{n}$ or $2^{n-1}$), so not all quantum states mapped
by $Q$ yield orthogonal matrices. In Section Fermionic representation of the
construction we characterize the elements of $\mathrm{Pin}(n)$ and
$\mathrm{Spin}(n)$ as a class of well-studied quantum states, namely, pure
fermionic Gaussian states.

#### Fermionic representation of the construction

##### Notation

First we establish some notation. A system of $n$ fermionic modes, described
by the creation operators $a_{1}^{\dagger},\ldots,a_{n}^{\dagger}$, can be
equivalently represented by the $2n$ Majorana operators

$\displaystyle\gamma_{i}$ $\displaystyle=a_{i}+a_{i}^{\dagger},$ (823)
$\displaystyle\widetilde{\gamma}_{i}$
$\displaystyle=-\mathrm{i}(a_{i}-a_{i}^{\dagger}),$ (824)

for all $i\in[n]$. These operators form a representation for the Clifford
algebra $\mathrm{Cl}(2n)$, as they satisfy494949Note that we adopt the
physicist’s convention here, which takes the generators to be Hermitian, as
opposed to Eq. (893) wherein they square to $-\mathbb{I}$.

$\displaystyle\gamma_{i}\gamma_{j}+\gamma_{j}\gamma_{i}=\widetilde{\gamma}_{i}\widetilde{\gamma}_{j}+\widetilde{\gamma}_{j}\widetilde{\gamma}_{i}$
$\displaystyle=2\delta_{ij}\mathbb{I},$ (825)
$\displaystyle\gamma_{i}\widetilde{\gamma}_{j}+\widetilde{\gamma}_{j}\gamma_{i}$
$\displaystyle=0.$ (826)

The Jordan–Wigner mapping allows us to identify this fermionic system with an
$n$-qubit system via the relations

$\displaystyle\gamma_{i}$ $\displaystyle=Z^{\otimes(i-1)}\otimes
X\otimes\mathbb{I}_{2}^{\otimes(n-i)},$ (827)
$\displaystyle\widetilde{\gamma}_{i}$ $\displaystyle=Z^{\otimes(i-1)}\otimes
Y\otimes\mathbb{I}_{2}^{\otimes(n-i)}.$ (828)

We will work with the two representations interchangeably.

A central tool for describing noninteracting fermions is the Bogoliubov
transformation $\bm{\gamma}\mapsto O\bm{\gamma}$, where $O\in\mathrm{O}(2n)$
and

$\bm{\gamma}\coloneqq\begin{bmatrix}\widetilde{\gamma}_{1}&\cdots&\widetilde{\gamma}_{n}&\gamma_{1}&\cdots&\gamma_{n}\end{bmatrix}^{\mathsf{T}}.$
(829)

This transformation is achieved by fermionic Gaussian unitaries, which are
equivalent to matchgate circuits on qubits under the Jordan–Wigner mapping
[knill2001fermionic, terhal2002classical, jozsa2008matchgates]. In particular,
we will make use of a subgroup of such unitaries corresponding to
$\mathrm{O}(n)\times\mathrm{O}(n)\subset\mathrm{O}(2n)$. For any
$U,V\in\mathrm{O}(n)$, let $\mathcal{U}_{(U,V)}$ be the fermionic Gaussian
unitary with the adjoint action

$\displaystyle\mathcal{U}_{(U,V)}\widetilde{\gamma}_{i}\mathcal{U}_{(U,V)}^{\dagger}$
$\displaystyle=\sum_{j\in[n]}U_{ij}\widetilde{\gamma}_{j},$ (830)
$\displaystyle\mathcal{U}_{(U,V)}\gamma_{i}\mathcal{U}_{(U,V)}^{\dagger}$
$\displaystyle=\sum_{j\in[n]}V_{ij}\gamma_{j}.$ (831)

In contrast to arbitrary $\mathrm{O}(2n)$ transformations, these unitaries do
not mix between the $\gamma$\- and $\widetilde{\gamma}$-type Majorana
operators.

##### Linear optimization as free-fermion models

Applying the representation of Majorana operators under the Jordan–Wigner
transformation, Eqs. (827) and (828), to the Clifford algebra automorphisms,
Eqs. (811) to (813), we see that
$\lambda_{i}^{\dagger}=\mathrm{i}\widetilde{\gamma}_{i}$ and
$\rho_{j}\alpha=\gamma_{j}$. Therefore the Pauli operators $P_{ij}$ defining
the quadratic map $Q$ are equivalent to fermionic one-body operators,

$P_{ij}=\mathrm{i}\widetilde{\gamma}_{i}\gamma_{j}.$ (832)

Consider now a linear objective function $\ell(X)\coloneqq\langle C,X\rangle$
for some fixed $C\in\mathbb{R}^{n\times n}$, which we wish to optimize over
$\mathrm{O}(n)$:

$\max_{X\in\mathrm{O}(n)}\ell(X)=\max_{X\in\mathrm{O}(n)}\langle C,X\rangle.$
(833)

Because we require $X\in\mathrm{O}(n)$, it is equivalent to search over all
$x\in\mathrm{Pin}(n)$ through $Q$:

$\max_{X\in\mathrm{O}(n)}\langle C,X\rangle=\max_{x\in\mathrm{Pin}(n)}\langle
C,Q(x)\rangle.$ (834)

Writing out the matrix elements explicitly, we see that the objective takes
the form

$\begin{split}\ell(X)&=\sum_{i,j\in[n]}C_{ij}[Q(x)]_{ij}\\\
&=\sum_{i,j\in[n]}C_{ij}\langle x|P_{ij}|x\rangle\\\ &=\langle
x|F(C)|x\rangle,\end{split}$ (835)

where we have defined the noninteracting fermionic Hamiltonian

$F(C)\coloneqq\sum_{i,j\in[n]}C_{ij}P_{ij}=\mathrm{i}\sum_{i,j\in[n]}C_{ij}\widetilde{\gamma}_{i}\gamma_{j}.$
(836)

The linear optimization problem is therefore equivalent to solving a free-
fermion model,

$\max_{X\in\mathrm{O}(n)}\langle C,X\rangle=\max_{x\in\mathrm{Pin}(n)}\langle
x|F(C)|x\rangle,$ (837)

the eigenvectors of which are fermionic Gaussian states. As such, this problem
can be solved efficiently by a classical algorithm. In fact, the known
classical algorithm for solving the optimization problem is exactly the same
as that used for diagonalizing $F(C)$.

We now review the standard method to diagonalize $F(C)$. Consider the singular
value decomposition of $C=U\Sigma V^{\mathsf{T}}$, which is computable in time
$\mathcal{O}(n^{3})$. This decomposition immediately reveals the diagonal form
of the Hamiltonian:

$\begin{split}F(C)&=\mathrm{i}\sum_{i,j\in[n]}[U\Sigma
V^{\mathsf{T}}]_{ij}\widetilde{\gamma}_{i}\gamma_{j}\\\
&=\mathrm{i}\sum_{k\in[n]}\sigma_{k}(C)\mathopen{}\left(\sum_{i\in[n]}[U^{\mathsf{T}}]_{ki}\widetilde{\gamma}_{i}\right)\mathclose{}\mathopen{}\left(\sum_{j\in[n]}[V^{\mathsf{T}}]_{kj}\gamma_{j}\right)\mathclose{}\\\
&={\mathcal{U}}_{(U,V)}^{\dagger}\mathopen{}\left(\sum_{k\in[n]}\sigma_{k}(C)\mathrm{i}\widetilde{\gamma}_{k}\gamma_{k}\right)\mathclose{}{\mathcal{U}}_{(U,V)}.\end{split}$
(838)

Because $\mathrm{i}\widetilde{\gamma}_{k}\gamma_{k}=Z_{k}$, it follows that
the eigenvectors of $F(C)$ are the fermionic Gaussian states

$|\phi_{b}\rangle={\mathcal{U}}_{(U,V)}^{\dagger}|b\rangle,\quad
b\in\\{0,1\\}^{n},$ (839)

with eigenvalues

$E_{b}=\sum_{k\in[n]}(-1)^{b_{k}}\sigma_{k}(C).$ (840)

The maximum energy is $E_{0^{n}}=\operatorname{tr}\Sigma$ since all singular
values are nonnegative. The corresponding eigenstate $|\phi_{0^{n}}\rangle$ is
the maximizer of Eq. (837), so it corresponds to an element
$\phi_{0^{n}}\in\mathrm{Pin}(n)$. It is straightforward to see this by
recognizing that
$[Q(\phi_{0^{n}})]_{ij}=\langle\phi_{0^{n}}|\mathrm{i}\widetilde{\gamma}_{i}\gamma_{j}|\phi_{0^{n}}\rangle=[UV^{\mathsf{T}}]_{ij}$.
The fact that $Q(\phi_{0^{n}})\in\mathrm{O}(n)$ if and only if
$\phi_{0^{n}}\in\mathrm{Pin}(n)$ concludes the argument.

Indeed, the standard classical algorithm [schonemann1966generalized] for
solving Eq. (833) uses precisely the same decomposition. From the cyclic
property of the trace and the fact that $\mathrm{O}(n)$ is a group, we have

$\max_{X\in\mathrm{O}(n)}\langle U\Sigma
V^{\mathsf{T}},X\rangle=\max_{X^{\prime}\in\mathrm{O}(n)}\langle\Sigma,X^{\prime}\rangle,$
(841)

where we have employed the change of variables $X^{\prime}\coloneqq
U^{\mathsf{T}}XV$. Again, because $\Sigma$ has only nonnegative entries,
$\langle\Sigma,X^{\prime}\rangle$ achieves its maximum,
$\operatorname{tr}\Sigma$, when $X^{\prime}=\mathbb{I}_{n}$. This implies that
the optimal solution is $X=UV^{\mathsf{T}}$. Note that this problem is
equivalent to minimizing the Frobenius-norm distance, since

$\begin{split}\operatorname*{arg\,min}_{X\in\mathrm{O}(n)}\|C-X\|_{F}^{2}&=\operatorname*{arg\,min}_{X\in\mathrm{O}(n)}\mathopen{}\left(\|C\|_{F}^{2}+\|X\|_{F}^{2}-2\langle
C,X\rangle\right)\mathclose{}\\\
&=\operatorname*{arg\,max}_{X\in\mathrm{O}(n)}\,\langle
C,X\rangle.\end{split}$ (842)

Now suppose we wish to optimize $\ell$ over $\mathrm{SO}(n)$. In this setting,
one instead computes $X=U\widetilde{V}^{\mathsf{T}}$ from the special singular
value decomposition of $C=U\widetilde{\Sigma}\widetilde{V}^{\mathsf{T}}$. This
ensures that $\det(X)=1$ while maximizing $\ell(X)$, as only the smallest
singular value $\sigma_{n}(C)$ has its sign potentially flipped to guarantee
the positive determinant constraint. This sign flip also has a direct analogue
within the free-fermion perspective. Recall that the determinant of
$Q(x)\in\mathrm{O}(n)$ is given by the parity of $x\in\mathrm{Pin}(n)$, or
equivalently the parity of the state $|x\rangle$ in the computational basis.
Note also that all fermionic states are eigenstates of the parity operator. To
optimize over $\mathrm{SO}(n)$, we therefore seek the maximal eigenstate
$|\phi_{b}\rangle$ of $F(C)$ which has even parity. If
$\langle\phi_{0^{n}}|Z^{\otimes n}|\phi_{0^{n}}\rangle=1$ then we are done. On
the other hand, if $\langle\phi_{0^{n}}|Z^{\otimes n}|\phi_{0^{n}}\rangle=-1$
then we need to flip only a single bit in $0^{n}$ to reach an even-parity
state. The smallest change in energy by such a flip is achieved from changing
the occupation of the mode corresponding to the smallest singular value of
$C$. The resulting eigenstate $|\phi_{0^{n-1}1}\rangle$ is then the even-
parity state with the largest energy,
$E_{0^{n-1}1}=\operatorname{tr}\Sigma-2\sigma_{n}(C)$.

Finally, we point out that all elements of $\mathrm{Pin}(n)$ are free-fermion
states. To see this, observe that $C$ is arbitrary. We can therefore construct
the family of Hamiltonians $\\{F(C)\mid C\in\mathrm{O}(n)\\}$. Clearly, the
maximum $\langle C,X\rangle=n$ within this family is achieved when $X=C$, each
of which corresponds to a fermionic Gaussian state $|\phi\rangle$ satisfying
$F(C)|\phi\rangle=n|\phi\rangle$ and $Q(\phi)=C$. We note that this argument
generalizes the mathematical one presented in Ref.
[saunderson2015semidefinite], which only considered the eigenvectors lying in
$\mathrm{Spin}(n)$.

##### Mixed states and the convex hull

First we review descriptions of the convex hull of orthogonal and rotation
matrices, the latter of which was characterized by Saunderson _et al._
[saunderson2015semidefinite]. The convex hull of $\mathrm{O}(n)$ is the set of
all matrices with operator norm bounded by 1,

$\operatorname{conv}\mathrm{O}(n)=\mathopen{}\left\\{X\in\mathbb{R}^{n\times
n}\mid\sigma_{1}(X)\leq 1\right\\}\mathclose{}.$ (843)

On the other hand, the convex hull of $\mathrm{SO}(n)$ has a more complicated
description in terms of special singular values:

$\operatorname{conv}\mathrm{SO}(n)=\mathopen{}\left\\{X\in\mathbb{R}^{n\times
n}\mathrel{\Bigg{|}}\sum_{i\in[n]\setminus
I}\widetilde{\sigma}_{i}(X)-\sum_{i\in I}\widetilde{\sigma}_{i}(X)\leq
n-2\quad\forall I\subseteq[n],|I|\text{ odd}\right\\}\mathclose{}.$ (844)

Saunderson _et al._ [saunderson2015semidefinite] establish that this convex
body is a spectrahedron, the feasible region of a semidefinite program. The
representation that we will be interested in is called a PSD lift:

$\operatorname{conv}\mathrm{SO}(n)=\mathopen{}\left\\{\begin{bmatrix}\langle\widetilde{P}_{11},\rho\rangle&\cdots&\langle\widetilde{P}_{1n},\rho\rangle\\\
\vdots&\ddots&\vdots\\\
\langle\widetilde{P}_{n1},\rho\rangle&\cdots&\langle\widetilde{P}_{nn},\rho\rangle\end{bmatrix}\mathrel{\Bigg{|}}\rho\succeq
0,\operatorname{tr}\rho=1\right\\}\mathclose{},$ (845)

where the $2^{n-1}\times 2^{n-1}$ matrices $\widetilde{P}_{ij}$ are defined in
Eq. (821).505050Technically, Saunderson _et al._ [saunderson2015semidefinite]
use the definition
$\widetilde{P}_{ij}=-\Pi_{0}\lambda_{i}\rho_{j}\Pi_{0}^{\mathsf{T}}$ because
they employ the standard adjoint representation, which differs from our use of
the twisted adjoint representation which includes the parity automorphism
$\alpha$. However since $\alpha(x)=x$ for all $x\in\mathrm{Cl}^{0}(n)$, both
definitions of $\widetilde{P}_{ij}$ coincide.

Recall that the density operators on a Hilbert space $\mathcal{H}$ form the
convex hull of its pure states:

$\mathcal{D}(\mathcal{H})\coloneqq\operatorname{conv}\\{|\psi\rangle\\!\langle\psi|\mid|\psi\rangle\in\mathcal{H},\langle\psi|\psi\rangle=1\\}=\\{\rho\in\mathcal{L}(\mathcal{H})\mid\rho\succeq
0,\operatorname{tr}\rho=1\\}.$ (846)

From Eq. (845) one immediately recognizes that the PSD lift of
$\operatorname{conv}\mathrm{SO}(n)$ corresponds to
$\mathcal{D}(\mathcal{H}_{2^{n-1}})$, where we recognize that
$\mathcal{H}_{2^{n-1}}\cong\mathrm{Cl}^{0}(n)$. Furthermore, the projection of
the lift is achieved through the convexification of the map
$Q:\mathrm{Cl}(n)\to\mathbb{R}^{n\times n}$, where the fact that $Q$ is
quadratic in $\mathrm{Cl}(n)$ translates to being linear in
$\mathcal{D}(\mathrm{Cl}(n))$. Specifically, by a slight abuse of notation we
shall extend the definition of $Q$ to act on density operators as

$Q(\rho)=\sum_{\mu}p_{\mu}Q(x_{\mu}),\quad\text{where
}\rho=\sum_{\mu}p_{\mu}|x_{\mu}\rangle\\!\langle x_{\mu}|.$ (847)

Then Eq. (845) is the statement that
$Q(\mathcal{D}(\mathrm{Cl}^{0}(n)))=\operatorname{conv}\mathrm{SO}(n)$.

In Appendix PSD lift of $\operatorname{conv}\mathrm{O}(n)$ and
$\operatorname{conv}\mathrm{SO}(n)$ we show that this statement
straightforwardly generalizes for
$Q(\mathcal{D}(\mathrm{Cl}(n)))=\operatorname{conv}\mathrm{O}(n)$. We prove
this using the fermionic representation developed in Section Linear
optimization as free-fermion models, and furthermore use these techniques to
provide an alternative derivation for the PSD lift of
$\operatorname{conv}\mathrm{SO}(n)$. The core of our argument is showing that
the singular-value conditions of Eqs. (843) and (844) translate into bounds on
the largest eigenvalue of corresponding $n$-qubit observables:

$\displaystyle\sigma_{i}(X)$
$\displaystyle=\operatorname{tr}(\mathrm{i}\widetilde{\gamma}_{i}\gamma_{i}\rho)\leq
1,$ (848) $\displaystyle\sum_{i\in[n]\setminus
I}\widetilde{\sigma}_{i}(X)-\sum_{i\in I}\widetilde{\sigma}_{i}(X)$
$\displaystyle=\operatorname{tr}\mathopen{}\left[\rho_{0}\mathopen{}\left(\sum_{i\in[n]\setminus
I}\mathrm{i}\widetilde{\gamma}_{i}\gamma_{i}-\sum_{i\in
I}\mathrm{i}\widetilde{\gamma}_{i}\gamma_{i}\right)\mathclose{}\right]\mathclose{}\leq
n-2,$ (849)

where $\rho\in\mathcal{D}(\mathrm{Cl}(n))$ and
$\rho_{0}\in\mathcal{D}(\mathrm{Cl}^{0}(n))$. The physical interpretation here
is that not all pure quantum states map onto to orthogonal or rotation
matrices (which is clear from the fact that fermionic Gaussian states are only
a subset of quantum states). However, all density operators _do_ map onto to
their convex hulls, and the distinction between
$\operatorname{conv}\mathrm{O}(n)$ and $\operatorname{conv}\mathrm{SO}(n)$ can
be automatically specified by restricting the support of $\rho$ to the even-
parity subspace.

### Quantum relaxation for the quadratic problem

We now arrive at the primary problem of interest in this work, the little
noncommutative Grothendieck problem over the (special) orthogonal group. While
the linear problem of Eq. (833) can be solved classically in polynomial time,
quadratic programs are considerably more difficult. Here, we use the quantum
formalism of the Pin and Spin groups developed above to construct a quantum
relaxation of this problem. Then in Section Rounding algorithms we describe
rounding procedures to recover a collection of orthogonal matrices from the
quantum solution to this relaxation.

Recall the description of the input to Problem (781). Let $(V,E)$ be a graph,
and associate to each edge $(u,v)\in E$ a matrix $C_{uv}\in\mathbb{R}^{n\times
n}$. We label the vertices as $V=[m]$. We wish to maximize the objective

$f(R_{1},\ldots,R_{m})=\sum_{(u,v)\in E}\langle
C_{uv},R_{u}R_{v}^{\mathsf{T}}\rangle.$ (850)

over $(R_{1},\ldots,R_{m})\in\mathrm{O}(n)^{m}$. First, expand this expression
in terms of matrix elements:

$\sum_{(u,v)\in E}\langle C_{uv},R_{u}R_{v}^{\mathsf{T}}\rangle=\sum_{(u,v)\in
E}\sum_{i,j\in[n]}[C_{uv}]_{ij}\sum_{k\in[n]}[R_{u}]_{ik}[R_{v}^{\mathsf{T}}]_{kj}.$
(851)

From the quadratic mapping $Q:\mathrm{Cl}(n)\to\mathbb{R}^{n\times n}$, we
know that for each $R\in G$ there exists some $\phi\in\mathrm{Pin}(n)$ such
that $R_{ij}=\langle\phi|P_{ij}|\phi\rangle$. Hence we can express the matrix
product as

$\begin{split}[R_{u}]_{ik}[R_{v}^{\mathsf{T}}]_{kj}&=\langle\phi_{u}|P_{ik}|\phi_{u}\rangle\langle\phi_{v}|P_{jk}|\phi_{v}\rangle\\\
&=\langle\phi_{u}\otimes\phi_{v}|P_{ik}\otimes
P_{jk}|\phi_{u}\otimes\phi_{v}\rangle,\end{split}$ (852)

which is now the expectation value of a $2n$-qubit Pauli operator with respect
to a product state of two Gaussian states $|\phi_{u}\rangle$,
$|\phi_{v}\rangle$. To extend this over the entire graph, we define a Hilbert
space of $m$ registers of $n$ qubits each. For each edge $(u,v)\in E$ we
introduce the Hamiltonian terms

$H_{uv}\coloneqq\sum_{i,j\in[n]}[C_{uv}]_{ij}\Gamma_{ij}^{(u,v)},$ (853)

where

$\Gamma_{ij}^{(u,v)}\coloneqq\mathopen{}\left(\sum_{k\in[n]}P_{ik}^{(u)}\otimes
P_{jk}^{(v)}\right)\mathclose{}\bigotimes_{w\in
V\setminus\\{u,v\\}}\mathbb{I}_{2^{n}}^{(w)}.$ (854)

To simplify notation, we shall omit the trivial support $\bigotimes_{w\in
V\setminus\\{u,v\\}}\mathbb{I}_{2^{n}}^{(w)}$ when the context is clear.

The problem is now reformulated as optimizing the $mn$-qubit Hamiltonian

$H\coloneqq\sum_{(u,v)\in E}H_{uv}=\sum_{(u,v)\in
E}\sum_{i,j\in[n]}[C_{uv}]_{ij}\sum_{k\in[n]}P_{ik}^{(u)}\otimes
P_{jk}^{(v)}.$ (855)

The exact LNCG problem over $\mathrm{O}(n)$ then corresponds to

$\max_{\mathbf{R}\in\mathrm{O}(n)^{m}}f(\mathbf{R})=\max_{|\psi\rangle\in\mathcal{H}_{2^{n}}^{\otimes
m}}\langle\psi|H|\psi\rangle\quad\text{subject to
}\begin{cases}\langle\psi|\psi\rangle=1,\\\
|\psi\rangle=\bigotimes_{v\in[m]}|\phi_{v}\rangle,\\\
|\phi_{v}\rangle={\mathcal{U}}_{(R_{v},\mathbb{I}_{n})}|0^{n}\rangle,&R_{v}\in\mathrm{O}(n)\
\forall v\in[m].\end{cases}$ (856)

The hardness of this problem is therefore related to finding the optimal
separable state for local Hamiltonians, which is $\mathsf{NP}$-hard in general
[gurvits2004classical, ioannou2006computational, gharibian2008strong].
Dropping these constraints on the state provides a relaxation of the problem,
since

$\max_{\begin{subarray}{c}|\psi\rangle\in\mathcal{H}_{2^{n}}^{\otimes m},\\\
\langle\psi|\psi\rangle=1\end{subarray}}\langle\psi|H|\psi\rangle\geq\max_{\mathbf{R}\in
G^{m}}f(\mathbf{R}).$ (857)

We point out here that the Hamiltonian terms $H_{uv}$ can be interpreted as
two-body fermionic interactions. Note that there is an important distinction
between two-body fermionic operators (Clifford-algebra products of four
Majorana operators) and two-body qudit operators (tensor products of two qudit
Pauli operators). Recall that
$P_{ij}=\mathrm{i}\widetilde{\gamma}_{i}\gamma_{j}$ is one-body in the
fermionic sense. While the operators $P_{ik}^{(u)}\otimes P_{jk}^{(v)}$ appear
to mix both notions, here they in fact coincide. To see this, we consider a
global algebra of Majorana operators
$\\{\gamma_{i+(v-1)n},\widetilde{\gamma}_{i+(v-1)n}\mid i\in[n],v\in[m]\\}$
acting on a Hilbert space of $mn$ fermionic modes. While it is not true that
the local single-mode Majorana operators map onto the global single-mode
operators, i.e.,

$\displaystyle\gamma_{i}^{(v)}\bigotimes_{w\in
V\setminus\\{v\\}}\mathbb{I}_{2^{n}}^{(w)}\neq\gamma_{i+(v-1)n},$ (858)

the local two-mode Majorana operators in fact do correspond to global two-mode
operators:

$\widetilde{\gamma}_{i}^{(v)}\gamma_{j}^{(v)}\bigotimes_{w\in
V\setminus\\{v\\}}\mathbb{I}_{2^{n}}^{(w)}=\widetilde{\gamma}_{i+(v-1)n}\gamma_{j+(v-1)n}.$
(859)

Thus, taking the tensor product of two local two-mode Majorana operators on
different vertices is equivalent to taking the product of two global two-mode
Majorana operators:

$\widetilde{\gamma}_{i}^{(u)}\gamma_{j}^{(u)}\otimes\widetilde{\gamma}_{k}^{(v)}\gamma_{l}^{(v)}\bigotimes_{w\in
V\setminus\\{u,v\\}}\mathbb{I}_{2^{n}}^{(w)}=\widetilde{\gamma}_{i+(u-1)n}\gamma_{j+(u-1)n}\widetilde{\gamma}_{k+(v-1)n}\gamma_{l+(v-1)n}.$
(860)

Therefore Eq. (855) can be equivalently expressed as a Hamiltonian with two-
body fermionic interactions.

Finally, when we wish to optimize over
$(R_{1},\ldots,R_{m})\in\mathrm{SO}(n)^{m}$, it is straightforward to see that
we can simply replace the terms $P_{ij}$ with $\widetilde{P}_{ij}$. Defining

$\displaystyle\widetilde{\Gamma}_{ij}^{(u,v)}$
$\displaystyle\coloneqq\mathopen{}\left(\sum_{k\in[n]}\widetilde{P}_{ik}^{(u)}\otimes\widetilde{P}_{jk}^{(v)}\right)\mathclose{}\bigotimes_{w\in
V\setminus\\{u,v\\}}\mathbb{I}_{2^{n-1}}^{(w)},$ (861)
$\displaystyle\widetilde{H}_{uv}$
$\displaystyle\coloneqq\sum_{i,j\in[n]}[C_{uv}]_{ij}\widetilde{\Gamma}_{ij}^{(u,v)},$
(862)

the quantum relaxation for the $\mathrm{SO}(n)$ problem is given by the
$m(n-1)$-qubit Hamiltonian

$\widetilde{H}\coloneqq\sum_{(u,v)\in E}\widetilde{H}_{uv}.$ (863)

### Rounding algorithms

Optimizing the energy of a local Hamiltonian is a well-studied problem, both
from the perspective of quantum and classical algorithms. In this section we
will assume that such an algorithm has been used to produce the state
$\rho\in\mathcal{D}(\mathcal{H}_{d}^{\otimes m})$ which (approximately)
maximizes the energy $\operatorname{tr}(H\rho)$. We wish to round this state
into the feasible space, namely the set of product states of Gaussian states.
We do so by rounding the expectation values of $\rho$ appropriately, such that
we return some valid approximation $R_{1},\ldots,R_{m}\in G$. In this section
we propose two approaches to perform this quantum rounding.

The first uses insight from the fact that our quantum relaxation is equivalent
to a classical semidefinite relaxation with additional constraints based on
the convex hull of the orthogonal group. This is approach is particularly
advantageous when optimizing over $\mathrm{SO}(n)$, as
$\operatorname{conv}\mathrm{SO}(n)$ has a matrix representation exponential in
$n$ (its PSD lift). To build the semidefinite variable from the quantum state,
we require measurements of the expectation values of the two-vertex operators
$\Gamma_{ij}^{(u,v)}=\sum_{k\in[n]}P_{ik}^{(u)}\otimes P_{jk}^{(v)}$ for each
pair of vertices $(u,v)$. We refer this procedure as
_$\operatorname{conv}\mathrm{SO}(n)$ -based rounding_.515151This rounding can
also be applied to the optimization problem over $\mathrm{O}(n)$ as well, but
we are particularly interested in the $\operatorname{conv}\mathrm{SO}(n)$
constraints due to their exponentially large classical representation. Our
second rounding protocol uses the expectation values of $P_{ij}^{(v)}$ of each
vertex $v$ directly. In this case, rather than expectation values of two-
vertex operators as before, we only require the information of single-vertex
marginals $\rho_{v}\coloneqq\operatorname{tr}_{\neg v}(\rho)$. Therefore we
call this approach _vertex-marginal rounding_.

If $\rho$ is produced by a deterministic classical algorithm, then the
relevant expectation values can be exactly computed (to machine precision).
However if the state is produced by a randomized algorithm, or is otherwise
prepared by a quantum computer, then we can only estimate the expectation
values to within statistical error by some form of sampling. In the quantum
setting, this can be achieved either by partial state tomography
[bonet2020nearly, zhao2021fermionic] or a more sophisticated measurement
protocol [huggins2021nearly].525252For the present discussion we do not
consider the effects of finite sampling, although we expect that rounding is
fairly robust to such errors since it will always return a solution in the
feasible space. See Appendix Measurement schemes for further comments on this
quantum measurement aspect. The rounding algorithms then operate entirely as
classical postprocessing after estimating the necessary expectation values.

#### Approximation ratio for rounding the classical SDP

Before describing our quantum rounding protocols, we first review classical
relaxations and rounding procedures for Problem (781). The standard
semidefinite relaxation can be expressed as the SDP

$\max_{M\in\mathbb{R}^{mn\times mn}}\langle C,M\rangle\quad\text{subject to
}\begin{cases}M\succeq 0,\\\ M_{vv}=\mathbb{I}_{n}&\forall
v\in[m],\end{cases}$ (864)

where $C\in\mathbb{R}^{mn\times mn}$ is the matrix with $n\times n$ blocks
$C_{uv}$. If an additional nonconvex constraint $\mathrm{rank}(M)=n$ is
imposed, then the solution would be exact:

$M=\mathbf{R}\mathbf{R}^{\mathsf{T}}=\begin{bmatrix}\mathbb{I}_{n}&R_{1}R_{2}^{\mathsf{T}}&\cdots&R_{1}R_{m}^{\mathsf{T}}\\\
R_{2}R_{1}^{\mathsf{T}}&\mathbb{I}_{n}&\cdots&R_{2}R_{m}^{\mathsf{T}}\\\
\vdots&\vdots&\ddots&\vdots\\\
R_{m}R_{1}^{\mathsf{T}}&R_{m}R_{2}^{\mathsf{T}}&\cdots&\mathbb{I}_{n}\end{bmatrix}.$
(865)

Problem (864) is therefore a relaxation of the original problem. However, the
solution $M\in\mathbb{R}^{mn\times mn}$ is still PSD, so it can be decomposed
as $M=\mathbf{X}\mathbf{X}^{\mathsf{T}}$, where

$\mathbf{X}=\begin{bmatrix}X_{1}\\\ \vdots\\\ X_{m}\end{bmatrix},\quad
X_{v}\in\mathbb{R}^{n\times mn}.$ (866)

The rounding algorithm of Bandeira _et al._ [bandeira2016approximating] then
computes, for each $v\in[m]$,

$O_{v}=\mathcal{P}(X_{v}Z)\coloneqq\operatorname*{arg\,min}_{Y\in\mathrm{O}(n)}\|Y-X_{v}Z\|_{F},$
(867)

where $Z$ is an $mn\times n$ Gaussian random matrix whose entries are drawn
i.i.d. from $\mathcal{N}(0,1/n)$. When optimizing over $G=\mathrm{O}(n)$, this
rounded solution guarantees (in expectation) an approximation ratio of

$\alpha_{\mathrm{O}(n)}^{2}=\operatorname*{\mathbb{E}}\mathopen{}\left[\frac{1}{n}\sum_{i\in[n]}\sigma_{i}(Z_{1})\right]^{2}\mathclose{},$
(868)

where $Z_{1}$ is a random $n\times n$ matrix with i.i.d. entries from
$\mathcal{N}(0,1/n)$.

In Appendix Classical approximation ratio for $\mathrm{SO}(n)$ we extend the
argument used to obtain this result for the optimization problem over
$G=\mathrm{SO}(n)$, and we show a corresponding approximation ratio of

$\alpha_{\mathrm{SO}(n)}^{2}=\operatorname*{\mathbb{E}}\mathopen{}\left[\frac{1}{n}\sum_{i\in[n-1]}\sigma_{i}(Z_{1})\right]^{2}\mathclose{}$
(869)

where the only change to the rounding algorithm is that we project to the
nearest $\mathrm{SO}(n)$ element via $\widetilde{\mathcal{P}}$, which is
defined as

$\widetilde{\mathcal{P}}(X)\coloneqq\operatorname*{arg\,min}_{Y\in\mathrm{SO}(n)}\|Y-X\|_{F}.$
(870)

Note that singular values are nonnegative, and in particular we show that
$\operatorname*{\mathbb{E}}[\sigma_{n}(Z_{1})]>0$ for all finite $n$. Hence it
follows that $\alpha_{\mathrm{SO}(n)}^{2}<\alpha_{\mathrm{O}(n)}^{2}$, which
provides evidence for the claim that solving for rotations is generally a more
difficult problem (see Ref. [pumir2021generalized, Section 4.3] for a brief
discussion). For small values of $n$, the numerical values of these
approximation ratios are (computed using Mathematica):

$\displaystyle\alpha_{\mathrm{O}(2)}^{2}\approx 0.6564,$
$\displaystyle\qquad\alpha_{\mathrm{SO}(2)}^{2}\approx 0.3927,$ (871)
$\displaystyle\alpha_{\mathrm{O}(3)}^{2}\approx 0.6704,$
$\displaystyle\qquad\alpha_{\mathrm{SO}(3)}^{2}\approx 0.5476,$ (872)
$\displaystyle\alpha_{\mathrm{O}(4)}^{2}\approx 0.6795,$
$\displaystyle\qquad\alpha_{\mathrm{SO}(4)}^{2}\approx 0.6096.$ (873)

In Appendix Classical approximation ratio for $\mathrm{SO}(n)$ we provide an
integral expression for $\alpha_{\mathrm{SO}(n)}$ which can be evaluated for
arbitrary $n$.

For the problem over $\mathrm{SO}(n)$, Saunderson _et al._
[saunderson2014semidefinite] propose augmenting this SDP by adding the
constraints that each block of $M$ lies in
$\operatorname{conv}\mathrm{SO}(n)$:

$\max_{M\in\mathbb{R}^{mn\times mn}}\langle C,M\rangle\quad\text{subject to
}\begin{cases}M\succeq 0,\\\ M_{vv}=\mathbb{I}_{n}&\forall v\in[m],\\\
M_{uv}\in\operatorname{conv}\mathrm{SO}(n)&\forall u,v\in[m].\end{cases}$
(874)

Although they do not prove approximation guarantees for this enhanced SDP,
they first show that, if one reintroduces the rank constraint on $M$, then the
convex constraint $M_{uv}\in\operatorname{conv}\mathrm{SO}(n)$ in fact
suffices to guarantee the much stronger condition $M_{uv}\in\mathrm{SO}(n)$.
Then, when dropping the rank constraint (but leaving the
$\operatorname{conv}\mathrm{SO}(n)$ constraint) they show that the relaxed
problem is still exact over certain types of graphs, such as tree graphs.
Finally, they provide numerical evidence that even when the relaxation is not
exact, it returns substantially more accurate approximations than the standard
SDP (864).

#### Quantum Gram matrix

Analogous to the classical SDP solution $M$, we can form a matrix
$\mathcal{M}\in\mathbb{R}^{mn\times mn}$ from the expectation values of $\rho$
as

$\mathcal{M}\coloneqq\begin{bmatrix}\mathbb{I}_{n}&T_{12}&\cdots&T_{1m}\\\
T_{21}&\mathbb{I}_{n}&\cdots&T_{2m}\\\ \vdots&\vdots&\ddots&\vdots\\\
T_{m1}&T_{m2}&\cdots&\mathbb{I}_{n}\end{bmatrix},$ (875)

where

$T_{uv}\coloneqq\frac{1}{n}\begin{bmatrix}\operatorname{tr}(\Gamma_{11}^{(u,v)}\rho)&\cdots&\operatorname{tr}(\Gamma_{1n}^{(u,v)}\rho)\\\
\vdots&\ddots&\vdots\\\
\operatorname{tr}(\Gamma_{n1}^{(u,v)}\rho)&\cdots&\operatorname{tr}(\Gamma_{nn}^{(u,v)}\rho)\end{bmatrix}$
(876)

and $T_{vu}=T_{uv}^{\mathsf{T}}$ for $u<v$. Just as $\langle C,M\rangle$ gives
the relaxed objective value (up to rescaling and constant shifts), here we
have that $\langle
C,\mathcal{M}\rangle=\frac{2}{n}\operatorname{tr}(H\rho)+\operatorname{tr}(C)$.
In Appendix Convex hull of orthogonal matrices and quantum states we show that
for any quantum state, $\mathcal{M}$ satisfies the following properties:

$\begin{cases}\mathcal{M}\succeq 0,\\\ \mathcal{M}_{vv}=\mathbb{I}_{n}&\forall
v\in[m],\\\ \mathcal{M}_{uv}\in\operatorname{conv}\mathrm{O}(n)&\forall
u,v\in[m].\end{cases}$ (877)

Furthermore, we show that when $\rho$ is supported only on the even subspace
of each single-vertex Hilbert space (or equivalently, if we replace
$\Gamma_{ij}^{(u,v)}$ with $\widetilde{\Gamma}_{ij}^{(u,v)}$ in Eq. (876)),
then

$\frac{1}{n}\begin{bmatrix}\operatorname{tr}(\widetilde{\Gamma}_{11}^{(u,v)}\rho)&\cdots&\operatorname{tr}(\widetilde{\Gamma}_{1n}^{(u,v)}\rho)\\\
\vdots&\ddots&\vdots\\\
\operatorname{tr}(\widetilde{\Gamma}_{n1}^{(u,v)}\rho)&\cdots&\operatorname{tr}(\widetilde{\Gamma}_{nn}^{(u,v)}\rho)\end{bmatrix}\in\operatorname{conv}\mathrm{SO}(n)\quad\forall\rho\in\mathcal{D}(\mathcal{H}_{2^{n-1}}^{\otimes
m}).$ (878)

Therefore when optimizing the relaxed Hamiltonian $\widetilde{H}$ for the
$\mathrm{SO}(n)$ setting, we are guaranteed to automatically satisfy the
$\operatorname{conv}\mathrm{SO}(n)$ constraints.

####  $\operatorname{conv}\mathrm{SO}(n)$-based rounding

Given the construction of the $\mathcal{M}$ from quantum expectation values,
we proceed to round the Gram matrix as in the classical SDP with
$\operatorname{conv}\mathrm{SO}(n)$ constraints [saunderson2014semidefinite].
This consists of computing the matrices

$R_{v}=\widetilde{\mathcal{P}}(\mathcal{M}_{1v}),$ (879)

where the projection to $\mathrm{SO}(n)$ can be efficiently computed from the
special singular value decomposition, i.e.,

$\widetilde{\mathcal{P}}(X)=U\widetilde{V}^{\mathsf{T}}$ (880)

(recall Eq. (807)). Our choice of rounding using the first $n\times mn$ “row”
of $\mathcal{M}$ amounts to fixing $R_{1}=\mathbb{I}_{n}$. We note that the
same rounding procedure can naturally be applied to the $\mathrm{O}(n)$
setting as well, replacing $\widetilde{\mathcal{P}}$ with $\mathcal{P}$.

#### Vertex-marginal rounding

The single-vertex marginals are obtained by tracing out the qudits associated
to all but one vertex $v\in[m]$,

$\rho_{v}=\operatorname{tr}_{\neg v}(\rho).$ (881)

As $\rho_{v}\in\mathcal{D}(\mathcal{H}_{d}^{\otimes m})$, from Section Mixed
states and the convex hull we have that $Q(\rho_{v})\in\operatorname{conv}G$,
where we linearly extend the definition of $Q$ to

$Q(\rho_{v})\coloneqq\begin{bmatrix}\operatorname{tr}(P_{11}^{(v)}\rho)&\cdots&\operatorname{tr}(P_{1n}^{(v)}\rho)\\\
\vdots&\ddots&\vdots\\\
\operatorname{tr}(P_{n1}^{(v)}\rho)&\cdots&\operatorname{tr}(P_{nn}^{(v)}\rho)\end{bmatrix}.$
(882)

The rounding scheme we propose here then projects $Q(\rho_{v})$ to $G$ using
either $\mathcal{P}$ or $\widetilde{\mathcal{P}}$:

$R_{v}=\operatorname*{arg\,min}_{Y\in G}\|Y-Q(\rho_{v})\|_{F}.$ (883)

We point out that the relaxed Hamiltonian only has two-vertex terms which we
seek to maximize. In Appendix LNCG Hamiltonian symmetries we show that $H$
commutes with $\mathcal{U}_{(\mathbb{I}_{n},V)}^{\otimes m}$ for all
$V\in\mathrm{O}(n)$, which we further show implies that $H$ may possess
eigenstates whose single-vertex marginals obey $Q(\sigma_{v})=0$. This
indicates that there may exist eigenstates of $H$ whose single-vertex
marginals yield no information, despite the fact that their two-vertex
marginals are nontrivial. In our numerical studies, we observe that breaking
this symmetry resolves this issue. We accomplish this by including small
perturbative one-body terms which correspond to the trace of $Q(\sigma_{v})$:

$H_{1}=\sum_{v\in[m]}\sum_{i\in[n]}P_{ii}^{(v)},$ (884)

Note that this trace quantity is importantly invariant with respect to the
choice of basis for $\mathbb{R}^{n}$. We then augment the objective
Hamiltonian with $H_{1}$, defining

$H^{\prime}(\zeta)\coloneqq H+\zeta H_{1}$ (885)

where $\zeta>0$ is a small regularizing parameter. While this one-body
perturbation does not correspond to any terms in the original quadratic
objective function, any arbitrarily small $\zeta>0$ suffices to break the
$\mathrm{O}(n)$ symmetry. Furthermore, the rounding procedure always
guarantees that the solution is projected back into the feasible space
$G^{m}$. When $G=\mathrm{SO}(n)$ we define $\widetilde{H}^{\prime}(\zeta)$
analogously.

### Numerical experiments

To explore the potential of our quantum relaxation and rounding procedures, we
performed numerical experiments on randomly generated instances of the group
synchronization problem. Because the Hilbert-space dimension grows
exponentially in both $m$ and $n$, our classical simulations here are limited
to small problem sizes. However, optimizing over rotations in $\mathbb{R}^{3}$
(requiring only two qubits per vertex) is highly relevant to many practical
applications, so here we focus on the problem of $\mathrm{SO}(3)$ group
synchronization. For example, this problem appears in the context of cryo-EM
as described in Section Group synchronization. To model the problem, we
generated random instances by selecting random three-regular graphs $([m],E)$,
uniformly randomly sampling $m$ rotations
$g_{1},\ldots,g_{m}\in\mathrm{SO}(n)$, and then constructing
$C_{uv}=g_{u}g_{v}^{\mathsf{T}}+\sigma W_{uv}$ for each $(u,v)\in E$, where
the Gaussian noise matrix $W_{uv}\in\mathbb{R}^{n\times n}$ has i.i.d.
elements drawn from $\mathcal{N}(0,1)$ and $\sigma\geq 0$ represents the
strength (standard deviation) of this noise.

While the classical $\operatorname{conv}\mathrm{SO}(n)$-based SDP is not
guaranteed to find the optimal solution, the problems studied here were
selected for such that this enhanced SDP in fact does solve the exact problem.
We verify this property by confirming that $\mathrm{rank}(M)=n$ before
rounding on each problem instance. In this way we are able to calculate an
approximation ratio for the other methods (as it is not clear how to solve for
the globally optimal solution in general, even with an exponential-time
classical algorithm). The methods compared here include our quantum relaxation
with $\operatorname{conv}\mathrm{SO}(n)$-based rounding (denoted CR), vertex-
marginal rounding (VR), and the classical SDP (without
$\operatorname{conv}\mathrm{SO}(n)$ constraints but using the
$\widetilde{\mathcal{P}}$ projection to guarantee that the rounded solutions
are elements of $\mathrm{SO}(n)$). When using the vertex-rounding method, we
employ $\widetilde{H}^{\prime}(\zeta)$ as the objective Hamiltonian with
$\zeta=10^{-6}$.

#### Exact eigenvectors

Figure 24: Approximation ratios for solutions obtained from rounding the
maximum eigenvector of the relaxed Hamiltonian $\widetilde{H}$. Violin plots
show the distribution of approximation ratios over 50 randomly generated
instance, and with the median being indicated by the center marker. CR refers
to rounding according to the $\operatorname{conv}\mathrm{SO}(n)$-based scheme
(Section $\operatorname{conv}\mathrm{SO}(n)$-based rounding), while VR denotes
the vertex-marginal rounding scheme (Section Vertex-marginal rounding). The
classical SDP solution was rounded by the standard randomized algorithm
[bandeira2016approximating], and we report the best solution over 1000
rounding trials. (Left) Varying the number of vertices $m$ in the graph
(random 3-regular graphs). Note that the number of qubits required here is
$2m$. (Right) Varying the noise strength parameter $\sigma$ which defines the
problem via $C_{uv}=g_{u}g_{v}^{\mathsf{T}}+\sigma W_{uv}$.

First, we consider the solution obtained by rounding the maximum eigenvector
of $\widetilde{H}$. Although the hardness of preparing such a state is
equivalent that of the ground-state problem, this nonetheless provides us with
a benchmark for the ultimate approximation quality of our quantum relaxation.
In Figure 24 we plot the approximation ratio of the rounded quantum states and
compare to that of the classical SDP on the same problem instances. Each
violin plot was constructed from the results of 50 random instances.

The results here demonstrate that, while the approximation quality of the
classical SDP quickly falls off with larger graph sizes, our rounded quantum
solutions maintain high approximation ratios, at least for the problem sizes
probed here. Notably, the $\operatorname{conv}\mathrm{SO}(n)$-based rounding
on the quantum state is significantly more powerful and consistent than the
vertex-marginal rounding. This feature is not unexpected since, as discussed
in Section Vertex-marginal rounding, we are maximizing an objective
Hamiltonian with only two-body terms, whereas the single-vertex rounding uses
strictly one-body expectation values. Furthermore, as demonstrated in previous
works [saunderson2014semidefinite, matni2014convex] the
$\operatorname{conv}\mathrm{SO}(n)$ constraints are powerful in practice, and
so we expect that the quantum rounding protocol which makes use of this
structure enjoys the same advantages.

Meanwhile, when varying the noise parameter $\sigma$, we observe that all
methods are fairly consistent. In particular, the
$\operatorname{conv}\mathrm{SO}(n)$-based rounding only shows an appreciable
decrease in approximation quality when the noise is considerable (note that
$\sigma=0.5\approx 10^{-0.3}$ is a relatively large amount of noise, since
$g_{u}g_{v}^{\mathsf{T}}$ is an orthogonal matrix and therefore has matrix
elements bounded in magnitude by 1).

#### Quasi-adiabatic state preparation

Because it may be unrealistic to prepare the maximum eigenvector of
$\widetilde{H}$, here we consider preparing states using ideas from adiabatic
quantum computation [RevModPhys.90.015002]. Specifically, we wish to
demonstrate that states whose relaxed energy may be far from the maximum
eigenvalue can still provide high-quality approximations after rounding. If
this is the case then we do not need to prepare very close approximations to
the maximum eigenstate of $\widetilde{H}$, so the rigorous conditions of
adiabatic state preparation may not be required in this context. Hence we
consider “quasi-adiabatic” state preparation, wherein we explore how time-
evolution speeds far from the adiabatic limit may still return high-quality
approximations. Our numerical experiments here provide a preliminary
investigation into this conjecture.

For simplicity of the demonstration, we consider a linear annealing schedule
according to the time-dependent Hamiltonian

$H(t)=\mathopen{}\left(1-\frac{t}{T}\right)\mathclose{}H_{i}+\frac{t}{T}H_{f},$
(886)

which prepares the state

$|\psi(T)\rangle=\mathcal{T}\exp\mathopen{}\left(-\mathrm{i}\int_{0}^{T}\mathrm{d}t\,H(t)\right)\mathclose{}|\psi(0)\rangle$
(887)

for some $T>0$, where $\mathcal{T}$ is the time-ordering operator. The final
Hamiltonian $H_{f}$ is the desired objective LNCG Hamiltonian,

$H_{f}=\widetilde{H}.$ (888)

The initial Hamiltonian $H_{i}$ is the parent Hamiltonian of the initial
state, which we choose to be the approximation obtained from the classical
SDP, as it can be obtained classically in polynomial time. Let
$R_{1},\ldots,R_{m}\in\mathrm{SO}(n)$ be the SDP solution. Our initial state
is then the product of Gaussian states

$|\psi(0)\rangle=\bigotimes_{v\in[m]}|\phi(R_{v})\rangle,$ (889)

where each $|\phi(R_{v})\rangle$ is the maximum eigenvector of the free-
fermion Hamiltonian

$F(R_{v})=\mathrm{i}\sum_{i,j\in[n]}[R_{v}]_{ij}\widetilde{\gamma}_{i}\gamma_{j}.$
(890)

Therefore the initial Hamiltonian $H_{i}$ is a sum of such free-fermion
Hamiltonians (here we include the even-subspace projection since we are
working with $\mathrm{SO}(n)$):

$H_{i}=\sum_{v\in[m]}\Pi_{0}F(R_{v})\Pi_{0}^{\mathsf{T}}.$ (891)

As a Gaussian state, $|\phi(R_{v})\rangle$ can be prepared exactly from a
quantum circuit of $\mathcal{O}(n^{2})$ gates [jiang2018quantum]. Note that
since we are working directly in the even subspace of $n-1$ qubits here, this
$n$-qubit circuit must be projected appropriately using $\Pi_{0}$. We discuss
how to perform this circuit recompilation in Appendix Details for working in
the even-parity subspace. We comment that this choice of initial state is that
of a mean-field state for non-number-preserving fermionic systems, for
instance as obtained from Hartree–Fock–Bogoliubov theory. Suitably, the final
Hamiltonian we evolve into is non-number-preserving two-body fermionic
Hamiltonian.

Figure 25: Demonstration of a typical instance of adiabatic state preparation
for preparing relaxed quantum solutions. The initial state is the product of
Gaussian states corresponding to the rounded solution of the classical SDP. As
the total evolution time $T$ increases, the evolution becomes more adiabatic,
indicated by the convergence of the relaxed value to the maximum eigenvalue
(in units of the original problem’s optimal value). The rounded solutions of
course can never exceed the original problem’s optimal value.

In adiabatic state preparation, the total evolution time $T$ controls how
close the final state $|\psi(T)\rangle$ is to the maximum eigenstate535353We
remind the reader that we are starting in the maximum eigenstate of the
initial Hamiltonian, whereas in the physics literature, adiabatic theorems are
typically stated in terms of ground states. Of course, the two perspectives
are equivalent by simply an overall sign change (note that all Hamiltonians
here are traceless). of the final Hamiltonian $H_{f}$. One metric of closeness
is how the energy of the prepared state,
$\langle\psi(T)|H_{f}|\psi(T)\rangle$, compares to the maximum eigenvalue of
$H_{f}$. On the other hand, as a relaxation, this maximum energy is already
larger than the optimal objective value of the original problem. We showcase
this in Figure 25, using one random problem instance as a demonstrative
(typical) example on a graph of $m=6$ vertices (12 qubits). For each total
evolution time point $T$, we computed $|\psi(T)\rangle$ by numerically
integrating the time-dependent Schrödinger equation, and we plot its relaxed
energy as well as its rounded objective values. For large $T$ we approach the
maximum eigenstate of $\widetilde{H}$ as expected (thereby also demonstrating
that the initial “mean-field” state $|\psi(0)\rangle$ has appreciable
overlap). Particularly interesting is the behavior for relatively small total
evolution times $T$, wherein the energy of $|\psi(T)\rangle$ is far from the
maximum eigenenergy. Despite this, the approximation quality after rounding
the state using $\mathcal{M}$ is nearly exact around $T\approx 1$. On the
other hand, the approximation quality of vertex-marginal rounding is highly
inconsistent, which again we attribute to the fact that the single-vertex
information is not directly seen by the final Hamiltonian $H_{f}$.

Figure 26: Approximation ratios for solutions obtained from rounding the
“adiabatically” evolved state $|\psi(T)\rangle$ with fixed $T=1$ for all
$m,\sigma$. Problem instances and visualization is the same as in Figure 24.
We also include the maximum eigenvalue of the relaxed Hamiltonian and the
energy of the prepared unrounded state, to demonstrate how far
$|\psi(T)\rangle$ is from the exact maximum eigenvector. (Left) Varying the
number of vertices $m$ in the graph. Note that the number of qubits required
here is $2m$. (Right) Varying the noise strength parameter $\sigma$.

Then in Figure 26 we plot the same 50 problem instances (per graph size/noise
level) as in Figure 24, but using the quasi-adibatically prepared state
$|\psi(T)\rangle$ where we have fixed $T=1$ for all graph sizes. The classical
SDP results are the same as in Figure 26, and for reference we include the
energy of the unrounded quantum state and the maximum eigenvalue of the
relaxed Hamiltonian (normalized with respect to the optimal objective value).
Qualitatively, we observe features similar to those seen in Figure 25. Namely,
although the annealing schedule is too fast to prepare a close approximation
to the maximum eigenstate, the rounded solutions (using the
$\operatorname{conv}\mathrm{SO}(n)$-based protocol) consistently have high
approximation ratios. Meanwhile, the vertex-rounded solutions are highly
inconsistent, which reflects the highly fluctuating behavior seen in Figure
25.

### Discussion and future work

In this chapter we have developed a quantum relaxation for a quadratic program
over orthogonal and rotation matrices, known as an instance of the little
noncommutative Grothendieck problem. The embedding of the classical objective
is achieved by recognizing an intimate connection between the geometric-
algebra construction of the orthogonal group and the structure of quantum
mechanics, in particular the formalism of fermions in second quantization.
From this perspective, the determinant condition of $\mathrm{SO}(n)$ is
succinctly captured by a simple linear property of the state—its parity—and
the convex bodies $\operatorname{conv}\mathrm{O}(n)$ and
$\operatorname{conv}\mathrm{SO}(n)$ (relevant to convex relaxations of
optimization over orthogonal matrices) are completely characterized by density
operators on $n$ and $n-1$ qubits, respectively. Recognizing that the reduced
state on each vertex therefore corresponds to an element of this convex hull,
we proposed vertex-marginal rounding which classically rounds the measured
one-body reduced density matrix of each vertex.

We additionally showed that these convex hulls are characterized by density
operators on $2n$ and $2(n-1)$ qubits as well, where the linear functionals
defining this PSD lift are the Hamiltonian terms appearing in our quantum
relaxation. This insight enables our second proposed rounding scheme,
$\operatorname{conv}G$-based edge rounding, which is inspired by the fact that
the a quantum Gram matrix $\mathcal{M}$ can be constructed from the
expectation values of the quantum state which obeys the same properties as the
classical SDP of Saunderson  _et al._ [saunderson2014semidefinite].
Numerically we observe that this approach to quantum rounding is significantly
more accurate and consistent than vertex rounding, and it consistently
achieves larger approximation ratios than the basic SDP relaxation. However,
we are severely limited by the exponential scaling of classically simulating
quantum states; further investigations would be valuable to ascertain the
empirical performance of these ideas at larger scales.

The primary goal of this work was to formulate the problem of orthogonal-
matrix optimization into a familiar quantum Hamiltonian problem, and to
establish the notion of a quantum relaxation for such optimization problems
over continuous-valued decision variables. A clear next step is to prove
nontrivial approximation ratios from our quantum relaxation. If such
approximation ratios exceed known guarantees by classical algorithms, for
example on certain types of graphs, then this would potentially provide a
quantum advantage for a class of applications not previously considered in the
quantum literature. We have proposed one standard, realistically preparable
class of states—quasi-adiabatic time evolution—but a variety of energy-
optimizing ansatze exist in the literature, especially considering that the
constructed Hamiltonian is an interacting-fermion model. From this
perspective, it would also be interesting to see if a classical many-body
method can produce states which round down to high-quality approximations,
even heuristically. Such an approach would constitute a potential example of a
quantum-inspired classical algorithm.

From a broader perspective, the quantum formalism described here may also
provide new insights into the computational hardness of the classical problem.
First, the $\mathsf{NP}$-hard thresholds for Problem (781) are not currently
known. However, by establishing the classical problem as an instance of
Gaussian product state optimization on the many-body Hamiltonian, it may be
possible to import tools from quantum computational complexity to study the
classical problem. This idea also applies to the more general instances of
noncommutative Grothendieck problems,

$\max_{U,V\in\mathrm{O}(N)}\sum_{i,j,k,l\in[N]}T_{ijkl}U_{ij}V_{kl},$ (892)

where the $N\times N\times N\times N$ tensor $T$ specifies the problem input.
It is straightforward to apply our quantum relaxation construction to this
problem, yielding a $2N$-qubit Hamiltonian whose terms are of the form
$P_{ij}\otimes P_{kl}$. While Briët _et al._ [briet2017tight] showed that the
$\mathsf{NP}$-hardness threshold of approximating this problem is $1/2$, it
remains an open problem to construct an algorithm which is guaranteed to
achieve this approximation ratio.

Although we have provided new approximation ratios for the instance of Problem
(781) over $\mathrm{SO}(n)$, it is unclear precisely how much harder the
$\mathrm{SO}(n)$ problem is compared to the $\mathrm{O}(n)$ problem. The work
by Saunderson _et al._ [saunderson2015semidefinite] establishes a clear
distinction between the representation sizes required for
$\operatorname{conv}\mathrm{O}(n)$ and $\operatorname{conv}\mathrm{SO}(n)$,
and this work has connected this structure to properties of quantum states on
$n$ qubits. However this does not yet establish a difference of hardness for
the corresponding quadratic programs. Again it would be interesting to see if
the tools of quantum information theory can be used to further understand this
classical problem. For example, one might study the $\mathsf{NP}$-hardness
threshold of Problem (892) where instead $U,V\in\mathrm{SO}(N)$ and leverage
the quantum (or equivalently, Clifford-algebraic) representation of
$\mathrm{SO}(N)$. In such a setting, the size of the problem is given by a
single parameter $N$ and so the exponentially large parametrization of
$\operatorname{conv}\mathrm{SO}(N)$ appears to signify a central difficulty of
this problem.

We note that it is straightforward to extend our quantum relaxation to the
unitary groups $\mathrm{U}(n)$ and $\mathrm{SU}(n)$, essentially by doubling
the number of qubits per vertex via the inclusions
$\mathrm{U}(n)\subset\mathrm{O}(2n)$ and
$\mathrm{SU}(n)\subset\mathrm{SO}(2n)$. However this is likely an inefficient
embedding, since the $n$-qubit Majorana operators already form a
representation of $\mathrm{Cl}(2n)$. It may therefore be possible to encode
complex-valued matrices via a complexification of $Q$, using the same amount
of quantum space. It is interesting to note that Briët _et al._
[briet2017tight] in fact utilize a “complex extension” of Clifford algebras
when considering Problem (892) over the unitary group, although the usage is
different from ours.

### Appendix

### Clifford algebras and the orthogonal group

In this appendix we review the key components for constructing the orthogonal
and special orthogonal groups from a Clifford algebra. Our presentation of
this material broadly follows Refs. [atiyah1964clifford,
saunderson2015semidefinite].

The Clifford algebra $\mathrm{Cl}(n)$ of $\mathbb{R}^{n}$ is a
$2^{n}$-dimensional real vector space, equipped with an inner product
$\langle\cdot,\cdot\rangle:\mathrm{Cl}(n)\times\mathrm{Cl}(n)\to\mathbb{R}$
and a multiplication operation satisfying the anticommutation relation

$e_{i}e_{j}+e_{j}e_{i}=-2\delta_{ij}\mathbb{I},$ (893)

where $\\{e_{1},\ldots,e_{n}\\}$ is an orthonormal basis of $\mathbb{R}^{n}$
and $\mathbb{I}$ is the multiplicative identity of the algebra. The basis
elements $e_{i}$ are called the generators of the Clifford algebra, in the
sense that they generate all other basis vectors of $\mathrm{Cl}(n)$ as

$e_{I}\coloneqq e_{i_{1}}\cdots e_{i_{k}},\quad
I=\\{i_{1},\ldots,i_{k}\\}\subseteq[n].$ (894)

By convention we order the indices $i_{1}<\cdots<i_{k}$, and the empty set
corresponds to the identity, $e_{\varnothing}=\mathbb{I}$. Taking all subsets
$I\subseteq[n]$ and extending the inner product definition from
$\mathbb{R}^{n}$ to $\mathrm{Cl}(n)$, it follows that $\\{e_{I}\mid
I\subseteq[n]\\}$ is an orthonormal basis with $2^{n}$ elements. Specifically,
we can write any element $x\in\mathrm{Cl}(n)$ as

$x=\sum_{I\subseteq[n]}x_{I}e_{I}$ (895)

with each $x_{I}\in\mathbb{R}$, and the inner product on $\mathrm{Cl}(n)$
is545454Equipping an inner product to the vector representation of
$\mathrm{Cl}(n)$ elements is achieved using the fact that algebra elements
square to a multiple of the identity.

$\langle x,y\rangle=\sum_{I\subseteq[n]}x_{I}y_{I},$ (896)

where $y=\sum_{I\subseteq[n]}y_{I}e_{I}$. Hence $\mathrm{Cl}(n)$ is isomorphic
as a Hilbert space to $\mathbb{R}^{2^{n}}$.

Now we show how to realize the orthogonal group $\mathrm{O}(n)$ from this
algebra. First observe the inclusion
$\mathbb{R}^{n}=\operatorname{span}\\{e_{i}\mid
i\in[n]\\}\subset\mathrm{Cl}(n)$. We shall identify the sphere
$S^{n-1}\subset\mathbb{R}^{n}$ as all $u\in\mathbb{R}^{n}$ satisfying $\langle
u,u\rangle=1$. We then define the Pin group as all possible products of
$S^{n-1}$ elements:

$\mathrm{Pin}(n)\coloneqq\\{u_{1}\cdots u_{k}\mid u_{1},\ldots,u_{k}\in
S^{n-1},0\leq k\leq n\\}.$ (897)

It is straightforward to check that this is indeed a group. Each
$x\in\mathrm{Pin}(n)$ is also normalized, $\langle x,x\rangle=1$. In fact, an
equivalent definition of this group is all elements $x\in\mathrm{Cl}(n)$
satisfying $x\overline{x}=\mathbb{I}$, where conjugation $\overline{x}$ is
defined from the linear extension of

$\overline{e_{I}}\coloneqq(-1)^{|I|}e_{i_{k}}\cdots e_{i_{1}}.$ (898)

The Pin group is a double cover of $\mathrm{O}(n)$, which can be seen from
defining a quadratic map $Q:\mathrm{Cl}(n)\to\mathbb{R}^{n\times n}$. This map
arises from the so-called twisted adjoint action, introduced by Atiyah _et
al._ [atiyah1964clifford]:555555Saunderson _et al._
[saunderson2015semidefinite] consider the standard adjoint action, which is
sufficient for describing rotations. However, the “twist” due to $\alpha$ is
necessary to construct arbitrary orthogonal transformations.

$v\mapsto\alpha(x)v\overline{x},\quad x,v\in\mathrm{Cl}(n),$ (899)

where the linear map $\alpha:\mathrm{Cl}(n)\to\mathrm{Cl}(n)$ is the parity
automorphism, defined by linearly extending

$\alpha(e_{I})\coloneqq(-1)^{|I|}e_{I}.$ (900)

Then for any $x\in\mathrm{Cl}(n)$, the linear map
$Q(x):\mathbb{R}^{n}\to\mathbb{R}^{n}$ is defined as

$Q(x)(v)\coloneqq\pi_{\mathbb{R}^{n}}(\alpha(x)v\overline{x})\quad\forall
v\in\mathbb{R}^{n},$ (901)

where $\pi_{\mathbb{R}^{n}}$ is the projection from $\mathrm{Cl}(n)$ onto
$\mathbb{R}^{n}$. To show that $Q(\mathrm{Pin}(n))=\mathrm{O}(n)$, it suffices
to recognize that, for any $u\in S^{n-1}$,
$\alpha(u)v\overline{u}\in\mathbb{R}^{n}$ is the reflection of the vector
$v\in\mathbb{R}^{n}$ across the hyperplane normal to $u$. To see this, first
observe that $uv+vu=-2\langle u,v\rangle\mathbb{I}$, which follows from Eq.
(893) by linearity. Then

$\begin{split}\alpha(u)v\overline{u}&=uvu\\\ &=(-vu-2\langle
u,v\rangle\mathbb{I})u\\\ &=v-2\langle u,v\rangle u,\end{split}$ (902)

which is precisely the elementary reflection as claimed. By the
Cartan–Dieudonné theorem, one can implement any orthogonal transformation on
$\mathbb{R}^{n}$ by composing $k\leq n$ such reflections about arbitrary
hyperplanes $u_{1},\ldots,u_{k}$ [gallier2011geometric]. This characterization
coincides precisely with the definition of the Pin group provided in Eq.
(897), through the composition of the linear maps $Q(u_{1}),\ldots,Q(u_{k})$
on $\mathbb{R}^{n}$. Hence for all $x\in\mathrm{Pin}(n)$, $Q(x)$ is an
orthogonal transformation on $\mathbb{R}^{n}$. The double cover property
follows from the fact that $Q$ is quadratic in $x$, so $Q(x)=Q(-x)$.

The special orthogonal group arises from the subgroup
$\mathrm{Spin}(n)\subset\mathrm{Pin}(n)$ containing only even-parity Clifford
elements. First observe that $\mathrm{Cl}(n)$ is a $\mathbb{Z}_{2}$-graded
algebra:

$\mathrm{Cl}(n)=\mathrm{Cl}^{0}(n)\oplus\mathrm{Cl}^{1}(n),$ (903)

where

$\displaystyle\mathrm{Cl}^{0}(n)$
$\displaystyle\coloneqq\operatorname{span}\\{e_{I}\mid|I|\text{ even}\\},$
(904) $\displaystyle\mathrm{Cl}^{1}(n)$
$\displaystyle\coloneqq\operatorname{span}\\{e_{I}\mid|I|\text{ odd}\\}.$
(905)

By a $\mathbb{Z}_{2}$ grading we mean that for each $x\in\mathrm{Cl}^{a}(n)$
and $y\in\mathrm{Cl}^{b}(n)$, their product $xy$ lies in
$\mathrm{Cl}^{a+b\>\mathrm{mod}\>2}(n)$. We say that elements in
$\mathrm{Cl}^{0}(n)$ (resp., $\mathrm{Cl}^{1}(n)$) have even (resp., odd)
parity. In particular, this grading implies that $\mathrm{Cl}^{0}(n)$ is a
subalgebra, hence its intersection with the Pin group is also a group, which
defines

$\mathrm{Spin}(n)\coloneqq\mathrm{Pin}(n)\cap\mathrm{Cl}^{0}(n)=\\{u_{1}\cdots
u_{2k}\mid u_{1},\ldots,u_{2k}\in S^{n-1},0\leq k\leq\lfloor n/2\rfloor\\}.$
(906)

Just as the Pin group double covers $\mathrm{O}(n)$, so does the Spin group
double cover $\mathrm{SO}(n)$. This is again a consequence of the
Cartan–Dieudonné theorem, wherein all rotations on $\mathbb{R}^{n}$ can be
decomposed into an even number of (at most $n$) arbitrary reflections.

### Convex hull of orthogonal matrices and quantum states

Recall the following characterizations of the convex hulls:

$\displaystyle\operatorname{conv}\mathrm{O}(n)$
$\displaystyle=\mathopen{}\left\\{X\in\mathbb{R}^{n\times
n}\mid\sigma_{1}(X)\leq 1\right\\}\mathclose{}$ (907)
$\displaystyle\operatorname{conv}\mathrm{SO}(n)$
$\displaystyle=\mathopen{}\left\\{X\in\mathbb{R}^{n\times
n}\mathrel{\Bigg{|}}\sum_{i\in[n]\setminus
I}\widetilde{\sigma}_{i}(X)-\sum_{i\in I}\widetilde{\sigma}_{i}(X)\leq
n-2\quad\forall I\subseteq[n],|I|\text{ odd}\right\\}\mathclose{},$ (908)

where $\\{\sigma_{i}(X)\\}_{i\in[n]}$ and
$\\{\widetilde{\sigma}_{i}(X\\}_{i\in[n]}$ are the singular values and special
singular values of $X$ in descending order, respectively. Note that
$\sigma_{i}(X)=\widetilde{\sigma}_{i}(X)$ for all $i\leq n-1$ and
$\sigma_{n}(X)=\operatorname{sign}(\det(X))\sigma_{n}(X)$.

#### PSD lift of $\operatorname{conv}\mathrm{O}(n)$ and
$\operatorname{conv}\mathrm{SO}(n)$

In this section we show that
$Q(\mathcal{D}(\mathrm{Cl}(n)))=\operatorname{conv}\mathrm{O}(n)$ and
$Q(\mathcal{D}(\mathrm{Cl}^{0}(n)))=\operatorname{conv}\mathrm{SO}(n)$, using
the quantum formalism described in Section Quantum formalism for optimization
over orthogonal matrices.

First we show that for all $X\in\operatorname{conv}\mathrm{O}(n)$, there
exists some $\rho\in\mathcal{D}(\mathcal{H}_{2^{n}})$ which generates $X$,
essentially by the convex extension of $Q$. Every
$X\in\operatorname{conv}\mathrm{O}(n)$ can be expressed as a convex
combination ($\sum_{\mu}p_{\mu}=1$, $p_{\mu}\geq 0$) of orthogonal matrices
$R_{\mu}\in\mathrm{O}(n)$:

$X=\sum_{\mu}p_{\mu}R_{\mu}.$ (909)

For each $R_{\mu}$ there exists some $x_{\mu}\in\mathrm{Pin}(n)$ such that
$[R_{\mu}]_{ij}=\langle x_{\mu}|P_{ij}|x_{\mu}\rangle$. Therefore the matrix
elements of $X$ can be expressed as

$X_{ij}=\operatorname{tr}\mathopen{}\left(P_{ij}\sum_{\mu}p_{\mu}|x_{\mu}\rangle\\!\langle
x_{\mu}|\right)\mathclose{}=\operatorname{tr}(P_{ij}\rho),$ (910)

where $\rho\coloneqq\sum_{\mu}p_{\mu}|x_{\mu}\rangle\\!\langle
x_{\mu}|\in\mathcal{D}(\mathcal{H}_{2^{n}})$.

Next we show the reverse direction, that for all
$\rho\in\mathcal{D}(\mathcal{H}_{2^{n}})$, the matrix
$X\coloneqq[\operatorname{tr}(P_{ij}\rho)]_{i,j\in[n]}$ is an element of
$\operatorname{conv}\mathrm{O}(n)$. Recall that
$X\in\operatorname{conv}\mathrm{O}(n)$ if and only if $\sigma_{1}(X)\leq 1$.
Therefore we take the singular value decomposition of $X=U\Sigma
V^{\mathsf{T}}$ and, using
$P_{ij}=\mathrm{i}\widetilde{\gamma}_{i}\gamma_{j}$, each singular value is
equal to

$\begin{split}\sigma_{k}(X)&=[U^{\mathsf{T}}XV]_{kk}\\\
&=\sum_{i,j\in[n]}U_{ik}\operatorname{tr}(\mathrm{i}\widetilde{\gamma}_{i}\gamma_{j}\rho)V_{jk}\\\
&=\operatorname{tr}(\mathrm{i}\mathcal{U}_{(U,\mathbb{I}_{n})}^{\dagger}\widetilde{\gamma}_{k}\mathcal{U}_{(U,\mathbb{I}_{n})}\mathcal{U}_{(\mathbb{I}_{n},V)}^{\dagger}\gamma_{k}\mathcal{U}_{(\mathbb{I}_{n},V)}\rho)\\\
&=\operatorname{tr}(\mathrm{i}\widetilde{\gamma}_{k}\gamma_{k}\rho^{\prime}),\end{split}$
(911)

where
$\rho^{\prime}\coloneqq\mathcal{U}_{(U,V)}\rho\mathcal{U}_{(U,V)}^{\dagger}$.
Because $\mathrm{i}\widetilde{\gamma}_{k}\gamma_{k}$ has eigenvalues $\pm 1$,
we see that $\sigma_{k}(X)\leq 1$ for all $k\in[n]$.

For the restriction to $\operatorname{conv}\mathrm{SO}(n)$, the first argument
is essentially the same. One merely replaces $P_{ij}$ with
$\widetilde{P}_{ij}$, hence $\sum_{\mu}p_{\mu}|x_{\mu}\rangle\\!\langle
x_{\mu}|\in\mathcal{D}(\mathcal{H}_{2^{n-1}})$. For the reverse direction, we
instead employ the special singular value decomposition which yields

$\widetilde{\sigma}_{k}(X)=\operatorname{tr}(\mathrm{i}\widetilde{\gamma}_{k}\gamma_{k}\rho^{\prime}),$
(912)

where now
$\rho^{\prime}\coloneqq\mathcal{U}_{(U,\widetilde{V})}\rho\mathcal{U}_{(U,\widetilde{V})}^{\dagger}\in\mathcal{D}(\mathcal{H}_{2^{n-1}})$.
Note that we have not projected to the even subspace this time, as it is more
convenient to work in the full $n$-qubit space when handling the Gaussian
unitaries. Instead, we will impose the constraint that $\rho$ only has support
on the even-parity subspace, so $\operatorname{tr}(Z^{\otimes n}\rho)=1$.
Furthermore, because the special singular value decomposition guarantees that
$\det(U)\det(\widetilde{V})=1$, $\mathcal{U}_{(U,\widetilde{V})}$ is parity
preserving so that $\operatorname{tr}(Z^{\otimes n}\rho^{\prime})=1$ as well.
Now recall that $X\in\operatorname{conv}\mathrm{SO}(n)$ if and only if

$\sum_{k\in[n]\setminus I}\widetilde{\sigma}_{k}(X)-\sum_{k\in
I}\widetilde{\sigma}_{k}(X)\leq n-2$ (913)

for all subsets $I\subseteq[n]$ of odd size. By linearity,

$\sum_{k\in[n]\setminus I}\widetilde{\sigma}_{k}(X)-\sum_{k\in
I}\widetilde{\sigma}_{k}(X)=\operatorname{tr}\mathopen{}\left(\rho^{\prime}\sum_{k\in[n]}(-1)^{z_{k}}Z_{k}\right)\mathclose{},$
(914)

where $z=z_{1}\cdots z_{n}\in\\{0,1\\}^{n}$ is defined as $z_{k}=1$ if $k\in
I$ and $z_{k}=0$ otherwise, and we have used the fact that
$\mathrm{i}\widetilde{\gamma}_{k}\gamma_{k}=Z_{k}$. It therefore suffices to
examine the spectrum of $A_{z}\coloneqq\sum_{k\in[n]}(-1)^{z_{k}}Z_{k}$:

$A_{z}|b\rangle=\mathopen{}\left(\sum_{k\in[n]}(-1)^{[z\oplus
b]_{k}}\right)\mathclose{}|b\rangle,\quad b\in\\{0,1\\}^{n},$ (915)

where $\oplus$ denotes addition modulo 2. As we are only interested in the
subspace spanned by even-parity states, we restrict attention to the
eigenvalues for which $|b|\>\mathrm{mod}\>2=0$. Because $|I|$ is odd, so too
is $|z|$, hence $|z\oplus b|\>\mathrm{mod}\>2=1$. This implies that there must
be at least one term in the sum of Eq. (915) which is negative, so it can only
take integer values at most $n-2$. This establishes Eq. (913), hence
$X\in\operatorname{conv}\mathrm{SO}(n)$.

#### Relation to $\operatorname{conv}\mathrm{SO}(n)$-based semidefinite
relaxation

Here we provide details for our claim that the relaxed quantum solution obeys
the same constraints as the classical SDP which uses the exponentially large
representation of $\operatorname{conv}\mathrm{SO}(n)$. Recall that this
relaxation can be formulated as

$\max_{M\in\mathbb{R}^{mn\times mn}}\sum_{(u,v)\in E}\langle
C_{uv},M_{uv}\rangle\quad\text{subject to }\begin{cases}M\succeq 0,\\\
M_{vv}=\mathbb{I}_{n}&\forall v\in[m],\\\
M_{uv}\in\operatorname{conv}\mathrm{SO}(n)&\forall u,v\in[m].\end{cases}$
(916)

We will show that the Gram matrix $\mathcal{M}$ constructed from the
measurements of a quantum state $\rho$, defined in Section Quantum Gram matrix
as

$[\mathcal{M}_{uv}]_{ij}=\begin{cases}\delta_{ij}&u=v,\\\
\frac{1}{n}\operatorname{tr}(\Gamma_{ij}^{(u,v)}\rho)&u<v,\\\
[M_{vu}]_{ji}&u>v,\end{cases}$ (917)

obeys the constraints of Eq. (916). Specifically, when the marginals of $\rho$
on each vertex are even-parity states (recall this is equivalent to replacing
$\Gamma_{ij}$ with $\widetilde{\Gamma}_{ij}$), we obtain the
$\operatorname{conv}\mathrm{SO}(n)$ condition, whereas when the parity of
$\rho$ is not fixed then $M_{uv}\in\operatorname{conv}\mathrm{O}(n)$.

First, we show that $\mathcal{M}$ is positive semidefinite for all quantum
states.

######  Lemma 0.33.

Let $\mathcal{M}\in\mathbb{R}^{mn\times mn}$ be defined as in Eq. (917). For
all $\rho\in\mathcal{D}(\mathcal{H}_{2^{n}}^{\otimes m})$, $\mathcal{M}\succeq
0$.

###### Proof.

We prove the statement by a sum-of-squares argument. To see where the fact of
$1/n$ appears in the quantum definition of $\mathcal{M}$ above, we first
construct a matrix $\mathcal{M}^{\prime}\succeq 0$ which turns out to simply
be $\mathcal{M}^{\prime}=n\mathcal{M}$.

For each $k\in[n]$ define the Hermitian operator

$A_{k}=\sum_{v\in V}\sum_{i\in[n]}c_{i}^{(v)}P_{ik}^{(v)},$ (918)

where $c_{i}^{(v)}\in\mathbb{R}$ are arbitrary coefficients. Consider its
square,

$\begin{split}A_{k}^{2}&=\mathopen{}\left(\sum_{v\in
V}\sum_{i\in[n]}c_{i}^{(v)}P_{ik}^{(v)}\right)^{2}\mathclose{}\\\ &=\sum_{v\in
V}\sum_{i,j\in[n]}c_{i}^{(v)}c_{j}^{(v)}P_{ik}^{(v)}P_{jk}^{(v)}+\sum_{\begin{subarray}{c}u,v\in
V\\\ u\neq
v\end{subarray}}\sum_{i,j\in[n]}c_{i}^{(u)}c_{j}^{(v)}P_{ik}^{(u)}\otimes
P_{jk}^{(v)}.\end{split}$ (919)

Note that the terms with $u\neq v$ feature the two-vertex operators as
desired, while the diagonal terms of the sum contain products of the Pauli
operators acting on the same vertex. Because
$P_{ik}=\mathrm{i}\widetilde{\gamma}_{i}\gamma_{k}$, the diagonal terms reduce
to (suppressing superscripts here)

$\begin{split}\sum_{i,j\in[n]}c_{i}c_{j}P_{ik}P_{jk}&=-\sum_{i,j\in[n]}c_{i}c_{j}\widetilde{\gamma}_{i}\gamma_{k}\widetilde{\gamma}_{j}\gamma_{k}\\\
&=\sum_{i,j\in[n]}c_{i}c_{j}\widetilde{\gamma}_{i}\widetilde{\gamma}_{j}\\\
&=\sum_{i\in[n]}c_{i}^{2}\mathbb{I}_{2^{n}}+\sum_{1\leq i<j\leq
n}c_{i}c_{j}(\widetilde{\gamma}_{i}\widetilde{\gamma}_{j}+\widetilde{\gamma}_{j}\widetilde{\gamma}_{i})\\\
&=\mathopen{}\left(\sum_{i\in[n]}c_{i}^{2}\right)\mathclose{}\mathbb{I}_{2^{n}}.\end{split}$
(920)

Plugging this result into Eq. (919) and summing over all $k\in[n]$, we obtain

$\begin{split}\sum_{k\in[n]}A_{k}^{2}&=\mathopen{}\left(\sum_{v\in
V}\sum_{i,k\in[n]}|c_{i}^{(v)}|^{2}\right)\mathclose{}\mathbb{I}_{2^{n}}+\sum_{\begin{subarray}{c}u,v\in
V\\\ u\neq
v\end{subarray}}\sum_{i,j,k\in[n]}c_{i}^{(u)}c_{j}^{(v)}P_{ik}^{(u)}\otimes
P_{jk}^{(v)}\\\ &=n\langle
c,c\rangle\mathbb{I}_{2^{n}}+\sum_{\begin{subarray}{c}u,v\in V\\\ u\neq
v\end{subarray}}\sum_{i,j\in[n]}c_{i}^{(u)}c_{j}^{(v)}\Gamma_{ij}^{(u,v)},\end{split}$
(921)

where we have collected the coefficients $c_{i}^{(v)}$ into a vector
$c\in\mathbb{R}^{mn}$. Similarly, if we arrange the expectation values
$\operatorname{tr}(\Gamma_{ij}^{(u,v)}\rho)$ into a matrix
$T\in\mathbb{R}^{mn\times mn}$ (where the $n\times n$ blocks on the diagonal
are $0$), then the expectation value of the sum-of-squares operator is

$\begin{split}\operatorname{tr}\mathopen{}\left(\sum_{k\in[n]}A_{k}^{2}\rho\right)\mathclose{}&=n\langle
c,c\rangle+\langle c,Tc\rangle\\\ &=\langle
c,\mathcal{M}^{\prime}c\rangle\end{split}$ (922)

where we have defined $\mathcal{M}^{\prime}\coloneqq n\mathbb{I}_{mn}+T$.
Because $\sum_{k\in[n]}A_{k}^{2}$ is a sum of PSD operators, its expectation
value is always nonnegative, hence $\langle c,\mathcal{M}^{\prime}c\rangle\geq
0$. This inequality holds for all vectors $c\in\mathbb{R}^{mn}$, so
$\mathcal{M}^{\prime}\succeq 0$ and hence
$\mathcal{M}=\mathcal{M}^{\prime}/n\succeq 0$ as claimed. ∎

The fact that the diagonal blocks $\mathcal{M}_{vv}=\mathbb{I}_{n}$ holds by
definition. Finally, we need to show that each block of $\mathcal{M}$ lies in
$\operatorname{conv}(\mathrm{SO}(n))$ when $\rho$ has even parity. A
straightforward corollary of this result is that the blocks lie in
$\operatorname{conv}\mathrm{O}(n)$ when $\rho$ does not have fixed parity.
Note that in Section Mixed states and the convex hull we showed that the
matrix of expectation values $\operatorname{tr}(P_{ij}\rho_{1})$ lies in
$\operatorname{conv}\mathrm{O}(n)$ for any $n$-qubit density operator
$\rho_{1}$, a straightforward extension of the PSD-lift representation of
$\operatorname{conv}\mathrm{SO}(n)$ presented in Ref.
[saunderson2015semidefinite]. Here we instead show that the matrix of
expectation values $\frac{1}{n}\operatorname{tr}(\Gamma_{ij}\rho_{2})$ for any
$2n$-qubit density operator $\rho_{2}$ also lies in
$\operatorname{conv}\mathrm{O}(n)$, and is an element of
$\operatorname{conv}\mathrm{SO}(n)$ when $\rho_{2}$ has support only in the
even-parity sector.

Because
$\mathbb{I}_{n}\in\operatorname{conv}\mathrm{SO}(n)\subset\operatorname{conv}\mathrm{O}(n)$,
and because these convex hulls are closed under transposition, all that
remains is to prove the statement for the blocks $\mathcal{M}_{uv}$ when
$u<v$. We show this by considering all density matrices on the reduced two-
vertex Hilbert space.

######  Lemma 0.34.

Let $\rho_{2}\in\mathcal{D}(\mathcal{H}_{2^{n}}^{\otimes 2})$. Define the
$n\times n$ matrix $T$ by

$T_{ij}\coloneqq\operatorname{tr}(\Gamma_{ij}\rho_{2})=\operatorname{tr}\mathopen{}\left(\sum_{k\in[n]}P_{ik}\otimes
P_{jk}\rho_{2}\right)\mathclose{},\quad i,j\in[n].$ (923)

Then $T/n\in\operatorname{conv}\mathrm{O}(n)$. Furthermore, if
$\operatorname{tr}(Z^{\otimes 2n}\rho_{2})=1$, then
$M_{uv}\in\operatorname{conv}\mathrm{SO}(n)$.

###### Proof.

Consider the special singular value decomposition of
$T=U\widetilde{\Sigma}\widetilde{V}^{\mathsf{T}}$. We can express the special
singular values as

$\begin{split}\widetilde{\sigma}_{i}(T)&=[U^{\mathsf{T}}T\widetilde{V}]_{ii}\\\
&=\sum_{j,\ell\in[n]}U_{ji}T_{j\ell}\widetilde{V}_{\ell i}\\\
&=\sum_{j,\ell\in[n]}U_{ji}\widetilde{V}_{\ell
i}\sum_{k\in[n]}\operatorname{tr}\mathopen{}\left(\mathrm{i}\widetilde{\gamma}_{j}\gamma_{k}\otimes\mathrm{i}\widetilde{\gamma}_{\ell}\gamma_{k}\rho_{2}\right)\mathclose{}\\\
&=\sum_{k\in[n]}\operatorname{tr}\mathopen{}\left({\mathcal{U}}_{(U,\mathbb{I}_{n})}^{\dagger}\mathrm{i}\widetilde{\gamma}_{i}\gamma_{k}{\mathcal{U}}_{(U,\mathbb{I}_{n})}\otimes{\mathcal{U}}_{(\widetilde{V},\mathbb{I}_{n})}^{\dagger}\mathrm{i}\widetilde{\gamma}_{i}\gamma_{k}{\mathcal{U}}_{(\widetilde{V},\mathbb{I}_{n})}\rho_{2}\right)\mathclose{}\\\
&=\sum_{k\in[n]}\operatorname{tr}\mathopen{}\left(\mathrm{i}\widetilde{\gamma}_{i}\gamma_{k}\otimes\mathrm{i}\widetilde{\gamma}_{i}\gamma_{k}\rho^{\prime}_{2}\right)\mathclose{},\end{split}$
(924)

where
$\rho^{\prime}_{2}\coloneqq\mathopen{}\left({\mathcal{U}}_{(U,\mathbb{I}_{n})}\otimes{\mathcal{U}}_{(\widetilde{V},\mathbb{I}_{n})}\right)\mathclose{}\rho_{2}\mathopen{}\left({\mathcal{U}}_{(U,\mathbb{I}_{n})}\otimes{\mathcal{U}}_{(\widetilde{V},\mathbb{I}_{n})}\right)^{\dagger}\mathclose{}$.
The fact that $T/n\in\operatorname{conv}\mathrm{O}(n)$ follows immediately
from the fact that the spectrum of
$\mathrm{i}\widetilde{\gamma}_{i}\gamma_{k}\otimes\mathrm{i}\widetilde{\gamma}_{i}\gamma_{k}$
is $\\{\pm 1\\}$:

$\begin{split}\sigma_{1}(T/n)&=\frac{1}{n}\mathopen{}\left|\widetilde{\sigma}_{1}(T)\right|\mathclose{}\\\
&=\frac{1}{n}\mathopen{}\left|\sum_{k\in[n]}\operatorname{tr}\mathopen{}\left(\mathrm{i}\widetilde{\gamma}_{1}\gamma_{k}\otimes\mathrm{i}\widetilde{\gamma}_{1}\gamma_{k}\rho^{\prime}_{2}\right)\mathclose{}\right|\mathclose{}\\\
&\leq\frac{1}{n}\sum_{k\in[n]}\mathopen{}\left|\operatorname{tr}\mathopen{}\left(\mathrm{i}\widetilde{\gamma}_{1}\gamma_{k}\otimes\mathrm{i}\widetilde{\gamma}_{1}\gamma_{k}\rho^{\prime}_{2}\right)\mathclose{}\right|\mathclose{}\\\
&\leq 1.\end{split}$ (925)

Now we examine the inclusion in $\operatorname{conv}\mathrm{SO}(n)$. Suppose
that $\rho_{2}$ is an even-parity state. By virtue of the special singular
value decomposition, we have that $\det(U\widetilde{V}^{\mathsf{T}})=1$ which
implies that the Gaussian unitary
${\mathcal{U}}_{(U,\mathbb{I}_{n})}\otimes{\mathcal{U}}_{(\widetilde{V},\mathbb{I}_{n})}$
preserves the parity of $\rho_{2}$: $\operatorname{tr}(Z^{\otimes
2n}\rho_{2})=\operatorname{tr}(Z^{\otimes 2n}\rho_{2}^{\prime})=1$. For $T/n$
to lie in $\operatorname{conv}\mathrm{SO}(n)$, the following inequality from
Eq. (908) must hold:

$\sum_{i\in[n]\setminus I}\widetilde{\sigma}_{i}(T)-\sum_{i\in
I}\widetilde{\sigma}_{i}(T)\leq n(n-2)$ (926)

for all subsets $I\subseteq[n]$ of odd size. To show this, first we write the
left-hand side in terms of the result derived from Eq. (924):

$\begin{split}\sum_{i\in[n]\setminus I}\widetilde{\sigma}_{i}(T)-\sum_{i\in
I}\widetilde{\sigma}_{i}(T)&=\sum_{i\in[n]\setminus
I}\sum_{k\in[n]}\operatorname{tr}\mathopen{}\left(\mathrm{i}\widetilde{\gamma}_{i}\gamma_{k}\otimes\mathrm{i}\widetilde{\gamma}_{i}\gamma_{k}\rho^{\prime}_{2}\right)\mathclose{}-\sum_{i\in
I}\sum_{k\in[n]}\operatorname{tr}\mathopen{}\left(\mathrm{i}\widetilde{\gamma}_{i}\gamma_{k}\otimes\mathrm{i}\widetilde{\gamma}_{i}\gamma_{k}\rho^{\prime}_{2}\right)\mathclose{}\\\
&=\operatorname{tr}\mathopen{}\left(\rho_{2}^{\prime}\sum_{i,k\in[n]}(-1)^{z_{i}}\mathrm{i}\widetilde{\gamma}_{i}\gamma_{k}\otimes\mathrm{i}\widetilde{\gamma}_{i}\gamma_{k}\right)\mathclose{},\end{split}$
(927)

where the string $z=z_{1}\cdots z_{n}\in\\{0,1\\}^{n}$ is defined as

$z_{i}=\begin{cases}1&i\in I,\\\ 0&i\notin I.\end{cases}$ (928)

Note that $|I|$ being odd implies that the Hamming weight of $z$ is also odd.
To bound Eq. (927) we shall seek a bound on the largest eigenvalue of the
Hermitian operator

$\sum_{i,k\in[n]}(-1)^{z_{i}}\mathrm{i}\widetilde{\gamma}_{i}\gamma_{k}\otimes\mathrm{i}\widetilde{\gamma}_{i}\gamma_{k}=\mathopen{}\left(\sum_{i\in[n]}(-1)^{z_{i}\oplus
1}\widetilde{\gamma}_{i}\otimes\widetilde{\gamma}_{i}\right)\mathclose{}\mathopen{}\left(\sum_{k\in[n]}\gamma_{k}\otimes\gamma_{k}\right)\mathclose{}\eqqcolon
A_{z}B$ (929)

over the space of even-parity states. Here we use $\oplus$ to denote addition
modulo 2.

The operator $A_{z}B$ can in fact be exactly diagonalized in the basis of Bell
states. First, observe that $[A_{z},B]=0$, which follows from the fact that
$\widetilde{\gamma}_{i}\gamma_{k}=-\gamma_{k}\widetilde{\gamma}_{i}$, hence
$[\widetilde{\gamma}_{i}\otimes\widetilde{\gamma}_{i},\gamma_{k}\otimes\gamma_{k}]=0$
for all $i,k\in[n]$. We can therefore seek their simultaneous eigenvectors,
which can be determined from looking at the Jordan–Wigner representation of
the Majorana operators:

$\displaystyle\widetilde{\gamma}_{i}\otimes\widetilde{\gamma}_{i}$
$\displaystyle=(Z_{1}\cdots Z_{i-1}Y_{i})\otimes(Z_{1}\cdots Z_{i-1}Y_{i}),$
(930) $\displaystyle\gamma_{k}\otimes\gamma_{k}$ $\displaystyle=(Z_{1}\cdots
Z_{k-1}X_{k})\otimes(Z_{1}\cdots Z_{k-1}X_{k}).$ (931)

These operators are diagonalized by the $2n$-qubit state

$|\beta(x,y)\rangle=\bigotimes_{j=1}^{n}|\beta(x_{j},y_{j})\rangle,$ (932)

where $x,y\in\\{0,1\\}^{n}$, and the Bell state $|\beta(x_{j},y_{j})\rangle$
between the $j$th qubits across the two subsystems is defined as

$|\beta(x_{j},y_{j})\rangle\coloneqq\frac{|0\rangle\otimes|y_{j}\rangle+(-1)^{x_{j}}|1\rangle\otimes|y_{j}\oplus
1\rangle}{\sqrt{2}}.$ (933)

Indeed there are $2^{2n}$ such states $|\beta(x,y)\rangle$, so they form an
orthonormal basis for the $2n$ qubits. The eigenvalues of Eqs. (930) and (931)
can be determined by a standard computation,

$\displaystyle(X_{j}\otimes X_{j})|\beta(x_{j},y_{j})\rangle$
$\displaystyle=(-1)^{x_{j}}|\beta(x_{j},y_{j})\rangle,$ (934)
$\displaystyle(Y_{j}\otimes Y_{j})|\beta(x_{j},y_{j})\rangle$
$\displaystyle=(-1)^{x_{j}\oplus y_{j}\oplus 1}|\beta(x_{j},y_{j})\rangle,$
(935) $\displaystyle(Z_{j}\otimes Z_{j})|\beta(x_{j},y_{j})\rangle$
$\displaystyle=(-1)^{y_{j}}|\beta(x_{j},y_{j})\rangle.$ (936)

Taking the appropriate products furnishes the eigenvalues of the $B$ and
$A_{z}$ as

$\displaystyle\begin{split}B|\beta(x,y)\rangle&=\sum_{k\in[n]}(Z_{1}\otimes
Z_{1})\cdots(Z_{k-1}\otimes Z_{k-1})(X_{k}\otimes X_{k})|\beta(x,y)\rangle\\\
&=\sum_{k\in[n]}(-1)^{y_{1}\oplus\cdots\oplus
y_{k-1}}(-1)^{x_{k}}|\beta(x,y)\rangle,\end{split}$ (937)
$\displaystyle\begin{split}A_{z}|\beta(x,y)\rangle&=\sum_{i\in[n]}(-1)^{z_{i}\oplus
1}(Z_{1}\otimes Z_{1})\cdots(Z_{k-1}\otimes Z_{k-1})(Y_{k}\otimes
Y_{k})|\beta(x,y)\rangle\\\ &=\sum_{i\in[n]}(-1)^{z_{i}\oplus
1}(-1)^{y_{1}\oplus\cdots\oplus y_{i-1}}(-1)^{x_{i}\oplus y_{i}\oplus
1}|\beta(x,y)\rangle\\\
&=\sum_{i\in[n]}(-1)^{z_{i}}(-1)^{y_{1}\oplus\cdots\oplus
y_{i}}(-1)^{x_{i}}|\beta(x,y)\rangle\end{split}$ (938)

Altogether we arrive at the expression for the eigenvalues of $A_{z}B$,

$\langle\beta(x,y)|A_{z}B|\beta(x,y)\rangle=\mathopen{}\left(\sum_{i\in[n]}(-1)^{z_{i}}(-1)^{y_{1}\oplus\cdots\oplus
y_{i}}(-1)^{x_{i}}\right)\mathclose{}\mathopen{}\left(\sum_{k\in[n]}(-1)^{y_{1}\oplus\cdots\oplus
y_{k-1}}(-1)^{x_{k}}\right)\mathclose{}.$ (939)

We wish to find the largest value this can take over even-parity states.
First, observe that the eigenstates $|\beta(x,y)\rangle$ have fixed parity
according to

$\langle\beta(x,y)|Z^{\otimes 2n}|\beta(x,y)\rangle=(-1)^{|y|},$ (940)

which follows from Eq. (936). Hence we shall only consider $y$ to have even
Hamming weight. Additionally recall that $z$ has odd Hamming weight, while the
Hamming weight of $x$ is unrestricted.

Let us denote the sums in Eq. (939) by

$\displaystyle a_{z}(x,y)$
$\displaystyle\coloneqq\sum_{i\in[n]}(-1)^{z_{i}}(-1)^{y_{1}\oplus\cdots\oplus
y_{i-1}\oplus y_{i}}(-1)^{x_{i}},$ (941) $\displaystyle b(x,y)$
$\displaystyle\coloneqq\sum_{k\in[n]}(-1)^{y_{1}\oplus\cdots\oplus
y_{k-1}}(-1)^{x_{k}}.$ (942)

Clearly, they can be at most $n$, and this occurs whenever all the terms in
their sum are positive. For $b(x,y)$, this is possible if and only if
$x=p(y)$, where $p:\\{0,1\\}^{n}\to\\{0,1\\}^{n}$ stores the parity
information of the $(k-1)$-length substring of its input into the $k$th bit of
its output:

$[p(y)]_{k}\coloneqq y_{1}\oplus\cdots\oplus y_{k-1}.$ (943)

For notation clarity we point out that $[p(y)]_{1}=0$ and $[p(y)]_{2}=y_{1}$.
The forward direction, $b(p(y),y)=n$, is clear by construction. The reverse
direction, that $b(x,y)=n$ implies $x=p(y)$, follows from the bijectivity of
modular addition.

Plugging this value of $x=p(y)$ into $a_{z}(x,y)$ yields

$a_{z}(p(y),y)=\sum_{i\in[n]}(-1)^{z_{i}}(-1)^{y_{i}}=\sum_{i\in[n]}(-1)^{[z\oplus
y]_{i}}.$ (944)

Because $y$ has even and $z$ has odd Hamming weight, their sum $y\oplus z$
must have odd Hamming weight. Therefore at least one term in Eq. (944) must be
negative, implying that $a_{z}(p(y),y)\leq n-2$. It follows that

$\langle\beta(p(y),y)|A_{z}B|\beta(p(y),y)\rangle=a_{z}(p(y),y)b(p(y),y)\leq(n-2)n.$
(945)

We now show that no other assignment of $(x,y)$ can exceed this bound. Recall
that $b(x,y)=n$ if and only if $x=p(y)$. Thus any other choice of $x$
necessarily returns a smaller value of $b(x,y)$. Because sums of $\pm 1$
cannot yield $n-1$, the next largest value would be $b(x,y)=n-2$. However we
can always trivially bound $a_{z}(x,y)\leq n$ for all $x,y$. This implies that
such a choice of $x$ for which $b(x,y)=n-2$ (whatever it is) also cannot
provide a value of $a_{z}(x,y)b(x,y)$ exceeding $n(n-2)$.

Note that there is another assignment that saturates the upper bound, which is
simply considering a global negative sign in front of both products.
Specifically, let $x=p(y)\oplus 1^{n}$. In this case $b(x,y)=-n$ and
$a_{z}(x,y)\geq-(n-2)$, yielding the same bound $a_{z}(x,y)b(x,y)\leq n(n-2)$.

To conclude, we use the fact that the maximum eigenvalue of $A_{z}B$ in the
even-parity subspace is $n(n-2)$ for any odd-weight $z$ (equiv., any odd-size
$I\subseteq[n]$). This bounds the value of Eq. (927) by $n(n-2)$, hence
validating the inequality of Eq. (926). Thus
$T/n\in\operatorname{conv}\mathrm{SO}(n)$ whenever $\rho_{2}$ has even parity.
∎

As usual, replacing the operators $\Gamma_{ij}$ with $\widetilde{\Gamma}_{ij}$
is equivalent to enforcing the even-parity constraint. In fact, since
$\widetilde{\Gamma}_{ij}=\sum_{k\in[n]}\widetilde{P}_{ik}\otimes\widetilde{P}_{jk}$,
the equivalent constraint involves the reduced single-vertex marginals,
$\operatorname{tr}(Z^{\otimes
n}\otimes\mathbb{I}_{2^{n}}\rho_{2})=\operatorname{tr}(\mathbb{I}_{2^{n}}\otimes
Z^{\otimes n}\rho_{2})=1$, rather than the entire $2n$-qubit Hilbert space. Of
course, if both single-vertex parity constraints are satisfied, then the two-
vertex constraint automatically follows.

### Details for working in the even-parity subspace

First we provide an expression for $n$-qubit Pauli operators projected to
$\mathrm{Cl}^{0}(n)$. Let

$A\coloneqq W_{1}\otimes\cdots\otimes W_{n},\quad
W_{i}\in\\{\mathbb{I},X,Y,Z\\}.$ (946)

A straightforward calculation yields the conditional expression:

$\widetilde{A}=\Pi_{0}A\Pi_{0}^{\mathsf{T}}=\begin{cases}0&\text{if
}[A,Z^{\otimes n}]\neq 0\\\ \begin{cases}W_{2}\otimes\cdots\otimes
W_{n}&\text{if }W_{1}=\mathbb{I},X\\\
\mathrm{i}(W_{2}Z)\otimes\cdots\otimes(W_{n}Z)&\text{if }W_{1}=Y\\\
(W_{2}Z)\otimes\cdots\otimes(W_{n}Z)&\text{if }W_{1}=Z.\end{cases}&\text{if
}[A,Z^{\otimes n}]=0.\end{cases}$ (947)

Notably, if $A$ does not commute with the parity operator then
$\widetilde{A}=0$.

Now we generalize from the main text, defining the operator

$\Pi_{k}\coloneqq\frac{1}{\sqrt{2}}\mathopen{}\left(\langle+|\otimes\mathbb{I}_{2}^{\otimes(n-1)}+(-1)^{k}\langle-|\otimes
Z^{\otimes(n-1)}\right)\mathclose{},\quad k\in\\{0,1\\},$ (948)

which is the projector $\mathrm{Cl}(n)\to\mathrm{Cl}^{k}(n)$. These operators
obey

$\displaystyle\Pi_{k}\Pi_{k}^{\mathsf{T}}$
$\displaystyle=\mathbb{I}_{2^{n-1}},$ (949)
$\displaystyle\Pi_{k}^{\mathsf{T}}\Pi_{k}$
$\displaystyle=\frac{\mathbb{I}_{2^{n}}+(-1)^{k}Z^{\otimes n}}{2}.$ (950)

Given a state $\widetilde{\rho}\in\mathcal{D}(\mathcal{H}_{2^{n-1}})$, the
expectation values of $\widetilde{P}_{ij}$ satisfy

$\begin{split}\operatorname{tr}(\widetilde{P}_{ij}\widetilde{\rho})&=\operatorname{tr}(P_{ij}\Pi_{0}^{\mathsf{T}}\widetilde{\rho}\Pi_{0})\\\
&=\operatorname{tr}(P_{ij}\rho_{0}),\end{split}$ (951)

where
$\rho_{0}\coloneqq\Pi_{0}^{\mathsf{T}}\widetilde{\rho}\Pi_{0}\in\mathcal{D}(\mathcal{H}_{2^{n}})$
has only support on even-parity computational basis states. By Eq. (949) we
can “invert” this relation in the following sense: given a state
$\rho\in\mathcal{D}(\mathcal{H}_{2^{n}})$ which only has support on the even
subspace, its $(n-1)$-qubit representation is
$\Pi_{0}\rho\Pi_{0}^{\mathsf{T}}\in\mathcal{D}(\mathcal{H}_{2^{n-1}})$.

This translation is useful when using the fermionic interpretation of the
$P_{ij}$ but we wish to work directly in the $(n-1)$-qubit subspace. For
example, suppose we wish to prepare a product of Gaussian states in a quantum
computer (as is done in Section Quasi-adiabatic state preparation to prepare
the initial state of the quasi-adiabatic evolution). It is well known how to
compile linear-depth circuits for this task [jiang2018quantum], however this
is within the standard $n$-qubit representation. Here we show how to translate
those circuits into the $(n-1)$-qubit representation under $\Pi_{0}$. In fact,
these techniques apply to any sequence of gates which commute with the parity
operator $Z^{\otimes n}$.

Let $|\psi\rangle\coloneqq V|0^{n}\rangle\in\mathcal{H}_{2^{n}}$, where the
circuit $V$ is constructed from $L$ gates, $V=V_{L}\cdots V_{1}$. We can
assume without loss of generality that the initial state is the vacuum
$|0^{n}\rangle$ and that all gates $V_{\ell}$ commute with the parity
operator, as otherwise $|\psi\rangle$ would not lie in the even subspace of
$\mathcal{H}_{2^{n}}$.565656While it is possible to have an even number of
gates which anticommute with the parity operator, for simplicity we assume
that the circuit has been compiled such that each $V_{\ell}$ preserves parity.
For example, parity-preserving Gaussian unitaries can be decomposed into
single- and two-qubit gates of the form $e^{-\mathrm{i}\theta Z_{i}}$,
$e^{-\mathrm{i}\theta X_{i}X_{i+1}}$. We wish to obtain a circuit description
for preparing the state $\Pi_{0}|\psi\rangle\in\mathcal{H}_{2^{n-1}}$. By Eq.
(950), $\Pi_{0}^{\mathsf{T}}\Pi_{0}|0^{n}\rangle=|0^{n}\rangle$, and
furthermore $\Pi_{0}|0^{n}\rangle=|0^{n-1}\rangle$ is the $(n-1)$-qubit
representation of the vacuum. Therefore

$\begin{split}\Pi_{0}|\psi\rangle&=\Pi_{0}V|0^{n}\rangle\\\
&=\Pi_{0}V_{L}\cdots V_{1}\Pi_{0}^{\mathsf{T}}\Pi_{0}|0^{n}\rangle\\\
&=\Pi_{0}V_{L}\cdots V_{1}\Pi_{0}^{\mathsf{T}}|0^{n-1}\rangle.\end{split}$
(952)

Because $\Pi_{0}$ is not unitary ($\Pi_{0}^{\mathsf{T}}$ is merely an
isometry), we cannot simply insert terms like $\Pi_{0}^{\mathsf{T}}\Pi_{0}$ in
between each gate. However, observe that if each $V_{\ell}$ preserves parity,
then they can be block diagonalized into the even and odd subspaces of
$\mathcal{H}_{2^{n}}$,

$V_{\ell}=\begin{bmatrix}V_{\ell,0}&0\\\ 0&V_{\ell,1}\end{bmatrix},$ (953)

where $V_{\ell,k}\coloneqq\Pi_{k}V_{\ell}\Pi_{k}^{\mathsf{T}}$ are
$2^{n-1}$-dimensional unitary matrices. Thus

$V=\begin{bmatrix}V_{L,0}\cdots V_{1,0}&0\\\ 0&V_{L,1}\cdots
V_{1,1}\end{bmatrix},$ (954)

and conjugation by the projector $\Pi_{0}$ precisely extracts the first block
of this matrix:

$\Pi_{0}V\Pi_{0}^{\mathsf{T}}=V_{L,0}\cdots V_{1,0}.$ (955)

This sequence of gates is what we wish to implement on the physical
$(n-1)$-qubit register. When the gates $V_{\ell}$ take the form

$V_{\ell}=e^{-\mathrm{i}\theta A_{\ell}}$ (956)

for some $n$-qubit Pauli operator $A_{\ell}$, then

$\begin{split}V_{\ell,0}&=\Pi_{0}V_{\ell}\Pi_{0}^{\mathsf{T}}\\\
&=\Pi_{0}\mathopen{}\left(\mathbb{I}_{2^{n}}\cos\theta-\mathrm{i}A_{\ell}\sin\theta\right)\mathclose{}\Pi_{0}^{\mathsf{T}}\\\
&=\mathbb{I}_{2^{n-1}}\cos\theta-\mathrm{i}\widetilde{A}_{\ell}\sin\theta\\\
&=e^{-\mathrm{i}\theta\widetilde{A}_{\ell}},\end{split}$ (957)

where $\widetilde{A}_{\ell}\coloneqq\Pi_{0}A_{\ell}\Pi_{0}^{\mathsf{T}}$. Note
that this calculation assumes that $V_{\ell}$ commutes with parity, hence
$[A_{\ell},Z^{\otimes n}]=0$, which guarantees that $\widetilde{A}_{\ell}$ is
unitary and Hermitian and hences furnishes the final line of Eq. (957). (If
they did not commute then the decomposition of Eq. (953) would not be valid to
begin with.)

### Small $n$ examples

In this appendix we write down the Pauli operators $P_{ij}$ and
$\widetilde{P}_{ij}$ for small but relevant values of $n$, to provide the
reader with some concerte examples of how to construct the LNCG Hamiltonian.

#### The Ising model from the $\mathrm{O}(1)$ setting

As a warm-up we first demonstrate that the LNCG Hamiltonian reduces to the
Ising formulation of the commutative combinatorial optimization problem when
considering $G=\mathrm{O}(1)$. Each local Hilbert space $\mathcal{H}_{d}$ with
$d=2$ is simply a qubit, and we only have to consider a single Pauli operator
on each local qubit,

$P_{11}=\mathrm{i}\widetilde{\gamma}_{1}\gamma_{1}=Z.$ (958)

It then follows that the LNCG interaction terms are

$\Gamma_{11}^{(u,v)}=Z^{(u)}\otimes Z^{(v)},$ (959)

and so the full Hamiltonian acting on $\mathcal{H}_{2}^{\otimes|V|}$ is indeed
the classical Ising Hamiltonian,

$H=\sum_{(u,v)\in E}C_{uv}Z^{(u)}\otimes Z^{(v)},$ (960)

with weights $C_{uv}\in\mathbb{R}$. It is also instructive to write down the
elements of $\mathrm{Pin}(1)$ as quantum states. The Clifford algebra
$\mathrm{Cl}(1)$ is spanned by $e_{1}$ and $e_{\varnothing}=\mathbb{I}$, with
the only elements of $S^{0}$ being $\pm e_{1}$. Therefore, taking all possible
products of elements in $S^{0}$ (including the empty product), we arrive at

$\mathrm{Pin}(1)=\\{\pm e_{\varnothing},\pm e_{1}\\},$ (961)

which corresponds to the qubit computational basis states
$\\{|0\rangle,|1\rangle\\}$ (up to global phases), as expected. Indeed, one
sees that the mappings $Q(e_{\varnothing})=\langle 0|Z|0\rangle=1$ and
$Q(e_{1})=\langle 1|Z|1\rangle=-1$ fully cover $\mathrm{O}(1)$.

#### The projected operators of the $\mathrm{SO}(3)$ setting

As $\mathrm{SO}(3)$ is arguably the most ubiquitous group for physical
applications, for reference we explicitly write down its Pauli operators under
the projection to $\mathrm{Cl}^{0}(3)\cong\mathcal{H}_{4}$. As seen by the
dimension of this Hilbert space, only two qubits per variable (for a total of
$2m$ qubits) are required to represent this problem. Using Eq. (947) we have

$\begin{bmatrix}\widetilde{P}_{11}&\widetilde{P}_{12}&\widetilde{P}_{13}\\\
\widetilde{P}_{21}&\widetilde{P}_{22}&\widetilde{P}_{23}\\\
\widetilde{P}_{31}&\widetilde{P}_{32}&\widetilde{P}_{33}\end{bmatrix}=\begin{bmatrix}Z_{1}Z_{2}&-X_{1}&-Z_{1}X_{2}\\\
X_{1}Z_{2}&Z_{1}&-X_{1}X_{2}\\\ X_{2}&-Y_{1}Y_{2}&Z_{2}\end{bmatrix}.$ (962)

### Classical approximation ratio for $\mathrm{SO}(n)$

Approximation ratios for the rounded solution of the classical semidefinite
relaxation of Problem (781) were obtained in Ref. [bandeira2016approximating]
for the cases $G=\mathrm{O}(n)$ and $\mathrm{U}(n)$. However, no such
approximation ratios were derived for the case of $\mathrm{SO}(n)$. Here we
adapt the argument of Ref. [bandeira2016approximating] to this setting,
wherein the rounding algorithm performs the special singular value
decomposition to guarantee that the rounded solutions have unit determinant.
As we will see, this feature results in a approximation ratio for the
classical semidefinite program over $\mathrm{SO}(n)$ that is strictly worse
than the previously studied $\mathrm{O}(n)$ case.

Recall that the semidefinite relaxation of Problem (781) can be formulated as

$\max_{X_{1},\ldots,X_{m}\in\mathbb{R}^{n\times mn}}\sum_{(u,v)\in E}\langle
C_{uv},X_{u}X_{v}^{\mathsf{T}}\rangle\quad\text{subject to}\quad
X_{v}X_{v}^{\mathsf{T}}=\mathbb{I}_{n}.$ (963)

To round the relaxed solution back into the feasible space of orthogonal
matrices, Ref. [bandeira2016approximating] proposes the following randomized
algorithm with a guarantee on the approximation ratio.

######  Theorem 0.35 ([bandeira2016approximating, Theorem 4]).

Let $X_{1},\ldots,X_{n}\in\mathbb{R}^{n\times mn}$ be a solution to Problem
(963). Let $Z$ be an $mn\times n$ Gaussian random matrix whose entries are
drawn i.i.d. from $\mathcal{N}(0,1/n)$. Compute the orthogonal matrices

$Q_{v}=\mathcal{P}(X_{v}Z),$ (964)

where
$\mathcal{P}(X)=\operatorname*{arg\,min}_{Y\in\mathrm{O}(n)}\|Y-X\|_{F}$. The
expected value of this approximate solution (averaged over $Z$) obeys

$\operatorname*{\mathbb{E}}\mathopen{}\left[f(Q_{1},\ldots,Q_{m})\right]\mathclose{}\geq\alpha_{\mathrm{O}(n)}^{2}\max_{R_{1},\ldots,R_{m}\in\mathrm{O}(n)}f(R_{1},\ldots,R_{m}).$
(965)

The approximation ratio $\alpha_{\mathrm{O}(n)}^{2}$ is defined by the average
singular value of random Gaussian $n\times n$ matrices
$Z_{1}\sim\mathcal{N}(0,\mathbb{I}_{n}/n)$,

$\alpha_{\mathrm{O}(n)}\coloneqq\operatorname*{\mathbb{E}}\mathopen{}\left[\frac{1}{n}\sum_{i\in[n]}\sigma_{i}(Z_{1})\right]\mathclose{},$
(966)

where $\sigma_{i}(Z_{1})$ is the $i$th singular value of $Z_{1}$.

Our adaptation to the $\mathrm{SO}(n)$ setting simply replaces the rounding
operator $\mathcal{P}$ with
$\widetilde{\mathcal{P}}(X)\coloneqq\operatorname*{arg\,min}_{Y\in\mathrm{SO}(n)}\|Y-X\|_{F}$.
With this change we obtain an analogous result for optimizing over
$\mathrm{SO}(n)$ elements with the same classical semidefinite program:

######  Theorem 0.36.

Let $X_{1},\ldots,X_{n}$ and $Z$ be as in Theorem 0.35. Compute the rotation
matrices

$Q_{v}=\widetilde{\mathcal{P}}(X_{v}Z),$ (967)

where
$\widetilde{\mathcal{P}}(X)=\operatorname*{arg\,min}_{Y\in\mathrm{SO}(n)}\|Y-X\|_{F}$.
The expected value of this approximate solution (averaged over $Z$) obeys

$\operatorname*{\mathbb{E}}\mathopen{}\left[f(Q_{1},\ldots,Q_{m})\right]\mathclose{}\geq\alpha_{\mathrm{SO}(n)}^{2}\max_{R_{1},\ldots,R_{m}\in\mathrm{SO}(n)}f(R_{1},\ldots,R_{m}).$
(968)

The approximation ratio $\alpha_{\mathrm{SO}(n)}^{2}$ is defined by the
average $n-1$ largest singular values of random Gaussian $n\times n$ matrices
$Z_{1}\sim\mathcal{N}(0,\mathbb{I}_{n}/n)$,

$\alpha_{\mathrm{SO}(n)}\coloneqq\operatorname*{\mathbb{E}}\mathopen{}\left[\frac{1}{n}\sum_{i\in[n-1]}\sigma_{i}(Z_{1})\right]\mathclose{},$
(969)

where $\sigma_{i}(Z_{1})$ is the $i$th singular value of $Z_{1}$, in
descending order $\sigma_{1}(Z_{1})\geq\cdots\geq\sigma_{n}(Z_{1})\geq 0$.

Because singular values are nonnegative, it is clear that

$\alpha_{\mathrm{O}(n)}-\alpha_{\mathrm{SO}(n)}=\frac{1}{n}\operatorname*{\mathbb{E}}\mathopen{}\left[\sigma_{n}(Z_{1})\right]\mathclose{}\geq
0.$ (970)

In particular we will see that
$\operatorname*{\mathbb{E}}\mathopen{}\left[\sigma_{n}(Z_{1})\right]\mathclose{}>0$
for all finite $n$, so the rounding algorithm guarantees a strictly smaller
approximation ratio for the problem over $\mathrm{SO}(n)$ than over
$\mathrm{O}(n)$.

The proof of Theorem 0.35 requires two lemmas regarding the expected value of
random Gaussian matrices under the rounding operator $\mathcal{P}$.
Analogously, our proof of Theorem 0.36 requires a modification of those lemmas
when $\mathcal{P}$ is replaced by $\widetilde{\mathcal{P}}$.

######  Lemma 0.37 (Adapted from [bandeira2016approximating, Lemma 5]).

Let $M,N\in\mathbb{R}^{n\times mn}$ obey
$MM^{\mathsf{T}}=NN^{\mathsf{T}}=\mathbb{I}_{n}$. For
$Z\in\mathbb{R}^{mn\times n}$ with i.i.d. entries drawn from
$\mathcal{N}(0,n^{-1})$, we have

$\operatorname*{\mathbb{E}}\mathopen{}\left[\widetilde{\mathcal{P}}(MZ)(NZ)^{\mathsf{T}}\right]\mathclose{}=\operatorname*{\mathbb{E}}\mathopen{}\left[(MZ)\widetilde{\mathcal{P}}(NZ)^{\mathsf{T}}\right]\mathclose{}=\alpha_{\mathrm{SO}(n)}MN^{\mathsf{T}}.$
(971)

This lemma is proved with the help of the following lemma.

######  Lemma 0.38 (Adapted from [bandeira2016approximating, Lemma 6]).

Let $Z_{1}\in\mathbb{R}^{n\times n}$ with i.i.d. entries drawn from
$\mathcal{N}(0,1/n)$. Then

$\operatorname*{\mathbb{E}}\mathopen{}\left[\widetilde{\mathcal{P}}(Z_{1})Z_{1}^{\mathsf{T}}\right]\mathclose{}=\operatorname*{\mathbb{E}}\mathopen{}\left[Z_{1}\widetilde{\mathcal{P}}(Z_{1})^{\mathsf{T}}\right]\mathclose{}=\alpha_{\mathrm{SO}(n)}\mathbb{I}_{n}.$
(972)

Before we prove these two lemmas, we will use them to prove Theorem 0.36. The
proof idea here is entirely analogous to the original argument of Theorem 0.35
from Ref. [bandeira2016approximating], but with the appropriate replacements
of $\mathcal{P}$ by $\widetilde{\mathcal{P}}$. Nonetheless we sketch the proof
below for completeness.

###### Proof (of Theorem 0.36).

We wish to lower bound the average rounded value

$\operatorname*{\mathbb{E}}[f(Q_{1},\ldots,Q_{m})]=\operatorname*{\mathbb{E}}\mathopen{}\left[\sum_{(u,v)\in
E}\langle
C_{uv},\widetilde{\mathcal{P}}(X_{u}Z)\widetilde{\mathcal{P}}(X_{v}Z)^{\mathsf{T}}\rangle\right]\mathclose{}$
(973)

in terms of the relaxed value $\sum_{(u,v)\in E}\langle
C_{uv},X_{u}X_{v}^{\mathsf{T}}\rangle$. Assuming we have such a lower bound
with ratio $0<\alpha^{2}\leq 1$, this leads to a chain of inequalities
establishing the desired approximation ratio to the original problem:

$\begin{split}\operatorname*{\mathbb{E}}\mathopen{}\left[\sum_{(u,v)\in
E}\langle
C_{uv},\widetilde{\mathcal{P}}(X_{u}Z)\widetilde{\mathcal{P}}(X_{v}Z)^{\mathsf{T}}\rangle\right]\mathclose{}&\geq\alpha^{2}\sum_{(u,v)\in
E}\langle C_{uv},X_{u}X_{v}^{\mathsf{T}}\rangle\\\
&\geq\alpha^{2}\max_{R_{1},\ldots,R_{m}\in\mathrm{O}(n)}\sum_{(u,v)\in
E}\langle C_{uv},R_{u}R_{v}^{\mathsf{T}}\rangle\\\
&\geq\alpha^{2}\max_{R_{1},\ldots,R_{m}\in\mathrm{SO}(n)}\sum_{(u,v)\in
E}\langle C_{uv},R_{u}R_{v}^{\mathsf{T}}\rangle,\end{split}$ (974)

where the second inequality follows from the fact that the relaxation provides