Robustness of Magic and Symmetries of the Stabiliser Polytope
Markus Heinrich David Gross





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. <cit.>), 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_restrictions_2009 states that there are no quantum codes that allow for a universal set of transversal gates.

In the magic state model <cit.>, 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. <ref>).
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=(1,e^iπ/4), which is realised by the following magic state

    |H⟩ := T|+⟩ = 1/√(2) (|0⟩ + e^iπ/4|1⟩).

Moreover, there is a second magic state, |T⟩, which realises the non-Clifford gate (1,e^iπ/6). Their Bloch representation is shown in Fig. <ref>. 

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. <ref>) can turn many copies of some mixed state ρ into a pure magic state using Clifford unitaries and computational basis measurements <cit.>.



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 <cit.>, 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. <cit.>).
However, as qubits are the paradigmatic systems for quantum computation, quantitative resource theories for multi-qubit magic states have since been developed <cit.>.

The starting point of these theories is the Gottesman-Knill Theorem <cit.>. 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 _n, i.e. 

    ρ = ∑_i p_i s_i,

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 <cit.>, 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. <ref>).



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 non-negative. 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 <cit.> 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. <cit.> and applied to the magic state model in Ref. <cit.>.

More precisely, they define the robustness of magic (RoM) as

    (ρ) := min{x_1  |   x∈^N:  ρ = ∑_i=1^N x_i s_i },

where the sum ranges over stabiliser states {s_1,…,s_N} and  the ℓ_1-norm

    x_1 = ∑_i=1^N |x_i| = 1+ 2∑_i:  x_i ≤ 0 |x_i|

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(ℛ(ρ)^2) <cit.>.

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 <cit.>.
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 <ref>.

Interestingly, the RoM is submultiplicative, i.e. ℛ(ρ^⊗ 2) ≤ℛ(ρ)^2, where the inequality is usually strict <cit.>.
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:

    _reg(ρ):= lim_n→∞(ρ^⊗ n)^1/n.

Unfortunately, computing ℛ(ρ^⊗ n) seems to be a difficult task. For ρ being a single-qubit state, the tensor power ρ^⊗ n lives in an 4^n-dimensional space, and the sum over the s_i in the definition (<ref>) 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 (<ref>) will thus quickly became computationally intractable—and, indeed, howard_application_2017 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. <ref>) 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. (<ref>) 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. <cit.>. 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 <cit.>. 

This paper is organised as follows. Section <ref> 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. <ref> and certain classes of states are singled out in Sec. <ref> 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. <ref>, followed by enumerating equivalence classes of stabiliser states up to symmetry in Sec. <ref>. The numerical solutions for the constructed problems are presented and discussed in Section <ref>. Based on this, we prove a polytime relaxation of the RoM problem in Sec. <ref>. Our results are summarised in Sec. <ref>.







§ 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 <cit.> studied in this context. 
There, the robustness of a state can be interpreted as a measure for the worst-case 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

    R(a||b) := inf{ s ≥ 0  |  1/1+s( a + s b ) ∈ S }.

Depending on the choice of b, the robustness might be infinite. If it is finite, we can express a as a pseudo-mixture 

    a = (1+s)b^+ - s b^-,     with b^±∈ S.

Following vidal_robustness_1999, one can define the so-called total robustness by minimising over the set of free resources:

    R(a) := inf_b ∈ SR(a||b).


In the following, we choose S=_n to be the convex polytope spanned by the n-qubit stabiliser states.
More precisely, _n = (n), where (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 _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 ρ=1. Within this affine hyperplane, _n is full-dimensional and we usually consider it as the the ambient space of _n.


howard_application_2017 work with an equivalent robustness measure: the robustness of magic (RoM) introduced in Eq. (<ref>).




A straightforward calculation (c.f. Appendix <ref>) shows that the two measures are related by a simple affine transformation:

    (ρ) = 1 + 2 R(ρ).

The robustness of magic provides a proper resource monotone with the following  properties:

 The robustness of magic has the following properties:
 
  
  * Faithfulness: (ρ)=1 iff ρ∈_n
  
  * Monotonicity: (𝒳(ρ)) ≤(ρ) for all stabiliser operations 𝒳 with equality if 𝒳 is unitary.
  
  * Convexity: ((1-t)ρ+tσ)≤ (1-t)(ρ)+t (σ) for 0≤ t ≤ 1.
  
  * Submultiplicativity: (ρ⊗σ) ≤(ρ) (σ).
 






§ EXPLOITING STABILISER SYMMETRIES






 §.§ Definition of the RoM problem.



The Robustness of Magic is defined as the following optimisation problem.
[Robustness of Magic]

Let  (n)={s_1,…,s_N} be the set of stabiliser states.
Given a state ρ, solve the following problem:

    min     x_1   over x ∈^N
    s. t.     ρ = ∑_i=1^N x_i s_i.


Using standard techniques, this problem can be reformulated as a linear program (LP) with D^2+2N constraints and 2N variables <cit.>. 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 super-exponential 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.[Already the storage of (7) would require around 77 TiB of memory.  For n=8, this number increases to around 76 PiB which exceeds the state-of-the-art by a factor of 7.]

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 <ref>, 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. <cit.>.




 §.§ 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   <cit.>. 
Here, we will explain the basic ideas and refer the interested reader to App. <ref> for a mathematical review.

By stabiliser symmetries (_n), we mean the linear symmetry group of the stabiliser polytope.
This is the group of linear maps H_D → H_D that leave _n invariant. 
These maps necessarily have to preserve the set of vertices,  the set of stabiliser states (n). 
Clearly, the group of n-qubit Clifford unitaries _n induces such  symmetry transformations by conjugation. 
Another obvious symmetry of the set of stabilisers is the transposition:

    s_i = |ψ_i⟩⟨ψ_i|↦ s_i^T = 𝒞|ψ_i⟩⟨ψ_i|𝒞,

where 𝒞 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 _n <cit.>.
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. <ref> on symmetries of 3-designs and is proven in App. <ref>.

corollarystabsym

 The group of stabiliser symmetries (_n) is given by the adjoint representation of the extended Clifford group _n.


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 <cit.>. 

Note that anti-unitary symmetries in _n act in the adjoint representation as (C)∘ T, where 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 _ρ < _n be a (not necessarily maximal) subgroup fixing ρ.
The projection onto the subspace of _ρ-fixed points V^_ρ⊂ H_D, see App. <ref>, is given by

    Π__ρ(σ) = 1/|_ρ|∑_U∈_ρ Uσ U^†.

Note that Π__ρ is trace-preserving, hence the image of quantum states will again lie in the affine subspace ^-1({1}) ∩ V^G_ρ.

Recall that we can express the robustness of ρ as a minimisation over t≥ 0 and (mixed) stabiliser states σ^±∈_n such that

    ρ = (1+t)σ^+ - t σ^-.

Since Π__ρ preserves _n, every such decomposition yields a decomposition in terms of _ρ-invariant mixed stabiliser states:

    ρ = Π__ρ(ρ) =  (1+t) Π__ρ(σ^+) - t  Π__ρ(σ^-),

In particular, if the decomposition was optimal in the first place, the projected decomposition is also optimal. 

This shows that there is always _ρ-invariant optimal solution for the problem. Hence, instead of optimising over the whole set of stabiliser states, we only have to optimise over _ρ-invariant mixed stabiliser states _n:=_n∩ V^_ρ. By Lemma <ref> in App. <ref>, these are exactly given by _n=Π__ρ(_n) and can thus be computed by evaluating the projections (n):=Π__ρ((n)). Since Π__ρ(U s U^†) = Π__ρ(s) for all U∈_ρ and s∈(n), it is sufficient to compute the projections on representatives of (n)/_ρ. Finally, we remark that a majority of the projected states (n) are not extremal points of the projected polytope _n. Given an extremal subset 𝒱_n={v_1,…,v_M}⊂(n), the symmetry-reduced version of Prob. <ref> is given by substituting (n)↦𝒱_n and N↦ M.






 §.§ Identification of symmetries



The first step towards the explicit symmetry-reduced problem is to identify the group _ρ 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 <ref>, this class contains only local Clifford operations.
 Let (_m)_ψ be the stabiliser of |ψ⟩ within the m-qubit Clifford group _m, then the local symmetry group is given by (_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 _ρ that stabilises ρ within _n. Here, we focus on the subgroup of _ρ which is given by local symmetries and permutations:

    _ρ :=  (_m)^⊗ n_ψ⋊ S_n.

The following analysis suggests that for our choices of ρ, _ρ actually coincides with _ρ, meaning that there are no further global symmetries. However, since the study of symmetries in _n can be quite involved <cit.>, 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 _1 acts by rotating about the symmetry axes of the stabiliser polytope. It is easy to see that states |ψ⟩ with non-trivial stabilisers (_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:

    |H⟩⟨H| = 1/2(  + 1/√(2) ( X + Y ) ),   |T⟩⟨T| = 1/2(  + 1/√(3) ( X + Y + Z ) ).

Figure <ref> 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

    (_1)_H = ⟨ SX ⟩,      (_1)_T = ⟨ SH ⟩.

Recall that these should be understood in the adjoint representation and thus the order of these groups is indeed |(_1)_H|=2 and |(_1)_T|=3. 

Furthermore, there are antiunitary stabiliser symmetries 

    𝒜:      X ↦ X,     ℬ:      X ↦ Y,        𝒞:      X ↦ Z, 
        Y ↦ Y,         Y ↦ X,        Y ↦ Y, 
        Z ↦ -Z,        Z ↦ Z,        Z ↦ X,

such that |H⟩ is fixed by 𝒜 and ℬ and |T⟩ is fixed by ℬ and 𝒞. Recall that these can only contribute global symmetries such as 𝒜^⊗ n. However, the common +1 eigenspace of 𝒜^⊗ n and ℬ^⊗ 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:

    _H := ⟨ SX ⟩^⊗ n⋊ S_n ,     _T := ⟨ SH ⟩^⊗ n⋊ S_n.


Since the symmetric group S_n is always a subgroup of the symmetry group, the  fixed point subspace V^_ρ is always a subspace of the totally symmetric subspace (H_D). Let us first consider a generic state ρ with no further symmetries. Then, V^_ρ coincides with (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 , X, Y, Z as follows

    N_i,j,k = ( X^⊗ i⊗ Y^⊗ j⊗ Z^⊗ k⊗^⊗(n-i-j-k)), 
    for  i,j,k∈{0,…,n}  such that   i+j+k ≤ n.

Here, the symmetrisation operator ≡Π_S_n 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= 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 that π(g)=X^⊗ i⊗ Y^⊗ j⊗ Z^⊗ k⊗^⊗(n-i-j-k). The appearing exponents i=_X(g), j=_Y(g) and k=_Z(g) are exactly the weights of g,  the number of X, Y, Z factors, respectively. By the invariance of  under permutations, we thus get (g)=(π(g))=N_i,j,k. We define weight indicator functions,

    A_i,j,k(g) := 
                        1      if _X = i,  _Y = j,  _Z = k, 
    
                        0      else,

such that we can write the S_n-projection of a Pauli operator g as 

    (g) = ∑_i,j,k A_i,j,k(g) N_i,j,k.

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<_n be a stabiliser group stabilising a state s. The projection of this state is

    (s)    = 1/2^n∑_g∈ S(g) (g) 
       = 1/2^n∑_i,j,k( ∑_g∈ S(g) A_i,j,k(g) ) N_i,j,k
       = 1/2^n∑_i,j,k A^±_i,j,k(S)  N_i,j,k.

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⊂_d^n, the complete weight enumerator is the degree-n polynomial in d variables given by

    ∑_c∈ C x_0^_0(c)… x_d-1^_d-1(c)
    	=:
    	∑_i_1, …, i_d-1 
    	A_i_1, …, i_d-1(C) 
    	x_0^n-(i_1+…+i_d-1) x_1^i_1… x_d-1^i_d-1,

where _i(c) gives the number of times i∈_d appears in c <cit.>.
The analogy should be clear.
Unsigned weight enumerators for quantum codes have been studied since the early days of quantum coding theory <cit.>.
Much less seems to be known about their signed counterparts, with Refs. <cit.> 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^_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:

    E_1^H    := 1/√(2)( X + Y ),     E_1^T    := 1/√(3)( X + Y + Z ), 
    
      E_2^H    := 1/√(2)( X - Y ),     E_2^T    := 1/√(6)( X - 2 Y + Z ) 
    
      E_3^H    := Z,     E_3^T    := 1/√(2)( X - Z ).

Note that this choice of basis is such that the orthogonal decompositions of state space H_2 = ⟨⟩⊕⟨ E_1^H ⟩⊕⟨ E_2^H ⟩⊕⟨ E_3^H ⟩ = ⟨⟩⊕⟨ E_1^T ⟩⊕⟨ E_2^T,E_3^T ⟩ correspond to (real) irreps of the respective Clifford stabilisers (_1)_H,T, as can be seen from the matrix representation of the generators in the rotated basis:

    SX ≃[  1  0  0  0;  0  1  0  0;  0  0 -1  0;  0  0  0 -1 ],    
     SH ≃1/2[     2     0     0     0;     0     2     0     0;     0     0    -1  √(3);     0     0 -√(3)    -1 ].

In general, a basis for the trivial representation of (_1)_H,T^⊗ n in the n-qubit state space H_2^n is given by ℬ^H,T:={, E_1^H,T}^⊗ n. To construct a basis for the full invariant subspace, we have to symmetrise ℬ^H,T resulting in {N^H,T_i,0,0=:N^H,T_i | i=0,…,n}. Here, N^H,T_i,j,k  is the occupation number basis associated to the rotated basis {,E^H,T_1,E^H,T_2,E^H,T_3} 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 (H_D) from N_i,j,k to N^H,T_i,j,k. 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.

    X^⊗ i⊗ Y^⊗ j   = (1/√(2))^i+j(E_1^H+E_2^H)^⊗ i⊗(E_1^H-E_2^H)^⊗ j
       = (1/√(2))^i+j(E_1^H)^⊗ (i+j) + orth. terms

Here, we left out possible identity factors and all orthogonal terms on the RHS,  those containing E_2^H. This result implies that we can write the projection of a stabiliser state s as

    Π_H(s) = 1/2^n∑_i=0^n(∑_j=0^i A^±_i-j,j,0(S)) N^H_i/2^i/2 =: 1/2^n∑_i=0^n B^±_i(S)  N^H_i/2^i/2.

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:

    X^⊗ i⊗ Y^⊗ j⊗ Z^⊗ k   = (E_1^T/√(3)+E_2^T/√(6)+E_3^T/√(2))^⊗ i⊗(E_1^T/√(3)-√(2/3)E_2^T)^⊗ j⊗
               ⊗(E_1^T/√(3)+E_2^T/√(6)-E_3^T/√(2))^⊗ k
       = (1/√(3))^i+j+k(E_1^T)^⊗ (i+j+k) + orth. terms

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:

    Π_T(s) = 1/2^n∑_i=0^n(∑_j=0^i∑_k=0^i-j A^±_i-j-k,j,k(S)) N^T_i/3^i/2 =: 1/2^n∑_i=0^n C^±_i(S) N^T_i/3^i/2.


Note that all projections ≡Π_S_n, Π_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:

    |H⟩⟨H|^⊗ n   = 1/2^n(  + E_1^H )^⊗ n = 1/2^n∑_i=0^n ni N^H_i, 
    |T⟩⟨T|^⊗ n   = 1/2^n(  + E_1^T )^⊗ n = 1/2^n∑_i=0^n ni N^T_i.



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.



The following is true for Π=,Π_H,Π_T: If the projection Π(s) of an arbitrary stabiliser state s is non-extremal, so is Π(s ⊗ s') for any other stabiliser state s'.


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 Π= 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,  _X(g⊗ g') = _X(g) + _X(g'). This implies that we can write the indicator function as A_i,j,k(g⊗ g')=A_i',j',k'(g)A_i-i',j-j',k-k'(g') for i', j', k' being the weights of g. However, since A_i',j',k'(g) is zero if i', j', k' 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

    A^±_i,j,k(S× S')    = ∑_g∈ S,g'∈ S'(g⊗ g') A_i,j,k(g⊗ g') 
       = ∑_i'=0^i∑_j'=0^j∑_k'=0^k∑_g∈ S∑_g'∈ S'(g)(g') A_i',j',k'(g) A_i-i',j-j',k-k'(g') 
       = ∑_i'=0^i∑_j'=0^j∑_k'=0^k A^±_i',j',k'(S) A^±_i-i',j-j',k-k'(S').

Suppose S is the stabiliser of a state s and (s) can be written as convex combination,
 
    (s) = ∑_l=1^M λ_l (s_l)   ⇔   A^±_i,j,k(S) = ∑_l=1^M λ_l A^±_i,j,k(S_l),

 with stabiliser states s_l, stabilised by the groups S_l. Let s' be stabilised by S', then we find by Eq. (<ref>),
 
    A^±_i,j,k(S × S') = ∑_i',j',k'∑_l=1^M λ_l A^±_i',j',k'(S_l) A^±_i-i',j-j',k-k'(S')
       = ∑_l=1^M λ_l A^±_i,j,k(S_l × S'),

 and hence the projection of the product state s⊗ s' is non-extremal.


Note that Eq. (<ref>) allows us to compute the projection of products Π(s⊗ s') from Π(s) and Π(s') via the signed quantum weight enumerators using (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 (n)/_ρ with respect to the symmetry group _ρ. 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 (n). For every simple,  self-loop free, graph G of n vertices, there is a state vector |G⟩ that is stabilised by operators of the form

    K_j = X_j ∏_k=1^n Z_k^θ_jk,      (j=1,…,n),

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_stabilizer_2001 proved that every stabiliser state is equivalent to a graph state under the action of the local Clifford group L_n = _1^⊗ n:

    (n) = L_n ·(n).

This result can be used to label every stabiliser state vector |C,G⟩ by a local Clifford unitary C∈ L_n and a graph state |G⟩∈(n) such that |C,G⟩= C|G⟩. However, L_n-equivalent graph states generate the same L_n-orbit and are equally well suited to represent a stabiliser state. nest_graphical_2004,hein_multiparty_2004 discovered that that two graph states are L_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 _ρ 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,  graphs up to permutation of their vertices, since for any graph state |G⟩ and a permuted version |π G⟩≡π|G⟩ the L_n-orbits are isomorphic: π C|G⟩ = C_π|π G⟩ with the permuted local Clifford unitary C_π = π C π^†∈ L_n. Moreover, it is straightforward to show that the composition of graph isomorphism and local complementation is symmetric and thus a equivalence relation ∼_LC,S_n on graphs whose equivalence classes are isomorphic to  (n)/∼_L_n,S_n. These equivalence classes have been studied in the context of graph codes and entanglement in graph states <cit.> and were enumerated by danielsen_website. However, different local Clifford unitaries can still result in equivalent states. To see this, pick some symmetry π∈(G) of the graph,  π G = G, then the actions of C and C_π yield isomorphic states. Hence, it is enough to act with L_n/(G) on the graph state |G⟩.

For the computation of the L_n-orbits it is enough to consider L_n/_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 _n/_n is isomorphic to the binary symplectic group (2n,_2) which is the foundation of the phase space formalism. We make use of this formalism to compute the L_n-orbits of graph states G by evaluating the orbits of the local symplectic group (2,_2)^× n up to the stabiliser of G and (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 (2,_2)/⟨Ŝ⟩≃{, Ĥ, ĤŜ} and (2,_2)/⟨ŜĤ⟩≃{, Ŝ}, 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 of deciding whether a state will be extremal after projection or not. However, Lemma <ref> 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,S_n and compute the projections of product states directly from lower-dimensional vertices using the appropriate version of Eq. (<ref>). 









§ 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 ^H/T_c(n)=Π_H/T(_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 𝒱^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. <cit.>: 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 <ref> shows the vertex sets 𝒱^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. <ref> in Sec. <ref>. 

The database of vertices and the program code can be found on the arXiv <cit.>. For a discussion of the algorithmic details see App. <ref>.






Table <ref> 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. <ref>, the required ℓ_1-minimisation 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)(2N)=4N^2+2Nd. The runtime is thus reduced as

    2^O(n^2)⟶ 2^O(n),

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. 





 §.§ Robustness of the |H⟩^⊗ n and |T⟩^⊗ n states






Figure <ref> shows the Robustness of Magic of |H⟩^⊗ n for n=1,…,9, computed from the vertices 𝒱_n^H of the projected stabiliser polytope. Note that the data for n≤ 5 is in perfect agreement with the so-far computed values in Ref. <cit.>. We are particularly interested in the submultiplicative behaviour of . 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. <ref>) resulting in (1.059± 0.015)× (1.283±0.002)^n. From previous works it is known that the regularised robustness _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). 




The previously known time complexity for simulating a circuit with t T gates using the RoM algorithm is O(2^0.753t) <cit.>. Our findings improve this to O((|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 (|T⟩^⊗ n) for up to 10 qubits, see Fig. <ref>. 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.988 m) 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 σ^±∈^H,T_n, cp. Eq. (<ref>):

    ρ = 1/2[ ((ρ)+1) σ^+ - ((ρ)-1) σ^- ].

The states σ^± can be obtained from the optimal solution of the ℓ_1-minimisation ρ=∑_i x^*_i v_i as follows:

    σ^+ = 2/(ρ)+1∑_i:  x^*_i > 0 x^*_i v_i,     σ^- = -2/(ρ)-1∑_i:  x^*_i < 0  x^*_i v_i.

Recall from the discussion in Sec. <ref> that replacing every vertex v_i in the optimal solution by a stabiliser representative in its preimage Π_H,T^-1(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. <ref>. 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 <ref> 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. <ref>, 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:

    ^H_n = Π_H( {all  n  qubits states that are products of |±⟩ and |Ψ^±⟩}).

By Eq. (<ref>), 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 ^H_n. Let 𝒲^H_n be the set of vertices of ^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 

    K := |𝒲^H_n| = 2 ∑_i=0^m (m+1-i) = (m+1)(m+2).

We define the approximate robustness of |H⟩^⊗ n as the robustness with respect to the polytope ^H_n:

    r^H_n := min{x_1  |   x∈^K:  ρ = ∑_i=1^K x_i w_i  with  w_i∈𝒲^H_n}.



Since the optimisation is over a subset of all projected stabiliser states, r^H_n is an upper bound for (|H⟩^⊗ n). Moreover, it can be efficiently evaluated since both the complexity of computing 𝒲^H_n and of the ℓ_1-minimisation is O(n^4). Figure <ref> 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 𝒲^H_n. This indicates that the dimension of the subpolytope ^H_n becomes too small. A solution to these infeasibility problems will be discussed in Sec. <ref>.






  
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 {X_1Z_2,Z_1X_2} and {-X_1Z_2,-Z_1Y_2}, 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. <ref>. As in the case of |H⟩, the σ^- state shows only partial structure, see Fig. <ref> . 











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

    ^T_n = Π_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 𝒲^T_n can be efficiently computed using the same procedure as in the |H⟩ case and the approximation is exact for n ≤ 3. Figure <ref> 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 (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 (ρ) 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 <ref> and Section <ref>, the set of S_n-projections of k-partite entangled n-qubit stabiliser states can be constructed from the vertices of the projected polytopes _i=(_i) for 1≤ i≤ k which have fully entangled representatives. Let us denote the sets of representatives by V_i⊂(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 := ( ⋃_λ⊢_k n⊗_i∈λ V_i ),

and the k-th level of the RoM hierarchy by the relaxation of Prob. <ref> to the subpolytope Q_n,k. Clearly, this defines an upper bound r_n,k(ρ) to the exact RoM (ρ).

To bound the runtime of the k-th level of the hierarchy, we have to count the vertices 𝒲_n,k of Q_n,k. An upper bound to this number is given by the number of tensor products appearing in Eq. (<ref>) up to permutations. Thus, let λ be a (descending) partition of n into r parts, with no part larger than k:

    n = λ_1 + … + λ_r,      k≥λ_1 ≥λ_2 ≥…≥λ_r > 0.

This can be rewritten as

    n = ∑_i=1^k m_i   i,

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

    ∏_i=1^k | V_i^m_i/S_m_i| = ∏_i=1^k m_i + L_i - 1L_i - 1,      L_i := |V_i|.

Using that the number of fully entangled vertices is increasing with i, we can bound this number by

    ∏_i=1^k m_i + L_i - 1L_i - 1≤∏_i=1^k m_i^L_i≤ n^k  L_k.

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

    p_k(n) = n^k-1/k!(k-1)! + O(n^k-2).

Thus, we can bound the number of vertices  𝒲_n,k to be

    |𝒲_n,k| ≤ p_k(n)  n^k  L_k = n^k (L_k+1)-1/k!(k-1)! + O(n^k (L_k+1)-2).

Since the dimension is O(n^3), this implies that the runtime of the relaxation of Problem <ref> 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 |𝒲_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+1/2⌋ and L_1=L_2=L_3=2, we find 

    |𝒲_n,2|    ≤( ⌊n/2⌋+1 ) n+2-12-1^2 = O(n^3), 
    
     |𝒲_n,3|    ≤⌊(n+3)^3/12+1/2⌋ n^3 = O(n^6).

Note that we derived |𝒲_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 (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 <cit.>.
The violations to additivity in <cit.> 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).












§ EQUIVALENCE OF THE TWO ROBUSTNESS MEASURES



The equivalence given in Eq. (<ref>) is stated implicitly in <cit.>. 
Here, we give an explicit proof.

vidal_robustness_1999 defined the so-called total robustness which is given by

    R(a) := inf_b ∈ SR(a||b).

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

    b^+ := 1/1+s^*( b^* + s^* a ) ∈ S.

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:

    a = (1+s^*)b^+ - s^* b^* = ∑_i=1^N ((1+s^*)λ_i - s^* μ_i )_=: x(s^*)_i v_i.

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:

    (a) := min{x_1  |   x ∈^N  :   a = ∑_i=1^N x_i v_i  and  1 = ∑_i=1^N x_i}.

Moreover, the optimal value can be bounded as follows:

    (a) ≤∑_i=1^N |x(s^*)_i| ≤ (1+s^*) ∑_i=1^N λ_i + s^* ∑_i=1^N μ_i = 1 + 2s^* = 1 + 2 R(a).


Assume x^* is the optimal solution for (a). Then, we can rewrite (a), using ∑_i x_i =1, as follows:

    (a) = x^*_1 = ∑_i: x^*_i ≥ 0 x^*_i - ∑_i: x^*_i < 0 x^*_i = 1  + 2 s(x^*),     with  s(x^*) := - ∑_i: x^*_i < 0 x^*_i.

Hence, the optimal affine combination for a becomes

    a    = ∑_i: x^*_i ≥ 0 x^*_i v_i - ∑_i: x^*_i < 0 |x^*_i| v_i 
       = (1+s(x^*)) ∑_i: x^*_i ≥ 0x^*_i/1+s(x^*) v_i_=:β^+  -   s(x^*) ∑_i: x^*_i < 0|x^*_i|/s(x^*)_=:β^- v_i.

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:

    R(a) ≤ s(x^*)   ⇔  (a) ≥ 1 + 2 R(a).

Combined with Eq. (<ref>), this shows that the two measures are equivalent:

    (a) = 1 + 2 R(a).


Finally, let us remark that β^- constructed from the optimal affine combination for a is such that 

    R(a) = R(a||β^-).






§ 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:
[Dualised Robustness of Magic]

Let  (n)={s_1,…,s_N} be the set of stabiliser states.
Given a state ρ, solve the following problem:

    max     (ρ Y)   over Y ∈ H_D, 
    s. t.      | (Y s_i) | ≤ 1.


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:

    Y = 1/2^n∑_i=1^4^n((ρ w_i)) w_i.

Here, {w_1,…,w_4^n} denote the n-qubit Pauli operators which generate the n-qubit Pauli group _n.
Feasibility follows from the following calculation for a stabiliser state s with stabiliser group S<_n:

    |(Ys)|    = |1/2^n∑_i=1^4^n((ρ w_i)) (w_i s)| 
       ≤1/4^n∑_i=1^4^n|((ρ w_i))| ∑_g∈ S|(w_i g)| 
       = 1/2^n∑_i=1^4^n ( δ(w_i ∈ S) + δ(-w_i ∈ S) ) 
       ≤1/2^n |S| = 1.

The corresponding objective value is

    (ρ Y) = 1/2^n∑_i=1^4^n((ρ w_i)) (ρ w_i) = ∑_i=1^4^n| (ρ w_i)/2^n| = p(ρ)_1,

where p(ρ) ∈^D^2 is the coefficient vector of ρ in the Pauli basis, i.e. p(ρ)_i = 2^-n(ρ 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 <cit.> with different techniques and gives the following lower bound on the RoM of |H⟩^⊗ n and |T⟩^⊗ n:

    1.207^n ≤(|H⟩^⊗ n),      1.366^n ≤(|T⟩^⊗ n).







§ 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)_i=1^N of unit vectors in ^d such that 

    1/N∑_i=1^N |ψ_i⟩⟨ψ_i|^⊗ t = _[t]/D_[t],

where 

    _[t]:   (^d)^⊗ t⟶ (^d)^⊗ t,     _[t]:=
      1/t!∑_π∈ S_tπ

is the orthogonal projection onto the totally symmetric subspace ((^d)^⊗ t). Furthermore, D_[t]=d+t-1t is its dimension and π∈ S_t acts by permuting the factors of the tensor product (^d)^⊗ t.



Taking a partical trace of Eq. (<ref>) 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):=(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.


 
 
  Let (ψ_i)_i=1^N ⊂^d be a set of unit vectors.
 Let L∈(H_d) be a linear map on Hermitian operators that permutes the projectors (|ψ_i⟩⟨ψ_i|)_i=1^N. 
 
 	
  * If (ψ_i)_i=1^N is a 1-design, then L is unital (i.e. L()=).
 	
  * If (ψ_i)_i=1^N is a 2-design, then L is orthogonal and trace-preserving.
 	
  * If (ψ_i)_i=1^N is a 3-design, then L is of the form L= U· U^†, where U is either a unitary or an antiunitary operator on ^d.




  Define ρ_i:=|ψ_i⟩⟨ψ_i|.

  1.—Using _[1]=, we get
  
    L() = d/N∑_i=1^N L(ρ_i) = d/N∑_i=1^N ρ_i = ,

  hence L is unital.

  2.—Let us define the traceless operators operators f_i:=ρ_i - /d. Using the fact that {ψ_i}_i forms a 2-design, Eq. (<ref>), and the “swap trick”
  
    (A B)=(A⊗ B π)
 
  valid for the non-trivial element π of S_2, one verifies the following for any traceless Hermitian operator A∈ H_d^0:
  
    1/N∑_i=1^N (f_i, A)^2    = 1/N∑_i=1^N [ (ρ_i A) - 1/d(A) ]^2 
       = 1/N∑_i=1^N (ρ_i^⊗ 2 A^⊗ 2) 
       = 1/D_[2](_[2] A^⊗ 2) 
       = 1/2D_[2]( (A)^2 + (A^2) ) 
       = A_2^2/2D_[2].

  In other words, the operators (f_i) form a tight frame for the subspace H_d^0⊂ H_d of traceless Hermitian matrices.
  Moreover, setting f_0 = c with c^2=N(1-d)/d^2+d^3, a similar calculation shows that the set {f_0, …, f_N} forms a tight frame for all of H_d.
  
  By 1., L is unital and thus permutes the tight frame {f_0, …, f_N}. However, any map permuting the elements of a tight frame is orthogonal. Finally, orthogonal and unital maps preserve the trace:
 
    L(A) =  L(A) =  L^†() A =  L^-1() A =   A =  A.

 
 3.—Consider the following trilinear function on H_d:
  

    F(A,B,C) := 1/N∑_i=1^N(A⊗ B⊗ C ρ_i^⊗ 3).

F is invariant under L since L^† = L^-1 is also a symmetry of the projectors ρ_i:

    F(L(A),L(B),L(C))    = 1/N∑_i=1^N( L(A)⊗ L(B)⊗ L(C) ρ_i^⊗ 3) 
       = 1/N∑_i=1^N A⊗ B⊗ C (L^†(ρ_i))^⊗ 3
       = 1/N∑_i=1^N A⊗ B⊗ C ρ_i^⊗ 3
       = F(A,B,C).

We can explicitely evaluate F by expanding _[3] in terms of permutations and arguing as in Eq. (<ref>).
This yields

    F(A,B,C)    = (A⊗ B⊗ C 1/N∑_i=1^Nρ_i^⊗ 3) 
       = 1/D_[3](A⊗ B⊗ C _[3]) 
       =  1/6 D_[3]( (A)(B)(C) + (A)(BC) + 
            + (AB)(C) + (AC)(B) + (ABC) + (BAC) ).

The first four terms are individually L-invariant since L preserves the trace and is orthogonal. 
Hence, the L-invariance of F implies

    (ABC) + (BAC) 
         =  L(A)L(B)L(C) +  L(B) L(A) L(C) 
       =  L^†(L(A)L(B))C +  L^†( L(B) L(A) ) C.






Since this holds ∀ C∈ H_d, we get

    AB + BA    = L^†(L(A)L(B)) + L^†( L(B) L(A) ) 
    ⇔ L(AB + BA)    = L(A)L(B) + L(B)L(A) 
    ⇔ L({A,B})    = {L(A),L(B)}.

A linear automorphism on a matrix algebra fulfilling (<ref>) is called a Jordan automorphism. 
Our goal is to apply a known structure theorem that restricts that form of such maps <cit.>.
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() of all matrices.
To this end, we use that every A∈ M_d() can be written uniquely as A = A_1 + i A_2 where A_1,2∈ H_d, and set

    L̂(A) := L(A_1) + i L(A_2) ∈ M_d().

Clearly, this continuation yields a linear automorphism on M_d(). Morover, since the anticommutator {·,·} is bilinear, we get ∀ A,B∈ M_d():

    L̂({A,B}) = {L̂(A),L̂(B)},

 the continuation L̂ to M_d() is a Jordan automorphism. 
It is also straightforward to check that orthogonality of L implies that L̂ 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 <cit.>. Since every algebra automorphism is inner and L̂ is unitary, L̂ (and thus also L≡L̂|_H_d) can in the first case be written as L̂=U· U^† for some U∈ U(d). In the second case, we can write L̂ as a composition L̂=L̂' ∘ T, where 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𝒞·𝒞U^†, where U∈ U(d) and 𝒞 is complex conjugation on ^d. Hence, L is in this case given by conjugation with the anti-unitary operator U𝒞.


Since the qubit stabiliser state vectors in Hilbert space form a complex projective 3-design <cit.>, we get the following corollary:

*



  Theorem <ref> implies that every qubit stabiliser symmetry is given by conjugation with either an unitary or anti-unitary operator on the Hilbert space ^2^n. 
 
  Theorem 2 in <cit.> 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 𝒞 preserves the set of stabiliser states. Thus, if A is an anti-unitary operator preserving this set, 𝒞A is a perserving unitary operator. Hence, up to a phase, 𝒞A is Clifford and thus A is anti-Clifford. Finally, this implies our claim that (_n)=(_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. <ref> is not in general true for 2-designs.
Concretely, take (ψ_i)_i to be the set of stabiliser state vectors for ^d, with d a prime number larger than or equal to 5.
Then (ψ_i)_i=1^N 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 ^d.


We sketch the proof of this claim in the language of <cit.>.
With each a∈_d^2, 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(_d^2) 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 _d^2 <cit.>.
Clearly, the GL(_d^2)-action introduced above preserves the set of affine lines and thus permutes the ρ_i.
As argued in the proof of Corollary <ref>, the group of (anti-)linear operators acting on the state vectors ψ_i is the extended Clifford group _n.
To each U in _n, one can associate a g∈_d^2 such that U A(a) U^-1 = A(g a). 
But g's that arise this way have determinant g=1 d (if U is unitary) or g = -1 d (if U is anti-unitary) <cit.>. 
The claim follows,

as for d≥ 5, there are elements g∈(_d^2) with determinant different from ± 1.





§ NUMERICAL IMPLEMENTATION






Based on the discussion in Sec. <ref>, we can construct a generic algorithm for generating projected stabiliser states by calling various oracles. GraphRepresentativesn generates suitable representatives of graph states. Here, these are given by connected representatives of (n)/∼_LC,S_n which were classified by danielsen_classification_2006 up to 12 qubits and can by found in Ref. <cit.>. GeneratorMatrixG computes the binary generator matrix of the graph state |G⟩. Furthermore, LocalSymplecticn,G returns the set of local symplectic matrices, ideally up to the considered symmetry group. For the discussed cases in Sec. <ref>, this is either the set of direct sums of {, Ĥ, ĤŜ} or {, Ŝ} up to the symmetry of the graph G. Finally, ProjectStateM',s and ProductStatev_1,…,v_k basically evaluate the weight enumerator formulas (<ref>) and (<ref>).



Furthermore, we use a output-sensitive algorithm by dula_new_1996 to compute the extremal points of the projected stabiliser polytope. This algorithm has time complexity O(d N m) where d is the dimension, N the input size and m the output size,  the number of extremal points. It performs way better than a naive approach since the input size N=O(2^n^2) is much larger than the number of extremal points m=O(2^n). 






§ SYMMETRY REDUCTION OF CONVEX OPTIMISATION PROBLEMS






 §.§ Group projections


The central tool for performing a symmetry reduction with respect to some (finite) group G is the so-called G-projection. Suppose g is represented by ρ: G→(V) on a (real or complex) vector space V, we define a linear map Π_G: V → V, the G-projection, by

    Π_G := 1/|G|∑_g∈ Gρ(g).

The G-projection is well known in the representation theory of finite groups. In the physics literature, it is often called a twirl or twirling operation. Thus, we will also sometimes refer to it as G-twirl. 
For reference, we state some elementary properties of these maps without proof.



 Let Π_G: V → V be a G-projection. Then the following holds:
 
  
  * Π_G is a projection operator. 
  Its image is the subspace V^G of fixed points of G. 
  
  * If V is an inner product space and ρ is an orthogonal/unitary representation, then the projection is orthogonal/unitary.
  
  * For all x∈ V, it holds that
    
    Π_G(x) = 1/|G· x|∑_y ∈ G· x y.

  
  * Π_G is constant on every orbit in V/G
  
  * If G=N ⋊ H, then Π_G = Π_N ∘Π_H =  Π_H  ∘Π_N.
 








 §.§ Symmetries in convex optimisation


A convex optimisation is the problem of minimising a convex function F over a convex set 𝒳. It can always be rewritten in standard form as follows: Let F: ^N→ be a convex function and C: ^N→^K be a (generalised) convex function with respect to the component-wise partial order ≼ on ^K,  every component of C is convex. Furthermore, let A: ^N→^M be an affine function. The problem is defined as <cit.>

    Minimise      F(x),   for  x ∈^N 
    subject to      A(x) = 0, 
          C(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 (<ref>), if it acts on ^N such that the feasible set

    𝒳 = { x ∈^N   |   A(x) = 0,   C(x) ≼ 0 },

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,  for all x∈^N and g∈ G it holds

    F(g· x)    = F(x), 
    
     A(g· x)    = g · A(x), 
    
     C(g· x)    = g · C(x).

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 ^K is given by order automorphisms,  

    p ≼ q   ⟺   g· p ≼ g· q   ∀ p,q∈^K, g∈ G

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∈^N is a feasible solution of Eq. (<ref>) 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

    Π_G = 1/|G|∑_g∈ G g,

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:

    A∘Π_G(x) = 1/|G|∑_g∈ G A∘ g(x) = 1/|G|∑_g∈ G g ∘ A(x) = Π_G ∘ A(x).

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:

    C∘Π_G(x) ≼1/|G|∑_g∈ G C∘ g(x) = 1/|G|∑_g∈ G g ∘ C(x) ≼Π_G ∘ C(x).

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:


 Every G-symmetric convex optimisation problem has a G-invariant optimal solution.



 Be x^* a optimal solution, then Π_G(x^*) is G-invariant, feasible and F(Π_G(x^*)) ≤ F(x^*). Hence, Π_G(x^*) is optimal, too.





 The convex optimisation problem (<ref>) with symmetry group G is equivalent to the following, symmetry-reduced convex optimisation problem: 
 
    Minimise      F^G(x),   for  x ∈ X^G 
    subject to      A^G(x) = 0, 
          C^G(x) ≼ 0.

 With F^G: X^G→, A^G: X^G→ Y^G and C^G: X^G→ Z^G being functions such that
 
    F^G∘Π_G = F,      A^G∘Π_G = Π_G ∘ A,      C^G∘Π_G = Π_G ∘ C,

 and X^G, Y^G, Z^G being the G-invariant subspace of X=^N, Y=^M and Z=^K.



 First, it should be clear that the functions F^G, A^G and C^G exist and are well-defined by Eq. (<ref>). Moreover, we compute for x,y∈ X^G and t∈[0,1],s∈:
 
    F^G( t x + (1-t) y )    = F( t x + (1-t) y ) 
       ≤ t F(x) + (1-t) F(y) 
       = t F^G(x) + (1-t) F^G(y), 
    
       A^G( s x + (1-s) y )    = Π_G∘ A( s x + (1-s) y ) 
       =  s Π_G ∘ A(x) + (1-s) Π_G∘ A(y) 
       = s A^G(x) + (1-s) A^G(y), 
    
       C^G( t x + (1-t) y )    = Π_G∘ C( t x + (1-t) y ) 
       ≼  t Π_G ∘ C(x) + (1-t) Π_G∘ C(y) 
       = t C^G(x) + (1-t) C^G(y).

 Hence, F^G and C^G are convex and A^G as an affine function.

 Suppose x∈^N is a feasible solution of the original problem (<ref>) which can be assumed to be G-invariant,  x∈ X^G. It will be feasible for the reduced problem since
 
    A^G(x)    = A^G∘Π_G(x) = Π_G∘ A(x) = 0, 
    
       C^G(x)    = C^G∘Π_G(x) = Π_G∘ C(x) ≼ 0,
 
 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. (<ref>):
 
    0 = A^G(x^G) = A^G∘Π_G(x) = Π_G∘ A(x) = A∘Π_G(x) = A(x).
 
 In the same fashion, we compute using Eq. (<ref>):
 
    0 ≽ C^G(x^G) = C^G∘Π_G(x) = Π_G∘ C(x) ≽ C∘Π_G(x) = C(x).
 
 Hence, x^G is feasible for the original problem and F^G(x^G)=F^G∘Π_G(x^G) = F(x^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.

 








 §.§ 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 𝒱 that span a polytope 𝒫. The symmetry group G leaves 𝒫 invariant and hence introduces permutations on 𝒱. 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

    A(x) = ∑_i=1^N x_i v_i + v_0

Here, 𝒱:={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 𝒱 invariant and fixes v_0[In general, G-equivariance requires that the action preserves the range of A which is a weaker condition.]. 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_π_g^-1(i). This action is clearly linear and such that for all g∈ G:

    ∑_i=1^N x_i  (g· v_i) = ∑_i=1^N x_i   v_π_g(i) = ∑_i=1^N x_π_g^-1(i)  v_i = ∑_i=1^N (x· g)_i  v_i.

In particular, the function A is G-equivariant:

    g· A(x) = ∑_i=1^N x_i (g· v_i) + g· v_0 =∑_i=1^N (x· g)_i  v_i + v_0 = A(x· g).

To make use of Thm. <ref>, 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:

    Π_G∘ A(x)    = ∑_i=1^N x_i  Π_G(v_i) + v_0 
       = ∑_O∈ [N]/G∑_j∈ O x_j  Π_G(v_j) + v_0 
       = ∑_O∈ [N]/G(∑_j∈ O x_j) w_O + v_0 
       = ∑_O∈ [N]/G y_O  w_O + v_0,

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=^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 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∈ Oe_j. Note that if we normalise that basis as e_O = 1/|O|ẽ_O, then the new components are x̅_O = 1/|O| y_O, which are exactly the components of Π_G(x). Hence, the induced map on X^G is

    A^G(x) = ∑_j=1^L x_j w_j + v_0.


As stated in the beginning of this subsection, the points 𝒱 are the extremal points of a polytope  and G is as subgroup of the polytope symmetries (). 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:



 Be G < () a subgroup. Then, the G-projection of  is contained in , Π_G() ⊂. More precisely, Π_G() = ∩ X^G.



 For all x∈, we have G· x⊂𝒫 and Π_G(x) is a convex combination of points in , hence in  itself. Moreover, it holds ∩ X^G = Π_G(∩ X^G) ⊂Π_G(). The converse direction follows since  Π_G() ⊂ X^G and Π_G() ⊂, thus Π_G()⊂∩ X^G, which shows Π_G() = ∩ 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∈𝒱/G since the projection only depends on the orbit.





