Vernon Austel
email: austel@us.ibm.com

Sanjeeb Dash
email: sanjeebd@us.ibm.com

Oktay Gunluk
email: gunluk@us.ibm.com

Lior Horesh
email: lhoresh@us.ibm.com

Leo Liberti
email: liberti@polytechnique.fr

Giacomo Nannicini
email: nannicini@us.ibm.com

Baruch Schieber

Globally Optimal Symbolic Regression

In this study we introduce a new technique for symbolic regression that guarantees global optimality. This is achieved by formulating a mixed integer non-linear program (MINLP) whose solution is a symbolic mathematical expression of minimum complexity that explains the observations. We demonstrate our approach by rediscovering Kepler's law on planetary motion using exoplanet data and Galileo's pendulum periodicity equation using experimental data.

Introduction

Discovering mathematical models that explain the behavior of a system has a broad range of applications. Symbolic regression is an important field in machine learning whose goal is to find a symbolic mathematical expression that explains a dependent variable in terms of a number of independent variables for a given data set, most commonly as an explicit function of the independent variables. Unlike traditional (numerical) regression schemes, the functional form of the expression is not assumed to be known a priori [START_REF] Kroll | Workflow to set up substantial target-oriented mechanistic process models in bioprocess engineering[END_REF]. The utility of the approach has been established for a broad range of applications [CT77, WHM + 97, DSW99, DSW03], including the discovery of equations [START_REF] Todorovski | Declarative bias in equation discovery[END_REF][START_REF] Langley | Data-driven discovery of physical laws[END_REF][START_REF] Schmidt | Distilling free-form natural laws from experimental data[END_REF][START_REF] Schmidt | Symbolic regression of implicit equations[END_REF].

Starting with [START_REF] Koza | Genetic Programming: On the Programming of Computers by Means of Natural Selection[END_REF], symbolic regression problems have been typically solved with genetic programming [START_REF] Banzhaf | Genetic programming: an introduction[END_REF][START_REF] Adriano | Symbolic regression via genetic programming[END_REF], an evolutionary metaheuristic related to genetic algorithms; another such technique is grammatical evolution [?]. Even before [START_REF] Koza | Genetic Programming: On the Programming of Computers by Means of Natural Selection[END_REF], heuristics to find explicit functional relationships were developed as part of the BACON system, see [START_REF] Simon | Heuristics for empirical discovery[END_REF]. There has been much research into improving symbolic regression techniques. In [START_REF] Jiang | A divide and conquer method for symbolic regression[END_REF], the separability of the desired functional form is exploited to speed up symbolic regression, whereas a set of candidate basis functions is used in [?] for this purpose. Other recent work has focused on finding accurate constants in the derived symbolic mathematical expression [?]. To discover meaningful functions in the context of physical systems [START_REF] Schmidt | Distilling free-form natural laws from experimental data[END_REF] search for functions that not only match the data, but also have the property that partial derivatives also match the empirically computed partial derivatives. The same authors search for implicit functional relationships in [START_REF] Schmidt | Symbolic regression of implicit equations[END_REF]. Another approach populates a large hypothesis space of functional forms, on which sparse selection is applied [START_REF] Steven L Brunton | Discovering governing equations from data by sparse identification of nonlinear dynamical systems[END_REF].

In this study we present a Mixed-Integer Non-Linear Programming (MINLP) formulation that produces the simplest symbolic mathematical expression that satisfies a prescribed upper bound on the prediction error. Our MINLP formulation can be solved to optimality using existing state-of-theart global optimization solvers. The key advantage of our approach is that it produces a globally optimal mathematical expression while avoiding exhaustive search of the solution space. Another advantage is that it produces correct real-valued constants (within a tolerance) directly; most other methods use specialized algorithms to refine constants [?, ?] and cannot guarantee global optimality. In addition, our formulation can, in principle, seamlessly incorporate additional constraints and objectives to capture domain knowledge and user preferences.

A MINLP Formulation for Free-Form Model Discovery

A symbolic regression scheme consists of a space of valid mathematical expressions together with a mechanism for its exploration. A mathematical expression can be represented by a (rooted) expression tree where each node is labeled by one of the following entities: operators (such as +, -, ×, log), variables (i.e., the independent variables), and constants. The edges of the tree link these entities in a way that is consistent with a prescribed grammar. For example, an expression tree for 2(w 1 +w 2 ) 3 +1 is:

+ root 1 × pow(•, •) 3 + w1 w2 2
Here +, × and the power function are the operators, w 1 , w 2 are the variables, and numbers 1, 2 and 3 are the constants.

In our MINLP based approach we model the grammar of valid expression trees by a set of constraints. The discrete variables of the formulation are used to define the structure of the expression tree, and the continuous variables to evaluate the resulting symbolic expression for specific numerical values associated with the data [START_REF] Horesh | Globally convergent system and method for automated model discovery[END_REF]?]. To choose among the many possible expressions that fit the data (within an error bound), we set the objective to minimize the description complexity of the expression. Consequently, the MINLP formulation has the form

min C(f cqwz ) (Complexity) s.t.
(c, q, w, z) ∈ T (Grammar)

v i = f cqwz (x i ), ∀i ∈ I (Prediction) D(v, y) ≤ (Error)
where c, q, w and z are decision variables, f cqwz is the expression tree defined by these variables, T is the set of values of the decision variables that define valid expression trees of bounded size, C measures the description complexity of the expression tree, and, D measures error of the predicted values v and (x, y) are the dataset observations. In practice, MINLP solvers (e.g., BARON [START_REF] Tawarmalani | A polyhedral branch-and-cut approach to global optimization[END_REF], COUENNE[BLL + 09], SCIP[MFG + 17]) employ various convex relaxation schemes to obtain lower bounds on the objective function value and use these bounds in divide-and-conquer strategies to obtain globally optimal solutions [START_REF] Tawarmalani | A polyhedral branch-and-cut approach to global optimization[END_REF].

We next give the details of our formulation. The input to the formulation consists of a rooted binary tree, a set of candidate operators and an observed dataset. For each observation i ∈ I, we denote the value of the independent variables by x i ∈ R m and the dependent variable values by y i ∈ R. Let D = {1, . . . , m} denote the indices of the independent variables. Let T = (N, E) be the input binary tree, let s ∈ N denote root node and L ⊂ N denote the leaf nodes. Each node n ∈ N \ {s} in the tree has exactly one predecessor and each node n ∈ N \ L has exactly two successors. Let O = U ∪ B denote the set of candidate operators where U contains the unary operators (such as squareroot) and B contains the binary operators (such as addition). Our formulation has two parts:

The first part consists of constraints and decision variables that construct the expression tree and the second part is used to evaluate the difference between the estimated and actual label of each observation using the constructed expression tree.

Grammar. The first part of the formulation chooses a subtree of T and assigns either an operator, a constant, or an input variable to each node of the chosen subtree in a consistent way. For each n ∈ N , we have a decision variable u n ∈ {0, 1} to indicate if the node is used (active) in the expression tree. We also have a decision variable z n,o ∈ {0, 1} to denote if the operator o ∈ O is assigned to the node, and variables w n,d ∈ {0, 1} to denote if the variable d ∈ D is assigned to the node. Finally we have a decision variable q n ∈ {0, 1} if a constant value is assigned the node. The constraint

q n + d∈D w n,d + o∈O z n,o = u n , ∀n ∈ N, (1)
enforces that each active node is assigned an operator, a constant, or one of the variables. Let n ∈ N \ L be a non-leaf node with successors l, r ∈ N where l has a smaller index than r. Note that if node n is assigned a binary operator then both nodes l, r have to be active in the tree; and if n is assigned a unary operator then exactly one of the nodes (say l) has to be active. Consequently, we have constraints u r = o∈B z n,o and u l = o∈B∪U z n,o . Notice that these constraints also enforce that if a node is inactive, or is assigned a constant or is a variable, then its successor nodes cannot be active. As the leaf nodes cannot be assigned operators, we also set z n,o = 0 for all o ∈ O, l ∈ L.

Finally, we define a continuous variable c n for n ∈ N to denote the value of the constant when q n = 1. It is possible to show that any c ∈ R |N | together with binary vectors q, w, z, and u (of appropriate dimensions) that satisfy the constraints above define an expression tree and vice versa.

In addition to the basic model described above, we impose some additional constraints to improve computational performance (without compromising optimality). For example we make sure that if a node is assigned a binary operator, then both its successor nodes cannot be assigned constants at the same time. Similarly the left successor of a unary operator cannot be a constant. We also have several constraints to deal with symmetry as expression trees are inherently symmetric objects in the sense that the same function can be represented with different trees, for example by flipping two branches succeeding a commutative binary operator. To avoid numerical problems, we also set bounds on the absolute value of the continuous variables.

Prediction. Using the first set of variables, the second part of the formulation computes the value of the function (defined by the expression tree) for each observation. More precisely, we define a variable v n,i to denote the value of node n ∈ N for observation i ∈ I. Value of a node clearly depends on the operator assigned to it and the values of its successor nodes. For example, if the addition operator a ∈ B assigned to node n ∈ N that has successor nodes r, l ∈ N , then v n,i would be the sum of v r,i and v l,i . Therefore we have a constraint:

v n,i = q n c n + d∈D w n,d x i d + o∈B z n,o f o (v r,i , v l,i ) + o∈U z n,o f o (v l,i ) ∀i ∈ I, n ∈ N, (2)
where f o (•) denotes the function on the argument(s) of the operator o ∈ O. Also note that v n,i = 0 only if u n = 1. Leaf nodes cannot be assigned operators, therefore the last two summations in (2) are zero for these nodes.

Error. Recall that s ∈ N denotes the root node of the expression tree and therefore v s,i corresponds to the estimated output variable for observation i ∈ I. We define the error for observation i ∈ I as (v s,i -y i ) and relative error as (v s,i /y i -1). Using these expressions, we can define a discrepancy function such as D(v s , y) = i∈I (v s,i -y i ) 2 . In our computational experiments, we control the relative error by a constraint added to the model i∈I (v s,i /y i -1) 2 ≤ for a given > 0.

Complexity. We can define a function to measure description complexity of the model in various ways, for example C(f cqwz ) = n∈N o∈O c o z n,o , where c o ∈ R is a complexity weight assigned to operator o ∈ O by the user. In our computational experiments, we use the expression n∈N u n to minimize the total number of nodes in the expression tree.

Computational Experiments

We report our computational experience on the discovery of physical laws using real-world datasets. We solve the optimization problem described in Section 2 with the MINLP solver BARON [START_REF] Tawarmalani | A polyhedral branch-and-cut approach to global optimization[END_REF] v17.8.9, using IBM ILOG CPLEX and IPOPT [START_REF] Wächter | On the implementation of a primal-dual interior point filter line search algorithm for large-scale nonlinear programming[END_REF] as subsolvers. We can only use operators that are handled by BARON. Furthermore, when multiple globally optimal solution exists, we do not have control over which one is returned. Experiments are run on the cloud, so CPU speed cannot be precisely quantified.

Exoplanet Data. The first set of experiments concerns the discovery of Kepler's law on planetary motion on a dataset taken from NASA [START_REF]NASA exoplanet archive[END_REF] reporting information on planetary systems. The dataset consists of tuples of the form (planet, star, τ, M, m, d) where τ is the orbital period of the planet, M the mass of the star, m the mass of the planet, d the major semi-axis of the orbit. Planet features are normalized to Earth's, star features to Sol's. τ is further normalized since its range is very large. We build three test problems from this dataset, including the systems listed in brackets: EP1 (Sol and Trappist), EP2 (HD 219134 and Kepler), EP3 (GJ 667C, HD 147018, HD 154857, HD 159243 and HD 159868).

Table 1: Results for the exoplanet dataset. Time limit is 6 hours. Grey cells exceed the time limit. c represents constants (we do not report its value for brevity).

Maximum relative error

Dataset 2% 5% 10% 20% 30% 50% EP1 3 √ cτ 2 M 3 √ cτ 2 M 3 √ τ 2 (M + c) 3 √ τ ( √ τ + M ) √ cτ 2 √ τ M EP2 3 √ cτ 2 M c 3 √ τ 2 c 3 √ τ 2 3 √ cτ 2 √ τ √ τ EP3 3 √ cτ 2 M 3 √ cτ 2 M √ τ M + τ 3 √ cτ 2 M √ cτ + c √ τ
Table 2: Results for the pendulum dataset. Time limit is 6 hours. All formulas are certified optimal.

Maximum relative error Dataset 0.5% 1% 2% 5% DS1 j( -c) -

j √ + c j √ cj DS2 j √ j √ cj cj DS3 j √ j √ -cj cj DS4 j( -c) -c j √ cj cj
We consider the set of operators +, * , √ , 3 √ , and full binary expression trees with depth at most 3 (the root has depth 0) in which each operator can appear at most 3 times. We aim to predict d in terms of the other input variables. Results are reported in Table 1, for different upper bounds applied to the relative model error. For small model errors, we always find a refined expression of Kepler's third law: d = 3 √ cτ 2 M . Given the data, the more comprehensive formula would be d = 3 cτ 2 (M + m), but the dependency on m is not picked up because m is negligible compared to M . As model error upper bound increases, we find simpler formulas that do not accurately reflect Kepler's third law.

Most of the formulas returned are certified globally optimal in a relatively short period of time: on average, 50 minutes for the instances solved to global optimality (two instances hit the time limit and thus are not certified globally optimal within the time limit), standard deviation 75. The number of nodes explored by the Branch-and-Bound algorithm varies greatly: from 147 for the simplest formula √ τ , to over 600k for one instance that is not solved to optimality within the allotted runtime. The geometric mean is 19917. By contrast, the number of binary expression trees of depth 3 with 4 candidate operations, 3 input variables, and numerical constants is ≈ 2 • 10 7 without accounting for symmetry. We remark that the number of nodes explored by Branch-and-Bound does not correspond to the number of candidate expression trees that have been examined: it is merely a measure of difficulty of the search.

Pendulum Data. The second set of experiments concerns ten pendulums with different link lengths, and timestamps for the time at which each pendulum crosses the midpoint of its arc. The data is obtained from a high-resolution video recording of the pendulums, yielding hundreds of timestamps. From the raw data, we randomly extract ten tuples (one for each link length) of the form ( , i, t i , j, t j ), where is the link length, i and j are integers, t i and t j are the timestamps for the i-th and j-th crossing of the midpoint, and j > i. We aim to predict t j .

The period of the pendulum is given by τ = 2π /g. Hence, we aim to obtain the equation t j = πj /g. The experimental setup is similar to the exoplanet dataset, but we add "-" to the list of operators. Many experiments return the correct formula except the multiplicative constant, since π/ √ g ≈ 1.002 and its discrepancy from 1 is too small to be significant, up to experimental error.

With error bound 0.5%, we observe overfitting on DS1 and DS4: the simpler formula j √ does not satisfy the error bound due to experimental error. Because the link lengths of the ten pendulum are close to each other (ranging from 0.28m to 0.307m), for larger error bounds a simpler expression of the form cj suffices. The average runtime is 2 minutes, standard deviation 2.3. The number of nodes ranges from 57 to 8134, geometric mean 299. Despite the presence of redundant variables, the algorithm quickly identifies the simplest formula explaining the data.

In this extended abstract we focused on small datasets and did not address scalability, which is a known limitation of Branch-and-Bound-based methods for MINLP. This is left for future research.


Sonia Toubaline
email: sonia.toubaline@dauphine.fr

Claudia D'ambrosio

Leo Liberti

Pierre-Louis L Poirion

Baruch Schieber

Hadas Shachnai

Claudia D'

Complexity and inapproximability results for the Power Edge Set problem

Keywords: PMU placement problem, Power Edge Set, NP-hardness, inapproximability

Introduction

Monitoring an electrical network is an important and challenging task. To ensure ongoing reliability and quality of electricity supply to customers, the state of the electrical network must be monitored continuously. The state of such a network is usually dened as the values of all voltages on its nodes and the branch currents. Phasor measurement units (PMUs) are monitoring devices that can be used for this purpose. PMUs are designed to be placed at (sub)stations and can measure their voltage and the current on all their outgoing transmission lines [START_REF] Manousakis | Taxonomy of PMU placement methodologies[END_REF]. Note that it is not necessary to place PMUs at all stations, as some of the currents and voltages can be deduced using Ohm and Kirchho Laws. Due to their high cost, nding a placement of PMUs that minimizes their number while still ensuring monitoring of the whole network is an important problem.

Let the electrical network be modelled by a graph where the vertices represent electrical nodes and the edges correspond to transmission lines joining two nodes. The (multi channel) PMU placement problem, also called the power dominating set (PDS) problem, consists of nding a minimum number of PMUs to install on the vertices such that all the graph is observed, that is, all voltages and currents are measured. [START_REF] Brueni | The PMU placement problem[END_REF] showed that the observability of a graph by (multi channel) PMUs can be dened by two rules: (R1) if a PMU is placed at a vertex then this vertex and all its neighbours are observed; (R2) if all the neighbours of an observed vertex except one are observed, then this latter is observed.

Several complexity results have been shown for PDS: NP-completeness proofs for bipartite, cographs [START_REF] Haynes | Domination in graphs applied to electric power networks[END_REF], and planar bipartite graphs [START_REF] Brueni | The PMU placement problem[END_REF], and polynomial-time algorithms for trees, grids [START_REF] Doring | A note on power domination in grid graphs[END_REF], block graphs [START_REF] Xu | Power domination in block graphs[END_REF], and bounded treewidth [START_REF] Guo | Improved algorithms and complexity results for power domination in graphs[END_REF]. Approximation algorithms and hardness results are presented in [START_REF] Aazami | Approximation algorithms and hardness for domination with propagation[END_REF]: O( √ n)-approximation for planar graphs and NP-hardness of approximability within a factor 2 log 1-n . Also, various solution methods have been proposed to solve the PMU placement problem [START_REF] Manousakis | Taxonomy of PMU placement methodologies[END_REF]. Some of the PMUs available in the market have a limited number of channels. A PMU with k channels installed at a vertex v can observe only v and k of its neighbors (and the edges connecting v to these k neighbors). The identity of these k neighbors is determined at the time of installation. We consider PMUs with a single channel that when placed at vertex v observe only v and one of its neighbors. Suppose that the single neighbor of v observed by the PMU is u. Since both v and u are observed, we can equivalently view this as placing the PMU on the edge (v, u), where a PMU placed on an edge is assumed to observe both its endpoints.

The observability of a graph using single channel PMUs can also be determined by two rules similar to the multi channel case. The second rule (R2) is the same as above while (R1) needs to be modied to account for the single channel to: (R1 E ) if a PMU is installed on an edge then its two endpoint vertices are observed. The objective is to nd a minimum placement of single channel PMUs that ensures the observability of the whole graph. We call this problem the Power Edge Set (PES) problem.

Emami et al. proposed a binary linear program for this PMU placement problem [START_REF] Emami | Robust measurement design by placing synchronized phasor measurements on network branches[END_REF], which considers only rule R1 E , and hence turns out to be equivalent to the minimum edge cover problem, which has polynomial-time solution. The authors discussed the consideration of a restricted version of R2 in [START_REF] Emami | Optimal placement of phasor measurements for enhanced state estimation: A case study[END_REF]. In [START_REF] Poirion | The power edge set problem[END_REF], we studied this problem from a practical point of view. We rst proposed a naturally iterative index binary linear model that turns out to be too large for practical purposes. Using a xed point argument, we removed the iteration indices and obtained a bilevel formulation. We then reformulated the latter to a single-level mixed-integer linear program, which performs better than the natural formulation. We then provided a cutting plane algorithm that solves the bilevel program much faster than an o-the-shelf solver can solve the previous models. In this paper, we study the complexity of this problem when both rules of observability are considered. We show that the PES problem is NP-hard to approximate within a factor (1.12)-, for any > 0. We also prove that it is polynomial-time solvable for trees and grids.

Problem statements and preliminaries

Let G = (V, E) be an undirected graph with |V | = n. Denote by N (v) the set of neighbors of a vertex v ∈ V . G is k-regular if all vertices have a degree k.

Given a subset S ⊆ V , dene the set B(S) ⊆ V of "observed" vertices starting at S. Initially, B(S) is set to be S. Then, as long as there exists a vertex v / ∈ B(S) such that v has a neighbor u ∈ B(S) and N (u) \ {v} ⊆ B(S), i.e., v is the only neighbor of u that is not in B(S), then v is added to B(S) (rule R2). In this case, we say that u is the parent of v and v is the child of u. Note that each parent has at most one child and each child has exactly one parent. We extend the parent/child relations to ancestor/descendant relations in a natural way.

In the Power Edge Set (PES) problem, we are given a subset of edges F ⊆ E and the initial set of observed vertices S(F ) is the set of all endpoints of the edges in F (rule R1 E ). B(S(F )) is constructed as describe above (rule R2). We say that F is a Power Edge Set (PES) if B(S(F )) = V . The Power Edge Set problem is to select a PES F of minimum cardinality.

Power Edge Set Problem

Input: A graph G = (V, E). Output: A minimum cardinality set F ⊂ E such that B(S(F )) = V .

To establish the NP-hardness approximation result, we use the notion of an E-reduction described in the following. For a function f , an algorithm is an f (n)-approximation, if for every instance I of the problem, it returns a solution S such that r(I, S) ≤ f (|I|).

E-reduction

Khanna et al. introduced the notion of an E-reduction (error-preserving reduction) in [START_REF] Khanna | On syntactic versus computational views of approximability[END_REF]. A problem Π is called E-reducible to a problem Π , if there exist polynomial time computable functions f , g and a constant β such that f maps an instance I of Π to an instance I of Π such that opt(I) and opt(I ) are related by a polynomial factor, i.e. there exists a polynomial p such that opt(I ) ≤ p(|I|)opt(I), g maps any solution S of I to one solution S of I such that ε(I, S) ≤ βε(I , S ).

An important property of an E-reduction is that it can be applied uniformly to all levels of approximability; that is, if Π is E-reducible to Π and Π belongs to C then Π belongs to C as well, where C is a class of optimization problems with any kind of approximation guarantee (see also [START_REF] Khanna | On syntactic versus computational views of approximability[END_REF]).

NP-Hardness of Approximation for PES

We prove that PES is hard to approximate within some constant, unless P=NP. To this end, we dene an E-reduction from Min Vertex Cover restricted to 3-regular graphs. This problem is NP-hard to approximate within a factor of 1.36 [START_REF] Dinur | On the hardness of approximating minimum vertex cover[END_REF][START_REF] Feige | Vertex cover is hardest to approximate on regular graphs[END_REF]. We rst describe our reduction.

Let I be an instance of Min Vertex Cover formed by a 3-regular graph G = (V, E). We construct an instance I of PES consisting of a graph G = (V , E ) as follows (see Fig. 2). We associate to each vertex v ∈ V a gadget G v in G (see Fig. 1). For v ∈ V , G v is composed of 10 vertices {v 0 , . . . , v 9 }. For each edge (v, u) ∈ E, we add the edge (v i , u j ) to E for one i ∈ {1, 2, 3}, and u j ∈ G u for one j ∈ {1, 2, 3}. The vertex subsets {v 1 , . . . , v 5 } and {v 6 , . . . , v 9 } form two cliques linked to each other by the edges {v i , v i+5 }, i = 1, . . . , 4. We also add the edge (v 0 , v 5 ) to E . Let u, w, t be the three neighbors of v ∈ V . Without lose of generality, we assume that (v 1 , u 1 ), (v 2 , w 2 ), and

(v 3 , t 3 ) are in E , for v i ∈ G v , i = 1, 2, 3, u 1 ∈ G u , w 2 ∈ G w and t 3 ∈ G t .
The vertices u 1 , w 2 , and t 3 are called "neighbor vertices" of G v and v 1 , v 2 , v 3 its junction vertices.

v 1 v 2 v 3 v 4 v 6 v 7 v 8 v 9 u 1 w 2 t 3 v 5 v 0 Fig. 1 Gadget Gv associated to a vertex v ∈ V . u 1 ,
v in G. v u t w ; v 1 v 2 v 3 u 1 u 3 u 2 w 2 w 3 w 1 t 3 t 2 t 1 G v G u G w G t
Fig. 2 Gadget associated to a vertex, u ∈ V with three neighbors v, w and t

In the following, we present four results on the number of PMUs to install according to the degree of vertices and the observability status of the junction and/or neighbor vertices.

Lemma 1 Let G = (V, E) be a connected graph with at least four vertices, such that only one PMU is placed at a given edge {u, v} of the graph. If no vertex of G has degree 2, then it is impossible that all the vertices of G are observed. Proof By the PMU placed at {u, v}, the vertices u and v are observed using R1 E . Since no vertex in G has degree 2, d(u) or d(v) is at least 3 (either both have degree at least 3 or one has degree 1 and the other degree at least 3). Let us assume that d(u) ≥ 3. The vertex u has then at least 2 adjacent vertices that are not observed. The same holds if d(v) ≥ 3. Then R2 cannot be used at either u or v. Then G is not observed.

Denition 1 Given a gadget G v in G , a junction vertex v i of G v , i ∈ {1, 2, 3}, is "observed externally" if v i can be observed by applying R2 on the neighbor vertex of G v adjacent to v i .
Lemma 2 Given a gadget G v of G , the placement of two PMUs is sucient to observe all the vertices of G v and its neighbor vertices. Furthermore, unless the neighbor vertices are observed externally, two PMUs are also necessary to observe all the vertices of the gadget. Proof We place one PMU on {v 6 , v 7 } and one on {v 8 , v 9 }. The vertices v 6 , . . . , v 9 are observed by R1 E . Using R2, the vertices v i-5 , for i = 6, . . . , 9, are then observed. By applying R2 to v 4 (which is now observed), the vertices v 5 and then v 0 are observed. Finally applying R2 to the junction vertices allow us to observe the neighbor vertices.

Now assume that all neighbor vertices of G v are not observed and that only one PMU installed in G v is sucient to observe all the vertices of G v . By Lemma 1, G v cannot be observed by only one PMU since no vertex in G v has a degree 2, contradiction. Then, at least two PMUs need to be installed.

Lemma 3 Given a gadget G v of G where v 1 , v 2 and v 3 are observed externally, one PMU is necessary and sucient to observe all the remaining vertices of G v . Proof For each i ∈ {1, 2, 3}, v i has three non-observed adjacent vertices, with v 4 and v 5 in common. It is then impossible to observe all the vertices of G v with no PMU. Let us place one PMU on {v 4 , v 5 }. Using R2 on v 5 , v 0 is observed. Furthermore, since all the vertices of the clique {v 1 , ..., v 5 } are now observed, for i ∈ {6, . . . , 9}, the vertices v i are observed by applying R2 on v i-5 .

Lemma 4 Given a gadget G v of G , such that at most two vertices of G v are observed externally, then the remaining vertices of G v cannot be observed if only one PMU is placed in G v . Proof Assume that v 1 and v 2 are observed externally (by symmetry we can prove the two other cases) and that only one PMU is installed in G v . The vertices v 1 and v 2 have each four non observed adjacent vertices, with vertices v 3 , v 4 and v 5 in common. Placing the PMU on one of the edges linking v 1 and v 2 to these non-observed vertices will not allow us to apply R2. Similarly, v 3 and v 4 having three non-observed adjacent vertices, v 5 having three nonobserved adjacent vertices, and v 6 , . . . , v 9 having at least three non-observed adjacent vertices, R2 cannot be used. Therefore, one PMU is not sucient to observe G v .

We now present our inapproximability result.

Theorem 1 PES is NP-hard to approximate within a factor (1.12)-, for any > 0.

Proof We prove rst that opt(I ) and opt(I) are polynomially related. Consider an optimal solution C * of I. Consider now a solution F ⊂ E of I . In F , PMUs can be installed on edges between gadgets and some gadgets can have more than two PMUs placed on them. We show that there exists another solution F of I that is at least as good as F and containing no PMUs on the edges linking gadgets and only one or two PMUs by gadget.

Let Π = {{v 6 , v 7 }, {v 8 , v 9 } : v ∈ C * } ∪ {{v 4 , v 5 } : v ∈ C * }
1. If a PMU is placed on {u 1 , v 1 } the edge linking G u and G v (the other cases are symmetric). If v i , for i ∈ {2, . . . , 6}, is the child of v 1 then F = F \ {{u 1 , v 1 }} ∪ {{v 1 , v i }}. Since all the vertices of G v are observed (F is a PES), u 1 can be observed using R2 on v 1 . If v 1 has no child then F = F \ {{u 1 , v 1 }}. Since all the other vertices of G v are observed (F is a PES), v 1 is observed by one of v i , for i ∈ {2, . . . , 6}. u 1 can be observed using R2 on v 1 . 2. There are at least 3 PMUs installed in a given gadget G v . Let F v be subset of edges of G v where a PMU is installed. The construction of F depends on the observability of neighbor vertices u 1 , w 2 and t 3 of G v . We distinguish 3 cases: a. All the neighbor vertices are observed: by Lemma 3, only one PMU is necessary and sucient to observe G v . Then considering the placement proposed in Lemma 3, Consider C = {u : G u has two PMUs installed on it} a subset of vertices in G. We prove that C is a cover by contradiction. We assume that there exists an edge {u, v} that is not covered by C, i.e. u, v ∈ C. Thus, G u and G v has only one PMU installed on them on G . By Lemmas 3 and 4, u 1 , u 2 , and u 3 are externally observed and so are v 1 , v 2 ,and v 3 . If we assume that {u 1 , v 1 } is the edge linking G u and G v , then u 1 is observed externally by v 1 using R2 and v 1 is observed externally by u 1 using R2, which is impossible, contradicting the assumption that F is a PES of I . Then C is a cover and val(I, C) = val(I , F ) -n. In particular, when F is an optimum solution, we have opt(I ) = val(I, C) + n ≥ opt(I) + n. It follows from the previous result that opt(I ) = opt(I) + n. Since I is an instance of a Minimum Vertex Cover dened on a 3-regular graphs we have that n 2 ≤ opt(I) ≤ 3n 4 [START_REF] Feige | Vertex cover is hardest to approximate on regular graphs[END_REF]. Then 3n 2 ≤ opt(I ) ≤ 7n 4 . We obtain that n ≤ 2opt(I ) 3 , and we get opt(I ) opt(I )-n ≤ 3. Therefore ε(I, C) ≤ 3 val(I ,F )-opt(I )

F = F \ F v ∪ {{v 4 , v 5 }}. b.
opt(I ) = 3 ε(I , F ).
Thus, r(I, C) -1 ≤ 3(r(I , F ) -1) and then r(I , F ) ≥ r(I,C)+2

3

. Since r(I, C) = ρ = 1.36, we have r(I , F ) ≥ ρ+2 3 = 1.12.

Polynomial-time cases for PES

We prove in this section that the PES problem has polynomial-time solutions for trees and grids.

Trees

We prove that, in trees, PES is equivalent to the Path Cover Problem, dened as follows.

Path Cover Problem

Input: A graph G = (V, E). Output: A minimum cardinality set of vertex disjoint paths, such that each vertex belongs to a path. (A singleton vertex is also considered a path.)

The Path Cover Problem is NP-hard for general graphs (as the Hamiltonian Path problem is easily reduced to it), but polynomial-time solvable on trees [START_REF] Moran | Optimal covering of cacti by vertex-disjoint paths[END_REF].

Theorem 2 The Power Edge Set problem is polynomial-time solvable on trees. On trees with n vertices, the algorithm runs in O(n) time.

Proof We prove that PES is equivalent to the Path cover problem.

Consider any solution of the PES problem of size k. We determine the parentchild paths starting from all the endpoints of these k edges. This results in a path cover of size at most 2k. (It will result in less than 2k paths if some edges in the PES share endpoints.) We now show how to reduce the size of the path cover to k. Let e 1 , . . . , e k be the edges in the PES. For an edge e i , let X i and Y i be the paths starting from the endpoints of e i . We construct the path cover P 1 , . . . , P k iteratively. In the i-th iteration, we consider X i and Y i : if both are not in the current partial cover, we add the path P i = (X i , e i , Y i ) to the partial cover. If only one of X i or Y i , say X i , is not in the partial cover (implying that e i shares an endpoint with one of the edges e 1 , . . . , e i-1 ), add X i to the partial cover. Note that due to the minimality of the PES, we cannot have a situation where both X i and Y i are not in the partial cover. It is easy to see that, after k iterations, we end up with a path cover of size k. Thus, any optimal solution for PES induces a solution of the same cardinality for the Path Cover Problem.

Consider now a path cover P 1 , . . . , P k . We construct a PES by putting a PMU on the extremity edge of each path. If for some i = 1, . . . , k, P i is a singleton vertex, we put a PMU on one of its incident edge. We claim that R2 can be applied to observe the rest of the vertices along the paths. We show this by contradiction: suppose that after the vertices in some prexes of P 1 , . . . , P k are observed, we reached a point where R2 cannot be applied anymore. Let i 1 , . . . , i be the indices of the paths that still have unobserved vertices, and let x ij be the last observed vertex in each such path. Consider the vertex y i1 that follows x i1 in path P i1 . Since x i1 cannot observe y i1 , it must have another unobserved neighbor. This neighbor cannot be on P i1 , since then we have a cycle. Assume that this neighbor is on P i2 and consider x i2 . Vertex x i2 also has an unobserved neighbor in addition to y i2 , but this neighbor cannot be on P i2 and P i1 , because in both cases we close a cycle. So assume it is on P i3 . We can repeat the process at most k times until we must have a cycle, a contradiction. Then any solution of Path Cover Problem induces a solution of the same cardinality for the PES problem.

Therefore, since the Path cover problem can be solved in O(n) time for trees [START_REF] Moran | Optimal covering of cacti by vertex-disjoint paths[END_REF], we deduce that the PES problem is also solvable in O(n) time for trees.

Grids

Let G m×n = (V, E) be a grid graph that is the graph Cartesian product of P m × P n of path graphs on m and n nodes respectively. We prove in the following that the PES problem is polynomial-time solvable for grids with opt(G m×n ) = 1 2 min{m, n} . Let = min{m, n}.

Lemma 5 For any optimal solution S ⊆ E of the PES problem and for each P i , for i = 1, . . . , , at least one node of P i is observed by an edge in S using rule R1 E . Proof Assume that = m. For contradiction, assume that Lemma 5 does not hold, i.e., there exists a solution S such that there exists a row path P i , for i = 1, . . . , m, with none of its nodes observed using R1 E . Consider P k such a path, k ∈ {1, . . . , m}. Thus its n nodes are observed from n distinct nodes of P k-1 or P k+1 using rule R2. By the same reasoning, these n nodes of P k-1 or P k+1 that are not observed by R1 E are observed from the n distinct nodes of P k-2 or P k+2 , and so on until we get to the border of the grid P 1 or P m .

We deduce that at least n nodes of G m×n are observed by rule R1 E using at least n 2 edges in S. However, G m×n can be observed by exactly m 2 PMUs. Let v 1 , v 2 , ..., v m/2 be the nodes of P 1 . By placing the PMUs on the edges {v 2i , v 2i-1 }, for i = 1, . . . , m/2 and one more PMU on {v m-1 , v m } if m is odd, we can observe the nodes of P 1 using R1 E and then all the nodes of P k from the nodes of P k-1 using R2, for k = 2, . . . , n. Hence, S is not optimal, contradiction.

The case = n is proved in a similar way.

Theorem 3 The Power Edge Set problem is polynomial-time solvable on grids. On grids with size m × n, an optimal solution is obtained in O( ) time and its optimal value in O(1) time, with = min{m, n}. Proof Let S be an optimal solution for the PES problem. By Lemma 5, since at least one node of each P i , for i = 1, . . . , , is observed by an edge in S using rule R1 E , then |S| ≥ 2 . We prove that there exists an optimal solution S * with size exactly 2 . Assume that = m. As proved in the proof of Lemma 5, if we note v 1 , v 2 , . . . , v m/2 the nodes of the column path P 1 , then S * consists of the set of edges {v 2i , v 2i-1 }, for i = 1, . . . , m/2 and the edge {v m-1 , v m } if m is odd. The nodes of P 1 are observed using R1 E and then all the nodes of P k are observed from the nodes of P k-1 using R2, for k = 2, . . . , n. We have then B(∪ m i=1 {v i }) = V . Now, if = n, we consider the row path P 1 and denote its nodes by v 1 , v 2 , . . . , v n/2 . As for the previous case we prove that for S * consisting of the set of edges {v 2i , v 2i-1 }, for i = 1, . . . , n/2 and the edge {v n-1 , v n } if n is odd, we have B(∪ n i=1 {v i }) = V . Therefore, S * is an optimal solution for the PES problem and is obtained in O( ) time. Its optimal value opt(G m×n ) = 1 2 min{m, n} and is obtained in O(1) time.

Remark 1 Dening the PES problem using R1 E and R2 gives rise to a natural variant of this problem that corresponds to PMUs with "zero" channels. In this variant PMUs are placed on a vertex and can observe only the vertex they are installed on. The observability of a graph using "zero" channel PMUs can also be determined by two rules. The second rule (R2) is the same as above while (R1) needs to be modied to: (R1 V ) if a PMU is installed on a vertex then only this vertex is observed. Again, the objective is a minimal placement of PMUs that ensures the observability of the whole graph. We call this problem the Power Vertex Set (PVS) problem. It is also known as the Zero Forcing Set problem (Minimum Rank-Special Graphs Work Group, 2008). The PVS problem is NP-hard for general graphs [START_REF] Aazami | Hardness results and approximation algorithms for some problems on graphs[END_REF] and polynomial-time solvable for trees (Minimum Rank-Special Graphs Work Group, 2008).

It is easy to observe that the size of the optimal PVS is lower bounded by the size of the optimal PES which in turn is lower bounded by the size of the optimal PDS. While the ratio of the optimal PVS size to the optimal PES size is at most 2, the ratio of the optimal PES size to the optimal PDS size may be as large as n -2. To see this consider a "star" graph with n vertices and n -1 edges. It is easy to see that the PDS size of this graph is 1 (placing the PMU in the center), while the size of the PVS and the PES is n -2.

Therefore, we can deduce the following bounds for the optimum value of an instance of the PES problem: 5 Conclusions

  Consider an NP Optimization problem and an instance I of this problem. We denote |I| the size of I, opt(I) the optimum value of I, and val(I, S) the value of a feasible solution S of I. The approximation factor r(I, S) of S is given by max val(I,S) opt(I) , opt(I) val(I,S) . The error of S, noted ε(I, S), is dened by ε(I, S) = r(I, S) -1.

  be the placement consisting of installing two PMUs in the gadgets associated to vertices in C * and one PMU in the other gadgets. According to Lemma 2, all the vertices of G v , v ∈ C * and their neighbor vertices are observed. For v ∈ C * , its neighbor vertices u, w, t ∈ C * are observed since otherwise the edges (v, u), (v, w) and (v, t) are not covered. Then two PMUs are installed in G u , G w and G t , and by Lemma 2, all their neighbor vertices that include v 1 , v 2 and v 3 , are observed. Hence, by Lemma 3, the remaining vertices of G v can be observed by placing only one PMU in G v . Therefore, opt(I ) ≤ |Π| = opt(I) + n. Since |C * | ≥ 2 27 (n-1), we have that n ≤ 27 2 opt(I)+1. Then opt(I ) ≤ 31 2 opt(I). Therefore, for n large enough, opt(I ) ≤ n opt(I).

  Therefore, we have opt(I ) ≤ n opt(I) and ε(I, C) = val(I,C) opt(I) -1 = val(I ,F )-n opt(I )-n -1 = val(I ,F )-opt(I ) opt(I )-n = val(I ,F )-opt(I ) opt(I ) × opt(I ) opt(I )-n .

  I D )} ≤ opt(I E ) ≤ min{opt(I V ), (n -2)opt(I D )},where I V , I E and I D are instances of the PVS, PES and PDS problems.

  w 2 and t 3 are vertices of Gu, Gw and Gt respectively where u, w and t are the adjacent vertices of

  All the neighbor vertices are not observed: By Lemma 2, only two PMUs are necessary and sucient to observe G v . Then considering the placement proposed in Lemma 2,F = F \ F v ∪ {{v 6 , v 7 }, {v 8 , v 9 }}.c. At most 2 are observed: By Lemma 4, one PMU is not enough to observe all the vertices of G v . By Lemma 2, two PMUs are sucient to observe all the vertices of G v . Then considering the placement proposed in Lemma 2, F = F \ F v ∪ {{v 6 , v 7 }, {v 8 , v 9 }}. Therefore, F is such that |F | ≤ |F | ≤ k, all the PMUs are only placed on edges of gadgets and each gadget has either one or two PMUs installed on it.

We presented NP-hardness of approximability for the PES problem in general graphs and polynomality results in trees and grids. An interesting avenue for future work would be to study the complexity of this problem on other classes of graphs, such as graphs with bounded treewidth, cographs, and regular and bipartite graphs. We conjecture that PVS and PES are easier than the PDS problem. This can be substantiated by nding a polynomial time algorithm for a class of graphs on which PDS is NP-hard. Another direction for further study is nding approximation algorithms for PES. Also, further work would be to nd interesting upper and lower bounds for the optimum value of an instance of the PES problem.

This work was carried out as part of the SOGRID project (www.so-grid.com), co-funded by the French agency for Environment and Energy Management (ADEME) and developed in collaboration between participating academic and industrial partners.


Leo Liberti
email: liberti@lix.polytechnique.fr

Mathematical programming bounds for kissing numbers

We give a short review of existing mathematical programming based bounds for kissing numbers. The kissing number in K dimensions is the maximum number of unit balls arranged around a central unit ball in such a way that the intersection of the interiors of any pair of balls in the configuration is empty. It is a cornerstone of the theory of spherical codes, a good way to find n equally spaced points on the surface of a hypersphere, and the object of a diatribe between Isaac Newton and David Gregory.

1 A brit and a scot went down the pub. . . "Kissing" is billiard jargon. British players would say two adjacent billiard balls on the table "kiss". The term found its way into mathematics thanks to Isaac Newton (whom everyone knows) and David Gregory (a professor of Mathematics at Edinburgh -without having ever obtained a degree -and then Savilian Professor of Astronomy at Oxford thanks to Newton's influence). In the 1690s, scared of the social unrest in Scotland, Gregory left and visited Newton in Cambridge. According to rumours and well-established British protocol, the brit and the scot went down the pub for a few pints of ale and a game of pool. There, among kissing balls and fumes of alcohol, they got into a brawl about the number of balls that could kiss a central ball on the billiard table. Still sober enough, they counted them, and came to agree on the number six. As the number of pints increased, the two started blabbering about gravity-defying floating balls passionately kissing in three dimensions, and disagreed: Newton, embracing the voice of Kepler, said no more than twelve balls could be arranged around a central one. Gregory, who had to gain his master's approval by attempting to be brilliant and surpris-

ing, said that perhaps thirteen could fit? George Szpiro [START_REF] Szpiro | Newton and the kissing problem[END_REF] recounts a different story in plus.maths.org/content/newton-and-kissing-problem (some nonsense about astronomy and planets), but since neither he nor I were present at the quabble, his word on the matter is just as good as mine.

Applications

Quite aside from the satisfaction I get out of spreading academic gossip on Isaac Newton, it turns out that arranging balls in a kissing configuration has applications.

If we only consider the points of contact of n surrounding balls of radius r with the central unit ball, we obtain a set of unit vectors x 1 , . . . , x n in R K that have pairwise distances at least 2r 1+r . This can be seen in the figure on the left (a two-dimensional section of part of some K-dimensional balls configuration), together with the proportion 2r : (1 + r) = d i j : 1, where d i j = x ix j 2 . This is useful if you ever want to send one of the vectors x i (for i ≤ n) over a noisy communication channel. You might receive a vector y different from all x i 's, but assuming the channel is not too noisy, you can simply assume that y is a corruption of the closest x i . This type of error correcting code is called a spherical code, and denoted by A(n, K, r). There is interest in maximizing r, since this corresponds to larger balls and consequently more errors being corrected by the code. Kissing number configurations correspond to spherical codes A(n, K, 1) where n is maximum for a given K. Finding a spherical code given n, K, r is the SPHERICAL CODE PROBLEM (SCP).

The other (less cited) application is finding n equally spaced points on the Kdimensional sphere S K-1 . On websites such as stackoverflow.com or math. stackexchange.com, it seems people expect this to be an easy problem (possibly because the circle is a simple case: place x i on the circle at an angle 2iπ/n). Some 3D solutions advise scattering equally spaced points on a spiral going from one pole to the opposite, warning readers that it does not guarantee equal spacing. The problem can be formulated for any K by means of spherical codes where the smallest r is maximum. Finding n equally spaced unit vectors on S K-1 is the EQUALLY SPACED SPHERICAL POINTS PROBLEM (ESSPP).

Contents

Most of the results in this paper are known. I propose a new SDP-based heuristic for finding lower bounds, and I give a practitioner's view of Delsarte's Linear Programming (LP) upper bound [START_REF] Delsarte | Spherical codes and designs[END_REF], which is useful for conducting experiments in view of trying to improve current bounds ( [START_REF] Antreicher | The thirteen spheres: a new proof[END_REF][START_REF] Mittelmann | High-accuracy semidefinite programming bounds for kissing numbers[END_REF] also discuss practical issues of computing these bounds). All the experiments reported in this paper are preliminary. Lastly, I enjoyed writing this paper more informally than is usual -I hope readers won't object!

The Kissing Number Problem

Finding kissing numbers and kissing configurations is known as the KISSING NUM-BER PROBLEM (KNP). Formally, this consists in finding the maximum number n of vectors x 1 , . . . , x n ∈ R K such that x i 2 2 = 1 for all i ≤ n and x i • x j ≤ 1 2 for all i < j ≤ n. If n is the kissing number in K dimensions, we write kn(K) = n. We know that kn(2) = 6, kn(3) = 12 (so Newton was right), kn(4) = 24, and we do not know kn(5) but it is at least 40. We also know kn(8) = 240 and kn(24) = 196560 [4, p. 510], and have bounds in many other dimensions (see https://en.wikipedia.org/wiki/Kissing_number_problem).

Computational complexity

The computational complexity of solving KNP, as an optimization or even a decision problem, is a prominent and embarassing question mark. Since the input is a pair of integers n, K, mapping an NP-complete problem (e.g. SAT) to a KNP, such that an instance of one problem is yes if and only if the corresponding instance of the other problem is also yes, really seems quite hard.

On the other hand, the KNP is certainly very hard to solve empirically, and noone working on this problem ever suggested that there might be a polynomial-time algorithm for solving it -even on a real RAM computational model.

Reductions between decision versions of KNP, SCP and ESSPP are as follows: the KNP is included in the SCP by definition (so it trivially reduces to the SCP), and, as pointed out in Sect. 3.2, the KNP can be decided by the decision version of the ESSPP (so, again, it reduces to the ESSPP). The decision versions of the SCP and ESSPP are really the same (though the optimization versions differ). The only reduction I cannot immediately prove is from SCP/ESSCP to the KNP, since only the latter fixes the angular separation at π/3.

I am not aware of any method for proving NP-hardness of problems whose input consists of a constant number of integers. For example, the complexity status of the well-known PACKING EQUAL CIRCLES IN A SQUARE (PECS) problem is as yet undetermined [START_REF] Costa | On the impact of symmetry-breaking constraints on spatial branch-and-bound for circle packing in a square[END_REF]. This is certainly an interesting open problem in computational complexity.

A referee pointed out an interesting link with another problem having undetermined complexity status: the PALLET LOADING PROBLEM (PLP), which asks whether n identical a × b rectangles can be packed in a given X × Y rectangle [12, §2.C]. Thus, PLP instances, like KNP ones, are described by a constant number of integers. Again, establishing reductions between KNP, SCP, ESSPP, PECS and PLP is an open question.

Lower bounds

Since the KNP is a maximization problem, the cardinality of any kissing configuration of balls yields a lower bound. Since any heuristic method might be able to find a good configuration, lower bounds for the KNP are considered "easy" to obtain.

The formulation of Maculan, Michelon and Smith

We start with the MMS95 formulation proposed in [START_REF] Maculan | Bounds on the kissing numbers in R n : Mathematical programming formulations[END_REF], a Mixed-Integer Nonlinear Program (MINLP) which correctly formulates the KNP:

max x∈[-1,1] τK ,α∈{0,1} τ τ ∑ i=1 α i ∀i ≤ τ x i 2 2 = α i ∀i < j ≤ τ x i -x j 2 2 ≥ α i α j ,        (1)
where τ is some (estimated) upper bound to the kissing number. The α i (binary) variables choose whether vector x i is part of the configuration or not. Note that the angular separation constraint x i • x j ≤ 1 2 is replaced by a Euclidean distance separation

x i -x j 2 2 ≥ 1, which is equivalent since 1 ≤ x i -x j 2 2 = x i 2 2 + x j 2 2 -2x i • x j = 2-2x i •x j (whence x i •x j ≤ 1/2)
as the vectors all have unit norm. Since the problem is formulated exactly, an exact MINLP solver would provide a feasible and optimal solution if it exists, given some guessed upper bound.

Unfortunately, the state of the art in MINLP solver technology cannot even provide an answer in K = 2 if τ = 7 (one more than kn(2) = 6) in "reasonable time" of a "reasonable laptop" using Eq. (1).

A feasibility formulation

Given K and n, the formulation below finds the configuration of n unit vectors where the minimal separation between closest vectors is maximum [START_REF] Kucherenko | New formulations for the kissing number problem[END_REF]:

max x∈[-1,1] nK ,α≥0 α ∀i ≤ n x i 2 2 = 1 ∀i < j ≤ n x i -x j 2 2 ≥ α.      (2)
Eq. ( 2) has a single scalar α variable, which is continuous and represents the minimum distance between pairs of vectors. Solving this nonconvex Nonlinear Program (NLP) to global optimality and obtaining an optimal α ≥ 1 yields a proof that kn(K) ≥ n; if α < 1 then kn(K) < n. The issue with this strategy for proving kissing numbers is the same as for Eq. ( 1): current solvers just cannot solve these instances to global optimality for interesting values of n, K.

On the other hand, Eq. ( 2) can be used heuristically to find KNP configurations given n, K. Moreover, these feasible solutions will in general spread points over the K-sphere quite evenly, each point being at roughly the same distance from its closest points. They will therefore provide a practical solution to the ESSPP.

Semidefinite Programming relaxations

Semidefinite Programming (SDP) relaxations of Eq. ( 1)-( 2) have not been looked at yet, as far as I know. I performed a few preliminary tests on both, and while the SDP from Eq. (1) seems extremely slack (trivially yielding the upper bound τ for whatever given τ, probably due to relaxed integrality), I found the SDP relaxation of Eq. ( 2) more interesting:

max X∈[-1,1] n 2 ,α≥0 α ∀i ≤ n X ii = 1 ∀i < j ≤ n X ii + X j j -2X i j ≥ α X 0.          (3)
Based on past experience with other problems involving Euclidean distance constraints, I tried the following heuristic strategy:

1. solve Eq. ( 3) and obtain an optimal solution ( X, ᾱ); 2. perform Principal Component Analysis (PCA) using the K largest eigenvalues of X to obtain an n × K matrix x such that the i-th row xi is a vector in R K ; 3. use x as a starting point for a local NLP solver deployed on Eq. ( 2), obtain a "good" feasible solution (x * , α * ).

Using Mosek [START_REF] Aps | The mosek manual[END_REF], a Python implementation of PCA, and IPOPT [START_REF] Coin-Or | Introduction to IPOPT: A tutorial for downloading, installing, and using IPOPT[END_REF], I was able to derive the following KNP configurations on my "resonable laptop" based on a 3.1 GHz Intel Core i7 with 16GB RAM.

(n, K) (6, 2) (12, 3) (24, 4) (40, 5) (72, 6) ε 0 0 0.04 0.05 0.07 CPU (s) 0.02 0.02 0.32 1.57 12.26 The first row measures the (additive) solution error with respect to a valid KNP configuration: α * = 1ε for all α < 1 (whereas α * > 1 whenever ε = 0). Unfortunately, the next interesting case, n = 12 and K = 7, made Mosek crash for lack of RAM (the computational bottleneck on solving SDPs is well known).

Upper bounds

In general, upper bounds to the KNP are considered hard to obtain. Each new upper bound either requires a completely new theoretical point of view, or a substantial body of theoretical work and some computation.

Direct LP and SDP bounds

A couple of easy upper bounding techniques, however, are readily available in the optimization literature: LP and SDP relaxation. Expanding the Euclidean distances in Eq. ( 1) and Eq. ( 2) yields a MINLP and, respectively, a nonconvex NLP involving products of decision variables as the only type of nonlinearity. Using McCormick's [START_REF] Mccormick | Computability of global solutions to factorable nonconvex programs: Part I -Convex underestimating problems[END_REF] for bilinear products and the secant relaxation for square terms one can obtain an LP. For having tested it in the past while I worked on [START_REF] Kucherenko | New formulations for the kissing number problem[END_REF], I know that this kind of LP bound is as slack as they come.

Based on some preliminary experiments using the Mosek SDP solver [START_REF] Aps | The mosek manual[END_REF] on Eq. ( 3) for K ∈ {2, 3}, I observed a regular decrease in the optimal value ᾱ of the objective function of Eq. ( 3) as n increases for a given fixed K, as shown in the figures below.

The decreasing sequence of optimal ᾱ values appears to be converging to 2 for both K = 2 and K = 3. I have not even started considering how to prove it (or disprove it), nor whether it would have any interesting consequence.

Delsarte's LP bound

This is an adaptation to the spherical context of Delsarte's upper bound technique for binary codes [START_REF] Ph | Bounds for unrestricted codes by linear programming[END_REF], based on LP.

Broadly speaking, Delsarte's idea is based on deciding the distance distribution of the code so as to maximize its cardinality. Let a t be the fraction of the vectors in a KNP configuration code A(n, K, 1) that have scalar product equal t. Since these codes contain n vectors, there are at most n 2 values of t (there could be fewer if many scalar products have the same value). In any case, summing over the a t we obtain n 2 /n = n, which is the cardinality of the code. Also, since these are fractions, a t ≥ 0. Moreover, there are exactly n scalar products having value 1, namely

x i • x i = x 2 2 = 1 for all i ≤ n, thus a 1 = n/n = 1. Hence, the LP max{ ∑ t a t | a 1 = 1 ∧ a ≥ 0} ( †)
is of interest to us, insofar as it maximizes the cardinality of the code it describes. The issue at this point is that this LP is unbounded -there is nothing that links the decision variables a t to the "code structure". Delsarte's idea consists in finding a family

F = {φ 1 , φ 2 , . . .} of functions φ : [-1, 1] → R such that ∀φ ∈ F ∑ t a t φ (t) ≥ 0 ( ‡)
is a valid constraint for the LP ( †). Delsarte's "main theorem" has been proved many times, and generalized in various ways [4, §2.2]. Its statement is also valid for the SPC and the ESSPP, as it considers an arbitrary separation angle ζ with cos ζ = z.

Theorem 1. Let c 0 > 0 and f : [-1, 1] → R such that: (i) ∑ i, j≤n f (x i • x j ) ≥ 0 (ii) ∀t ∈ [-1, z] f (t) + c 0 ≤ 0 (iii) f (1) + c 0 ≤ 1.
Then n ≤ 1 c 0 . My favorite proof is the one-liner given in [START_REF] Pfender | Improved delsarte bounds for spherical codes in small dimensions[END_REF]. Let g(t) = f (t) + c 0 , then:

n 2 c 0 ≤ n 2 c 0 + ∑ i, j≤n f (x i • x j ) = ∑ i, j≤n g(x i • x j ) ≤ ∑ i≤n g(x i • x i ) = ng(1) ≤ n,
whence n ≤ 1/c 0 . This suggests that the problem max{c | (i)-(iii)} ( * ) is relevant, as higher values for c 0 correspond to tighter bounds. We look for a function f written as a linear combination of functions φ h ∈ F , and introduce the coefficients c h so that f (t) = ∑ h c h φ h (t).

We now come to the family F : Delsarte's LP bounding techniques require an orthogonal family of polynomials. In the KNP case, this family consist of Gegenbauer polynomials C (λ ) h (t) [1, p. 776], that encode certain properties of S K-1 . For example, so far our description of Delsarte's LP has failed to include any information about K, which is provided by this choice of F , notably by setting λ = (K -2)/2 [19, §3]. Another crucial property is that if a function f is a conic combination of Gegenbauer polynomials, then ( f (x i • x j ) i j ) 0, whence condition (i) of Thm. 1 holds; moreover, conversely, any function satisfying ( f (x i • x j ) i j ) 0 can be written as a conic combination of Gegenbauer polynomials [START_REF] Schoenberg | Positive definite functions on spheres[END_REF]. We therefore adjoin the constraints ( ‡) quantified over G K = {C (K-2)/2 h (t) | h ∈ H} (for some set H) to the LP ( †).

We explicitly write the LP ( * ) by requiring that the f appearing in Thm. 1 should be a conic combination of elements of G K :

max c≥0 {c 0 | ∀t ( f (t) = ∑ φ h ∈G K c h φ h (t) ∧ c 0 + f (t) ≤ 0) ∧ c 0 + f (1) ≤ 1}. (4)
As observed in [START_REF] Odlyzko | New bounds on the number of unit spheres that can touch a unit sphere in n dimensions[END_REF], ( †) and ( * ) are a pair of dual LPs. You have to "massage" ( * ) somewhat before the duality relation with (4) becomes apparent, though (eliminate c 0 and minimize ∑ h c h , see [2, Eq. ( 5), p. 615]). By duality, we only focus on one LP, namely Eq. (4).

Getting your hands dirty

As stated, Eq. ( 4) is infinite in both dimensions: F is countably infinite (if we restrict λ to be integer) and t varies in the uncountably infinite set T = [-1, 1/2]∪{1}. Only a couple among the many works in the literature mention this as a difficulty (without providing a discussion, however). The only paper I found that provides a satisfactory discussion of this issue is [START_REF] Antreicher | The thirteen spheres: a new proof[END_REF].

1. For an upper bound to kn(K), start with polynomials in G K up to degree, say, 15 [START_REF] Odlyzko | New bounds on the number of unit spheres that can touch a unit sphere in n dimensions[END_REF], and gradually work your way up the degrees to see if the bounds improve (be wary of floating point errors arising from evaluating polynomials large degrees -they may invalidate the bound, see [2, Lemma 1]). You can also use any polynomial from G for ≥ K (though not for < K), since Gegenbauer polynomials having lower λ "dominate" -in the sense that they yield larger values of c 0 -those for higher λ . Essentially, the minimum for which you choose elements in G determines the dimensionality of the KNP instance you are going to compute a bound for. 2. A safe way to deal with the fact that T is infinite is to choose any finite T ⊆ T , as removing constraints yields a relaxation (and hence a valid bound). A few experiments will convince you that this discretization has a very small impact on the value of the bound -moreover, since these are "small LPs", you can still instantly obtain answers even when |T | is pretty large (e.g. O(10 5 )).

In order to obtain closed form expressions for Gegenbauer polynomials, I used

Mathematica's GegenbauerC[h,λ ,t]. However, Eq. ( 4) expects Gegenbauer polynomials to be normalized so that C λ h (1) = 1, whereas Mathematica's implementation normalizes them differently -you have to remember to evaluate GegenbauerC

Extensions

I would like to emphasize three among the extensions to Delsarte's LP bound:

(a) Oleg Musin's extension [START_REF] Musin | The kissing number in four dimensions[END_REF], which brought us the proof that kn(4) = 24; (b) Pfender's extension [START_REF] Pfender | Improved delsarte bounds for spherical codes in small dimensions[END_REF], yielding new upper bounds for K ∈ {10, 16, 17, 25, 26}; (c) Bachoc and Vallentin's extension [START_REF] Bachoc | New upper bounds for kissing numbers from semidefinite programming[END_REF], which brought us kn(5) ≤ 45 (further improved to kn(5) ≤ 44.99899685 = 44 in [START_REF] Mittelmann | High-accuracy semidefinite programming bounds for kissing numbers[END_REF], obtained using the GNU Multiple Precision (GMP) arithmetic library: if you trust the GMP library, you should also trust this bound.

The first two extensions are based on enriching the function family F : in case (a) by including some polynomials f such that condition (ii) of Thm. 1 is violated for a fixed x i and finitely many x j (the exceptions that ensue are dealt with combinatorially); and in case (b) by adding a new family of functions to F that are not convex combinations of Gegenbauer polynomials but satisfy the conditions of Thm. 1. The case (c) is even more interesting as it considers distance distributions on triplets of vectors (rather than pairs) and derives a semidefinite programming (SDP) formulation instead of an LP. See [START_REF] Boyvalenkov | A survey on the kissing numbers[END_REF] for further extensions.

  [h,λ ,t] / GegenbauerC[h,λ ,1]. I obtained the following (standard) LP upper bounds using Gegenbauer families C 0.5 h and C 1 h respectively, for h ≤ 10: kn(3) ≤ 13.1583 , kn(4) ≤ 25.5581 and kn(5) ≤ 46.3365 .