Better bounds on Grothendieck constants of finite orders
Sébastien Designolle Tamás Vértesi Sebastian Pokutta



sebastien.designolle@gmail.com
Zuse-Institut Berlin, 14195 Berlin, Germany

HUN-REN Institute for Nuclear Research, 4001 Debrecen, Hungary

Zuse-Institut Berlin, 14195 Berlin, Germany
Institute of Mathematics, Berlin Institute of Technology, 10587 Berlin, Germany




§ ABSTRACT 
  Grothendieck constants d bound the advantage of d-dimensional strategies over 1-dimensional ones in a specific optimisation task.
  They have applications ranging from approximation algorithms to quantum nonlocality.
  However, apart from d=2, their values are unknown.
  Here, we exploit a recent Frank-Wolfe approach to provide good candidates for lower bounding some of these constants.
  The complete proof relies on solving difficult binary quadratic optimisation problems.
  For d∈{3,4,5}, we construct specific rectangular instances that we can solve to certify better bounds than those previously known; by monotonicity, our lower bounds improve on the state of the art for d≤9.
  For d∈{4,7,8}, we exploit elegant structures to build highly symmetric instances achieving even greater bounds; however, we can only solve them heuristically.
  We also recall the standard relation with violations of Bell inequalities and elaborate on it to interpret generalised Grothendieck constants d2 as the advantage of complex quantum mechanics over real quantum mechanics.
  Motivated by this connection, we also improve the bounds on d2.


Better bounds on Grothendieck constants of finite orders
    Sebastian Pokutta
    5th September 2024
========================================================





§ INTRODUCTION


Published in French in a Brazilian journal, Grothendieck's pioneering work on Banach spaces from 1953 <cit.>, now informally known as his Résumé, has long remained unnoticed.
In 1968, Lindenstrauss and Pełczyński <cit.> discovered it and rephrased the main result, the Grothendieck inequality, which proves a relationship between three fundamental tensor norms through the so-called Grothendieck constant, denoted K_G.
Since then, this far-reaching theorem has found numerous applications <cit.>, in particular in combinatorial optimisation where it is at the heart of an algorithm to approximate the cut-norm of a matrix <cit.>.

Quantum information is another field where this result is influential: following early observations by Tsirelson <cit.>, an explicit connection has been established with the noise robustness in Bell experiments <cit.>.
These experiments aim at exhibiting a fascinating property of quantum mechanics in correlation scenarios: nonlocality <cit.>.
The link with Grothendieck's theorem has then raised a surge of interest for the value of 3, the Grothendieck constant of order three.
Many works have thus demonstrated increasingly precise lower bounds <cit.> and upper bounds <cit.> on its value.
More recently, a numerical method has also been developed to come up with an even more precise (but not provable) estimate of 3 <cit.>.

For Grothendieck constants of higher orders, the link with quantum nonlocality remains <cit.>.
However, the bounds on their values have been less studied and are less tight <cit.>.
There are quite a few difficulties that explain this relative scarcity of results, many of them being manifestations of the curse of dimensionality.
Finding suitable high-dimensional ansätze indeed becomes increasingly hard and resulting instances involve sizes that rapidly become intractable.

In this article, we combine the recent projection technique from <cit.> with the powerful solver developed in <cit.> to obtain better lower bounds on d for 3≤ d≤9.
Following <cit.>, we also consider symmetric structures in high dimensions emerging from highly symmetric line packings <cit.> to suggest even better bounds on 4, 7, and 8, that we unfortunately cannot prove as they involve optimisation problems that we only solve heuristically.
We also consider the generalised Grothendieck constants d2 and interpret them as the advantage of d-dimensional quantum mechanics over real quantum mechanics, a fact that was already studied in <cit.>, but recently received more attention through the powerful results of <cit.>.
Our bounds on d are analytical, while those on d2 strongly rely on numerical methods: we post-process the upper bound on 32 to convert it into an exact result, but the lower bounds remain inaccurate.

We first formally define the constants we want to bound in <ref> before presenting in <ref> our method and main results on lower bounds on d, summarised in <ref>.
Then we turn to generalised constants in <ref>, where we review existing bounds before deriving ours, which are particularly tight on the value of 32.
Finally, we recall the connection with quantum mechanics in <ref> and conclude in <ref>.



§ PRELIMINARIES



Given a real matrix M of size m_1× m_2, we define

    dM=max{∑_x=1^m_1∑_y=1^m_2 M_xya_xb_y | ∀ x∈[m_1], a_x∈ S^d-1, ∀ y∈[m_2], b_y∈ S^d-1},

where [m]={1,…,m} and where S^d-1 is the d-dimensional unit sphere, which reduces to S^0={-1,1} for d=1.
By introducing the set

    dm_1,m_2={X | X_xy=a_xb_y, ∀ x∈[m_1], a_x∈ S^d-1, ∀ y∈[m_2], b_y∈ S^d-1},

we get the equivalent definition

    dM=max{MX | X∈dm_1,m_2}.

The reason for this name comes from the interpretation of the Grothendieck constant that we are about to define in the context of rank-constrained semidefinite programming <cit.>.
In the following, when m_1=m_2, we simply write dm,m=dm; also, when the size m_1,m_2 is either clear from the context of irrelevant, we will use the shorthand notation dm_1,m_2=d.

In essence, the Grothendieck inequality states that there exists a (finite) constant independent of the size of the matrix M that bounds the ratio between the quantities dM for various d.
More formally, given n≤ d, for all real matrices M, we have <cit.>

    dM≤dnnM,

so that the exact definition of these generalised Grothendieck constants is

    dn=sup{dM/nM | M∈^m_1× m_2,   m_1,m_2∈}.

The standard Grothendieck constant of order d is obtained when n=1, in which case we use the shorthand notation d1=d.
The Grothendieck constant of infinite order K_G initially studied in <cit.> corresponds to the limit lim_d→∞d, but it is out of the scope of our work, although we briefly mention it in <ref>.
We refer the reader interested in general aspects and further generalisations of these constants to <cit.>.



§ LOWER BOUNDS ON KG(D)



Given the definition of d in <ref>, any matrix M automatically provides a valid lower bound.
However, there are two main difficulties when looking for good lower bounds.

  Given m_1 and m_2, how to find a matrix M such that the inequality in <ref> is as tight as possible?
  


  Given such a matrix M, how to compute the resulting 1M or at least a close upper bound?
  

<ref> is the most limiting one, as the computation of 1M is equivalent to MaxCut and therefore NP-hard, so that the size of the matrices to consider are limited by the methods and resources available to compute this number.
In general, solving <ref> exactly is out of reach, but good candidates can be found and this suffices to derive bounds.

In the rest of this section, we first recall the method from <cit.> to both problems above, in particular to <ref>.
To illustrate the idea of the method, we give up on solving <ref> for the sake of the elegance of the solution to <ref>, providing instances of remarkable symmetry for which future works may solve <ref> to improve on the bounds on 4, 7, and 8.
We then focus on <ref> and consider rectangular matrices allowing us to obtain certified bounds on 3, 4, and 5 that beat the literature up to d=9, by monotonicity.



 §.§ Obtaining facets of the symmetrised correlation polytope



When n=1, the set 1 is a polytope, called the correlation polytope <cit.>.
In <cit.> the method used to solve <ref> is to start from a point P∈d and to derive M as a hyperplane separating P from 1.
This hyperplane is obtained by solving the projection of P onto the correlation polytope 1 via Frank-Wolfe algorithms.
We refer to <cit.> for a gentle introduction to FW algorithms, to <cit.> for a complete review, and to <cit.> for the details of our implementation.
Note that some accelerations developed in <cit.> and used in this work are now part of the  package <cit.>.

In particular, the symmetrisation described in <cit.> is crucial, but as it depends on the underlying group and its action of the matrices we consider, exposing the full structure of the symmetrised correlation polytope for each group would be tedious.
Instead, we give an brief general definition of this polytope and refer to <cit.> for a detailed example.
Given a group G acting on [m_1] and [m_2] (by means of signed permutations), the symmetrised correlation polytope G(1m_1,m_2) is the convex hull of the averages of all orbits of the vertices of 1m_1,m_2 under the action of G.
Note that this polytope lives in a subspace of the space invariant under the action of G on matrices of size m_1× m_2.

In <ref> we illustrate the various geometrical cases that we encounter in this work.
Their understanding justifies the procedure that we describe in <ref> to derive facets of the symmetrised correlation polytope, that is, separating hyperplanes touching the polytope on a space of codimension one.
Importantly, for symmetric instances, the dimension of the ambient space is strictly smaller than m_1m_2.
In the following, we consistently denote these facets by A, while M will indicate separating hyperplanes without this extra property.
Putative facets will also be denoted by A, that is, in cases where we rely on heuristic methods to compute 1A.

These heuristic methods play an important role throughout the FW algorithm as they quickly give a good direction to make primal progress.
Here we present them in a generalised framework that will turn useful when generalising our algorithm to d2 in <ref>.
At a given step of the FW algorithm minimising the squared distance to the set n, finding the best direction with respect to the current gradient M is exactly the problem in <ref>.
This subroutine is called the Linear Minimisation Oracle (LMO) and, in the course of the algorithm, we usually use an alternating minimisation to obtain heuristic solutions that are enough to make progress <cit.>.

More formally, for all x∈[m_1], we pick a random a^h_x∈ S^n-1 and we compute, for all y∈[m_2],

    b^h_y = _b_y∈ S^n-1∑_y=1^m_2 b_y (∑_x=1^m_1 M_xya^h_x),    that is,     b^h_y = ∑_x=1^m_1 M_xya^h_x/∑_x=1^m_1 M_xya^h_x.

In the unlikely case where the denominator happens to be zero, we set b^h_y to be a predetermined vector, for instance, (1,0,…,0)∈ S^n-1.
Similarly, we then use these b^h_y to compute

    a^h = _a_x∈ S^n-1∑_x=1^m_1 a_x (∑_y=1^m_2 M_xyb^h_y),    that is,     a^h_x = ∑_y=1^m_2 M_xyb^h_y/∑_y=1^m_2 M_xyb^h_y,

and we repeat <ref> until the objective value in <ref> stops increasing, up to numerical precision when n>1.
Depending on the initial choice of a^h_x, the value attained will vary.
Therefore we restart the procedure a large number of times: from a few hundreds or thousands within the FW algorithm, to 10^5 in <ref> and 10^8 in <ref>.

In the following, we clearly mention when the last part of our method — the computation of 1M (or 2M in <ref>) where M is the last gradient return by our FW algorithm — is done with the heuristic procedure just described, which comes with no theoretical guarantee, or with a more elaborate solver for which the exact value can be certified or rigourously upper bounded.







 §.§ Heuristic results with highly symmetric line packings



Starting from a good point P∈d in the procedure described above is key to obtaining good bounds on d.
Following <cit.>, we observed that spanning the d-dimensional projective unit sphere in the most uniform way seems to be favourable, which served as our guiding principle to come up with structures in higher dimensions.
In this <ref>, we use the same distribution of points on both sides in <ref>, that is, m_1=m_2=m and a_1=b_1,… a_m=b_m; the matrix P is therefore a Gram matrix.

The intuitive but vague notion of uniform spreading on the sphere can be formally seen as the problem of finding good line packings <cit.>.
In <ref>, for various symmetric d-dimensional line packings, we give the values of dP/1P and dA/1A.
Since the Gram matrix P is positive semidefinite, the ratio dP/1P is upper bounded by K_G^≽(d), the positive semidefinite Grothendieck constant of order d <cit.>, whose value is known to be γ(d)π/2 <cit.>, where γ(d) is defined in <ref> below.
Note that, for d≥9, this bound is the current best known lower bound on d which is naturally greater than K_G^≽(d).

Among the structures that we present in <ref>, the d-dimensional configurations reaching the best known kissing number in dimension d are of particular interest.
In particular, in dimension d=3, different realisations of the kissing number give rise to different values of 3A/1A.
As the icosahedron is an optimal line packing, this indicates that giving good bounds on d does not, in general, boil down to finding good line packings.
Note that the unicity of the configuration in dimension d=4 was recently proven <cit.>.

Upon inspection of <ref> it is apparent that, in the small dimensions that we consider in this work (d≤8), root systems play an important role in obtaining good kissing configurations.
These structures are inherently symmetric and play an important role in the theory of Lie groups and Lie algebras.
We refer to <cit.> for definitions and elementary properties.
The family D_d has been studied in <cit.>; there, starting from the Gram matrix P of size d(d-1), the exact value of 1P is computed by invoking symmetry arguments.
Moreover, for d∈{3, 4, 5}, the optimal diagonal modification of λ=2/3 is derived for D_d; in other words, <cit.> shows that, among all matrices with shape P-λ𝕀, the matrix P-2/3𝕀 is giving the best lower bound on d.
Here we go a bit further by extending this result to d∈{6, 7, 8}, a fact that was suspected in <cit.>, but we also observe that the corresponding matrices A=P-2/3𝕀 are actually facets of the symmetrised polytope D_d(1d(d-1)).

Moreover, we computationnally establish optimal diagonal modifications for other configurations, see <ref>.
About the root system E_8, we note that the value 45/31 reached without diagonal modification is already provided in <cit.>, where it is attributed to Reeds and Sloane.
It indeed follows from <cit.> and from the transitivity of the Weyl group, see the acknowledgements in <ref>.
However, these arguments do not directly apply to solve the diagonally modified case.

Note that in all cases involving irrational numbers (except the icosahedron), the facet A cannot be obtained via a diagonal modification of P.
This is because the off-diagonal elements of P feature rational and irrational numbers.
The exception of the icosahedron is due to its extra property of being an equiangular tight frame (ETF), all off-diagonal elements being ±1/√(5), which can then lead to a facet by taking an irrational diagonal modification.



The 600-cell and the 120-cell have already been studied in <cit.>, where symmetry arguments are exploited to compute the value of 1P, even for the 300×300 matrix arising from the 120-cell.
Interestingly, although the difference of value for 4P/1P is almost negligible, the facets that we obtain here “activate” the advantage of the 120-cell.
Also, we emphasise the advantage of our geometric approach for these structures: in <cit.>, the optimal diagonal modification of 23/9 was indeed derived for the 600-cell, giving rise the the value of 35(-37+27√(5))/569≈1.4378, which is significantly smaller than our result of about 1.4740, see <ref>.
This is because we are not restricted to diagonal modifications in our algorithm, which gives a significant superiority for these configurations.

For some instances of 1P and 1A, the size of the matrix is too large for numerical solvers to handle it.
The values presented with an asterisk in <ref> and summarised in <ref> are then obtained heuristically and cannot be considered final results until these values are confirmed to be optimal.
For this, the high symmetry of these instances, detrimental for our branch-and-bound algorithm and not exploited in available QUBO/MaxCut solvers, could play a crucial role.
We leave this open for further research.



 §.§ Exact results in asymmetric scenarios



Now, we discuss how to address <ref>, namely, the computation of

    1M=max_a_x=±1b_y=±1∑_x=1^m_1∑_y=1^m_2 M_xya_xb_y

In <cit.> this problem is reformulated into a Quadratic Unconstrained Binary Optimisation (QUBO) instance, which is then given to the solver QuBowl from <cit.>.
The instance solved there involves m_1=m_2=97 and is bigger than the ones from <cit.> (of maximal size m_1=m_2=92).
However, it is worth noticing that the solution of the exact instance solved in <cit.> stands out in the landscape of all binary variables.
More precisely, the exact solution is likely to be unique and its value is far above the other feasible points, so that the solver could efficiently exclude vast portions of the search space.
This empirical observation is supported by two facts: instances from later stages of convergence of the Frank-Wolfe method cannot be solved by QuBowl (even within ten times as much time), and the symmetric instance A of size 63×63 obtained from the E_7 root system (see <ref>) cannot be solved by QuBowl either, although it is way smaller.

Given these difficulties, we consider other formulations exploiting the specificities of the problem.
In particular, the branch-and-bound algorithm from <cit.> works by breaking the symmetry between the binary variables a_x and b_y in <ref>, fixing the latter to obtain:

    1M=max_a_x=±1∑_y=1^m_2|∑_x=1^m_1 M_xya_x|,

which has half as many variables.
Importantly, the parameter m_2 plays a less critical role in the complexity of the resulting algorithm, which is indeed still exponential in m_1, but now linear in m_2.
This suggests to use rectangular matrices M in our procedure, with a large m_2 to increase the achievable lower bound on d and a relatively small m_1 to maintain the possibility to solve <ref>, that is, to compute 1M.

Such rectangular matrices can be obtained by following the procedure described in <ref> starting from rectangular matrices P∈dm_1,m_2.
Good starting points are still obtained by using well spread distribution on the sphere, but these distributions a_1,…,a_m_1 and b_1,…,b_m_2 are now different, so that the matrix P with entries a_xb_y is not a Gram matrix any more.
We give our best lower bounds on d in <ref> and present relevant details below.
All matrices can be found in the supplementary files [Supplementary files can be found on Zenodo: <https://doi.org/10.5281/zenodo.13693164>.].

In dimension d=3, despite our efforts, we could not find a rectangular instance beating the 97×97 from <cit.>.
However, the quantum point provided therein, namely, the polyhedron on the Bloch sphere consisting of 97 pairs of antipodal points, is not the best quantum strategy to violate the inequality.
By using the Lasserre hierarchy at its first level, we can indeed obtain a slightly better violation, hence a tighter lower bound on 3.
Summarising, with the matrix M from <cit.> that had been constructed starting from a Gram matrix P∈397, we have:

  
  * 3M≥MP≈2.0000×10^22,
  
  * 3M≥MP'≈2.0001×10^22 where P'∈397 is obtained at the first level of the Lasserre hierarchy (rational expression in the supplementary file <cit.>),
  
  * 1M=13921227005628453160441 (solved with QuBowl <cit.>).


In dimension d=4, we start with the 600-cell on one side and use the compound formed by the same 600-cell together with its dual (the 120-cell) on the other side.
This creates a matrix P∈460,360 we can run our separation procedure on.
Similarly to <cit.>, we exploit the symmetry of P throughout the algorithm to reduce the dimension of the space and accelerate the convergence.
The resulting matrix A already has integer coefficients and satisfies:

  
  * 4A≥AP=30(227668+322725√(2)+170064√(5)+182375√(10)),
  
  * the first level of the Lasserre hierarchy numerically confirms this value as optimal,
  
  * 1A=33135128 (solved with the branch-and-bound algorithm from <cit.>).


In dimension d=5, we make use of structures in http://neilsloane.com/grass/Sloane's database <cit.>.
On one side, we directly use their five-dimensional structure with 65 lines.
On the other side, we take their five-dimensional structure with 37 lines and augment it by adding the center of all edges, properly renormalised to lie on the sphere, which creates a five-dimensional structure with 385 lines.
The resulting matrix P∈565,385 goes through our separation algorithm until a satisfactory precision is attained.
After rounding we obtain a matrix M satisfying:

  
  * 5M≥MP≈2.1061×10^8,
  
  * 5M≥MP'≈2.1068×10^8 where P'∈565,385 is obtained at the first level of the Lasserre hierarchy (rational expression in the supplementary files <cit.>),
  
  * 1M=141074623 (solved with the branch-and-bound algorithm from <cit.>).






§ BOUNDS ON KG(D->2)



We now turn to the generalised Grothendieck constants d2, which have not received a lot of attention so far.
Motivated by the quantum interpretation of these constants (see <ref>), we extend the techniques presented above and in previous works to refine the bounds known on their values.



 §.§ Bounds arising from the literature


There is an explicit lower bound on dn for any d>n that uses a particular one-parameter family of linear functionals of infinite size.
This bound, initially proven in <cit.> (see also <cit.> for the proof), reads

    dn≥γ(d)/γ(n),  with  γ(d)=2/d(Γ(d+1/2)/Γ(d/2))^2,

where Γ(z)=∫_0^∞ t^z-1^-t t is the gamma function.
When n=2, γ(n)=π/4 and <ref> can be rewritten:

    d2≥dd-1k^2/2^2d-3  when   d=2k    and    d2≥2^2d+1/dd-1k^2π^2  when   d=2k+1.


The bounds in <ref> are better than indirect bounds derived from natural combinations of bounds on the Grothendieck constants d, even when considering the improved bounds proven in this article.
Starting with the Grothendieck inequality (<ref>), dividing both sides with 1M and lower and upper bounding the left and right sides, we indeed obtain the following relation

    d≤dnn.

Then, inserting 2=√(2) from <cit.> and using our lower bounds exposed in <ref> yields 32≥1.0159, 42≥1.0506, and 52≥1.0560, which does not improve on the values from <ref>.

As can be seen in <ref>, <ref> also outperforms all known bounds obtained with finite matrices.
Therefore, to the best of our knowledge, <ref> represents the state of the art.



As far as we authors know, there are no nontrivial upper bounds on dn readily available in the literature.
The higher complexity of finding such bounds is partly explained by the natural direction imposed by the definition of dn as a supremum, see <ref>.
In <ref> below, we give a general algorithmic method to overcome this difficulty.



 §.§ Numerical lower bounds on KG(d->2)



Similarly to <ref>, providing a matrix M together with an upper bound on 2M and a lower bound on dM suffices to obtain a lower bound on d2.
However, solving the optimisation problem 2M is even harder than 1M: it is also nonlinear and nonconvex, but with continuous variables instead of binary ones.
In this section we explain how we obtain numerical upper bounds on 2M and summarise our results in <ref>.

We first rewrite <ref> by introducing the (real) coordinates of the d-dimensional vectors a_x and b_y:

    dM=max{∑_x=1^m_1∑_y=1^m_2 M_xy[ a_x,1;     ⋮; a_x,d ][ b_y,1;     ⋮; b_y,d ] | ∀ x∈[m_1], ∑_i=1^da_x,i^2=1, ∀ y∈[m_2], ∑_i=1^db_y,i^2=1}.

With this reformulation, we see that obtaining upper bounds on dM amounts to solving the Lasserre hierarchy of this polynomial system with d(m_1+m_2) variables and m_1+m_2 constraints at a certain level <cit.>.
At the first level, this method is used in <ref> to numerically confirm the value of dA, as the upper bound computed matches, up to numerical accuracy, our analytical lower bound.
In the following, we use it at the second level and for n=2 to numerically obtain upper bounds on 2M.

Similarly to <ref>, the different values that we reach are given in <ref>.
The solution to the Lasserre hierarchy was obtained with the implementation in <cit.>.
Importantly, these results suffer from numerical imprecision.
This is all the more relevant here because of the huge size of most instances, making the use of a first-order solver, COSMO <cit.> in our case, almost mandatory, which is detrimental to the precision of the optimum returned.
This explains the deviations observed in <ref>: the fourth digit of the bound computed via the Lasserre hierarchy is not significant.

Interestingly, the icosahedron now provides us with a better bound than the cuboctahedron, whereas this was the opposite in <ref>.
The dodecahedron also gives no advantage over the icosahedron, up to numerical precision at least.
Similarly, the 600-cell and the 120-cell give very close bounds, which is in strong contrast with <ref>.

Note that the bounds presented in <ref> allow us to prove that <ref> is strict for n=2 and d=3, a fact that could not be verified with <ref>.
In this case, we indeed have that 3≤1.455 from <cit.> and that 232≥√(2)×1.103≈1.560 from <cit.> and with the bound given in <ref>.





 §.§ Exact upper bound on KG(3->2)



Here we reformulate the procedure first introduced in <cit.> and later used to obtain better upper bounds on 3 in <cit.> in a way that makes the generalisation to other Grothendieck constants more transparent.


  Let P∈dm_1,m_2 with underlying vectors a_1,…,a_m_1,b_1,…,b_m_2∈ S^d-1 such that there exist η_A,η_B∈[0,1] for which η_AS^d-1⊂{a_x}_x=1^m_1 and η_BS^d-1⊂{b_y}_y=1^m_2.
  If α is such that α P∈nm_1,m_2, then
  
    dn≤1/α η_Aη_B.

  



  Consider a matrix N∈dm_1,m_2 and observe that the assumptions on P imply that α η_Aη_B N∈nm_1,m_2.
  Then, for any matrix M,
  
    α η_Aη_B dM=max_N∈dm_1,m_2Mα η_Aη_B N≤nM,

  so that <ref> follows from the definition of dn in <ref>.


To obtain our bound on 32, we start from the Gram matrix P obtained from a centrosymmetric polyhedron with 912 vertices, that is, 406 lines.
Note that this polyhedron is different from the one used in <cit.> and features a slightly higher shrinking factor through the following construction: starting from an icosahedron, apply four times the transformation adding the normalised centers of each (triangular) facet.
Taking v_0=0.8962 and running our Frank-Wolfe algorithm, we obtain a decomposition using 7886 extreme points of 2406 and recovering v_0P up to a Euclidean distance of ε≈2.7×10^-4.
This algorithm is strictly identical to the one in <cit.> except for the LMO, which corresponds to the case decribed in <ref> when n=2.
Following <cit.> we can make this decomposition purely analytical by noting that the unit ball for the Euclidean norm is contained in 1 and therefore in 2.
Note that, contrary to <cit.> where the extremal points of 1 were automatically analytical as they only have ±1 elements, an extra step is needed here to make numerical extremal points of 2 analytical.
Here we follow a rationalisation procedure very similar to the one described in <cit.> and we obtain a value α=v_0/(1+ε) that can be used in <ref> to derive the bound

    32≤1.1233…

whose analytical expression is provided in the supplementary files <cit.>.
This upper bound is indeed entirely rigourous, contrary to the lower bounds presented above in <ref>, which rely on numerical methods that we did not convert into exact results.



 §.§ Implications in Bell nonlocality



The connection between d and Bell nonlocality has first been established by Tsirelson <cit.> and was later investigated more thoroughly in <cit.>.
We refer to the appendix of <cit.> for a detailed construction of the quantum realisation of the Bell inequalities constructed in our work.
Note that the new lower bound mentioned in <ref> slightly improves on the upper bound on v_c^Wer=1/3, namely, it reduces it from about 0.69606 to 0.69604; see also <cit.> for the definition of this number and its interpretation as a noise robustness.

Works anticipating the connection with d2 can be traced back to 2008 <cit.>, when Vértesi and Pál studied the realisation of Bell inequalities with real quantum mechanics and the advantage of d-dimensional quantum systems over real ones.
This continued with the work of Briët et al. <cit.> but gained traction when the powerful result by Renou et al. <cit.> showed that complex numbers are necessary for the complex formalism.

We refer to <cit.> for the explicit connection; for this, the examination of <ref> should be of some help as it gives the bridge between their notations and ours.



§ DISCUSSION



We generalised the method used in <cit.> to obtain lower bounds on 3 to Grothendieck constant of higher order d.
We obtained certified lower bounds on 3, 4, and 5 beating previous ones, and setting, by monotonicity, the best known lower bounds on d for 3≤ d≤9.
For these constants of higher orders, we also showed how our Frank-Wolfe algorithm can be used to derive putative facets, but solving the corresponding quadratic binary problems remains open.
The instances given in this article would allow immediate improvement on 4, 7, and 8.
We expect future work to exploit their symmetry to make the computation possible.

Beyond the perhaps anecdotal improvement brought to their respective constants, developing such tools could unlock access to higher-dimensional structures with a size way larger than anything that numerical methods could ever consider solving.
This could in turn give rise to lower bounds on Grothendieck constants of finite order strong enough to compete with the best known lower bound on the Grothendieck constant (of infinite order), namely, the one presented by Davie in 1984 <cit.> and independently rediscovered by Reeds in 1991 <cit.>, both works being unpublished.
This bound reads

    K_G≥sup_0<λ<11-ρ(λ)/max{ρ(λ),F_ρ(λ)(λ)}≈1.676956674…,
        attained for  λ≈0.255730213…

where

    ρ(λ)=√(2/π)λ^-λ^2/2  and  
      F_ρ(λ)=2/π^-λ^2+ρ(1-2√(2/π)∫_λ^∞^-x^2/2 x).


We also extended the procedure of <cit.> to derive bounds on the generalised constant d2.
These constants appear quite naturally when considering the advantage of quantum mechanics over real quantum mechanics.
In particular, our analytical upper bound on 32 delivers some quantitative insight on this question, while our numerical lower bounds on d2 could turn useful when looking for estimates of the benefit of considering higher-dimensional quantum systems.

We also note that our techniques naturally apply in the complex case, for which similar Grothendieck constants have also been defined.
However, since the complex Grothendieck constants of finite order do not have known interpretations, we refrained from deriving bounds on them, although our code is already capable of dealing with this case.



§ CODE AVAILABILITY


The code used to obtain the results presented in this article is part of the  package  introduced in <cit.> and based on the  package <cit.>.
While the algorithm used to derive bounds on d is directly available in the standard version of this package, the extension to generalised constants dn can be found on the branch  of this repository.
Installing this branch can be directly done within  package manager by typing .



§ ACKNOWLEDGEMENTS



The authors thank
Péter Diviánszky for fruitful discussions and his help to run the branch-and-bound algorithm from <cit.> on asymmetric instances;
Jim Reeds and Neil Sloane for sharing the details of the derivation of the value 45/31 for the E_8 root system mentioned in <cit.>;
Timotej Hrga, Nathan Krislock, Thi Thai Le, and Daniel Rehfeldt for trying to solve our symmetric instances with various solvers;
Mathieu Besançon, Daniel Brosch, Wojciech Bruzda, Dardo Goyeneche, Volker Kaibel, David de Laat, Frank Valentin, and Karol Życzkowski for inspiring discussions.
This research was partially supported by the DFG Cluster of Excellence MATH+ (EXC-2046/1, project id 390685689) funded by the Deutsche Forschungsgemeinschaft (DFG).
T. V. acknowledges the support of the European Union (QuantERA eDICT) and the National Research, Development and Innovation Office NKFIH (Grants No. 2019-2.1.7-ERA-NET-2020-00003 and No. K145927).


sd2

