Robustness of Magic and Symmetries of the Stabiliser Polytope

We give a new algorithm for computing the robustness of magic - a measure of the utility of quantum states as a computational resource. Our work is motivated by the magic state model of fault-tolerant quantum computation. In this model, all unitaries belong to the Clifford group. Non-Clifford operations are effected by injecting non-stabiliser states, which are referred to as magic states in this context. The robustness of magic measures the complexity of simulating such a circuit using a classical Monte Carlo algorithm. It is closely related to the degree negativity that slows down Monte Carlo simulations through the infamous sign problem. Surprisingly, the robustness of magic is submultiplicative. This implies that the classical simulation overhead scales subexponentially with the number of injected magic states - better than a naive analysis would suggest. However, determining the robustness of n copies of a magic state is difficult, as its definition involves a convex optimisation problem in a 4^n-dimensional space. In this paper, we make use of inherent symmetries to reduce the problem to n dimensions. The total run-time of our algorithm, while still exponential in n, is super-polynomially faster than previously published methods. We provide a computer implementation and give the robustness of up to 10 copies of the most commonly used magic states. Guided by the exact results, we find a finite hierarchy of approximate solutions where each level can be evaluated in polynomial time and yields rigorous upper bounds to the robustness. Technically, we use symmetries of the stabiliser polytope to connect the robustness of magic to the geometry of a low-dimensional convex polytope generated by certain signed quantum weight enumerators. As a by-product, we characterised the automorphism group of the stabiliser polytope, and, more generally, of projections onto complex projective 3-designs.


Introduction
In fault-tolerant quantum computation (for a recent review, see Ref. [8]), each logical qubit is encoded in a non-local subspace of a number of physical qubits. There are several ways of effecting a unitary transformation of logical qubits. In the simplest case, logical unitaries can be implemented transversally, i.e. by local gates acting on the physical qubits. Unfortunately, a no-go theorem by Eastin and Knill [13] states that there are no quantum codes that allow for a universal set of transversal gates.
In the magic state model [5], the logical gate set is chosen to be the Clifford group, which can be implemented transversally in various quantum codes using their physical counterparts. Any logical non-Clifford gate would promote the Clifford group to universality. This remaining problem is solved by providing an auxiliary qubit in a non-stabiliser state. Using a circuit gadget (which only requires Clifford operations), one can turn this auxiliary state into a non-Clifford gate (Fig. 1). The auxiliary qubit state is consumed in the process, so that one such input needs to be injected for each non-Clifford gate. These inputs are the magic states from which the protocol derives its name.  A common choice for a non-Clifford gate is the T -gate T = diag(1, e iπ/4 ), which is realised by the following magic state Moreover, there is a second magic state, |T , which realises the non-Clifford gate diag(1, e iπ/6 ). Their Bloch representation is shown in Fig. 5. Interestingly, it has been found that even certain mixed states can "supply the magic" to promote a Clifford circuit to universality. Indeed, a process called magic state distillation (Fig. 2) can turn many copies of some mixed state ρ into a pure magic state using Clifford unitaries and computational basis measurements [5,34]. Magic state distillation motivates the search for quantitative measures of the "computational utility" of auxiliary states. This analysis turns out to be slightly simpler for quantum systems with odd-dimensional Hilbert spaces [27,36,37], as the theory of stabiliser states is somewhat better-behaved in this case, and there is a better-developed toolbox of "phase space methods" available in this case (see e.g. Refs. [15,24,42]). However, as qubits are the paradigmatic systems for quantum computation, quantitative resource theories for multi-qubit magic states have since been developed [6,23].
The starting point of these theories is the Gottesman-Knill Theorem [30]. It states that quantum circuits consisting only of preparations of stabiliser states, Clifford unitaries, and computational basis measurements can be efficiently simulated on a classical computer. Therefore, if the auxiliary states are stabilisers, there can be no quantum computational advantage. Next, assume that an auxiliary n-qubit state ρ is an element of the stabiliser polytope SP n , i.e.
where (p i ) i is a probability distribution and the s i = |ψ i ψ i | are stabiliser states. This readily gives rise to an efficient classical randomised algorithm that will draw outcomes from the same distribution as a quantum computer would [39], provided that one can sample efficiently from the probability distribution (p i ) i : Indeed, draw s i with probability p i , and then continue to simulate the further time evolution using Gottesman-Knill. Thus, density matrices contained in the convex hull of stabiliser states are equally useless as computational resource states in the magic state model (Fig. 3).  Figure 3: Bloch representation of the the two most commonly considered magic states |H and |T . These states lie outside of the octahedron spanned by 1-qubit stabiliser states having a Bloch vector orthogonal to an edge (|H ) or a facet (|T ) of the stabiliser octahedron. The intersection of their Bloch vector with the facet or edge is marked with a blue dot. Certain mixed states can be used to distil these pure states using Clifford unitaries and measurements. However, states lying inside the stabiliser polytope are useless as a resource state.
Since the stabiliser states {s i } i span the space of Hermitian operators, any auxiliary state can be expanded as ρ = i x i s i , with coefficients x i that are not necessarily nonnegative. However, taking traces on both sides shows that the expansion is affine, i.e. i x i = 1. It is well-known in the theory of Quantum Monte Carlo methods [17] that the probabilistic algorithm sketched above can be extended to the more general scenario. However, the runtime will increase with the total amount of "negativity" in the expansion coefficients x i . This is the dreaded sign problem. A precise theory of the simulation runtime in the context of quantum computation has been developed in Ref. [31] and applied to the magic state model in Ref. [23]. More precisely, they define the robustness of magic (RoM) as where the sum ranges over stabiliser states {s 1 , . . . , s N } and the 1 -norm measures the "amount of negativity" in the affine combination. Then, the number of samples which have to be taken in the Monte Carlo simulation scales as O(R(ρ) 2 ) [23,31]. In addition to measuring the "computational utility" in the above precise sense, the RoM has further interpretations. For example, it can be used to systematically lower-bound the number of non-Clifford gates required to synthesise certain unitaries, namely those that allow for a magic state realisation [23]. Lastly, the RoM derives its name from the fact that it quantifies the robustness of a state's computational utility against noise processes. A precise account of this point of view is given in Section 2.
Interestingly, the RoM is submultiplicative, i.e. R(ρ ⊗2 ) ≤ R(ρ) 2 , where the inequality is usually strict [23]. That means that the simulation effort of a magic state circuit grows subexponentially with the number of injected magic states-an intriguing phenomenon. Therefore, a quantity of interest is the regularised RoM : Unfortunately, computing R(ρ ⊗n ) seems to be a difficult task. For ρ being a singlequbit state, the tensor power ρ ⊗n lives in an 4 n -dimensional space, and the sum over the s i in the definition (2) of the RoM has to range over the 2 O(n 2 ) stabiliser states defined for n-qubit systems. Any direct implementation of the optimisation problem (2) will thus quickly became computationally intractable-and, indeed, Howard and Campbell [23] could carry it out only up to n = 5.
The starting point of this work is the observation that there is a large symmetry group shared by ρ ⊗n and the stabiliser polytope. Thus, we formulate the optimisation in a space where the joint symmetries have been "modded out". The space of operators invariant under the joint symmetry group turns out to have a dimension mildly polynomial in n. For the especially interesting cases where the state is |H ⊗n or |T ⊗n , the dimension reduces further to exactly n. While the projection of the stabiliser polytope to this invariant space ( Fig. 4) still has exponentially many vertices, it turns out that formulating the optimisation problem in this symmetry-reduced way leads to a super-polynomially faster algorithm.
Equipped with the knowledge of the exact solution to Eq. (2) for the commonly used magic states |H ⊗n and |T ⊗n and n ≤ 10 qubits, we formulate a relaxation of the RoM problem for these states which yields an upper bound for the exact RoM. These approximations are in excellent agreement with the exact data for n ≤ 10 and can be carried out for up to 26 qubits. What is more, we can not only compute the RoM bounds for these approximations, but also find the corresponding affine decompositions ρ ⊗n = i x i s i , which can directly be used in Monte Carlo simulations. Furthermore, we find a hierarchy of such RoM approximations by restricting to k-partite entangled stabiliser states which converges to the exact RoM. Interestingly, every level of the hierarchy can be computed in polynomial time.
Finally, both the exact and approximate results imply a runtime of O(2 0.737t ) for simulating a circuit with t T gates using the RoM algorithm. Moreover, our analysis suggests that this runtime is the optimal one that can be achieved using a RoM algorithm. Our work improves on the previously known runtime of O(2 0.753t ) derived in Ref. [23]. Note that the RoM algorithm is able to simulate noisy circuits and mixed states. This is in contrast to simulation algorithms based on the so-called stabiliser rank which can achieve a runtime of O(2 0.48t ) for pure states [4,6,7]. This paper is organised as follows. Section 2 is devoted to a short discussion of the Robustness of Magic, giving an alternative definition to the one in the previous section and stating the properties of this resource monotone. Next, a series of techniques is presented which use the symmetries in the definition of the monotone to simplify the computation significantly. To this end, the symmetry group of the stabiliser polytope is characterised in Sec. 3.2 and certain classes of states are singled out in Sec. 3.3 which profit from a high degree of symmetry. For these states, we explicitly derive the symmetry-reduced problem by constructing a suitable basis for the invariant subspace in Sec. 3.3, followed by enumerating equivalence classes of stabiliser states up to symmetry in Sec. 3.4. The numerical solutions for the constructed problems are presented and discussed in Section 4. Based on this, we prove a polytime relaxation of the RoM problem in Sec. 4.3. Our results are summarised in Sec. 5.

Robustness of Magic
The resource theory of magic states can be developed in analogy to the more-established resource theory of entanglement and the robustness of entanglement [38] studied in this context. There, the robustness of a state can be interpreted as a measure for the worstcase separable noise that renders the state separable. However, its construction can be generalised to any resource theory as follows: Given a convex set S of free resources, the robustness of a relative to b ∈ S is defined as Depending on the choice of b, the robustness might be infinite. If it is finite, we can express a as a pseudo-mixture Following Vidal and Tarrach [38], one can define the so-called total robustness by minimising over the set of free resources: In the following, we choose S = SP n to be the convex polytope spanned by the nqubit stabiliser states. More precisely, SP n = conv stab(n), where stab(n) = {s 1 , . . . , s N } is the set of all n-qubit stabiliser states. Here, and in the following, by a "quantum state", we will always mean the density matrix representing it. In the case of pure states s i = |ψ i ψ i |, the associated vector |ψ i will be referred to as a state vector. The polytope SP n is a subset of the real vector space of (D × D)-dimensional Hermitian matrices H D where D = 2 n is the overall dimension of Hilbert space. More specifically, quantum states lie in the (D 2 − 1)-dimensional affine subspace given by tr ρ = 1. Within this affine hyperplane, SP n is full-dimensional and we usually consider it as the the ambient space of SP n .
Howard and Campbell [23] work with an equivalent robustness measure: the robustness of magic (RoM) introduced in Eq. (2). A straightforward calculation (c.f. Appendix A) shows that the two measures are related by a simple affine transformation: The robustness of magic provides a proper resource monotone with the following properties: Proposition 1 (Properties of Robustness of Magic [23]). The robustness of magic has the following properties: 3 Exploiting stabiliser symmetries 3.1 Definition of the RoM problem.
The Robustness of Magic is defined as the following optimisation problem.
Problem 1 (Robustness of Magic). Let stab(n) = {s 1 , . . . , s N } be the set of stabiliser states. Given a state ρ, solve the following problem: Using standard techniques, this problem can be reformulated as a linear program (LP) with D 2 + 2N constraints and 2N variables [3]. Although the time complexity of LPs is linear in the product of number of constraints and variables, these numbers themselves grow super-exponentially with the number of qubits n. Concretely, N = 2 O(n 2 ) and D 2 = 4 n . Moreover, the LP needs access to an oracle which provides the N stabiliser states. The implementation of such an oracle would necessarily have superexponential time complexity itself. However, even if an efficient oracle were provided, the storage of the states would quickly exceed the memory capacity of any computer. In practice, this limits the evaluation of the problem to n ≤ 5 on normal computers and renders it infeasible, even on supercomputers, for n ≥ 8. 1 A standard method in the analysis of optimisation problems is dualising the problem. Clearly, by Slater's condition, strong duality holds and thus the dual problem is an equivalent definition for the Robustness of Magic. In Appendix B, we state the dual problem and derive a lower bound from a feasible solution. However, this bound matches the one that was already found in Ref. [23].

Symmetry reduction
The complexity of the RoM problem can be significantly reduced by exploiting the symmetries of the problem, a procedure that we will call symmetry reduction and is well-known in convex optimisation theory, see e. g. [2]. Here, we will explain the basic ideas and refer the interested reader to App. E for a mathematical review.
By stabiliser symmetries Aut(SP n ), we mean the linear symmetry group of the stabiliser polytope. This is the group of linear maps H D → H D that leave SP n invariant. These maps necessarily have to preserve the set of vertices, i. e. the set of stabiliser states stab(n). Clearly, the group of n-qubit Clifford unitaries C n induces such symmetry transformations by conjugation. Another obvious symmetry of the set of stabilisers is the transposition: where C is the (anti-unitary) operation of complex conjugation in the computational basis. The group of unitary and anti-unitary operations generated by Clifford unitaries and complex conjugation is known as the extended Clifford group EC n [1]. Our first result states that any stabiliser symmetry is induced by the action of an element of the extended Clifford group on the Hilbert space. This is a corollary of the more general Thm. 1 on symmetries of 3-designs and is proven in App. C. We emphasise that this is a non-trivial result which is in general wrong for the case of odd-dimensional qudits where it is possible to construct explicit counter-examples. This turns out to be related to the fact that stabiliser states fail to form 3-designs in odd dimensions [25,40,41].
Note that anti-unitary symmetries in EC n act in the adjoint representation as Ad(C) • T , where C ∈ C n and T is the transposition map. Hence, there are only global antiunitary symmetries. Every tensor product of local antiunitary symmetries would involve a partial transposition and such a map could not preserve the set of entangled stabiliser states.
Let G ρ < EC n be a (not necessarily maximal) subgroup fixing ρ. The projection onto the subspace of G ρ -fixed points V Gρ ⊂ H D , see App. E, is given by Note that Π Gρ is trace-preserving, hence the image of quantum states will again lie in the affine subspace Recall that we can express the robustness of ρ as a minimisation over t ≥ 0 and (mixed) stabiliser states σ ± ∈ SP n such that Since Π Gρ preserves SP n , every such decomposition yields a decomposition in terms of G ρ -invariant mixed stabiliser states: In particular, if the decomposition was optimal in the first place, the projected decomposition is also optimal. This shows that there is always G ρ -invariant optimal solution for the problem. Hence, instead of optimising over the whole set of stabiliser states, we only have to optimise over G ρ -invariant mixed stabiliser states SP n := SP n ∩ V Gρ . By Lemma 3 in App. E, these are exactly given by SP n = Π Gρ (SP n ) and can thus be computed by evaluating the projections stab(n) := Π Gρ (stab(n)). Since Π Gρ (U sU † ) = Π Gρ (s) for all U ∈ G ρ and s ∈ stab(n), it is sufficient to compute the projections on representatives of stab(n)/G ρ . Finally, we remark that a majority of the projected states stab(n) are not extremal points of the projected polytope SP n . Given an extremal subset V n = {v 1 , . . . , v M } ⊂ stab(n), the symmetry-reduced version of Prob. 1 is given by substituting stab(n) → V n and N → M .

Identification of symmetries
The first step towards the explicit symmetry-reduced problem is to identify the group G ρ that fixes the state ρ of interest. Motivated by magic state distillation and the submultiplicativity problem, we are especially interested in the case ρ = |ψ ψ| ⊗n with |ψ being a m-qubit state. A large part of the analysis does not depend on the choice of |ψ , so we keep the discussion as general as possible and specialise later to m = 1 and particular choices of |ψ . The symmetries of |ψ ⊗n can be classified as follows: Permutation symmetry Clearly, |ψ ⊗n is invariant under permutations of the n tensor factors. Such permutations also preserve the stabiliser polytope. Thus, the symmetric group S n is contained in the symmetry group of the problem.
Local symmetries By local symmetries of |ψ ⊗n we mean products of m-qubit stabiliser symmetries of |ψ . By Corollary 1, this class contains only local Clifford operations. Let (C m ) ψ be the stabiliser of |ψ within the m-qubit Clifford group C m , then the local symmetry group is given by (C m ) ⊗n ψ . Global symmetries We refer to all other symmetries as global. The global symmetry group contains e.g. the transposition ρ → ρ T .
The maximal symmetry group for ρ = |ψ ψ| ⊗n is given by the subgroup C ρ that stabilises ρ within EC n . Here, we focus on the subgroup of C ρ which is given by local symmetries and permutations: The following analysis suggests that for our choices of ρ, G ρ actually coincides with C ρ , meaning that there are no further global symmetries. However, since the study of symmetries in EC n can be quite involved [16], we can not exclude the possibility that we missed some of the symmetries. For the rest of this paper, we will consider the case m = 1. Note that C 1 acts by rotating about the symmetry axes of the stabiliser polytope. It is easy to see that states |ψ with non-trivial stabilisers (C 1 ) ψ fall into three classes: Stabiliser states (with trivial robustness), and magic states that lie on the Clifford orbit of |H or |T . Since the RoM is Clifford-invariant, we can pick the following states for concreteness: Figure 5 shows the two states and their stabiliser symmetries. The respective unitary symmetries correspond to a two-fold rotation symmetry about the |H -axis and three-fold rotation symmetry about the |T -axis. In terms of Clifford operations, these stabiliser groups are represented by Recall that these should be understood in the adjoint representation and thus the order of these groups is indeed Furthermore, there are antiunitary stabiliser symmetries such that |H is fixed by A and B and |T is fixed by B and C. Recall that these can only contribute global symmetries such as A ⊗n . However, the common +1 eigenspace of A ⊗n and B ⊗n coincides with that of SX ⊗n and thus adding these symmetries to the symmetry group will not further reduce the invariant subspace. A similar argument holds also for the antiunitary symmetries of |T . Hence, the considered symmetry groups are as follows: Since the symmetric group S n is always a subgroup of the symmetry group, the fixed point subspace V Gρ is always a subspace of the totally symmetric subspace Sym(H D ). Let us first consider a generic state ρ with no further symmetries. Then, V Gρ coincides with Sym(H D ). Thus, the trace 1 subspace has dimension 1 6 (n+3)(n+2)(n+1)−1 and is thus exponentially smaller than the full space. A basis for the symmetric subspace is given by a Fock-style "occupation number basis" constructed from the Pauli basis 1, X, Y, Z as follows Here, the symmetrisation operator Sym ≡ Π Sn is given by averaging over all permutations of the tensor factors. The trace one subspace can be obtained as the span of all basis elements with the N 0,0,0 = 1 component set to 1/D.
Due to linearity, the symmetrisation map is completely determined by its action on the Pauli basis. Given a Pauli operator g, there is a permutation π ∈ S n such . The appearing exponents i = wt X (g), j = wt Y (g) and k = wt Z (g) are exactly the weights of g, i. e. the number of X, Y , Z factors, respectively. By the invariance of Sym under permutations, we thus get Sym(g) = Sym(π(g)) = N i,j,k . We define weight indicator functions, such that we can write the S n -projection of a Pauli operator g as By extending the functions A i,j,k linearly to H D , we thus get exactly the coefficients of the projection in the number basis. Let S < P n be a stabiliser group stabilising a state s. The projection of this state is The A ± i,j,k (S) are the coefficients of the complete signed quantum weight enumerators of the stabiliser code S. Recall that for a classical code C ⊂ F n d , the complete weight enumerator is the degree-n polynomial in d variables given by where wt i (c) gives the number of times i ∈ F d appears in c [26]. The analogy should be clear. Unsigned weight enumerators for quantum codes have been studied since the early days of quantum coding theory [28,Ch. 13]. Much less seems to be known about their signed counterparts, with Refs. [32,33] being the only related references we are aware of. There it is shown that, as their classical analogues, signed quantum weight enumerators are NP-hard to compute. Finally, we want to return to the cases |ψ = |H and |ψ = |T and discuss the invariant subspaces V H,T := V G H,T for these states. Let us rotate the Pauli basis such that the first basis vector corresponds to the Bloch representation of |H and |T , respectively: Note that this choice of basis is such that the orthogonal decompositions of state space irreps of the respective Clifford stabilisers (C 1 ) H,T , as can be seen from the matrix representation of the generators in the rotated basis: In general, a basis for the trivial representation of (C 1 ) ⊗n H,T in the n-qubit state space H } and is constructed analogously to before.
In general, the components of stabiliser states in the rotated bases can be written in terms of weight enumerators by computing the induced basis transformations on However, we are only interested in the projection onto j = k = 0 which simplifies this computation. First, let us rewrite the n-qubit Pauli operators in the H-basis. Note that every operator with non-vanishing Z-weight is already in the orthocomplement of V H .
Here, we left out possible identity factors and all orthogonal terms on the RHS, i. e. those containing E H 2 . This result implies that we can write the projection of a stabiliser state s as We call the numbers B ± i (S) the partial signed quantum weight enumerators of S. The analysis works the same way for the T -projection: In this case, the T -projection of a stabiliser state s with stabiliser group S involves total signed quantum weight enumerators C ± i (S) as follows: Note that all projections Sym ≡ Π Sn , Π H and Π T can be computed from the complete signed weight enumerators of the stabiliser codes which themselves are functions of the weight distributions. For numerical purposes, it is convenient to absorb all appearing factors in the bases such that the coefficients of stabiliser states are given by the integer weight enumerators.
Finally, we want to give expressions for the states |H ⊗n and |T ⊗n in the respective bases: |T T | ⊗n = 1 In general, we are not aware of any method which can predict whether the projection of a stabiliser state will be extremal within the projected polytope. However, the following lemma gives a necessary condition on the extremality of products s ⊗ s of stabiliser states which will be useful later. Proof. We prove the statement by showing it on the level of the complete signed weight enumerators A ± i,j,k . This proves the claim directly for Π = Sym and the other cases follow since the partial and total signed weight enumerators are linear functions of the complete ones.
Note that the Pauli X, Y , Z weights are additive under tensor products, e. g. wt X (g⊗ g ) = wt X (g) + wt X (g ). This implies that we can write the indicator function as are not the weights of g, we can instead sum over all possible decompositions on the right hand side. Hence, for any two stabiliser codes S, S we get Suppose S is the stabiliser of a state s and Sym(s) can be written as convex combination, with stabiliser states s l , stabilised by the groups S l . Let s be stabilised by S , then we find by Eq. (31), (33) and hence the projection of the product state s ⊗ s is non-extremal.
Note that Eq. (31) allows us to compute the projection of products Π(s ⊗ s ) from Π(s) and Π(s ) via the signed quantum weight enumerators using poly(n) operations. This is an important improvement over computing Π(s) for a general (fully entangled) stabiliser state s which requires O(2 n ) operations.

Representatives of inequivalent stabiliser states
Computing the projected polytope involves the computation of the signed quantum weight enumerators for all stabiliser states. However, from the previous discussions we know that we can restrict the computations to the orbits stab(n)/G ρ with respect to the symmetry group G ρ . In this section we will construct representatives for these orbits.
Our approach is based on a subset of the set of stabiliser states, the so-called graph states graph(n). For every simple, i. e. self-loop free, graph G of n vertices, there is a state vector |G that is stabilised by operators of the form where X j , Z j are the Pauli operators on the j-th qubit and θ is the adjacency matrix of the graph G. Graph states play a fundamental role in the studies of stabiliser states since Schlingemann [35] proved that every stabiliser state is equivalent to a graph state under the action of the local Clifford group LC n = C ⊗n 1 : This result can be used to label every stabiliser state vector |C, G by a local Clifford unitary C ∈ LC n and a graph state |G ∈ graph(n) such that |C, G = C |G . However, LC n -equivalent graph states generate the same LC n -orbit and are equally well suited to represent a stabiliser state. Hein, Eisert, and Briegel [19] and Nest, Dehaene, and De Moor [29] discovered that that two graph states are LC n -equivalent if and only if the underlying graphs are related by a graph theoretic transformation called local complementation (LC). Thus, it is sufficient to consider graphs up to local complementation. Furthermore, the symmetry group G ρ induces additional equivalence relations on the graph state representation. Let us again begin the discussion with the case of a generic state with S n -symmetry. This already allows us to restrict the representation to non-isomorphic graphs, i. e. graphs up to permutation of their vertices, since for any graph state |G and a permuted version |πG ≡ π |G the LC n -orbits are isomorphic: πC |G = C π |πG with the permuted local Clifford unitary C π = πCπ † ∈ LC n . Moreover, it is straightforward to show that the composition of graph isomorphism and local complementation is symmetric and thus a equivalence relation ∼ LC,Sn on graphs whose equivalence classes are isomorphic to graph(n)/ ∼ LCn,Sn . These equivalence classes have been studied in the context of graph codes and entanglement in graph states [11,20] and were enumerated by Danielsen [10]. However, different local Clifford unitaries can still result in equivalent states. To see this, pick some symmetry π ∈ Aut(G) of the graph, i. e. πG = G, then the actions of C and C π yield isomorphic states. Hence, it is enough to act with LC n / Aut(G) on the graph state |G .
For the computation of the LC n -orbits it is enough to consider LC n /P n , since Pauli operators will only change the possible 2 n signs of the final generators which are better added by hand. It is well known that the quotient C n /P n is isomorphic to the binary symplectic group Sp(2n, Z 2 ) which is the foundation of the phase space formalism. We make use of this formalism to compute the LC n -orbits of graph states G by evaluating the orbits of the local symplectic group Sp(2, Z 2 ) ×n up to the stabiliser of G and Aut(G).
The additional symmetries in the case of the |H and |T state can be taken into account by restricting the allowed symplectic transformations using the symplectic mapsŜ andŜĤ induced by the generators SX and SH, respectively. The corresponding cosets are given by the representatives Sp(2, Z 2 )/ Ŝ {1,Ĥ,ĤŜ} and Sp(2, Z 2 )/ ŜĤ {1,Ŝ}, respectively. However, the described generation procedure will quickly become computationally expensive. Moreover, most of the projected stabiliser states are non-extremal points for the projected polytope and thus redundant. Unfortunately, there is no simple way  [14,35], c.f. the description in the text. The convex hull of these vertices is shown in Fig. 4.
of deciding whether a state will be extremal after projection or not. However, Lemma 1 states at least a criterion for product states which allows us to restrict to projecting only fully entangled stabiliser states. To this end, we only have to iterate over connected graph representatives with respect to ∼ LC,Sn and compute the projections of product states directly from lower-dimensional vertices using the appropriate version of Eq. (31).

Computing the robustness of magic
Using the enumeration procedure of the last section, we generated the set of H-and T -projections of fully entangled stabiliser states stab H/T c (n) = Π H/T (stab c (n)) and the set of projected product states from lower-dimensional vertices. In an additional step, we removed non-extremal points from the set of projected states, resulting in vertex sets V H/T n of the projected stabiliser polytopes for n ≤ 9 and n ≤ 10, respectively. As described in the last section, we are labelling the vertices by certain stabiliser representatives. To this end, we use a notation in terms of "decorated graph states" compatible with Refs. [14,35]: A graph is decorated by symbols which indicate the action of local Clifford operations on the respective graph state. Nodes with signs indicate a sign change of the respective stabiliser generator, or alternatively, the action of Z on the respective qubit prior the any other gates. A hollow node in the graph denotes a Hadamard gate acting on the respective qubit and self-loops correspond to the action of phase gates (prior to possible Hadamard gates). Figure 6 shows the vertex sets V H n for n = 1, 2, 3. Since the dimension of the polytope is exactly n, it can be easily visualised for n ≤ 3, see also Fig. 4 in Sec. 2.
The database of vertices and the program code can be found on the arXiv [21]. For a discussion of the algorithmic details see App. D. Table 1 shows the number of vertices of the projected polytopes in comparison with the original number of stabiliser states. We see that the number of states N that have to be used in the 1 -minimisation is reduced drastically from 2 O(n 2 ) to a scaling which is approximately 2 n . Additionally, the dimension d of the ambient space is reduced exponentially from 4 n − 1 to exactly n. As discussed in Sec. 3.2, the required 1minimisation for RoM is computed via a linear program with 2N + d constraints and 2N variables and has a runtime that is linear in its size (2N + d)   The runtime is thus reduced as leading to a super-polynomial speed-up in the 1 -minimisation. Although both time and space complexity of the 1 -minimisation are exponential in n, it is in principle feasible for moderate n. Here, the limiting factor is the implementation of the oracle providing the projected states with runtime which is still super-exponential in n. Figure 7 shows the Robustness of Magic of |H ⊗n for n = 1, . . . , 9, computed from the vertices V H n of the projected stabiliser polytope. Note that the data for n ≤ 5 is in perfect agreement with the so-far computed values in Ref. [23]. We are particularly interested in the submultiplicative behaviour of R. Here, the new data for n > 5 turns out to be helpful: We can observe that the data points quickly approach an apparent exponential scaling with n. More precisely, submultiplicativity is clearly observable for 1 ≤ n ≤ 4, but the scaling becomes effectively multiplicative for larger n. We quantified this using an exponential fit of the data range 3 ≤ n ≤ 9 (shown in blue in Fig. 7) resulting in (1.059 ± 0.015) × (1.283 ± 0.002) n . From previous works it is known that the regularised robustness R reg (|H ) is bounded from below by 1.207. Our work, however, indicates that it converges from above to a constant which is given by the fit as (1.283 ± 0.002).

Robustness of the |H ⊗n and |T ⊗n states
The previously known time complexity for simulating a circuit with t T gates using the RoM algorithm is O(2 0.753t ) [23]. Our findings improve this to O(R(|H ⊗9 ) 2t 9 ) = O(1.667 t ) = O(2 0.737t ). Moreover, since we already explored an effectively multiplicative regime of the RoM, solving the problem for higher n > 9 will not much reduce the runtime. From our estimate for the asymptotic regularised robustness, we can estimate the best possible scaling to be 2 0.719t . Furthermore, we applied the same procedure to compute the robustness of the magic state |T ⊗n . Since the T -symmetry group is larger than in the previous case, we were able to compute R(|T ⊗n ) for up to 10 qubits, see Fig. 8. Qualitatively, the results agree very well with those of the last section. Quantitatively, the robustness of the T state is considerably higher than the one of the H state. Using again an exponential fit, we find the scaling (1.169 ± 0.011) × (1.3865 ± 0.0014) n which predicts a regularised robustness of (1.3865 ± 0.0014) n . By the RoM construction, the 10-qubit solution gives rise to a simulation algorithm with runtime O(1.984 m ) = O(2 0.988m ) where m is the total number of |T magic states used, or equivalently, the number of π/12 Z-rotation gates.

Analysis of the optimal solutions
Additionally, we studied the optimal solutions of the 1 -minimisation for the previously discussed cases of |H ⊗n and |T ⊗n . For this purpose, it is instructive to use the original formulation of the robustness of a state ρ in terms of an optimal affine combination of two (mixed) stabiliser states σ ± ∈ SP H,T n , cp. Eq. (5): The states σ ± can be obtained from the optimal solution of the 1 -minimisation ρ = i x * i v i as follows: Recall from the discussion in Sec. 3.2 that replacing every vertex v i in the optimal solution by a stabiliser representative in its preimage Π −1 H,T (v i ) yields an optimal solution for the original problem. Hence, we simply identify the vertices of the projected polytope by their stabiliser representatives constructed in Sec. 3.3. Surprisingly, these states seem to have a rather simple structure, especially the positive contributions σ + . We will discuss the solutions in the following for the H and T case separately. Optimal solutions for the |H ⊗n state The positive contributions σ + to the |H ⊗n state for n = 1, 2, 3 are simply given by the graph state |+ ⊗n . Figure 9 shows the remaining states for n = 4, . . . , 8. Note that these states have to lie on a facet of the polytope to minimise the robustness. But instead of the generic n contributions, they can be written using only log 2 n terms. The vertices themselves are products of |+ and the Bell state |Ψ + .
In contrast, the negative contributions σ − , shown in Fig. 10, have less structure and seem to be partially irregular. Of course, σ − has a non-unique convex combination and thus part of structure could be shadowed by the non-uniqueness. Nevertheless, since the dominant part of the contributions consists of products of |± and the Bell states |Ψ ± , it is reasonable to assume that the σ − can be approximated by Bell states. We suspect that this approximation is quite good, at least for a moderate number of qubits, due to the apparent suppression of vertices with more complex structure.
Motivated by these observations, we define the following polytope: Q H n = conv Π H all n qubits states that are products of |± and |Ψ ± .
By Eq. (31), we can compute these states efficiently from the signed weight enumerators of |± and |Ψ ± . Note that the projection of |+ ⊗ |− is a convex combination of the projected Bell states and thus only states with "all plus" or "all minus" contributions are extremal in Q H n . Let W H n be the set of vertices of Q H n and m = n/2 . We can explicitly enumerate its elements by tuples (i, j, k) ∈ {0, . . . , m} 3 such that i + j + k = m. Every such tuple corresponds to a product of i |Ψ + , j |Ψ − and (2k + n − 2m) |± states. Hence, the number of vertices is We define the approximate robustness of |H ⊗n as the robustness with respect to the polytope Q H n : Since the optimisation is over a subset of all projected stabiliser states, r H n is an upper bound for R(|H ⊗n ). Moreover, it can be efficiently evaluated since both the complexity of computing W H n and of the 1 -minimisation is O(n 4 ). Figure 11 shows a comparison of r H n with the exact robustness. From the previous analysis it is clear that the approximation is exact for n ≤ 4. The deviation from the exact data for 4 < n ≤ 9 is at most 0.06% and thus negligible. However, we expect that the deviation becomes larger the higher n is, since it is likely that the importance of multipartite entangled contributions increases. Nevertheless, the approximation seems to be surprisingly good. The approximate data again follows an exponential increase with n, predicting an asymptotic regularised robustness of about (1.2829 ± 0.0017) which is compatible with the prediction (1.283 ± 0.002) from the exact data.
However, this approach is limited to n ≤ 26. For larger n, the 1 -minimisation lacks a feasible solution, which can only be the case if |H ⊗n is not in the affine span of the product states W H n . This indicates that the dimension of the subpolytope Q H n becomes too small. A solution to these infeasibility problems will be discussed in Sec. 4.3.   Optimal solutions for the |T ⊗n state As in the previous case, the two connected vertices of the projected 2-qubit polytope constitute a dominant part in the optimal solutions. They are not projections of Bell states, so we will denote their representatives by |γ ± and define them to be the states stabilised by respectively. The analysis of the optimal solutions shows that the σ + states are convex combinations of products of |+ and the maximally entangled state |γ + . Moreover, they seem to be even more sparse than for the previous case, see Fig. 12. As in the case of |H , the σ − state shows only partial structure, see Fig. 14 .
The similarities suggest that the robustness for |T ⊗n can be well approximated using a similar procedure as in the last section. To this end, we define the polytope Q T n = conv Π T all n qubits states that are products of |± and |γ ± .
The approximate robustness r T n is again defined with respect to this polytope. The vertices W T n can be efficiently computed using the same procedure as in the |H case and the approximation is exact for n ≤ 3. Figure 13 shows the approximate robustness     compared to the exact results. The approximation is again surprisingly good with a maximum deviation from the exact data of around 0.8%. Although this error is still small, it is an order of magnitude larger than for the |H state. The approximation yields an asymptotic regularised robustness of (1.3916 ± 0.0014) which is slightly larger than the result from the exact data. Similar to the last section, the applicability of this approximation is limited to n ≤ 24 due to the infeasibility of the optimisation problem for larger n. In the next section, we will show how to generalise this approximation to overcome the feasibility problems.

Finite hierarchy of RoM approximations
In general, the idea of restricting to at most k-partite entangled stabiliser states leads to a hierarchy of approximations with levels 1 ≤ k ≤ n. Clearly, for k = n the exact problem is recovered. The set of at most k-partite entangled n-qubit stabiliser states can be constructed by taking all possible tensor products of states in stab(i) for 1 ≤ i ≤ k which result in n-qubit states. However, without the presence of additional symmetries, this will still result in an exponentially large set since already the set of fully separable stabiliser states (k = 1) has size 6 n . Hence, we assume that we want to compute approximations to R(ρ) where ρ is a symmetric n-qubit state (not necessarily pure) such that the stabiliser symmetry group contains at least the symmetric group S n . In particular, this applies to the magic states |H ⊗n and |T ⊗n . In this case, we are able to give poly(n) upper bounds on the runtime for every fixed level k < n.
Following Lemma 1 and Section 3.3, the set of S n -projections of k-partite entangled n-qubit stabiliser states can be constructed from the vertices of the projected polytopes SP i = Sym(SP i ) for 1 ≤ i ≤ k which have fully entangled representatives. Let us denote the sets of representatives by V i ⊂ stab(i). Since the order does not matter, the possible ways to take tensor products of these sets are exactly captured by (descending) partitions of n into parts with size at most k. We will denote such a partition by λ k n. Then, we define the subpolytope of projected k-partite entangled states as Q n,k := conv Sym and the k-th level of the RoM hierarchy by the relaxation of Prob. 1 to the subpolytope Q n,k . Clearly, this defines an upper bound r n,k (ρ) to the exact RoM R(ρ).
To bound the runtime of the k-th level of the hierarchy, we have to count the vertices W n,k of Q n,k . An upper bound to this number is given by the number of tensor products appearing in Eq. (43) up to permutations. Thus, let λ be a (descending) partition of n into r parts, with no part larger than k: This can be rewritten as where 0 ≤ m i ≤ n is the multiplicity of i in the partition λ. Since the permutations of the partition itself were already considered, the number of product states corresponding to the partition λ is given, up to permutations, by Using that the number of fully entangled vertices is increasing with i, we can bound this number by Finally, the number of partitions of n with parts no greater than k coincides with the number of partitions of n into at most k parts and is denoted by p k (n). A standard result in number theory is that Thus, we can bound the number of vertices W n,k to be Since the dimension is O(n 3 ), this implies that the runtime of the relaxation of Problem 1 is polynomial in n for a fixed k. Finally, we remark that one has to know the vertex sets V i up to k to run the k-th level of the hierarchy. Moreover, the bounds are very loose due to the fact we have not strictly bound the number of fully entangled vertices L i which is beyond the scope of this paper. However, by using the actual numbers for L i , one can obtain much better bounds on |W n,k | by evaluating the binomial coefficients. Let us illustrate this for the case of |H ⊗n and k = 2, 3: Using that p 2 (n) = n 2 + 1, p 3 (n) = (n+3) 3 12 Note that we derived |W n,2 | = O(n 2 ) in the previous section using further information about the extremality of products.

Conclusion & Outlook
In this work, we have studied the symmetries of the n-qubit stabiliser polytope and showed how to use these to greatly reduce the combinatorical complexity of computing the robustness of single-qubit magic states and to gain insight into the structure of the problem.
We have determined the symmetry groups for the two types of single-qubit magic states and have constructed explicit stabiliser state representatives of the symmetry orbits. This has allowed us to evaluate the robustness of |H ⊗n for n ≤ 9 and |T ⊗n for n ≤ 10 qubits. Using the structure of the solutions, we have proposed an approximation based on at most bipartite entangled states which is efficient in n and gives an upper bound on the exact robustness. Furthermore, the agreement with the exact data for n ≤ 10 qubits is excellent. Since the RoM becomes effectively multiplicative for larger n, we expect that the approximation is still very good in the regime n > 10. Moreover, by restricting to k-partite entangled stabiliser states, we obtained a finite hierarchy of approximations which recovers the exact RoM for k = n. We showed that a fixed level k < n of the hierarchy can be computed in poly(n) time.
We feel that the most interesting task left open in this work is to explain why even two-body entangled states are sufficient to produce excellent bounds on the RoM. This may be insightful in a wider context. Indeed, sub-additivity of resource costs occurs in several areas of quantum information theory, most famously for the entanglement of formation [18]. The violations to additivity in [18] can be proven to exist for randomised constructions in high dimensions. This makes it hard to study the structure of the optimal solutions, or their behavior in a limit of many copies. The combinatorial nature of the stabiliser polytope, and the observation that only few-body entanglement is enough to find almost-optimal solutions, suggest that RoM may provide an instance where understanding submultiplicativity is feasible.

Acknowledgements
We thank Earl Campbell, Mateus Araújo, Felipe M. Mora, Felix Huber, Frank Vallentin, Arne Heimendahl, and Huangjun Zhu for helpful discussions and comments. In particular, we want to thank Richard Kueng for productive discussions concerning the dual problem.
This work has been supported by the Excellence Initiative of the German Federal and State Governments (Grant ZUK 81), the German Research Foundation (DFG project B01 of CRC 183), and ARO under contract W911NF-14-1-0098 (Quantum Characterization, Verification, and Validation).

A Equivalence of the two robustness measures
The equivalence given in Eq. (6) is stated implicitly in [23]. Here, we give an explicit proof.
Vidal and Tarrach [38] defined the so-called total robustness which is given by For S being a (compact) polytope, this can be rewritten as follows. Since S is compact, the minimum b * is attained. Hence, R(a) = R(a||b * ) =: s * and Let {v 1 , . . . , v N } be the vertices of S and write b + , b * ∈ S as convex combinations with coefficients λ i and µ i . It follows: The last sum is an affine combination of the vertices since i x(s * ) i = 1. In other words, x(s * ) is a feasible solution for the following minimisation problem: Moreover, the optimal value can be bounded as follows: Assume x * is the optimal solution for R(a). Then, we can rewrite R(a), using i x i = 1, as follows: Hence, the optimal affine combination for a becomes Here, the renormalised modulus of the affine coefficients form a convex combination and hence β ± ∈ S. Thus, we found a pseudo-mixture for a and the parameter s(x * ) can not be smaller than the total robustness of a: Combined with Eq. (56), this shows that the two measures are equivalent: Finally, let us remark that β − constructed from the optimal affine combination for a is such that R(a) = R(a||β − ).

B On the dual RoM problem
At this point, any analytical insight could be helpful in simplifying the problem. A standard method is dualising the problem. Clearly, by Slater's condition, strong duality holds and thus the dual problem is an equivalent definition for the Robustness of Magic. The dual problem is straightforwardly obtained as follows: Problem 2 (Dualised Robustness of Magic). Let stab(n) = {s 1 , . . . , s N } be the set of stabiliser states. Given a state ρ, solve the following problem: This formulation of the RoM has a particularly nice form. Thus, it seems at first that the dual problem might be easier to solve. Indeed, one can guess the following feasible solution: Here, {w 1 , . . . , w 4 n } denote the n-qubit Pauli operators which generate the n-qubit Pauli group P n . Feasibility follows from the following calculation for a stabiliser state s with stabiliser group S < P n : The corresponding objective value is where p(ρ) ∈ R D 2 is the coefficient vector of ρ in the Pauli basis, i.e. p(ρ) i = 2 −n tr(ρw i ).
The objective value yields a lower bound to the RoM of ρ. Note that this bound, also called st-norm ρ st , was already found in [23] with different techniques and gives the following lower bound on the RoM of |H ⊗n and |T ⊗n :

C Symmetries of 3-designs
In this section, we characterise the symmery group associated with the projectors of certain t-designs.
A complex projective t-design is a finite family (ψ i ) N i=1 of unit vectors in C d such that where is the orthogonal projection onto the totally symmetric subspace Sym((C d ) ⊗t ). Furthermore, is its dimension and π ∈ S t acts by permuting the factors of the tensor product (C d ) ⊗t . Taking a partical trace of Eq. (67) shows that a t-design is also a t − 1 design.
As in the main part of this paper, we denote by H d the real vector space of Hermitian d×d matrices with the induced Hilbert-Schmidt inner product (A, B) := tr(AB). With respect to this inner product, we denote by L † the adjoint of a linear map L : H d → H d and call L orthogonal if it preserves the inner product, or equivalently, if L † = L −1 .
where U is either a unitary or an antiunitary operator on C d .
hence L is unital.
valid for the non-trivial element π of S 2 , one verifies the following for any traceless Hermitian operator A ∈ H 0 d : tr(Sym [2] A ⊗2 ) = 1 2D [2] tr(A) 2 + tr( 3.-Consider the following trilinear function on H d : F is invariant under L since L † = L −1 is also a symmetry of the projectors ρ i : We can explicitely evaluate F by expanding Sym [3] in terms of permutations and arguing as in Eq. (70). This yields = 1 6D [3] tr A linear automorphism on a matrix algebra fulfilling (77) is called a Jordan automorphism. Our goal is to apply a known structure theorem that restricts that form of such maps [22]. For the theorem to be applicable, we have to extend L from a map on the real vector space of Hermitian matrices, to a map on the algebra M d (C) of all matrices.
To this end, we use that every A ∈ M d (C) can be written uniquely as i. e. the continuationL to M d (C) is a Jordan automorphism. It is also straightforward to check that orthogonality of L implies thatL is unitary with respect to the trace inner product.
It is known that every Jordan automorphism is either an algebra automorphism or algebra anti-automorphism [22]. Since every algebra automorphism is inner andL is unitary,L (and thus also L ≡L| H d ) can in the first case be written asL = U · U † for some U ∈ U (d). In the second case, we can writeL as a compositionL =L • T , wherê L = U · U † is an algebra automorphism and T is the transposition map. For every Hermitian matrix, transposition coincides with complex conjugation as A T = (A † ) * = A * . Hence, we can write L = U C · CU † , where U ∈ U (d) and C is complex conjugation on C d . Hence, L is in this case given by conjugation with the anti-unitary operator U C.
Since the qubit stabiliser state vectors in Hilbert space form a complex projective 3-design [25,40,41], we get the following corollary: Corollary 1. The group of stabiliser symmetries Aut(SP n ) is given by the adjoint representation of the extended Clifford group EC n .
Proof. Theorem 1 implies that every qubit stabiliser symmetry is given by conjugation with either an unitary or anti-unitary operator on the Hilbert space C 2 n . Theorem 2 in [9] implies that every unitary operator that preserves the set of stabiliser states is an element of the Clifford group, up to a global phase.
Furthermore, note that complex conjugation C preserves the set of stabiliser states. Thus, if A is an anti-unitary operator preserving this set, CA is a perserving unitary operator. Hence, up to a phase, CA is Clifford and thus A is anti-Clifford. Finally, this implies our claim that Aut(SP n ) = Ad(EC n ) We note that the result is in general wrong for stabiliser states on odd-dimensional qudits. This also means that the third conclusion of Thm. 1 is not in general true for 2-designs. Concretely, take (ψ i ) i to be the set of stabiliser state vectors for C d , with d a prime number larger than or equal to 5. Then (ψ i ) N i=1 is a 2-design, but the group of linear symmetries of {|ψ i ψ i |} i contains maps that cannot be represented by a linear or anti-linear operator on C d .

Sketch of proof.
We sketch the proof of this claim in the language of [15]. With each a ∈ Z 2 d , one can associate a phase space point operator A(a). The {A(a)} a form a basis for H d . The finite general linear group GL(Z 2 d ) acts on this basis by permuting the indices g A(a) = A(g a). The expansion coefficients W ρ (a) of an operator ρ with respect to the phase space point basis are the Wigner function of the operator. The stabiliser state ρ i = |ψ i ψ i | are exactly the set of Hermitian operators whose Wigner function is the indicator function of an affine line in Z 2 d [15]. Clearly, the GL(Z 2 d )-action introduced above preserves the set of affine lines and thus permutes the ρ i . As argued in the proof of Corollary 1, the group of (anti-)linear operators acting on the state vectors ψ i is the extended Clifford group EC n . To each U in EC n , one can associate a g ∈ Z 2 d such that U A(a)U −1 = A(g a). But g's that arise this way have determinant det [1]. The claim follows, as for d ≥ 5, there are elements g ∈ GL(Z 2 d ) with determinant different from ±1.

D Numerical implementation
Based on the discussion in Sec. 3.3, we can construct a generic algorithm for generating projected stabiliser states by calling various oracles. GraphRepresentatives(n) generates suitable representatives of graph states. Here, these are given by connected representatives of graph(n)/ ∼ LC,Sn which were classified by Danielsen and Parker [11] up to 12 qubits and can by found in Ref. [10]. GeneratorMatrix(G) computes the binary generator matrix of the graph state |G . Furthermore, LocalSymplectic(n, G) returns the set of local symplectic matrices, ideally up to the considered symmetry group. For the discussed cases in Sec. 3.3, this is either the set of direct sums of {1,Ĥ,ĤŜ} or

E.2 Symmetries in convex optimisation
A convex optimisation is the problem of minimising a convex function F over a convex set X . It can always be rewritten in standard form as follows: Let F : R N → R be a convex function and C : R N → R K be a (generalised) convex function with respect to the component-wise partial order on R K , i. e. every component of C is convex. Furthermore, let A : R N → R M be an affine function. The problem is defined as [3] Minimise F (x), for x ∈ R N subject to A(x) = 0, Here, the function F is called the objective function and the functions C and A are the (in-)equality constraints. Depending on the convex set that is modelled, one distinguishes between many subclasses such as linear, conic or semi-definite programming.
We call G a symmetry of the problem (82), if it acts on R N such that the feasible set and the objective function F are left invariant. In particular, this will be the case if G acts linearly on all vector spaces such that the objective function is G-invariant and the constraints are G-equivariant, i. e. for all x ∈ R N and g ∈ G it holds Again, note that the G-action is different on the left and right hand side. Additionally, for G to be a proper symmetry, we require that its representation on R K is given by order automorphisms, i. e.
p q ⇐⇒ g · p g · q ∀p, q ∈ R K , g ∈ G (85) Consequently, both the inequality C(g · x) 0 and the equality constraint A(g · x) = 0 are fulfilled if and only if they hold for x. Hence, x ∈ R N is a feasible solution of Eq. (82) iff its orbit is feasible. Moreover, the objective function is constant on every orbit and thus any optimal solution x * will have an orbit of optimal solutions. The key point for the simplification of the problem is that all functions are convex (A is even affine). Let us again slightly abuse notation and denote with all G-projections on the respective spaces. Using this, we will derive two important consequences of G-equivariance of A and C. First, we evaluate the affine function A: Recall that C is convex w.r.t. to the component-wise order and that every g ∈ G preserves this order. Thus, Π G preserves order, too, and it follows: Suppose x is a feasible solution, then by these relations, its G-projection x G = Π G (x) is feasible, too. Following the same argument as above, we get F (x G ) ≤ F (x). Finally, we find the following results:
Theorem 2 (Symmetry reduction of convex optimisation problems). The convex optimisation problem (82) with symmetry group G is equivalent to the following, symmetryreduced convex optimisation problem: With F G : X G → R, A G : X G → Y G and C G : X G → Z G being functions such that and X G , Y G , Z G being the G-invariant subspace of X = R N , Y = R M and Z = R K .
Proof. First, it should be clear that the functions F G , A G and C G exist and are welldefined by Eq. (90). Moreover, we compute for x, y ∈ X G and t ∈ [0, 1], s ∈ R: Hence, F G and C G are convex and A G as an affine function. Suppose x ∈ R N is a feasible solution of the original problem (82) which can be assumed to be G-invariant, i. e. x ∈ X G . It will be feasible for the reduced problem since and F G (x) = F (x). Next, suppose x G is feasible for the reduced problem. By the same line of argumentation we get due to Eq. (87): In the same fashion, we compute using Eq. (88): Hence, x G is feasible for the original problem and F G ( Finally, this implies that the optimal objective values have to agree: Suppose x * and x G * are (G-invariant) optimal solutions for the original and the reduced problem, respectively. Then, F G (x * ) = F (x * ) and F (x G * ) = F G (x G * ). But since both x * and x G * are feasible for both problems, F (x G * ) = F (x * ) would be a contradiction to the optimality of the solutions.

E.3 Affine constraints and symmetries
In the remainder of this work, both A and C will be affine maps and originate from a set of points V that span a polytope P. The symmetry group G leaves P invariant and hence introduces permutations on V. This will lead naturally lead to G-equivariance of these functions, as we will see in the following.
To simplify the discussion, we will focus on the function A. We can write the affine function A as Here, V := {v 1 , . . . , v N } ⊂ Y are the column vectors of the matrix representing the linear part of A and v 0 is its affine part. Suppose G is represented on Y such that it leaves the set V invariant and fixes v 0 2 . Hence, it can by identified with the left action of some subgroup of the symmetric group S N on the index set [N ] = {1, . . . , N } via g · y i =: y πg(i) for some π g ∈ S N . We can associate a right action on X with this left action by (x · g) i := x π −1 g (i) . This action is clearly linear and such that for all g ∈ G: In particular, the function A is G-equivariant: To make use of Thm. 2, we have to compute the function A G . Note that Π G is constant on the every orbit O ∈ [N ]/G and hence Π G (v j ) =: w O for all j ∈ O: where in the last step we set y O = j∈O x j . Finally, we have to turn this into a map on X G . Note that the right permutation action of G on X = R N partitions the standard basis {e 1 , . . . , e N } into L orbits O 1 , . . . , O L corresponding to [N ]/G. Next, the linear spans X j = O j of these orbits provide a decomposition of X = j X j and G acts transitively on every orbit. Hence, Π G (X j ) is one-dimensional and Π G (X) = j Π G (X j ) due to linearity. This implies that dim X G = L = |[N ]/G|. Hence, the y O are the components of a vector y ∈ X G w.r.t. the basisẽ O = j∈O e j . Note that if we normalise that basis as e O = 1 |O|ẽ O , then the new components arex O = 1 |O| y O , which are exactly the components of Π G (x). Hence, the induced map on X G is As stated in the beginning of this subsection, the points V are the extremal points of a polytope P and G is as subgroup of the polytope symmetries Aut(P). We saw that the symmetry reduction corresponds to projecting the vertices of the polytope, and hence the polytope itself, onto the G-invariant subspace. This is equivalent to taking its intersection with this subspace as the following lemma states: Lemma 3 (Projection with Polytope Symmetries). Be G < Aut(P) a subgroup. Then, the G-projection of P is contained in P, Π G (P) ⊂ P. More precisely, Π G (P) = P ∩X G .
Proof. For all x ∈ P, we have G · x ⊂ P and Π G (x) is a convex combination of points in P, hence in P itself. Moreover, it holds P ∩ X G = Π G (P ∩ X G ) ⊂ Π G (P). The converse direction follows since Π G (P) ⊂ X G and Π G (P) ⊂ P, thus Π G (P) ⊂ P ∩ X G , which shows Π G (P) = P ∩ X G .
Finally, we want to remark that for computing the projection of the vertices {v 1 , . . . , v M }, it is sufficient to compute Π G (w O ) for some representatives w O of the orbits O ∈ V/G since the projection only depends on the orbit.