Expanding the reach of quantum optimization with fermionic embeddings

Quadratic programming over orthogonal matrices encompasses a broad class of hard optimization problems that do not have an efficient quantum representation. Such problems are instances of the little noncommutative Grothendieck problem (LNCG), a generalization of binary quadratic programs to continuous, noncommutative variables. In this work, we establish a natural embedding for this class of LNCG problems onto a fermionic Hamiltonian, thereby enabling the study of this classical problem with the tools of quantum information. This embedding is accomplished by a new representation of orthogonal matrices as fermionic quantum states, which we achieve through the well-known double covering of the orthogonal group. Correspondingly, the embedded LNCG Hamiltonian is a two-body fermion model. Determining extremal states of this Hamiltonian provides an outer approximation to the original problem, a quantum analogue to classical semidefinite relaxations. In particular, when optimizing over the \emph{special} orthogonal group our quantum relaxation obeys additional, powerful constraints based on the convex hull of rotation matrices. The classical size of this convex-hull representation is exponential in matrix dimension, whereas our quantum representation requires only a linear number of qubits. Finally, to project the relaxed solution back into the feasible space, we propose rounding procedures which return orthogonal matrices from appropriate measurements of the quantum state. Through numerical experiments we provide evidence that this rounded quantum relaxation can produce high-quality approximations.


Introduction
Finding computational tasks where a quantum computer could have a large speedup is a primary driver for the field of quantum algorithm development.While some examples of quantum advantage are known, such as quantum simulation [1,2], prime number factoring [3], and unstructured search [4], generally speaking computational advantages for industrially relevant calculations are scarce.Specifically in the field of optimization, which has attracted a large amount of attention from quantum algorithms researchers due to the ubiquity and relevance of the computational problems, substantial quantum speedups, even on model problems, are difficult to identify.This difficulty is in part because it is not obvious a priori how the unique features of quantum mechanics-e.g., entanglement, unitarity, and interference-can be leveraged towards a computational advantage [4][5][6].
In this work we take steps toward understanding how to apply quantum computers to optimization problems by demonstrating that the class of optimization problems involving rotation matrices as decision variables has a natural quantum formulation and efficient embedding.Examples of such problems include the joint alignment of points in Euclidean space by isometries, which has applications within the contexts of structural biology via cryogenic electron microscopy (cryo-EM) [7,8] and NMR spectroscopy [9], computer vision [10,11], robotics [12,13], and sensor network localization [14].The central difficulty in solving these problems is twofold: first, the set of orthogonal transformations O(n) is nonconvex, making the optimization landscape challenging to navigate in general.Second, the objectives of these problems are quadratic in the decision variables, making them examples of quadratic programming under orthogonality constraints [15].In this paper we specifically focus on the problem considered by Bandeira et al. [16], which is a special case of the real little noncommutative Grothendieck (LNCG) problem [17].While significant progress has been made in classical algorithms development for finding approximate solutions, for example by semidefinite relaxations [16,[18][19][20][21], guaranteeing highquality solutions remains difficult in general.This paper therefore provides a quantum formulation of the optimization problem, as a first step in exploring the potential use of a quantum computer to obtain more accurate solutions.
The difficulty of the LNCG problem becomes even more pronounced when restricting the decision variables to the group of rotation matrices SO(n) [22,23].One promising approach to resolving this issue is through the convex relaxation of the problem, studied by Saunderson et al. [21,24].They identified that the convex hull of rotation matrices, conv SO(n), is precisely the feasible region of a semidefinite program (SDP) [24].Therefore, standard semidefinite relaxations of the quadratic optimization problem can be straightforwardly augmented with this convex-hull description as an additional constraint [21].They prove that when the problem is defined over particular types of graphs, this enhanced SDP is exact, and for more general instances of the problem they numerically demonstrate that it yields significantly higher-quality approximations than the basic SDP.The use of this convex hull has since been explored in related optimization contexts [25][26][27].Notably however, the semidefinite description of conv SO(n) is exponentially large in n.Roughly speaking, this reflects the complexity of linearizing a nonlinear determinant constraint.One such representation is the socalled positive-semidefinite (PSD) lift of conv SO(n), which is defined through linear functionals on the trace-1, PSD matrices of size 2 n−1 × 2 n−1 .
One may immediately recognize this description as the set of density operators on n − 1 qubits.In this paper we investigate this statement in detail and make a number of connections between the optimization of orthogonal/rotation matrices and the optimization of quantum states, namely fermionic states in second quantization.The upshot is that these connections provide us with a relaxation of the quadratic program into a quantum Hamiltonian problem.Although this relaxation admits solutions (quantum states) which lie outside the feasible space of the original problem, we show that it retains much of the important orthogonal-group structure due to this natural embedding.The notion of quantum relaxations have been previously considered in the context of combinatorial optimization (such as the Max-Cut problem), wherein quantum rounding protocols were proposed to return binary decision variables from the relaxed quantum state [28].In a similar spirit, in this paper we consider rounding protocols which return orthogonal/rotation matrices from our quantum relaxation.
Within the broader context of quantum information theory, our work here also provides an alternative perspective to relaxations of quantum Hamiltonian problems.There is a growing interest in classical methods for approximating quantum many-body problems based on SDP relaxations [29][30][31][32][33][34][35][36][37][38].In that context, rounding procedures are more difficult to formulate because the space of quantum states is exponentially large.For instance, the algorithm may only round to a subset of quantum states with efficient classical descriptions such as product states [29-31, 33, 36] or low-entanglement states [32,34], effectively restricting the approximation from representing the true ground state.Nonetheless, these algorithms can still obtain meaningful approximation ratios of the optimal energy, indicating that such states can at least capture some qualitative properties of the generically entangled ground state.
Our quantum relaxation can be viewed as working in the opposite direction: we construct a manybody Hamiltonian where the optimal solution to the underlying classical quadratic program is essentially a product state.Therefore, we propose preparing an approximation to the ground state of the Hamiltonian, 1 which is then rounded to the nearest product state corresponding to the original classical solution space.This is not unlike quantum approaches to binary optimization such as quantum annealing or the quantum approximate optimization algorithm [6,[39][40][41][42], which explore a state space outside the classical feasible region before projectively measuring, or rounding, the quantum state to binary decision variables.We furthermore provide numerical evidence that the physical qualitative similarity between optimal product and entangled states may translate into quantitative accuracy for the classical optimization problem, in a context beyond discrete combinatorial optimization.
Finally, we remark that Grothendieck-type problems and inequalities have a considerable historical connection to quantum theory.Tsirelson [43] employed Clifford algebras to reformulate the commutative Grothendieck inequality into a statement about classical XOR games with entanglement.Regev and Vidick [44] 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 [45] 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. [30].Recasting noncommutative Grothendieck problems into a quantum Hamiltonian problem may therefore provide new insights into these connections.
The rest of this paper is organized as follows: Section 2 provides a formal description of the optimization problem that we study in this paper and reviews known complexity results of related problems.In Section 3 we describe two well-known applications of the problem: the group synchronization problem and the generalized orthogonal Procrustes problem.Section 4 provides a summary of our quantum relaxation which embeds the optimization problem into a Hamiltonian, and two accompanying rounding protocols.In Section 5 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 6 then establishes the quantum Hamiltonian relaxation of the quadratic optimization problem.Section 7 describes both classical and quantum rounding protocols for relaxations of the problem.Notably, for the classical SDP we derive an approximation ratio for SO(n) rounding.Finally, in Section 8 we demonstrate numerical experiments on random instances of the group synchronization problem for 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 9 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. [16]. 2 Let (V, E) be an undirected graph with m = |V | vertices and edge set E. For integer n ≥ 1, let C ∈ R mn×mn be a symmetric matrix, which for notation we partition into n × n blocks as The quadratic program we wish to solve is of the form max R1,...,Rm∈G (u,v)∈E where G is either the orthogonal group or the special orthogonal group on R n .Here, ⟨A, B⟩ = tr(A T B) denotes the Frobenius inner product on the space of real matrices and I n is the n × n identity matrix.Note that when where now C ∈ R m×m .This is sometimes referred to as the commutative instance of the little Grothendieck problem.Problem (2) 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 (5) is already 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. [46] proved that, assuming the Unique Games conjecture, it is NP-hard to approximate the optimal Max-Cut solution to better than a fraction of (2/π) min θ∈[0,π] θ 1−cos θ ≈ 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 [47].More generally, let (V, E) be fully connected and C ⪰ 0 arbitrary.In this setting, Nesterov [48] showed that GW rounding guarantees an approximation ratio of 2/π ≈ 0.636, which Alon and Naor [49] showed matches the integrality gap of the semidefinite program.Khot and Naor [50] later demonstrated that this approximation ratio is also Unique-Games-hard to exceed, and finally Briët et al. [17] strengthened this result to be unconditionally NP-hard.
For the noncommutative problem (2) 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. [20] provided a rounding procedure of its semidefinite relaxation.Their algorithm achieves an approximation ratio of at least 1/2 √ 2 ≈ 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. [17] for both the real-and complexvalued settings; specifically, they showed that this value is the NP-hardness threshold of a special case of the problem, which is the little noncommutative Grothendieck problem. 3 However, the threshold for Problem (2), which is an special case of the LNCG, is not known.Algorithmically, Bandeira et al. [16] demonstrated constant approximation ratios for Problem (2) when C ⪰ 0 and G = O(n) or U(n) via an (n × n)-dimensional generalization of GW rounding, along with matching integrality gaps.These approximation ratios exceed 1/2, indicating that this subclass is less difficult than the more general formulation of LNCG.However, although the optimization of rotation matrices is of central importance to many applications, we are unaware of any approximation ratio guarantees for the G = SO(n) setting.

Applications of the LNCG problem
Before describing our quantum relaxation, here we motivate the practical interest in Problem (2) by briefly discussing some applications.Throughout, let G = O(n) or 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 typically have low signal-to-noise ratio, and they feature the structure in different, unknown orientations.One approach to solving the group synchronization problem yields best-fit estimates for these orientations via least-squares minimization [51], from which one can produce a model of the desired 3D structure. 4See Ref. [8] for a further overview, and Ref. [11] for a survey of other applications of group synchronization.
The formal problem description is as follows.To each vertex v ∈ [m] := {1, . . ., m} there is an unknown g v ∈ G.An interaction between each pair of vertices connected by an edge (u, v) ∈ E is modeled as g uv = g u g T v .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 + σW uv , where σ ≥ 0 characterizes the strength of the noise and each W uv ∈ R n×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 where ∥A∥ F = ⟨A, A⟩ is the Frobenius norm and we employ the notation R ≡ (R 1 , . . ., R m ).It is straightforward to see that the minimzer of this problem is equivalent to the maximizer of which is precisely in the form of Problem (2).

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. [52] for a comprehensive review. Consider We wish to find an orthogonal transformation R v ∈ G for each S v that best aligns all sets of points simultaneously.That is, for each k ∈ [K] and u, v ∈ [m] we wish to minimize the Euclidean distance Taking least-squares minimization as our objective, we seek to solve From the relation between the vector 2-norm and matrix Frobenius norm, Eq. ( 8) can be formulated as max where each C uv ∈ R n×n is defined as

Summary of results
We now provide a high-level overview of the main contributions of this paper.We provide summary cartoon in Figure 1, 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

Quantum Hamiltonian relaxation
First, consider the setting in which R = (R 1 , . . ., R m ) ∈ O(n) m .We embed this problem into a Hamiltonian by placing n qubits on each vertex v ∈ [m], resulting in a total Hilbert space H ⊗m 2 n of mn qubits.Define the n-qubit Pauli operators where (similarly for X i , Y i ).The Hamiltonian defines our quantum relaxation of the objective f over 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 SO(n) m , we consider instead the (n − 1)-qubit Pauli operators where Π 0 : H 2 n → H 2 n−1 represents the projection onto the even-parity subspace of H 2 n .The construction of the relaxed Hamiltonian for SO(n) is then analogous to Eq. (13): where now the relaxed quantum problem is defined over m(n − 1) qubits.These Hamiltonians serve as relaxations to Problem (2) as follows.It is known that the operators P ij are a representation of the double cover of the orthogonal group, in the sense that for each R ∈ O(n) there exists a unit vector |ϕ(R)⟩ ∈ H 2 n such that all R ij = ⟨ϕ(R)|P ij |ϕ(R)⟩ [24,53].Additionally, R ∈ SO(n) if and only if such |ϕ(R)⟩ has even parity.We go beyond these classical results by identifying each P ij as a one-body fermion operator, which implies that each |ϕ(R)⟩ is a free-fermionic Gaussian state.Hence solving Problem (2) is equivalent to optimizing the Hamiltonian over a subset of quantum states: Dropping these product-state constraints on |ψ⟩ implies the inequalities where D(H) denotes the set of density operators on a Hilbert space H.This establishes our quantum Hamiltonian relaxation.

Quantum rounding
In order to recover orthogonal matrices from a relaxed quantum solution ρ, 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 ρ 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 (2).Consider the mn × mn matrix of expectation values where the off-diagonal blocks are defined as when G = O(n), and we replace the operators P ij with ‹ P ij when G = SO(n).We show that M satisfies the following properties for all states ρ: where conv G is the convex hull of G. Notably, when G = SO(n), M obeys the same constraints as the conv SO(n)-based semidefinite relaxation proposed by Saunderson et al. [21].However, whereas the classical representation of the conv SO(n) constraints involves matrices of size 2 n−1 × 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 tr(P (v) ij ρ) of ρ, as opposed to the two-vertex information tr(Γ (u,v) ij ρ).We consider this rounding procedure due to the fact that, if ρ is a pure Gaussian state satisfying the constraint of Eq. (17), then the matrix of expectation values lies in O(n).On the other hand, for arbitrary density matrices we have the relaxation and again when we replace Both rounding procedures use the standard projection [54] of the matrices X = T uv or Q v to some R ∈ G by finding the nearest orthogonal matrix according to Frobenius-norm distance: This can be solved efficiently as a classical postprocessing step, essentially by computing the singular value decomposition of we instead use the so-called special singular value decomposition [55] of X = U Σ ‹ V T , where Σ = ΣJ and ‹ V = V J, with J being the diagonal matrix assuming that the singular values σ i (X) are in descending order, σ 1 (X) ≥ • • • ≥ σ n (X).Then the solution to Eq. (26 5 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 [53].We review this mathematical construction in Appendix A and only discuss the main aspects here.The Clifford algebra Cl(n) is a 2 n -dimensional real vector space equipped with an inner product and multiplication operation satisfying the anticommutation relation where e 1 , . . ., e n is an orthonormal basis for R n and 1 is the multiplicative identity of the algebra.The orthogonal group is then realized through a quadratic map Q : Cl(n) → R n×n and the identification of a subgroup Pin(n Notably, the elements of Pin(n) have unit norm (with respect to the inner product on Cl(n)).The special orthogonal group, meanwhile, is constructed by considering only the even-parity elements of Cl(n), denoted by Cl 0 (n).The group Because the Clifford algebra Cl(n) is a 2 n -dimensional vector space, we observe that it can be identified with a Hilbert space of n qubits. 5In 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 Cl(n) and H 2 n := (R 2 ) ⊗n as Hilbert spaces.We denote the standard basis of Cl(n) by {e By convention we assume that the elements of I are ordered as i 1 < • • • < i k .Each basis element e I maps onto to a computational basis state |b⟩, 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| mod 2 = |b| mod 2, is therefore preserved, so Cl 0 (n) corresponds to the subspace of H 2 n with even Hamming weight.
To represent the multiplication of algebra elements in this Hilbert space, we use the fact that leftand right-multiplication are linear automorphisms on Cl(n), which are denoted by The action of the algebra can therefore be represented on H 2 n as linear operators.We shall use the matrix representation provided in Ref. [24], 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 It will also be useful to write down the parity automorphism α(e I ) = (−1) |I| e I under this matrix representation.As the notion of parity is equivalent between Cl(n) and H 2 n , α is simply the n-qubit parity operator, It will also be useful to represent the subspace Cl 0 (n) explicitly as an (n − 1)-qubit Hilbert space.This is achieved by the projection from Cl(n) to Cl 0 (n), expressed in Ref. [24] as the 2 n−1 × 2 n matrix It is straightforward to check that Π 0 |b⟩ = 0 if |b| 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 : Cl(n) → R n×n is defined as where π R n is the projector from Cl(n and the conjugation operation x → x is defined as the linear extension of e I = (−1) |I| e i k • • • e i1 .This map associates Clifford algebra elements with orthogonal matrices via the relations A for a review of the construction).In the standard basis of R n , the linear map Q(x) : R n → R n has the matrix elements Using the linear maps λ i , ρ j of left-and right-multiplication by e i , e j , as well as the conjugation identity ⟨x, yz⟩ = ⟨xz, y⟩ in the Clifford algebra, these matrix elements of Q(x) can be rearranged as We now transfer this classical 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. (38): where the expressions in terms of Pauli matrices follow from Eqs. (31) to (33).Then we may rewrite Eq. (38) as where |x⟩ ∈ H 2 n is the quantum state identified with x ∈ Cl(n).Hence, the matrix elements of Q(x) ∈ R n×n possess the interpretation as expectation values of a collection of n 2 Pauli observables , one can work in the even-parity sector directly by projecting the operators as These are (n − 1)-qubit Pauli operators, and we provide explicit expressions in Appendix C. When necessary, we may specify another map ‹ for which ‹ Q(Spin(n)) = SO(n).In general, these double covers are only a subset of the unit sphere in H d (d = 2 n or 2 n−1 ), so not all quantum states mapped by Q yield orthogonal matrices.In Section 5.3 we characterize the elements of Pin(n) and Spin(n) as a class of well-studied quantum states, namely, pure fermionic Gaussian states.

Notation
First we establish some notation.A system of n fermionic modes, described by the creation operators a † 1 , . . ., a † n , can be equivalently represented by the 2n Majorana operators for all i ∈ [n].These operators form a representation for the Clifford algebra Cl(2n), as they satisfy6 The Jordan-Wigner mapping allows us to identify this fermionic system with an n-qubit system via the relations We will work with the two representations interchangeably.
A central tool for describing noninteracting fermions is the Bogoliubov transformation γ → Oγ, where O ∈ O(2n) and This transformation is achieved by fermionic Gaussian unitaries, which are equivalent to matchgate circuits on qubits under the Jordan-Wigner mapping [56][57][58].In particular, we will make use of a subgroup of such unitaries corresponding to be the fermionic Gaussian unitary with the adjoint action In contrast to arbitrary O(2n) transformations, these unitaries do not mix between the γ-and γ-type Majorana operators.

Linear optimization as free-fermion models
Applying the representation of Majorana operators under the Jordan-Wigner transformation, Eqs. ( 47) and (48), to the Clifford algebra automorphisms, Eqs.(31) to (33), we see that λ † i = i γ i and ρ j α = γ j .Therefore the Pauli operators P ij defining the quadratic map Q are equivalent to fermionic one-body operators, Consider now a linear objective function ℓ(X) := ⟨C, X⟩ for some fixed C ∈ R n×n , which we wish to optimize over O(n): max Because we require Writing out the matrix elements explicitly, we see that the objective takes the form where we have defined the noninteracting fermionic Hamiltonian The linear optimization problem is therefore equivalent to solving a free-fermion model, 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 ΣV T , which is computable in time O(n 3 ).This decomposition immediately reveals the diagonal form of the Hamiltonian: Because i γ k γ k = Z k , it follows that the eigenvectors of F (C) are the fermionic Gaussian states with eigenvalues The maximum energy is E 0 n = tr Σ since all singular values are nonnegative.The corresponding eigenstate |ϕ 0 n ⟩ is the maximizer of Eq. (57), so it corresponds to an element ϕ Indeed, the standard classical algorithm [54] for solving Eq. ( 53) uses precisely the same decomposition.From the cyclic property of the trace and the fact that O(n) is a group, we have max where we have employed the change of variables X ′ := U T XV .Again, because Σ has only nonnegative entries, ⟨Σ, X ′ ⟩ achieves its maximum, tr Σ, when X ′ = I n .This implies that the optimal solution is X = U V T .Note that this problem is equivalent to minimizing the Frobenius-norm distance, since arg min Now suppose we wish to optimize ℓ over SO(n).In this setting, one instead computes X = U ‹ V T from the special singular value decomposition of C = U Σ ‹ V T .This ensures that det(X) = 1 while maximizing ℓ(X), as only the smallest singular value σ 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) ∈ O(n) is given by the parity of x ∈ Pin(n), or equivalently the parity of the state |x⟩ in the computational basis.Note also that all fermionic states are eigenstates of the parity operator.To optimize over SO(n), we therefore seek the maximal eigenstate |ϕ b ⟩ of F (C) which has even parity.If ⟨ϕ 0 n |Z ⊗n |ϕ 0 n ⟩ = 1 then we are done.On the other hand, if ⟨ϕ 0 n |Z ⊗n |ϕ 0 n ⟩ = −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 |ϕ 0 n−1 1 ⟩ is then the even-parity state with the largest energy, Finally, we point out that all elements of Pin(n) are free-fermion states.To see this, observe that C is arbitrary.We can therefore construct the family of Hamiltonians {F Clearly, the maximum ⟨C, X⟩ = n within this family is achieved when X = C, each of which corresponds to a fermionic Gaussian state |ϕ⟩ satisfying F (C)|ϕ⟩ = n|ϕ⟩ and Q(ϕ) = C.We note that this argument generalizes the mathematical one presented in Ref. [24], which only considered the eigenvectors lying in 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. [24].The convex hull of O(n) is the set of all matrices with operator norm bounded by 1, On the other hand, the convex hull of SO(n) has a more complicated description in terms of special singular values: Saunderson et al. [24] 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: where the 2 n−1 × 2 n−1 matrices ‹ P ij are defined in Eq. (41). 7ecall that the density operators on a Hilbert space H form the convex hull of its pure states: From Eq. (65) one immediately recognizes that the PSD lift of conv SO(n) corresponds to D(H 2 n−1 ), where we recognize that H 2 n−1 ∼ = Cl 0 (n).Furthermore, the projection of the lift is achieved through the convexification of the map Q : Cl(n) → R n×n , where the fact that Q is quadratic in Cl(n) translates to being linear in D(Cl(n)).Specifically, by a slight abuse of notation we shall extend the definition of Q to act on density operators as Then Eq. ( 65) In Appendix B.1 we show that this statement straightforwardly generalizes for We prove this using the fermionic representation developed in Section 5.3.2, and furthermore use these techniques to provide an alternative derivation for the PSD lift of conv SO(n).The core of our argument is showing that the singular-value conditions of Eqs. ( 63) and (64) translate into bounds on the largest eigenvalue of corresponding n-qubit observables: where ρ ∈ D(Cl(n)) and ρ 0 ∈ D(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 conv O(n) and conv SO(n) can be automatically specified by restricting the support of ρ 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. ( 53) 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 7 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 (2).Let (V, E) be a graph, and associate to each edge (u, v) ∈ E a matrix C uv ∈ R n×n .We label the vertices as V = [m].We wish to maximize the objective First, expand this expression in terms of matrix elements: From the quadratic mapping Q : Cl(n) → R n×n , we know that for each R ∈ G there exists some ϕ ∈ Pin(n) such that R ij = ⟨ϕ|P ij |ϕ⟩.Hence we can express the matrix product as which is now the expectation value of a 2n-qubit Pauli operator with respect to a product state of two Gaussian states |ϕ u ⟩, |ϕ v ⟩.To extend this over the entire graph, we define a Hilbert space of m registers of n qubits each.For each edge (u, v) ∈ E we introduce the Hamiltonian terms where To simplify notation, we shall omit the trivial support w∈V \{u,v} I (w) 2 n when the context is clear.The problem is now reformulated as optimizing the mn-qubit Hamiltonian The exact LNCG problem over O(n) then corresponds to The hardness of this problem is therefore related to finding the optimal separable state for local Hamiltonians, which is NP-hard in general [59][60][61].Dropping these constraints on the state provides a relaxation of the problem, since 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 (Cliffordalgebra products of four Majorana operators) and two-body qudit operators (tensor products of two qudit Pauli operators).Recall that P ij = i γ i γ j is one-body in the fermionic sense.While the operators jk appear to mix both notions, here they in fact coincide.To see this, we consider a global algebra of Majorana operators } 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., the local two-mode Majorana operators in fact do correspond to global two-mode operators: 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: Therefore Eq. ( 75) can be equivalently expressed as a Hamiltonian with two-body fermionic interactions.
Finally, when we wish to optimize over (R 1 , . . ., R m ) ∈ SO(n) m , it is straightforward to see that we can simply replace the terms P ij with ‹ P ij .Defining the quantum relaxation for the SO(n) problem is given by the m(n − 1)-qubit Hamiltonian Clearly, as local (fermionic) Hamiltonians, both H and ‹ H lie in QMA [62,63].In fact, they are also QMA-hard, which can be shown by a reduction to an arbitrary instance of the XY model [64].We prove this reduction in Appendix E. Therefore rather than trying to solve the QMA-hard problem, our goal will be to produce an approximate ground state in polynomial time.This approximation is then classically rounded to the desired collection of orthogonal matrices, which we subsequently describe in Section 7. Afterwards, in Section 8 we will explore the effectiveness of this approximation with numerical experiments.

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 ρ ∈ D(H ⊗m d ) which (approximately) maximizes the energy tr(Hρ).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 ρ appropriately, such that we return some valid approximation R 1 , . . ., R m ∈ 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 SO(n), as conv 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 for each pair of vertices (u, v).We refer this procedure as conv SO(n)-based rounding. 8Our second rounding protocol uses the expectation values of P (v) ij 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 ρ v := tr ¬v (ρ).Therefore we call this approach vertex-marginal rounding.
If ρ 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 [65,66] or a more sophisticated measurement protocol [67]. 9See Appendix I for further comments on this quantum measurement aspect.The rounding algorithms then operate entirely as classical postprocessing after estimating the necessary expectation values.

Approximation ratios for rounding the classical SDP
Before describing our quantum rounding protocols, we first review classical relaxations and rounding procedures for Problem (2) in order to put our results in context.The standard semidefinite relaxation can be expressed as the SDP where C ∈ R mn×mn is the matrix with n × n blocks C uv .If an additional nonconvex constraint rank(M ) = n is imposed, then the solution would be exact: Problem (84), without the rank constraint, is therefore a relaxation of the original problem.However, the solution M ∈ R mn×mn is still PSD, so it can be decomposed as M = XX T , where The rounding algorithm of Bandeira et al. [16] then computes, for each v ∈ [m], where Z is an mn × n Gaussian random matrix whose entries are drawn i.i.d.from N (0, 1/n).When optimizing over G = O(n), this rounded solution guarantees (in expectation) an approximation ratio of where Z 1 is a random n × n matrix with i.i.d.entries from N (0, 1/n).This expression can be evaluated for any given n, and some numerical values are provided in Ref. [16].

The SO(n) setting
The approximation ratios when one demands rounding to SO(n) elements were not originally considered by Bandeira et al. [16].However, as the special orthogonal setting is one of the main draws of our quantum formulation, we first establish analogous classical approximation results here.In Appendix F, we show how the argument of Bandeira et al. [16] can be extended to the SO(n) setting, yielding an approximation ratio of The only algorithmic change is that the rounding procedure is modified to project to the nearest SO(n) element via ‹ P, defined as We also show that E[σ n (Z 1 )] > 0 for all finite n, hence it follows that α 2 SO(n) < α 2 O(n) .This provides evidence that solving for rotations is generally a more difficult problem than optimizing over all orthogonal matrices (see Ref. [23,Section 4.3] for a brief discussion).For small values of n, the numerical values of both the O(n) and SO(n) approximation ratios are (computed using Mathematica): In Appendix F we provide an integral expression for α SO(n) which can be evaluated for arbitrary n.
The problem over SO(n) can be augmented further by introducing constraints to the SDP based on the convex hull of SO(n).Saunderson et al. [21] apply such constraints on each block of M : 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 ∈ conv SO(n) is sufficient to guarantee the much stronger condition M uv ∈ SO(n).Then, when dropping the rank constraint (but leaving the conv 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.
However, the cost of this augmented SDP is the exponential size of representing conv SO(n) elements, as seen by the PSD lift in Eq. (65).In a different work by the same authors, Saunderson et al. [24] proved that, roughly speaking, semidefinite representations of conv SO(n) necessarily have size at least 1 n 2 Ω(n) .

Quantum Gram matrix
We now introduce the rounding procedures for our quantum relaxed solution.Analogous to the classical SDP solution M , we can form a matrix M ∈ R mn×mn from the expectation values of ρ as where and T vu = T T uv for u < v. Just as ⟨C, M ⟩ gives the relaxed objective value (up to rescaling and constant shifts), here we have that ⟨C, M⟩ = 2 n tr(Hρ) + tr(C).In Appendix B we show that for any quantum state, M satisfies the following properties: (97) Furthermore, we show that when ρ is supported only on the even subspace of each single-vertex Hilbert space (or equivalently, if we replace Therefore when optimizing the relaxed Hamiltonian ‹ H for the SO(n) setting, we are guaranteed to automatically satisfy the conv SO(n) constraints.

conv SO(n)-based rounding
Given the construction of the M from quantum expectation values, we proceed to round the Gram matrix as in the classical SDP with conv SO(n) constraints [21].This consists of computing the matrices where the projection to SO(n) can be efficiently computed from the special singular value decomposition, i.e., (recall Eq. ( 27)).Our choice of rounding using the first n × mn "row" of M amounts to fixing R 1 = I n .We note that the same rounding procedure can naturally be applied to the O(n) setting as well, replacing ‹ P with P.

Vertex-marginal rounding
The single-vertex marginals are obtained by tracing out the qudits associated to all but one vertex v ∈ [m], As ρ v ∈ D(H ⊗m d ), from Section 5.3.3we have that Q(ρ v ) ∈ conv G, where we linearly extend the definition of Q to The rounding scheme we propose here then projects Q(ρ v ) to G using either P or ‹ P: We point out that the relaxed Hamiltonian only has two-vertex terms which we seek to maximize.In Appendix G we show that H commutes with U ⊗m (In,V ) for all V ∈ O(n), which we further show implies that H may possess eigenstates whose single-vertex marginals obey Q(σ 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(σ v ): Note that this trace quantity is importantly invariant with respect to the choice of basis for R n .We then augment the objective Hamiltonian with H 1 , defining where ζ > 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 ζ > 0 suffices to break the O(n) symmetry.Furthermore, the rounding procedure always guarantees that the solution is projected back into the feasible space G m .When G = SO(n) we define ‹ H ′ (ζ) 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 R 3 (requiring only two qubits per vertex) is highly relevant to many practical applications, so here we focus on the problem of SO(3) group synchronization.For example, this problem appears in the context of cryo-EM as described in Section 3.1.
To model the problem, we generated random instances by selecting random three-regular graphs ([m], E), uniformly sampling m rotations g 1 , . . ., g m ∈ SO(n) according to the Haar measure, and then constructing C uv = g u g T v + σW uv for each (u, v) ∈ E, where the Gaussian noise matrix W uv ∈ R n×n has i.i.d.elements drawn from N (0, 1) and σ ≥ 0 represents the strength (standard deviation) of this noise.This model describes a typical instance of group synchronization wherein the data is inherently corrupted by a low signal-to-noise ratio [68].Therefore our numerics serve to demonstrate the averagecase performance of our quantum relaxation and rounding, as opposed to a worst-case theoretical analysis.
While the classical conv 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 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 conv SO(n)-based rounding (denoted CR), vertex-marginal rounding (VR), and the classical SDP (without conv SO(n) constraints but using the ‹ P projection to guarantee that the rounded solutions are elements of SO(n)).When using the vertex-rounding method, we employ ‹ H ′ (ζ) as the objective Hamiltonian with ζ = 10 −6 .We explore two classes of quantum states to probe the effectiveness of our relaxation.First, in Section 8.1 we study the maximal eigenstate of the Hamiltonian.This represents the optimal solution to the relaxation, which can be seen as the quantum analogue to solving the classical SDP before rounding.However, computing extremal eigenstates of local Hamiltonians is QMA-hard in general (which is harder than the underlying classical NP-hard problem).Therefore, we next consider in Section 8.2 quantum states which can be prepared in at most polynomial time on a quantum computer.This allows us to study the potential for developing an efficient quantum algorithm derived from the quantum relaxation.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 conv SO(n)based scheme (Section 7.3), while VR denotes the vertex-marginal rounding scheme (Section 7.4).The classical SDP solution was rounded by the standard randomized algorithm [16], 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 σ which defines the problem via Cuv = gug T v + σWuv.

Exact eigenvectors
First, we consider the solution obtained by rounding the maximum eigenvector of ‹ 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 2 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 conv 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 7.4, 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 [21,25] the conv 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 σ, we observe that all methods are fairly consistent.In particular, the conv SO(n)-based rounding only shows an appreciable decrease in approximation quality when the noise is considerable (note that σ = 0.5 ≈ 10 −0.3 is a relatively large amount of noise, since g u g T v 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 ‹ H, here we consider preparing states using ideas from adiabatic quantum computation [69].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 approxima-  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.
tions to the maximum eigenstate of ‹ 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 timedependent Hamiltonian which prepares the state for some T > 0, where T is the time-ordering operator.The final Hamiltonian H f is the desired objective LNCG Hamiltonian, 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 , . . ., R m ∈ SO(n) be the SDP solution.Our initial state is then the product of Gaussian states where each |ϕ(R v )⟩ is the maximum eigenvector of the free-fermion Hamiltonian 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 SO(n)):  As a Gaussian state, |ϕ(R v )⟩ can be prepared exactly from a quantum circuit of O(n 2 ) gates [70,71].
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 Π 0 .We discuss how to perform this circuit recompilation in Appendix C. We comment that this choice of initial state is that of a mean-field state for non-numberpreserving 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.
In adiabatic state preparation, the total evolution time T controls how close the final state |ψ(T )⟩ is to the maximum eigenstate10 of the final Hamiltonian H f .One metric of closeness is how the energy of the prepared state, ⟨ψ(T )|H f |ψ(T )⟩, 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 3, 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 |ψ(T )⟩ 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 ‹ H as expected (thereby also demonstrating that the initial "mean-field" state |ψ(0)⟩ has appreciable overlap).Particularly interesting is the behavior for relatively small total evolution times T , wherein the energy of |ψ(T )⟩ is far from the maximum eigenenergy.Despite this, the approximation quality after rounding the state using M is nearly exact around T ≈ 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 .
Then in Figure 4 we plot the same 50 problem instances (per graph size/noise level) as in Figure 2, but using the quasi-adibatically prepared state |ψ(T )⟩ where we have fixed T = 1 for all graph sizes.The classical SDP results are the same as in Figure 4, 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 3. Namely, although the annealing schedule is too fast to prepare a close approximation to the maximum eigenstate, the rounded solutions (using the conv 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 3.

Discussion and future work
In this paper 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 SO(n) is succinctly captured by a simple linear property of the state-its parity-and the convex bodies conv O(n) and conv 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, conv G-based edge rounding, which is inspired by the fact that the a quantum Gram matrix 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. [21].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-quasiadiabatic 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 NP-hard thresholds for Problem (2) 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, where the N × N × N × 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 ⊗ P kl .While Briët et al. [17] showed that the 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 (2) over SO(n), it is unclear precisely how much harder the SO(n) problem is compared to the O(n) problem.The work by Saunderson et al. [24] establishes a clear distinction between the representation sizes required for conv O(n) and conv SO(n), and this paper 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 NP-hardness threshold of Problem (112) where instead U, V ∈ SO(N ) and leverage the quantum (or equivalently, Clifford-algebraic) representation of 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 conv 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 U(n) and SU(n), essentially by doubling the number of qubits per vertex via the inclusions U(n) ⊂ O(2n) and SU(n) ⊂ SO(2n).However this is likely an inefficient embedding, since the n-qubit Majorana operators already form a representation of Cl(2n).It may therefore be possible to encode complexvalued matrices via a complexification of Q, using the same amount of quantum space.It is interesting to note that Briët et al. [17] in fact utilize a "complex extension" of Clifford algebras when considering Problem (112) over the unitary group, although the usage is different from ours.

A 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.[24,53].
The Clifford algebra Cl(n) of R n is a 2 n -dimensional real vector space, equipped with an inner product ⟨•, •⟩ : Cl(n) × Cl(n) → R and a multiplication operation satisfying the anticommutation relation where {e 1 , . . ., e n } is an orthonormal basis of R n and 1 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 Cl(n) as By convention we order the indices x I e I (115) with each x I ∈ R, and the inner product on Cl(n) is 11⟨x, y⟩ = where y = I⊆ [n] y I e I .Hence Cl(n) is isomorphic as a Hilbert space to R 2 n .Now we show how to realize the orthogonal group O(n) from this algebra.First observe the inclusion R n = span{e i | i ∈ [n]} ⊂ Cl(n).We shall identify the sphere S n−1 ⊂ R n as all u ∈ R n satisfying ⟨u, u⟩ = 1.We then define the Pin group as all possible products of S n−1 elements: It is straightforward to check that this is indeed a group.Each x ∈ Pin(n) is also normalized, ⟨x, x⟩ = 1.
In fact, an equivalent definition of this group is all elements x ∈ Cl(n) satisfying xx = 1, where conjugation x is defined from the linear extension of The Pin group is a double cover of O(n), which can be seen from defining a quadratic map Q : Cl(n) → R n×n .This map arises from the so-called twisted adjoint action, introduced by Atiyah et al. [53]: where the linear map α : Cl(n) → Cl(n) is the parity automorphism, defined by linearly extending Then for any x ∈ Cl(n), the linear map Q(x) : R n → R n is defined as where π R n is the projection from Cl(n) onto R n .To show that Q(Pin(n)) = O(n), it suffices to recognize that, for any u ∈ S n−1 , α(u)vu ∈ R n is the reflection of the vector v ∈ R n across the hyperplane normal to u.To see this, first observe that uv + vu = −2⟨u, v⟩1, which follows from Eq. (113) by linearity.Then which is precisely the elementary reflection as claimed.By the Cartan-Dieudonné theorem, one can implement any orthogonal transformation on R n by composing k ≤ n such reflections about arbitrary hyperplanes u 1 , . . ., u k [72].This characterization coincides precisely with the definition of the Pin group provided in Eq. (117), through the composition of the linear maps Q(u 1 ), . . ., Q(u k ) on R n .Hence for all x ∈ Pin(n), Q(x) is an orthogonal transformation on R n .The double cover property follows from the fact that Q is quadratic in x, so The special orthogonal group arises from the subgroup Spin(n) ⊂ Pin(n) containing only evenparity Clifford elements.First observe that Cl(n) is a Z 2 -graded algebra: where By a Z 2 grading we mean that for each x ∈ Cl a (n) and y ∈ Cl b (n), their product xy lies in Cl a+b mod 2 (n).We say that elements in Cl 0 (n) (resp., Cl 1 (n)) have even (resp., odd) parity.In particular, this grading implies that Cl 0 (n) is a subalgebra, hence its intersection with the Pin group is also a group, which defines Just as the Pin group double covers O(n), so does the Spin group double cover SO(n).This is again a consequence of the Cartan-Dieudonné theorem, wherein all rotations on R n can be decomposed into an even number of (at most n) arbitrary reflections.

B Convex hull of orthogonal matrices and quantum states
Recall the following characterizations of the convex hulls: where {σ i (X)} i∈ [n] and { σ i (X} i∈ [n] are the singular values and special singular values of X in descending order, respectively.Note that σ i (X) = σ i (X) for all i ≤ n − 1 and σ n (X) = sign(det(X))σ n (X).

B.1 PSD lift of conv O(n) and conv SO(n)
In this section we show that , using the quantum formalism described in Section 5.
First we show that for all X ∈ conv O(n), there exists some ρ ∈ D(H 2 n ) which generates X, essentially by the convex extension of Q.Every X ∈ conv O(n) can be expressed as a convex combination For each R µ there exists some Therefore the matrix elements of X can be expressed as where Next we show the reverse direction, that for all ρ ∈ D(H 2 n ), the matrix Therefore we take the singular value decomposition of X = U ΣV T and, using P ij = i γ i γ j , each singular value is equal to For the restriction to conv SO(n), the first argument is essentially the same.One merely replaces ).For the reverse direction, we instead employ the special singular value decomposition which yields where now . 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 ρ only has support on the even-parity subspace, so tr(Z ⊗n ρ) = 1.Furthermore, because the special singular value decomposition guarantees that for all subsets I ⊆ [n] of odd size.By linearity, otherwise, and we have used the fact that i γ k γ k = Z k .It therefore suffices to examine the spectrum of where ⊕ 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| mod 2 = 0.Because |I| is odd, so too is |z|, hence |z ⊕ b| mod 2 = 1.This implies that there must be at least one term in the sum of Eq. (135) which is negative, so it can only take integer values at most n − 2. This establishes Eq. (133), hence X ∈ conv SO(n).

B.2 Relation to conv 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 conv SO(n).Recall that this relaxation can be formulated as We will show that the Gram matrix M constructed from the measurements of a quantum state ρ, defined in Section 7.2 as obeys the constraints of Eq. (136).Specifically, when the marginals of ρ on each vertex are even-parity states (recall this is equivalent to replacing Γ ij with Γ ij ), we obtain the conv SO(n) condition, whereas when the parity of ρ is not fixed then M uv ∈ conv O(n).First, we show that M is positive semidefinite for all quantum states.
Lemma B.1.Let M ∈ R mn×mn be defined as in Eq. (137).For all ρ ∈ D(H ⊗m 2 n ), M ⪰ 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 M above, we first construct a matrix M ′ ⪰ 0 which turns out to simply be M ′ = nM.
For each k ∈ [n] define the Hermitian operator ∈ R are arbitrary coefficients.Consider its square, jk . (139) Note that the terms with u ̸ = 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 = i γ i γ k , the diagonal terms reduce to (suppressing superscripts here) Plugging this result into Eq.( 139) and summing over all k ∈ [n], we obtain where we have collected the coefficients c (v) i into a vector c ∈ R mn .Similarly, if we arrange the expectation values tr(Γ (u,v) ij ρ) into a matrix T ∈ R mn×mn (where the n × n blocks on the diagonal are 0), then the expectation value of the sum-of-squares operator is tr where we have defined k is a sum of PSD operators, its expectation value is always nonnegative, hence ⟨c, M ′ c⟩ ≥ 0. This inequality holds for all vectors c ∈ R mn , so M ′ ⪰ 0 and hence M = M ′ /n ⪰ 0 as claimed.
The fact that the diagonal blocks M vv = I n holds by definition.Finally, we need to show that each block of M lies in conv(SO(n)) when ρ has even parity.A straightforward corollary of this result is that the blocks lie in conv O(n) when ρ does not have fixed parity.Note that in Section 5.3.3we showed that the matrix of expectation values tr(P ij ρ 1 ) lies in conv O(n) for any n-qubit density operator ρ 1 , a straightforward extension of the PSD-lift representation of conv SO(n) presented in Ref. [24].Here we instead show that the matrix of expectation values 1 n tr(Γ ij ρ 2 ) for any 2n-qubit density operator ρ 2 also lies in conv O(n), and is an element of conv SO(n) when ρ 2 has support only in the even-parity sector. Because and because these convex hulls are closed under transposition, all that remains is to prove the statement for the blocks M uv when u < v.We show this by considering all density matrices on the reduced two-vertex Hilbert space.
Proof.Consider the special singular value decomposition of T = U Σ ‹ V T .We can express the special singular values as where Now we examine the inclusion in conv SO(n).Suppose that ρ 2 is an even-parity state.By virtue of the special singular value decomposition, we have that det(U ‹ V T ) = 1 which implies that the Gaussian unitary U (U,In) ⊗ U ( ‹ V ,In) preserves the parity of ρ 2 : tr(Z ⊗2n ρ 2 ) = tr(Z ⊗2n ρ ′ 2 ) = 1.For T /n to lie in conv SO(n), the following inequality from Eq. (128) must hold: for all subsets I ⊆ [n] of odd size.To show this, first we write the left-hand side in terms of the result derived from Eq. (144): where the string z = z 1 • • • z n ∈ {0, 1} n is defined as Note that |I| being odd implies that the Hamming weight of z is also odd.To bound Eq. ( 147) we shall seek a bound on the largest eigenvalue of the Hermitian operator over the space of even-parity states.Here we use ⊕ 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 γ We can therefore seek their simultaneous eigenvectors, which can be determined from looking at the Jordan-Wigner representation of the Majorana operators: These operators are diagonalized by the 2n-qubit state where x, y ∈ {0, 1} n , and the Bell state |β(x j , y j )⟩ between the jth qubits across the two subsystems is defined as Indeed there are 2 2n such states |β(x, y)⟩, so they form an orthonormal basis for the 2n qubits.The eigenvalues of Eqs. ( 150) and ( 151) can be determined by a standard computation, Taking the appropriate products furnishes the eigenvalues of the B and A z as Altogether we arrive at the expression for the eigenvalues of A z B, We wish to find the largest value this can take over even-parity states.
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 → {0, 1} n stores the parity information of the (k − 1)-length substring of its input into the kth bit of its output: 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 ±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) ≤ 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)⊕1 n .In this case b(x, y) = −n and a z (x, y) ≥ −(n − 2), yielding the same bound a z (x, y)b(x, y) ≤ 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 ⊆ [n]).This bounds the value of Eq. (147) by n(n − 2), hence validating the inequality of Eq. ( 146).Thus T /n ∈ conv SO(n) whenever ρ 2 has even parity.
As usual, replacing the operators Γ ij with Γ ij is equivalent to enforcing the even-parity constraint.In fact, since P jk , the equivalent constraint involves the reduced single-vertex marginals, tr(Z ⊗n ⊗ I 2 n ρ 2 ) = tr(I 2 n ⊗ Z ⊗n ρ 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.
C Details for working in the even-parity subspace First we provide an expression for n-qubit Pauli operators projected to Cl 0 (n).Let A straightforward calculation yields the conditional expression: Notably, if A does not commute with the parity operator then A = 0. Now we generalize from the main text, defining the operator which is the projector Cl(n) → Cl k (n).These operators obey Given a state ρ ∈ D(H 2 n−1 ), the expectation values of ‹ P ij satisfy tr( ‹ where ρ 0 := Π T 0 ρΠ 0 ∈ D(H 2 n ) has only support on even-parity computational basis states.By Eq. (169) we can "invert" this relation in the following sense: given a state ρ ∈ D(H 2 n ) which only has support on the even subspace, its (n − 1)-qubit representation is Π 0 ρΠ T 0 ∈ D(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 8.2 to prepare the initial state of the quasiadiabatic evolution).It is well known how to compile linear-depth circuits for this task [70], however this is within the standard n-qubit representation.Here we show how to translate those circuits into the (n − 1)-qubit representation under Π 0 .In fact, these techniques apply to any sequence of gates which commute with the parity operator Z ⊗n . Let We can assume without loss of generality that the initial state is the vacuum |0 n ⟩ and that all gates V ℓ commute with the parity operator, as otherwise |ψ⟩ would not lie in the even subspace of H 2 n . 13For example, parity-preserving Gaussian unitaries can be decomposed into single-and two-qubit gates of the form e −iθZi , e −iθXiXi+1 .We wish to obtain a circuit description for preparing the state Π 0 |ψ⟩ ∈ H 2 n−1 .By Eq. (170), Π T 0 Π 0 |0 n ⟩ = |0 n ⟩, and furthermore Π 0 |0 n ⟩ = |0 n−1 ⟩ is the (n − 1)-qubit representation of the vacuum.Therefore Because Π 0 is not unitary (Π T 0 is merely an isometry), we cannot simply insert terms like Π T 0 Π 0 in between each gate.However, observe that if each V ℓ preserves parity, then they can be block diagonalized into the even and odd subspaces of where and conjugation by the projector Π 0 precisely extracts the first block of this matrix: This sequence of gates is what we wish to implement on the physical (n − 1)-qubit register.When the gates V ℓ take the form for some n-qubit Pauli operator A ℓ , then where A ℓ := Π 0 A ℓ Π T 0 .Note that this calculation assumes that V ℓ commutes with parity, hence [A ℓ , Z ⊗n ] = 0, which guarantees that A ℓ is unitary and Hermitian and hences furnishes the final line of Eq. (177).(If they did not commute then the decomposition of Eq. (173) would not be valid to begin with.)

D Small n examples
In this appendix we write down the Pauli operators P ij and ‹ P ij for small but relevant values of n, to provide the reader with some concerte examples of how to construct the LNCG Hamiltonian.
It then follows that the LNCG interaction terms are and so the full Hamiltonian acting on is indeed the classical Ising Hamiltonian, with weights C uv ∈ R. It is also instructive to write down the elements of Pin(1) as quantum states.
The Clifford algebra Cl( 1) is spanned by e 1 and e ∅ = 1, with the only elements of S 0 being ±e 1 .
E QMA-hardness of the relaxed Hamiltonian Here we show that the relaxed LNCG Hamiltonian is QMA-hard for both O(n) and SO(n), for all n ≥ 2. Note that this hardness result does not hold for n = 1, since O(1) corresponds precisely to the classic NP-complete problem of combinatorial optimization, while SO(1) is the trivial group.We prove this by reducing to an instance of the XY model.Placing m qubits on a graph G = (V, E), this model is defined as where α uv ∈ R are arbitrary.Deciding the ground-state energy of this Hamiltonian to within inversepolynomial precision is QMA-complete [64].We will show that the LNCG Hamiltonian for both O(2) and SO(2) can be reduced to H XY .Because O(2) and SO(2) can be trivially embedded into O(n) and SO(n) respectively, for any n ≥ 2, this proves the claim.
To begin, recall the definitions of the relaxed Hamiltonians.
are the two-qubit operators defined on each vertex v ∈ V .Defining an instance of C uv as we get To see that this Hamiltonian is equivalent to H XY up to a degeneracy of the spectrum (because H XY is over m qubits, while H O( 2) is over 2m qubits), observe that on each vertex we can apply the two-qubit Clifford unitary CNOT • (U ⊗ I 2 ), where U := 1 √ 2 (Y + Z).This transformation on each vertex achieves hence demonstrating the reduction.
Next we show the SO(2) reduction.The projected single-qubit operators per vertex can be determined via Eq.(167): Setting C uv as before, Eq. (187), we get14 Relabeling Z → Y via the local Clifford U furnishes the claim.

F Classical approximation ratio for SO(n)
Approximation ratios for the rounded solution of the classical semidefinite relaxation of Problem (2) were obtained in Ref. [16] for the cases G = O(n) and U(n).However, no such approximation ratios were derived for the case of SO(n).Here we adapt the argument of Ref. [16] 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 SO(n) that is strictly worse than the previously studied O(n) case.
Recall that the semidefinite relaxation of Problem (2) can be formulated as To round the relaxed solution back into the feasible space of orthogonal matrices, Ref. [16] proposes the following randomized algorithm with a guarantee on the approximation ratio.
Theorem F.1 ([16,Theorem 4]).Let X 1 , . . ., X n ∈ R n×mn be a solution to Problem (192).Let Z be an mn × n Gaussian random matrix whose entries are drawn i.i.d.from N (0, 1/n).Compute the orthogonal matrices where P(X) = arg min Y ∈O(n) ∥Y − X∥ F .The expected value of this approximate solution (averaged over Z) obeys The approximation ratio α 2 O(n) is defined by the average singular value of random Gaussian n where σ i (Z 1 ) is the ith singular value of Z 1 .
Our adaptation to the SO(n) setting simply replaces the rounding operator P with ‹ P(X) := arg min Y ∈SO(n) ∥Y − X∥ F .With this change we obtain an analogous result for optimizing over SO(n) elements with the same classical semidefinite program: Theorem F.2. Let X 1 , . . ., X n and Z be as in Theorem F.1.Compute the rotation matrices where ‹ P(X) = arg min Y ∈SO(n) ∥Y − X∥ F .The expected value of this approximate solution (averaged over Z) obeys The approximation ratio α 2 SO(n) is defined by the average n − 1 largest singular values of random Gaussian n × n matrices Z 1 ∼ N (0, Because singular values are nonnegative, it is clear that In particular we will see that E[σ n (Z 1 )] > 0 for all finite n, so the rounding algorithm guarantees a strictly smaller approximation ratio for the problem over SO(n) than over O(n).
The proof of Theorem F.1 requires two lemmas regarding the expected value of random Gaussian matrices under the rounding operator P. Analogously, our proof of Theorem F.2 requires a modification of those lemmas when P is replaced by ‹ P.
Lemma F.3 (Adapted from [16,Lemma 5]).Let M, N ∈ R n×mn obey M M T = N N T = I n .For Z ∈ R mn×n with i.i.d.entries drawn from N (0, n −1 ), we have This lemma is proved with the help of the following lemma.
Before we prove these two lemmas, we will use them to prove Theorem F.2.The proof idea here is entirely analogous to the original argument of Theorem F.1 from Ref. [16], but with the appropriate replacements of P by ‹ P. Nonetheless we sketch the proof below for completeness.
Proof (of Theorem F.2). We wish to lower bound the average rounded value in terms of the relaxed value (u,v)∈E ⟨C uv , X u X T v ⟩.Assuming we have such a lower bound with ratio 0 < α 2 ≤ 1, this leads to a chain of inequalities establishing the desired approximation ratio to the original problem: where the second inequality follows from the fact that the relaxation provides an upper bound to the original problem, and the third inequality is a consequence of SO(n) ⊂ O(n).The task is then to determine such an α which satisfies the first inequality of Eq. ( 203).The core argument is a generalization of the Rietz method [49], which proceeds by constructing a positive semidefinite matrix S ∈ R mn×mn whose (u, v)th block is defined as The expected value of this matrix is is the quantity we wish to bound.To compute the expected values of the two cross terms, we invoke Lemma F.3 which holds because Thus, setting α = α SO(n) , we obtain (207) Finally, using the fact that C, S ⪰ 0, we have that ⟨C, S⟩ ≥ 0 and so E⟨C, S⟩ ≥ 0, which implies that Then by Eq. ( 203) the claim follows.
We now establish the value of from Lemmas F.3 and F.4.Because Lemma F.3 is somewhat technical and the argument is virtually unchanged by replacing P with ‹ P, we refer the reader to Ref. [16] for proof details.Instead, we simply note that the only part of the proof for Lemma F.3 which does depend on the change to ‹ P is the final result, wherein it is established that where Z 1 ∈ R n×n has entries i.i.d.from N (0, 1/n).Thus proving Lemma F.4 is the key component in establishing the value of the approximation ratio α 2 SO(n) .
Proof (of Lemma F.4). Consider the singular value decomposition of Z 1 = U ΣV T ∈ R n×n .Its special singular value decomposition can be written as Note that J U V T = J U J V .Using the fact that the (special) rounding operator returns we have Because Z 1 is a random Gaussian matrix with i.i.d.entries, its singular values and left-and rightsingular vectors are distributed independently [73].In particular, both U and V are distributed according to the Haar measure on O(n).The expected value of Eq. ( 213) can therefore be split into three independent averages: The Haar average over U ∼ O(n) in Eq. ( 215) is well-known [74] to be proportional to the identity, E[U ΣU T ] = λI n , and the constant of proportionality can be determined by considering its trace: Hence λ = α SO(n) and so The corresponding statement for E î Z 1 ‹ P(Z 1 ) T ó follows completely analogously, essentially by interchanging the roles of U and V .The entire argument is equivalent because U and V are i.i.d.
To numerically evaluate α SO(n) we can use the linearity of expectation, The distribution of singular values of random Gaussian matrices can be analyzed from the theory of Wishart matrices.In particular, Z 1 Z T 1 = U Σ 2 U T is a Wishart matrix with covariance matrix I n /n, so the distribution of singular values Σ is the square root of the Wishart distribution of eigenvalues.The quantity α O(n) in terms of the marginal distribution p (avg) n (x) of Wishart eigenvalues x ∈ (0, ∞) was studied in Ref. [16], yielding the expression Note the factor of n −1/2 , which is introduced because the distribution p (avg) n (x) is normalized to have unit variance.An explicit expression of p (avg) n (x) can be found in Refs.[16,Lemma 21] and [75,Eq. (16)].
For our newly derived approximation ratio α SO(n) , we need to additionally evaluate the expected smallest singular value of this Wishart distribution.This minimum-eigenvalue distribution was studied in Ref. [76], wherein an analytical expression was derived (again assuming unit variance): Now suppose that |ψ⟩ is nondegenerate.Then we have that |ψ⟩⟨ψ| = |ψ(V )⟩⟨ψ(V )| for all V ∈ O(n), and in particular we can take the Haar integral over O(n) of this identity: where µ is the normalized Haar measure satisfying µ(O(n)) = 1.Because the Haar integral over linear functions vanishes, i.e., O(n) dµ(V ) V ij = 0 [74], it follows that Furthermore, because i γ

H The Pin group from quantum circuits
In the main text we showed that each x ∈ Pin(n) corresponds to the eigenstates of a family of freefermion Hamiltonians, which are (pure) fermionic Gaussian states.Here we provide an alternative perspective of this correspondence, using quantum circuits which prepare such states.
Recall that every x ∈ Pin(n) can be written as for some k ≤ n, where we may expand each u j ∈ S n−1 in the standard basis as for some unit vector v (j) ∈ R n .The product of these unit vectors can be expressed using the rightmultiplication operator ρ uj acting on the identity element, On the other hand, consider the so-called Clifford loader [77], a circuit primitive defined (in our notation) as for any unit vector v ∈ R n .It is straightforward to check that this operator is Hermitian and unitary, and Ref. [77] provides an explicit circuit constructions based on two-qubit Givens rotation primitives. 15sing the relation γ i = ρ i α and acting this circuit on the vacuum state |0 n ⟩ ≡ |e ∅ ⟩, we see that the state indeed is equivalent to x ∈ Pin(n), up to a global sign (recall that α 2 = 1 and α|e ∅ ⟩ = |e ∅ ⟩).In other words, the Clifford loader is precisely the quantum-circuit representation of generators of the Pin group.
It is worth noting that in Ref. [77] they construct "subspace states" from this composition of Clifford loaders.In the language of fermions, subspace states are Slater determinants: free-fermion states with fixed particle number.Preparing Slater determinants in this fashion requires that the unit vectors v (1) , . . ., v (k) be linearly independent (and thus, without loss of generality, they can be made orthonormal while preserving the subspace that they span, hence the alternative name).However, the definition of the Pin group demands all possible unit vectors in such products, not just those which are linearly independent.Indeed, one can see that if the state |x⟩ is a Slater determinant, then its trace is an integer, as where N = i∈[n] a † i a i is the total number operator.Clearly not all orthogonal matrices have integer trace, so Slater determinants are insufficient to cover all of Pin(n).To reach the remaining elements, we note that if the unit vectors are linearly dependent, then one can show that the |x⟩ = Γ(v (k) ) • • • Γ(v (1) )|0 n ⟩ does not have fixed particle number, so ⟨x|N |x⟩ is not necessarily an integer.
for each (u, v) ∈ E and i, j ∈ [n].When considering G = O(n), because each P (u) jk is a fermionic two-body operator, we can straightforwardly apply the partial tomography schemes developed for local fermionic systems, such as Majorana swap networks [65] or classical shadows [66].In either case, the measurement circuits required are fermionic Gaussian unitaries and the sample complexity is O(N 2 /ϵ 2 ), where N = n|V | is the total number of qubits and ϵ > 0 is the desired estimation precision of each expectation value.

I.2 Tomography of vertex marginals
The vertex-rounding procedure requires the expectation values of only single-qudit observables In this case the observables being measured commute across vertices, so it suffices to talk about the tomography of a single vertex, as the same process can be executed in parallel across all vertices.Again, because these operators are fermionic one-body observables, the same fermionic partial tomography technology [65,66] can be applied here, incurring a sampling cost of O(n/ϵ 2 ).In fact, further constant-factor savings can be achieved in the one-body setting by using the measurement scheme introduced in Ref. [78].This scheme requires only particleconserving fermionic Gaussian unitaries, which can be compiled with only half the depth of the more general Gaussian unitaries required of the previous two methods.Note that each operator P ij is of the form of either XX or Y Y when |i − j| = 1, and so they correspond precisely to the observables measured to reconstruct the real part of the fermionic one-body reduced density matrix [78].

I.3 Estimating observables via gradient method
Ref. [67] introduces a quantum algorithm for estimating a large collection of (generically noncommuting) M observables {O j | j ∈ [M ]} to precision ϵ by encoding their expectation values into the gradient of a function.This function is implemented as a quantum circuit which prepares the state of interest and applies ‹ O( √ M /ϵ) gates of the form c-e −iθOj , controlled on O(M log(1/ϵ)) ancilla qubits.Finally, using the algorithm of Ref. [79] for gradient estimation, one calls this circuit ‹ O( √ M /ϵ) times to estimate the encoded expectation values (the notation ‹ O(•) suppresses polylogarithmic factors).Although this approach demands additional qubits and more complicated circuitry, it has the striking advantage of a quadratically improved scaling in the number of state preparations with respect to estimation error ϵ, compared to the refinement of sampling error in tomographic approaches.In our context, we have either M = n 3 |E| or M = n 2 |V | observables of interest (satisfying a technical requirement of having their spectral norms bounded by 1), corresponding to the measurement of edge or vertex terms respectively.The gates required are then simply controlled Pauli rotations.

Figure 1 :
Figure 1: A cartoon description of our quantum encoding of the LNCG problem, compared to the classical formulation.(Left) The description of the problem that we consider, which is described by a graph ([m], E) and n × n matrices Cuv for each edge (u, v) ∈ E. We wish to assign elements of O(n) or SO(n) to each vertex such that Eq. (11) is maximized.(Center top) The standard classical relaxation optimizes an mn × mn PSD matrix M via semidefinite programming.(Right top) The classical rounding procedure, which returns a collection of orthogonal matrices from M .(Center bottom) Our quantum formulation of the problem as a two-body fermionic Hamiltonian.On each vertex we place a d-dimensional Hilbert space, with interaction terms Huv on the edges constructed from each Cuv.The classical solution to the LNCG problem lies as a subset of this full Hilbert space.(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.

Algorithm 2 :
Rounding vertex marginals Data: Quantum state ρ ∈ D(H ⊗m d ) over a graph of m vertices, each with local Hilbert space of dimension d

Figure 2 :
Figure 2: Approximation ratios for solutions obtained from rounding the maximum eigenvector of the relaxed Hamiltonian ‹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 conv SO(n)based scheme (Section 7.3), while VR denotes the vertex-marginal rounding scheme (Section 7.4).The classical SDP solution was rounded by the standard randomized algorithm[16], 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 σ which defines the problem via Cuv = gug T v + σWuv.

Figure 3 :
Figure 3: 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.

Figure 4 :
Figure4: Approximation ratios for solutions obtained from rounding the "adiabatically" evolved state |ψ(T )⟩ with fixed T = 1 for all m, σ.Problem instances and visualization is the same as in Figure2.We also include the maximum eigenvalue of the relaxed Hamiltonian and the energy of the prepared unrounded state, to demonstrate how far |ψ(T )⟩ 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 σ.

D. 1
The Ising model from the 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 = O(1).Each local Hilbert space H d with d = 2 is simply a qubit, and we only have to consider a single Pauli operator on each local qubit,

)(
We shall comment on the distribution D of singular values later.)Because O(n) is evenly divided into its unconnected (+1)-and (−1)-determinant components, the average determinant vanishes:E V ∼O(n) [det V ] = 0.This leaves us with E î ‹ P(Z 1 )Z T 1 ó = E U ΣU T ,
D.2 The projected operators of the SO(3) settingAs SO(3) is arguably the most ubiquitous group for physical applications, for reference we explicitly write down its Pauli operators under the projection to Cl 0 (3) ∼ = 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. (167) we have

)
(2)orem E.1.Let G = (V, E) be a graph with m vertices.For any collection of α uv ∈ R, where (u, v) ∈ E, there exists an instance of H O(2) and H SO(2)which are equivalent to H XY .