Esteban Salgado
email: esteban.salgado@polytechnique.edu

Andrea Scozzari
email: andrea.scozzari@unicusano.it

Fabio Tardella
email: fabio.tardella@uniroma1.it

Leo Liberti
email: liberti@lix.polytechnique.fr

Alternating Current Optimal Power Flow with Generator Selection

Keywords: smart grid, semidefinite programming, diagonal dominance, dimensionality reduction

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

The Alternating Current Optimal Power Flow (ACOPF) problem is as follows: given an electric power network consisting of nodes (called buses) and links (called lines) one seeks an optimal generation and distribution plan of active and reactive power under physical constraints (Ohm's and Kirchhoff's laws), and subject to power generation, voltage magnitude and current bounds on each line.

Not every bus can produce power. Those which can are called generators. There is often a planning issue related to their activation and deactivation. This is important because of two reasons. First, the minimum amount of power produced by a generator may be a (reasonably large) positive constant, so disabling a generator is not equivalent to keeping it active at minimum levels. Secondly, there can be a cost to keep it activated. Modelling this choice implies the addition of binary variables to the model, which yields a Mixed-Integer Quadratically Constrained Quadratic Programming (MIQCQP) problem.

Most of the recent ACOPF literature [START_REF] Bai | Semidefinite programming for optimal power flow problems[END_REF][START_REF] Chen | Bound tightening for the alternating current optimal power flow problem[END_REF][START_REF] Kuang | Alternative lp and socp hierarchies for acopf problems[END_REF][START_REF] Lavaei | Zero duality gap in optimal power flow problem[END_REF] ignores the generator activation/deactivation issue: every available generator is always active. Equivalently, This work was partially supported by the PGMO grant "Projet 2016-1751H", by a Siebel Energy Institute Seed Grant, and by Università di Roma "La Sapienza" under a visiting grant for the last author.

every binary variable is fixed to 1, which yields a continuous Quadratically Constrained Quadratic Programming (QCQP) formulation. We were drawn to the mixed-integer ACOPF variant by [START_REF] Ruiz | A progressive method to solve large-scale ac optimal power flow with discrete variables and control of the feasibility[END_REF], where the authors also consider the selection of Phase-Shifting Transformers (PST), used to control the flow of real power, as well as shunts, which are stabilizing devices. We do not consider PSTs nor shunts in this paper: they require technical data (specifically, more information about the admittance matrices) which we do not possess at this stage. From a theoretical point of view, however, the activation/deactivation of generators and the selection of PSTs yield formulations of the same MIQCQP class, that are likely to require similar solution approaches.

In this paper we study the ACOPF with selection of generators (ACOPFG). Based on the ideas in [START_REF] Ahmadi | Sum of squares basis pursuit with linear and second order cone programming[END_REF][START_REF] Ahmadi | Sum of squares basis pursuit with linear and second order cone programming[END_REF][START_REF] Ahmadi | Dsos and sdsos optimization: more tractable alternatives to sum of squares and semidefinite optimization[END_REF], we will derive Mixed-Integer Linear Programming (MILP) formulations using Diagonally Dominant Programming (DDP) for inner and outer approximations for the ACOPF with binary variables. More precisely, we propose the following relaxations:

1. a Mixed-Integer Semidefinite Programming (MISDP) relaxation of the original MIQCQP formulation of the ACOPFG; 2. an inner Diagonally Dominant Programming (DDP) [START_REF] Ahmadi | Sum of squares basis pursuit with linear and second order cone programming[END_REF][START_REF] Dias | Diagonally dominant programming in distance geometry[END_REF] approximation of the MISDP relaxation; 3. an outer DDP relaxation obtained by replacing the primal DD cone in the inner DDP approximation with its dual.

We shall exhibit some computational results showing that DDP yields tight upper and lower bounds on the optimal objective function value of the original MIQCQP.

One of the important reasons for the usefulness of a tight lower objective function bound is that ACOPF/ACOPFG problems are sometimes solved as lower-level subproblems of bilevel problems where the upper-level decisions concern unit electricity prices. A non-guaranteed heuristic ACOPF/ACOPFG solution might be detrimental to finding good solutions of the upper-level problem with cutting plane approaches. Since we cannot hope to find a guaranteed global optimum in reasonable times, a tight lower bound represents a good trade-off.

Notation

We remark first that alternating currents are commonly modelled by means of scalar and vector quantities over the complex field C. On the other hand, physicists have always denoted current by i. The ambiguity with the usual notation for √ -1 is resolved by denoting the latter by j. Accordingly, we shall refrain from using j as an index (as is common in Mathematical Programming). In agreement with standard notation in complex numbers, we denote complex conjugation by means of a bar over the corresponding quantity: ī, for example, is the complex conjugate of the current, rather than -√ -1. Complex conjugation is applied componentwise to tensors.

We often neglect to explicitly mention matrix sizes in formulations, in order to lighten the notation. Symbols denoting classes of matrices, such as S, D, etc., really stand for "the subset of matrices of appropriate size with the corresponding property".

Contents

The rest of this paper is organized as follows. In Sect. 2 we discuss formulations of the ACOPFG. In Sect. 3 we introduce a MISDP relaxation. In Sect. 4 we present our new inner and outer relaxations based on DDP. In Sect. 5 we briefly discuss dimensionality reduction issues in order to retrieve a good solution from the relaxation output. In Sect. 6 we report computational results.

ACOPF formulations

Consider an electric power network with set of buses N , set of generators G ⊆ N and set of lines L ⊆ N × N . The parameters of the problem are: We want to find the optimal values of active and reactive power at each generator g ∈ G that is switched on (S

G g ∆ = P G g + jQ G g
). We must also decide the current and the voltage for the nodes and lines of the system. Given that:

∀n ∈ N i n = m∈N Y nm v m ∀n ∈ N G v n īn = -S D n ∀n ∈ G v n īn = S G n -S D n
we can model the problem by only using voltage variables v n ∈ C (for n ∈ N ) [START_REF] Salgado | Fast relaxations for Alternating Current Optimal Power Flow[END_REF] and binary variables z g (for g ∈ G) representing the activation of a generator to produce power. The formulation of the ACOPFG is as follows:

min v∈C |N | z∈{0,1} |G| f (v, z) ∀n ∈ N G m∈N v n vm Ȳnm = -S D n ∀g ∈ G m∈N v g vm Ȳgm -S max g z g ≤ -S D g ∀g ∈ G m∈N v g vm Ȳgm -S min g z g ≥ -S D g ∀n ∈ N v n vn ≤ v max n ∀n ∈ N v n vn ≥ v min n ∀(n, m) ∈ L k, ∈N vk v Ȳ f (n,m)k Y f (n,m) ≤ i max nm ∀(n, m) ∈ L k, ∈N vk v Ȳ s (n,m)k Y s (n,m) ≤ i max nm ,                                              (ACOPFG C ) ( 1
)
where f is a real-valued function that is commonly considered as a polynomial on the total active power generated, i.e.:

f (v, z) ∆ = g∈G   cg,0zg + 1≤ ≤L c g, (P G g )   = g∈G   cg,0zg + ≤L c g,   Re   m∈N vg vm Ȳgm       .
We shall consider polynomials which are quadratic in v g and linear in P G g (e.g. L = 1). This makes (ACOPF C ) a complex-valued MIQCQP.

By doubling the number of the variables (which implies quadrupling the sizes of the matrices in each constraint), the problem can be exactly reformulated into a real-valued MIQCQP:

min v∈R 2|N | z∈{0,1} |G| v Cv + c z ∀k ∈ E v A k v = a k ∀ ∈ I v B v ≤ b ∀w ∈ Z v Q w v -q max w z w/2 ≤ q w ∀w ∈ Z v Q w v -q min w z w/2 ≥ q w ,                  (ACOPFG R ) ( 2
)
where E is the set indexing the 2|N G| equality constraints, Z is the set indexing the 2|G| inequalities with binary variables and I is the set indexing the other 2(2|G| + 2|N | + 2|L|) inequality constraints. We remark that this problem is non-convex. Eq. ( 2) therefore cannot generally be solved globally using local optimization methods (as a convex QCQP would). Indeed, it is well known that the decision problem associated to Eq. ( 2) is NP-hard [START_REF] Vavasis | Quadratic programming is in NP[END_REF] 3 MISDP relaxation Semidefinite Programming (SDP) [START_REF] Alizadeh | Interior point methods in semidefinite programming with applications to combinatorial optimization[END_REF][START_REF] Klerk | Aspects of Semidefinite Programming[END_REF] is widely employed in order to derive relaxations of the continuous version of the ACOPF [START_REF] Bai | Semidefinite programming for optimal power flow problems[END_REF][START_REF] Chen | Bound tightening for the alternating current optimal power flow problem[END_REF][START_REF] Lavaei | Zero duality gap in optimal power flow problem[END_REF]. We derive a corresponding MISDP relaxation for the ACOPFG. First we rewrite Eq. ( 2) by replacing the products of the form v M v as M, V = tr(M V ) and adding the rank constraint V = vv : min

V ∈S C, V + c z ∀k ∈ E A k , V = a k ∀ ∈ I B , V ≤ b ∀w ∈ Z Q w , V -q max w z w/2 ≤ q w ∀w ∈ Z Q w , V -q min w z w/2 ≥ q w V = vv ,                  (ACOPFG R ) ( 3
)
where S is the set of all n × n symmetric matrices. Given that V does not encode any integrality constraint, Eq. ( 3) is essentially a MILP with the additional (non-convex) constraint V = v v , which states that V must be positive semidefinite (PSD) and have rank 1. The standard MISDP relaxation is obtained by relaxing the rank 1 constraint to V -v v 0, i.e. by requiring that the Schur complement of V and v, defined as

S (V, v) = 1 v v V is PSD: min V ∈S,v∈R 2|N | z∈{0,1} |G| C, V + c z ∀k ∈ E A k , V = a k ∀ ∈ I B , V ≤ b ∀w ∈ Z Q w , V -q max w z w/2 ≤ q w ∀w ∈ Z Q w , V -q min w z w/2 ≥ q w S (V, v) 0                      (ACOPFG MISDP ) (4)
Most SDP solvers are unable to solve MISDP formulations directly (with some exceptions, e.g. Pajarito [START_REF] Lubin | Extended formulations in mixedinteger convex programming[END_REF] and PENLAB [START_REF] Fiala | PENLAB: A MATLAB solver for nonlinear semidefinite optimization[END_REF]). A possible workaround consists in reformulating the binary constraints on z by ∀g ∈ G (z 2 g = z g ). These can in turn be relaxed, within the SDP framework, with the constraints:

• diag(Z) = z • S (Z, z) 0,
This yields the following SDP relaxation:

min V ∈S,v∈R 2|N | Z∈S,z∈R |G| C, V + c z ∀k ∈ E A k , V = a k ∀ ∈ I B , V ≤ b ∀w ∈ Z Q w , V -q max w z w/2 ≤ q w ∀w ∈ Z Q w , V -q min w z w/2 ≥ q w diag(Z) -z = 0 S (V, v) 0 S (Z, z) 0.                              (ACOPFG SDP ) (5)
While Eq. ( 5) can be tackled by an SDP solver, there are two issues with it: (a) the fact that today, current SDP solver technology is far from allowing the systematic solution of ACOPF instances of even moderate sizes (into the hundreds or thousands of nodes and beyond); (b) the proposed SDP relaxation of the binary variables usually yields very poor bounds. Instead, we shall investigate below some MILP formulations derived from Eq. ( 5). While MILP is also NPhard, its state-of-the-art solvers are much faster than MIQCQP solvers.

Diagonal dominance

A real square symmetric matrix Y = (Y j ) is diagonally dominant (DD) if ∀j Y jj ≥ =j |Y j |. (6)
The main observation leading to our tight MILP relaxations of the MIQCQP is the well-known fact that every diagonally dominant matrix is PSD [START_REF] Gershgorin | Über die abgrenzung der eigenwerte einer matrix[END_REF]. This means that by replacing the constraint S (V, v) 0 with membership in the cone of DD matrices (also called the DD cone) we obtain an inner approximation of the SDP cone. We restrict our attention to two formulations of these Diagonally Dominant Programs (DDP) applied to (ACOPFG MISDP ): a formulation based on the extreme rays of the DD cone [START_REF] Barker | Cones of diagonally dominant matrices[END_REF], and the outer approximation obtained by replacing the extreme rays of the DD cone with the extreme rays of the dual DD cone. We devote particular attention to the latter, since it has two desirable properties: (i) it provides a relaxation of the MISDP, and (ii) we found empirically that it provides tight lower bounds. From Eq. ( 4) we derive the inner approximation:

min V ∈S , v∈R 2|N | , z∈{0,1} |G| C, V + c z ∀k ∈ E A k , V = a k ∀ ∈ I B , V ≤ b ∀w ∈ Z Q w , V -q max w z w ≤ q w ∀w ∈ Z Q w , V -q min w z w ≥ q w V ∈ D,                    (ACOPFG inner ) (7)
where D denotes the DD cone. Given that the decision variable vector v only appears in the conic constraint S (V, v) 0, we can replace the latter by V 0 without modifying the feasible region. We then inner-approximate it by means of the constraint V ∈ D, which is equivalent to Eq. ( 6) and can be written as follows using linear constraints:

∀j V jj ≥ =j Z j -Z ≤ V ≤ Z. (8)
By [START_REF] Barker | Cones of diagonally dominant matrices[END_REF] and [2, Lemma 3.2], we know that the extreme rays of D are square matrices defined as follows: ∀ j, such that j < ,

• E j = e j e j • E + j = (e j + e ) (e j + e ) • E - j = (e j -e ) (e j -e )

In other words, the rays either have one nonzero entry in the diagonal which is equal to 1 or have a single nonzero principal minor

1 ±1 ±1 1
where, for j < , the ones on the diagonal are in positions (j, j) and ( , ), and the ±1 components are in positions (j, ) and ( , j). The main relevant result is that every matrix M ∈ D can be written as a non-negative combination of these extreme rays:

∃ δ j , δ + j , δ - j ≥ 0 M = j δ j E j + j< (δ + j E + j + δ - j E - j ). ( 10
)
This allows us to write membership of a square symmetric matrix M in the dual DD cone D * by means of linear constraints as stated in the following result. Lemma 4.1 (From paper [START_REF] Salgado | Fast relaxations for Alternating Current Optimal Power Flow[END_REF]). The dual DD cone D * can be written as follows:

D * = {M ∈ S | ∀x ∈ X (x M x ≥ 0)}, ( 11
)
where X = {(e j ) j , (e j ± e

) j< } is a set of |N | + 2 |N | 2 = |N | 2 elements.
From Eq. ( 4) we derive the outer approximation:

min V ∈S,v∈R 2|N | z∈{0,1} |G| C, V + c z ∀k ∈ E A k , V = a k ∀ ∈ I B , V ≤ b ∀w ∈ Z Q w , V -q max w z w/2 ≤ q w ∀w ∈ Z Q w , V -q min w z w/2 ≥ q w V ∈ D * .                      (ACOPFG outer ) (12)
By our description of the extreme rays of D, the constraint V ∈ D * can be written using Lemma 4.1. Eq. ( 12) is an outer approximation of Eq. ( 4) because D ⊆ S + = (S + ) * ⊆ D * , where S + is the cone of PSD matrices.

Iterative inner approximation

There are two potential issues with the inner approximation Eq. ( 7): (a) given that we are modifying the conic constraints of Eq. ( 4) by an inner approximation of the original cone, the resulting problem may be infeasible; (b) even if the problem turns out to be feasible, the solution may not be a good approximation of the MISDP. The authors of [START_REF] Ahmadi | Sum of squares basis pursuit with linear and second order cone programming[END_REF][START_REF] Ahmadi | Sum of squares basis pursuit with linear and second order cone programming[END_REF] introduce a possible way around these issues by means of iteratively solving a sequence of auxiliary problems, until a cone C ⊆ S + is found on which Eq. ( 4) is feasible when replacing the conic constraint V ∈ D by V ∈ C. The auxiliary problems are obtained by varying the (square matrix) parameter U in the formulation below:

min V ∈S,z∈{0,1} |G| δ∈(R + ) 4|N | 2 C, V + c z ∀k ∈ E A k , V = a k ∀ ∈ I B , V ≤ b ∀w ∈ Z Q w , V -q max w z w/2 ≤ q w ∀w ∈ Z Q w , V -q min w z w/2 ≥ q w U x∈X δ x xx U = V.                          (ACOPFG inner (U )) (13)
We remark that, while U x∈X δ x xx U is an appropriate description of V ∈ D * , implementation performances improve if we use an additional auxiliary matrix variable:

∀i

W ii = δ i + j =i (δ + ij + δ - i,j ) ∀i, j W ij = δ + ij -δ - ij V = U W U (14)
The iterative procedure is the following:

1. Solve (ACOPFG inner ), obtaining the solution V * ; 2. Define U = chol(V * ); 3. Solve (ACOPFG inner (U )) to obtain V * . While the solution is improving, repeat from Step 2.

At the k-th iteration of this procedure, we consider the parametrization matrix U k and the optimum V * k of (ACOPFG inner (U k )). Then the following holds [START_REF] Ahmadi | Sum of squares basis pursuit with linear and second order cone programming[END_REF][START_REF] Ahmadi | Sum of squares basis pursuit with linear and second order cone programming[END_REF].

Proposition 4.2. V * k is feasible for (ACOPFG inner (U k+1 )) and cannot worsen the current objective value: val(ACOPFG inner (U k+1 )) ≤ val(ACOPFG inner (U k )). Moreover, if V * k is PSD and val(ACOPFG inner (U k )) is not optimal, then V * k is improving, i.e. val(ACOPFG inner (U k+1 )) < val(ACOPFG inner (U k )).

In particular, we know that if we start from a feasible problem, the iterated problems will always remain feasible. To find an initial U such that the parametrized problem is feasible, we solve a slightly different formulation, where we add a slack variable that we minimize in order to attempt to achieve feasibility:

min V ∈S,z∈{0,1} |G| δ∈(R + ) 4|N | 2 ,α∈R 4|N | 2 α ∀k ∈ E A k , V + αI = a k ∀ ∈ I B , V + αI ≤ b ∀w ∈ Z Q w , V + αI -q max w z w ≤ q w ∀w ∈ Z Q w , V + αI -q min w z w ≥ q w U x∈X δ x xx U = V + αI.                          (ACOPFG ph1 (U )) (15)
Algorithmically, we start with U = I. While α * > 0, we set U = chol(V * + α * I).

We solve (ACOPF ph1 (U )) using the procedure previously stated. The pseudocode for our inner approximation algorithm is given in Alg. 1.

Algorithm 1 Iterative inner approximation

1: procedure innerACOPFMI 2: ← 10 -5 , aux ← 1, U ← I 3:
# Achieve feasible U first 4:

while aux > 0 do 5:

(V * , z * , α * ) ← optimal solution of (ACOPFG ph1 (U )) 6:

aux ← α * , U ← chol(X + α * I) 7:

valf ← +∞, vali ← 0 8:

# Now improve the approximation 9:

while |vali -valf| > do 10:

vali ← valf 11:

((V * , z * ), valf ) ← optimal solution and value of (ACOPFG inner (U )) 12:

U ← chol(V * ) return ((V * , z * ), valf)

Negative rank of outer approximation

One of the issues with Eq. ( 12) is that, being an outer approximation of Eq. ( 4), it does not ensure that V 0. In fact, we empirically found that, although the bound was often tight, the solution V * had a considerable number of negative eigenvalues. We therefore propose to add the following cutting planes to the formulation: for each eigenvector p of V * corresponding to a strictly negative eigenvalue we add the following inequality to the formulation Eq. ( 12):

p V p ≥ 0. ( 16
)
These cuts make V * infeasible since, by definition of negative eigenvector, we have that p V * p < 0. Unfortunately, we found that in practice these cuts, by themselves, do not yield a PSD solution. Nevertheless, this procedure decreases the the number and/or the sum of the strictly negative eigenvalues.

The pseudocode for our outer approximation algorithm is given in Alg. 2. The function SpectralDecomposition() returns the spectral decomposition of its argument.

Algorithm 2 Iterative outer approximation 1: procedure outerACOPFMI 2:

P ← ∅ 3:

λ-← 1 4:

while λ-> 0 and time < maxtime do 5:

((V * , z * ), valf) ← optimal solution of (ACOPFGouter) 6:

(W, Λ) ← SpectralDecomposition(V * ) 7:

λ-← number of negative eigenvalues 8:

P-← eigenvectors associated to the negative eigenvalues 9:

P ← P ∪ Preturn (X * , z * , valf)

Dimensionality reduction

All our methods solve either inner or outer matrix approximations of the original MIQCQP Eq. ( 2), derived from Eq. ( 3) by relaxing the rank constraint. As such, they will not provide a solution vector v * for the voltage in the original problem Eq. ( 2), but rather a solution matrix V * , which will be very unlikely to have rank 1.

In order to heuristically extract a vector v * from a matrix V * with higher rank we considered two options for a dimensionality reduction algorithm: (a) Principal Component Analysis (PCA) [START_REF] Jolliffe | Principal Component Analysis[END_REF] and (b) Barvinok's naive algorithm [7, §5]. Both produce some estimate v of v * according to different analyses. We then used v as a starting point for a local descent carried out by a local Nonlinear Programming (NLP) solver.

Both options are simple and effective in different settings. PCA is well known and needs no introduction. Barvinok's algorithm is as follows:

1. factor V * into V * = F F ; 2. sample a random vector y componentwise from a standard Gaussian distribution N (0, 1); 3. let v = F y.

The analysis carried out by A. Barvinok shows that v has high probability of being "not too far" from the feasible set. We remark that we applied Barvinok's algorithm to a wrong setting (in general), since the analysis only holds whenever V * is a solution of an SDP, and some of our algorithms do not ensure PSD solutions. This being "just a heuristic", we proceeded nonetheless, based on the fact that the formulations that produce V * are themselves derived from SDP formulations.

  magnitude at bus n • ∀(n, m) ∈ L i max nm upper bound on current on line (n, m) • Y bus admittance matrix in C |N |×|N | • Y f , Y s line admittance matrices in C |L|×|N | .

Some empirical tests showed a promise of superiority of Barvinok's naive algorithm to PCA (also see [START_REF] Liberti | Barvinok's naive algorithm in distance geometry[END_REF]), which was therefore our dimensionality reduction method of choice in obtaining computational results (Sect. 6).

Numerical results

We tested our approach using Cplex-12.6.3 [START_REF]IBM: ILOG CPLEX 12.6 User's Manual[END_REF] on YALMIP [START_REF] Löfberg | YALMIP: A toolbox for modeling and optimization in MATLAB[END_REF] for the inner and outer approximations. The MIQCQP was solved with Bonmin-4.8.4 [START_REF] Bonami | BONMIN user's manual[END_REF] on AMPL [START_REF] Fourer | The AMPL Book[END_REF]. Our results were obtained on an Intel i7 dual-core CPU at 2.1GHz with 15GB RAM.

The time limits for decreasing the negative rank (outer approximation) and feasible solution improvement (inner approximation) were set to 300 seconds. The limit of time for solving MIQCQPs was set to 1200 seconds.

Inner-approximation

Outer As we can observe the bounds we obtained are tight with respect to the feasible solution we obtained for the MIQCQP formulation. In general, inner approximations seem to provide solutions with a smaller rank than outer ones. Unfortunately, we could not scale these experiments to larger sizes because of the slow convergence of the loop from Step 4 to Step 6 of Alg. 1. It is for this reason that we write "-" on some instances for the inner approximation. When we write "(x)" we mean that the local solver was not able to find a local optimum in the time limit.


Pierre-Olivier Bauguion
email: pierre-olivier.bauguion@irt-systemx.fr

Claudia D'ambrosio
email: dambrosio@lix.polytechnique.fr

Leo Liberti
email: liberti@lix.polytechnique.fr

Maximum concurrent flow with incomplete data

Keywords: Maximum concurrent flow, multi-commodity flow problems, incomplete data, unknown data, uncertainty, inverse optimization, transportation systems

The Maximum Concurrent Flow Problem (MCFP) is often used in the planning of transportation and communication networks. We discuss here the MCFP with incomplete data. We call this new problem the Incomplete Maximum Concurrent Flow Problem (IMCFP). The main objective of IMCFP is to complete the missing information assuming the known and unknown data form a MCFP and one of its optimal solutions. We propose a new solution technique to solve the IMCFP which is based on a linear programming formulation involving both primal and dual variables, which optimally decides values for the missing data so that they are compatible with a set of scenarios of different incomplete data sets. We prove the correctness of our formulation and benchmark it on many different instances.

Introduction

Network flows have been introduced long ago (see, e.g., [START_REF] Ford | Maximal flow through a network[END_REF], enhanced later by [START_REF] Edmonds | Theoretical improvements in algorithmic efficiency for network flow problems[END_REF]) to tackle single commodity flow problems, such as the max-flow problem. Since then, these models have been generalized for multiple commodities [START_REF] Shahrokhi | The maximum concurrent flow problem[END_REF] and grouped under the label of multi-commodity flow models. Nowadays, multicommodity flow formulations are extensively used in many contexts for their ability to capture the movements of different types of commodities in various real-world activities such as people in transportation models, data in telecommunication networks, water flows... (see, e.g., [START_REF] Ahuja | Network flows: theory, algorithms, and applications[END_REF]). These formulations are generally used to help make the best cost-effective solution for allocating resources; this leads to optimize a cost function.

In real-world applications, the available data are often uncertain or incomplete, and their actual values may only be revealed at a time when the overall decision strategy has already been chosen. This is often the case in transportation systems where the parameters are time-dependent and event-sensitive. Statistical inference and data mining represent convenient ways to deal with this uncertainty. One of the best known inference models in transportation systems is the Four Step Model [START_REF] Manheim | Fundamentals of transportation systems analysis[END_REF], which is an algorithm that iterates over time according to an equilibrium criterion. More recently, a lot of attention has been devoted to machine learning approaches, which generally performs better on large scale datasets. In this context, [START_REF] Shiliang | A bayesian network approach to traffic flow forecasting[END_REF] proposes bayesian networks and [START_REF] Lv | Traffic flow prediction with big data: A deep learning approach[END_REF] uses a deep learning approach to forecast flow in transportation systems.

However, optimization methods that deal with uncertainty actually do exist. To the best of our knowledge, [START_REF] Dantzig | Linear programming under uncertainty[END_REF] was the first to propose a stochastic approach to tackle incompleteness of input data. It assumes that uncertain data follow some given probability distribution, and that the objective of this approach is not to optimize a certain cost, but an expected cost instead. On the other hand, [START_REF] Soyster | Convex programming with set-inclusive constraints and applications to inexact linear programming[END_REF] proposed a complementary approach by optimizing a robust criterion such as the worst case or the maximum regret. This particular method received renewed attention from [START_REF] Ben-Tal | Adjustable robust solutions of uncertain linear programs[END_REF] and [START_REF] Bertsimas | The price of robustness[END_REF], while [START_REF] Ben-Ameur | Routing of uncertain traffic demands[END_REF] applies this approach to multi-commodity flow problems by considering a polyhedral uncertainty set of demands. Later, [START_REF] Atamturk | Two-stage robust network flow and design under demand uncertainty[END_REF] mixed the recourse variables introduced in [START_REF] Dantzig | Linear programming under uncertainty[END_REF] with the robust approach for a network flow and design problem.

Optimization methods can also be used to optimally fit experimental measurements. In [START_REF] Hooker | Inferring network flows from incomplete information with application to natural gas flows[END_REF], multi-commodity flow optimization is used to model a gas transportation network while retrieving missing data. The problem discussed in [START_REF] Hooker | Inferring network flows from incomplete information with application to natural gas flows[END_REF] consists in recomposing the flow on each arc, knowing only the global amount of incoming and outgoing flows for each node. The problem of finding a minimal adjustment of the cost function to ensure the optimality of a given solution generated a particular interest with [START_REF] Burton | On an instance of the inverse shortest paths problem[END_REF] under the label of inverse optimization. For example [START_REF] Ahuja | Inverse optimization, part I: Linear programming and general problem[END_REF][START_REF] Zhang | Calculating some inverse linear programming problems[END_REF] apply this concept to multi-commodity flow problems (especially min cost flow problem). The survey [START_REF] Heuberger | Inverse combinatorial optimization: A survey on problems, methods, and results[END_REF] on this subject includes situations where the inverse problem seeks parameters other than objective function coefficients.

The Maximum Concurrent Flow Problem (MCFP) has been extensively studied over time [START_REF] Shahrokhi | The maximum concurrent flow problem[END_REF][START_REF] Bienstock | Asymptotic analysis of the flow deviation method for the maximum concurrent flow problem[END_REF][START_REF] Bauguion | Efficient algorithms for the maximum concurrent flow problem[END_REF], but in this paper we present a new approach for finding optimal maximum concurrent flows using incomplete data. Our method seeks optimal solutions and completes the partial input. This problem typically arises when we have insights about the global behavior of a system while data are partially unknown [START_REF] Hooker | Inferring network flows from incomplete information with application to natural gas flows[END_REF]. Symmetrically it can validate/invalidate a hypothetical behavior by comparing it with the observed data. This is particularly relevant in transportation when the routing strategy of passengers is known while data are incomplete. We call this problem Incomplete Maximum Concurrent Flow Problem (IMCFP).

The rest of this paper is organized as follows: Section 2 recalls the MCFP and presents the IMCFP. In Section 3 we propose a formulation for the IMCFP by integrating both primal and dual formulations of the MCFP, and prove its correctness based on the complementary slackness conditions of Linear Programming (LP), which we recall for convenience [START_REF] Kuhn | Nonlinear programming[END_REF]. Then, in Section 4 we present preliminary experiments for these formulations and their practical interest. Finally, Section 5 concludes the paper.

∀i ∈ V {o, d} j∈N -(i) f ji = j∈N + (i) f ij . (1)
We also recall that N -(i) is the set of nodes j such that (j, i) ∈ A and N + (i) such that (i, j) ∈ A, for each i ∈ V .

The MCFP was introduced in [START_REF] Shahrokhi | The maximum concurrent flow problem[END_REF]. It can be formulated as follows using LP:

max γ≥0,f ≥0 γ ∀k = (o, d, D) ∈ K j∈N -(d) f k jd = γD + j∈N + (d) f k dj ∀k = (o, d, D) ∈ K, i ∈ V {o, d} j∈N -(i) f k ji = j∈N + (i) f k ij ∀a = (i, j) ∈ A k∈K f k a ≤ c a .                  (2)
Its dual is (one can choose appropriate inequalities instead of equalities for the first and second sets of constraints of (2)):

min p≥0,u≥0 a∈A u a c a ∀k = (o, d, D) ∈ K, ∀a = (i, j) ∈ A, p k i + u a ≥ p k j ∀k = (o, d, D) ∈ K, p k o = 0 k=(o,d,D)∈K p k d D ≥ 1            (3)
where u a , (for a ∈ A) are the dual variables associated to the capacity constraints (last set of constraints of (2)) and p k i (for k ∈ K, i ∈ V ) are the dual variables associated to the flow conservation constraints (first and second sets of constraints of ( 2)). This implies that the MCFP is polynomial-time solvable (for example with an interior point algorithm). The MCFP is also strongly polynomial-time solvable [START_REF] Norton | Using separation algorithms in fixed dimension[END_REF], but it appears to be common knowledge that, for instance sizes of current practical interest, it is empirically more efficient to use a good LP solver on Problem (2) rather than the algorithm in [START_REF] Norton | Using separation algorithms in fixed dimension[END_REF]. The MCFP is often used in real-life applications in order to design networks or evaluate the arcs with highest risk of becoming saturated [START_REF] Gerla | A cut saturation algorithm for topological design of packet-switched communication networks[END_REF]. The main applied interest in the MCFP is that, through the ratio variable γ, it ensures a fairness of arc capacity utilization over all flows.

Our motivation for studying this problem stems from transportation networks, be they road or rail-oriented. As a critical increase of the load can induce a decrease of the quality of services, an hypothesis consists of assuming that the passengers traffic tends naturally to balance itself to an equilibrium [START_REF] Manheim | Fundamentals of transportation systems analysis[END_REF]. One can model this problem by minimizing the maximum capacity utilization, and the latter can be reformulated as a MCFP [START_REF] Shahrokhi | The maximum concurrent flow problem[END_REF]. In our context, we have historical traffic data including a partial observation of the arc loads for a certain subset of arcs. For some networks, we are also given a subset A ′ ⊂ A of arcs with known capacities.

In general, however, we do not know the arc capacities. The problem we are interested in is the MCFP with incomplete arc capacities. The MCFP in LP formulation (2) without the capacity constraints is clearly an unbounded LP. To avoid this situation, we employ a given set S of scenarios from our historical arc load database. Each scenario s = (A s , ℓ s , K s ) ∈ S consists of a subset A s ⊂ A of arcs, a partial arc load function ℓ s : A s → R + , and a set of commodities K s . We require that: (i) missing capacities should be estimated so as to allow the maximum known arc loads over all scenarios, (ii) arc loads from computed flows should be as close as possible to the loads given in the scenarios, and (iii) each flow solution for a scenario should describe an optimal solution of the MCFP w.r.t. capacity and commodity values. We therefore define the following problem, which is new as far as we could ascertain.

IMCFP. Given a graph G = (V, A), a subset A ′ ⊂ A, a partial arc capacity function c : A ′ → R + , and a set S of scenarios (A s ⊂ A, ℓ s , K s ) where ℓ s : A s → R + and K s is a set of triplets k = (o, d, D) ∈ V 2 × R + , find: a threshold function γ : S → R + , a complementary arc capacity function c : A \ A ′ → R + , and a set of flows f sk (for s ∈ S and for k ∈ K s ) such that (i) for each arc a ∈ A and for each s ∈ S, the arc load of f on a is bounded above by c a ; (ii) for each s ∈ S and arc a ∈ A s , the arc load of f on a is as close as possible to the arc load γ s ℓ s a ; (iii) for each s ∈ S, the flows f s and γ s should be optimal with respect to an MCFP defined over the capacities c a over all a ∈ A and the commodities K s .

Although the IMCFP is natively cast in a multi-objective fashion (see condition (ii)), in practice we minimize a max norm over all arcs and all scenarios. We remark that condition (iii) is only apparently recursive: we want to decide f, c at the same time and also require that every f s should be optimal flows w.r.t. a putative MCFP instance defined over the values of the c variables and the K s parameters. We shall see below that the IMCFP can be formulated by means of a Mixed-Integer Linear Programming formulation that combines both primal and dual variables.

The IMCFP

In this section we shall first introduce a Mixed-Integer Linear formulation for the IMCFP, and then prove its correctness.

Formulation

-Sets:

• V : set of nodes • A: set of arcs • S: set of scenarios • K s : set of commodities for scenario s • A ′ : subset of arcs from which the capacity is known • A s : subset of arcs from which the load is known. -Decision variables:

-Parameters: • k = (o, d, D) for k ∈ K s , s ∈ S: commodity data (origin o ∈ V , destina- tion d ∈ V , demand value D ∈ R + ) • ℓ s : A s → R + :
• f sk a ≥ 0 for a ∈ A, s ∈ S, k ∈ K s : flow variable of arc a ∈ A, for scenario s ∈ S and demand k ∈ K s • γ s ≥ 0 for s ∈ S: threshold variable for scenario s ∈ S • p sk i ≥ 0 for i ∈ V, s ∈ S, k ∈ K s : potential variable (dual variable from MCFP's conservation constraint) for node i, scenario s, and demand k • u s a ≥ 0 for a ∈ A, s ∈ S: weight variable (dual variable from MCFP's capacity constraint) for arc a ∈ A and scenario s ∈ S • w s a ∈ {0, 1} for a ∈ A, s ∈ S: binary variable that allows for a corresponding weight u s a to be greater than 0 or not • x sk a ∈ {0, 1} for a ∈ A, s ∈ S, k ∈ K s : binary variable that allows for a corresponding flow f sk a to be greater than 0 or not • ∆ ≥ 0: maximal difference between the load parameters with the computed ones • c a ≥ 0 for a ∈ A \ A ′ : capacity variable of arc a.

-Objective function:

min f,γ,p,u,w,x,c ∆. (4)
-Constraints:

• flow conservation:

∀s ∈ S, ∀k = (o, d, D) ∈ K s , ∀i ∈ V \ {o, d} j∈N -(i) f sk ji = j∈N + (i) f sk ij (5)
• demand satisfaction:

∀s ∈ S, ∀k = (o, d, D) ∈ K s j∈N -(d) f sk jd - j∈N + (d) f sk dj = γ s D (6)
• min-cost node access:

∀s ∈ S, ∀k ∈ K s , ∀(i, j) ∈ A p sk i + u s ij -p sk j ≥ 0 (7)
• min-cost path condition:

∀s ∈ S, ∀k ∈ K s , ∀(i, j) ∈ A p sk i + u s ij -p sk j ≤ 1 -x sk ij (8)
• origin access:

∀s ∈ S, ∀k = (o, d, D) ∈ K s p sk o = 0 (9)
• complementary slackness condition on capacity:

∀s ∈ S, ∀a ∈ A c a ≤ k∈K s f sk a + M w (1 -w s a ) (10)
• binary flow constraint:

∀s ∈ S, ∀k ∈ K s , ∀a ∈ A f sk a ≤ M f x sk a (11)
• binary weight constraint:

∀s ∈ S, ∀a ∈ A u s a ≤ w s a (12)
• dual weights constraint:

∀s ∈ S a∈A u s a = 1 (13)
• threshold bound:

∀s ∈ S, ∀a ∈ A s c a ≥ γ s ℓ s a (14)
• feasibility bound:

∀s ∈ S, ∀a ∈ A k∈K s f sk a ≤ c a (15)
• max norm (i):

∀s ∈ S, ∀a ∈ A s ∆ ≥ γ s ℓ s a - k∈K s f sk a (16)
• max norm (ii):

∀s ∈ S, ∀a ∈ A s ∆ ≥ k∈K s f sk a -γ s ℓ s a . (17)
The model aims at minimizing the maximal error ∆, which is the absolute value of the difference for each scenario and each arc between the load of f and the load γ s ℓ s a (see constraints ( 16) and ( 17)) as stated in condition (ii) of the IMCFP. Constraints ( 5) and ( 6) are flow conservation constraints already mentioned in [START_REF] Ahuja | Network flows: theory, algorithms, and applications[END_REF]. Constraints ( 7), ( 8), [START_REF] Burton | On an instance of the inverse shortest paths problem[END_REF], and (11) provide the system that ensures the optimality of the routing in the sense of the weights u, which correspond to the dual variables in (3) of capacity constraint in MCFP [START_REF] Ahuja | Inverse optimization, part I: Linear programming and general problem[END_REF]. Constraints [START_REF] Dantzig | Linear programming under uncertainty[END_REF] and ( 12) describe the slackness conditions for the dual variables of the capacities. Constraint (13) ensures a feasible dual solution of (3) for each scenario. Finally, Constraints ( 15) and ( 14) guarantee the feasibility of (2) as condition (i) of IMCFP. Note that c a of ( 14), [START_REF] Hooker | Inferring network flows from incomplete information with application to natural gas flows[END_REF], and (10) can either be a parameter or a variable, depending whether a is in A ′ or not. The condition (iii) of IMCFP follows by Proposition 1 and Theorem 2 (below).

Correctness

Proposition 1 Let f = (f k a | a ∈ A, k ∈ K) ≥ 0 and γ ≥ 0 be a feasible flow solution of the MCFP (2) and p = (p k i | k ∈ K, i ∈ V ) ≥ 0 and u = (u a | a ∈ A)
≥ 0 be a feasible solution of its dual [START_REF] Atamturk | Two-stage robust network flow and design under demand uncertainty[END_REF]. These two assertions are equivalent: Proof. 2⇒1: If f and γ (resp. p and u) are optimal for (2) (resp. for (3)), the complementary slackness conditions [START_REF] Kuhn | Nonlinear programming[END_REF] state immediately the first two equations. The third complementary slackness condition of MCFP states k∈K p k d D = 1. Therefore, (3) ensures there exists a so that u a > 0, meaning a∈A u a > 0.

1. (a) ∀a ∈ A, we have u a (c a - k∈K f k a ) = 0 (b) ∀a = (i, j) ∈ A, ∀k ∈ K, we have f k a (p k i + u a -p k j ) = 0 (c)
1⇒2: If a∈A u a > 0 and ∀a ∈ A, u a ≥ 0, there exists a ∈ A so that u a > 0.

The first condition implies that for this a we have c a = k∈K f k a . We can assume the arcs have a non-zero capacity (otherwise the arc should not have existed), meaning there exists a k = (o, d, D) for which f k a > 0. The second condition implies that going through this arc satisfies the "min cost path condition", meaning

p k d ≥ u a > 0. We can then construct Q = k=(o,d,D)∈K p k d D > 0. Choose u ′ a so that u ′ a = ua Q ,
= (i, j) ∈ A, ∀k ∈ K f k a (p k i + u a -p k j ) = 0
still hold and (u, p) becomes a feasible solution of (3). The proposition follows from the LP complementary slackness conditions [START_REF] Kuhn | Nonlinear programming[END_REF].

⊓ ⊔

Theorem 2 Given a graph G = (V, A) and a feasible solution (f, c) of IMCFP, for each s ∈ S, f s corresponds to an optimal solution of an MCFP with respect to c and K s .

Proof. Let us choose a s ∈ S. f s is feasible for MCFP (2) with parameters c and K s due to ( 5), ( 6), [START_REF] Hooker | Inferring network flows from incomplete information with application to natural gas flows[END_REF]. p s and u s describes a feasible solution of (3) with parameters c and K s due to ( 13) and ( 7) (one can scale the solution to ensure

k=(o,d,D)∈K s p sk d D ≥ 1 )
. The combination of ( 10) and ( 12) ensure the first slackness condition, as u s a = 0 if the c a -k∈K s f sk a > 0. Constraints ( 7)-( 8) validate the second slackness condition. The third slackness condition is obtained by [START_REF] Gerla | A cut saturation algorithm for topological design of packet-switched communication networks[END_REF]. Proposition 1 validates the optimality of f s for MCFP [START_REF] Ahuja | Inverse optimization, part I: Linear programming and general problem[END_REF] with parameters c and K s . As it holds for each s ∈ S, it concludes the proof. ⊓ ⊔

Numerical results

We can solve IMCFP instances by simply formulating them as the Mixed-Integer Linear Program in Sect. 3.1 and solving them using an off-the-shelf solver. By the polynomial number of constraints and variables of our formulation, we know that the decision version of IMCFP is in NP. Although the MCFP is polynomialtime solvable, our IMCFP formulation introduces the need of binary variables to ensure optimality among the different scenarios, meaning it is still an open question whether IMCFP is NP-complete or not. To solve all the following instances, we solved the IMCFP model with the MILP solver IBM CPLEX 12.6 with default settings and a time limit of 1200 CPU seconds on a personal computer (Intel Core i7-6820HQ 2.70 GHz, 16 GB DDR3 RAM). All graphs used in the following experiments are based on the topology of the Paris subway network, often restricted to the left bank. Therefore, each node represents a connection between one or more different metro lines, and each arc represents a section of a line. The network is strongly connected, meaning we can generate complete sets of demands (|K| = (|V | -1)|V |). The demand value for each commodity k ∈ K is an integer uniformly chosen in the interval [START_REF] Ahuja | Network flows: theory, algorithms, and applications[END_REF][START_REF] Dantzig | Linear programming under uncertainty[END_REF].

Firstly, to evaluate the prediction performance of the proposed formulations, we generated integer capacities for the MCFP in the interval [START_REF] Ahuja | Network flows: theory, algorithms, and applications[END_REF][START_REF] Hooker | Inferring network flows from incomplete information with application to natural gas flows[END_REF] and then calculated an optimal solution of MCFP for each scenario of demands, keeping the same capacities C among them. It allows us to give the total configuration of loads ℓ s a and K s as input where an optimal solution with zero ∆ value exists. This let us compare the computed capacities in A \ A ′ , with those we chose to construct our instances of MCFP (the generated ones). Secondly, we generated ℓ s a , K s , and c a according to a feasible flow solution (in terms of conservation and capacity constraints) which is not MCFP compliant (i.e., whose input data do not follow an optimal MCFP solution pattern). The second set of instances aims at observing how our model deals with data based on a wrong hypothesis (structure of an optimal solution of MCFP) and, hence, how the objective value is impacted. Moreover, we studied the impact of the quantity of known and unknown data by giving a fixed percentage of capacities for MCFP ( |A ′ | |A| × 100) as an input.

The resulting tables are organized as follows. The number of nodes (|V |) and arcs (|A|) are reported in the caption of each table. The proportion of known capacities is specified in the first column "C(%)". The rest of the table is divided in 3 subsets of columns, 4 for each value of cardinality of S. The four columns report: (i) the amount of capacities "c(%)" that has been successfully predicted (the percentage of successfully predicted capacities may be lower than the given ones due to the truncation process of "C(%)|A ′ |"); (ii) the maximal absolute gap ("Gap") observed between the predicted capacities and the generated ones; (iii) the CPU time in seconds "T (s)" (we denote termination due to time limit by *); and (iv) objective value "∆". The most striking thing one can note is the effectiveness of our methodology in predicting arc capacities (see Table 1). Even if without insights of capacities the prediction correctness remains low, it skyrockets to more than 50% of additionnal correct predictions with only three scenarios. This can be explained by the fact that even with few capacities the γ s become strongly bounded. This proportion of corrected predictions tends naturally to arise when more and more capacity parameters are known, but the number of scenarios seems to remain a strong source of insights for the network, especially when the number of known capacities is low (see Table 2). However, the number of scenarios seems more and more important as the size of the instance grows. Indeed, the average per-centage of correct predictions is quite close compared from |S| = 4 to |S| = 5 (around 81% for both) in the instances of Table 1 as if this increase did not bring further information. But this average percentage rises from 70% to 79% between |S| = 4 and |S| = 5 in the instances of Table 2. This could mean that the larger the instance, the higher the number of required scenarios to reach a similar prediction performance. It is also interesting to note that some of these results are obtained without even reaching optimality, especially on the total metro network instances, meaning that the routing problems are quite hard to solve even when all capacities are known. This is confirmed by the amount of time consumed by the instances when all the capacities are set. This suggests that MCFPs are also quite hard to solve. This is a well-known result that motivated the design of approximation schemes to solve it (see e.g. [START_REF] Bienstock | Asymptotic analysis of the flow deviation method for the maximum concurrent flow problem[END_REF]). As one can expect, confronted with inconsistent data (input data that do not follow a structure of an optimal solution of the MCFP), the prediction effectiveness does not perform as well as the previous case (see Table 3). Nevertheless, our model often tends to bring a solution that fits significantly the input loads, especially when few capacities are known. Generally speaking, this means that a low objective value does not guarantee necessarily that the MCFP pattern hypothesis on our input data is right. However, when all capacities are set (A ′ = A), this seems to constraint drastically our model whose objective value skyrockets. Therefore, it suggests that our model can still be used to validate/invalidate this hypothesis on condition that all capacities are known.

|S| = 3 |S| = 4 |S| = 5 C(%) c(%) Gap T (s) ∆ c(%) Gap T (s) ∆ c(%) Gap T (s) ∆ 0

Conclusion and future work

In this paper we studied a new problem called the IMCFP and proposed a solution technique to tackle it. The purpose of the new formulation proposed is to find a routing of commodities that fits input data at best (namely loads and demand matrix) as [START_REF] Hooker | Inferring network flows from incomplete information with application to natural gas flows[END_REF]. But, in addition, the solution has to follow an optimal structure of MCFP regarding a set of known and unknown data. The practical interest of the problem we discussed arises in transportation systems when it comes to recompose unknown data assuming hypothesis of data structure (namely a MCFP and one of its optimal solutions). Moreover, this formulation can validate/invalidate hypothesis of a MCFP's optimal routing by confronting it to observed data. More theoretically, we showed that this problem can be tackled by embedding primal and dual formulations and complementary slackness conditions all together into a single Mixed Integer Linear Program. We proved the correctness of such models, and led experiments to emphasize empirical behaviours and computational hardness. Future work will focus on computational performance improvements so that larger graphs can be treated.

2

  The maximum concurrent flow Consider the following (well-known) problem. MCFP. Given a simple directed graph G = (V, A), an arc capacity function c : A → R + , a set K of triplets k = (o, d, D) ∈ V 2 × R + , find: a scalar γ ≥ 0 (called threshold ) and a set of flows f k on G for each k ∈ K such that (i) for each arc a ∈ A, the arc load of f on a (i.e., the sum of the flows on a) does not exceed the arc capacity c a ; (ii) each "od" flow f k has value γD for each k = (o, d, D) ∈ K; (iii) γ is maximum. We recall that, for a given k = (o, d, D) ∈ K, a flow f having value D in a graph G from node o to node d is a non-negative arc function f : A → R + such that j∈N -(d) f jd -j∈N + (d) f dj = D and the following flow balance equations hold (we omit k index for clarity):

  load function over the arcs A s for each scenario s ∈ S • c : A ′ → R + : capacity function over the arcs A ′ • M w : "Big M" parameter associated to the binary weights w • M f : "Big M" parameter associated to the binary flows x.

a∈A u a > 0 2 .

 2 f and γ are optimal for the MCF (2), and p and u (both scaled by k=(o,d,D)∈K p k d D) are optimal for its dual (3).

  and, for each k = (o, d, D) ∈ K, p ′ k d = p k d Q . Then we trivially have k=(o,d,D)∈K p ′ k d D = 1. Substituting u by u ′ and p by p ′ , one can observe that ∀a ∈ A, u a (c a -k∈K f k a ) = 0 and ∀a

Table 1 .

 1 IMCFP, Paris left bank subway network topology (|V | = 13, |A| = 38).

	42.11 13.00 2.95 0.00 42.11 14.00 4.31 0.00 63.16 7.00 6.11 0.00
	10 60.53 9.72 2.64 0.00 65.79 11.00 3.07 0.00 63.16 9.00 7.43 0.00
	20 68.42 11.00 2.12 0.00 71.05 10.00 118.00 0.00 73.68 10.00 4.82 0.00
	30 60.53	14 1.83 0.00 78.95 8.58 4.48 0.00 76.32 9.00	* 0.04
	40 44.74 12.83 2.18 0.00 81.58 6.00 3.00 0.00 73.68 13.00 5.40 0.00
	50 84.21 10.46 2.14 0.00 86.84 6.00 3.31 0.00 78.95 11.77 6.80 0.00
	60 81.58 10.00 2.59 0.00 81.58 7.00 4.12 0.00 89.47 13.00 14.09 0.00
	70 81.58 9.26 2.01 0.00 89.47 9.00 2.11 0.00 84.21 12.72 6.05 0.00
	80 86.84 12.00 2.56 0.00 97.37 6.00 3.09 0.00 94.74 11.25 6.63 0.00
	90 100.00 0.00 1.90 0.00 97.37 0.69 3.46 0.00 92.11 4.00 9.95 0.00
	100 100.00 0.00 1.42 0.00 100.00 0.00 2.47 0.00 100.00 0.00 3.74 0.00


Leo Liberti
email: liberti@lix.polytechnique.fr

Ky Vu

Barvinok's naive algorithm in Distance Geometry $

Keywords: distance geometry, concentration of measure, protein structure

In 1997, A. Barvinok gave a probabilistic algorithm to derive a near-feasible solution of a quadratically (equation) constrained problem from its semidefinite relaxation. We generalize this algorithm to handle matrix variables instead of vectors, and to twosided inequalities instead of equations. We derive a heuristic for the distance geometry problem, and showcase its computational performance on a set of instances related to protein conformation.

Introduction

We consider the following Distance Geometry Problem (DGP). Given an integer K > 0 and a simple graph G = (V, E, d) where d : E → R + , decide whether there is a realization x :

V → R K such that ∀{i, j} ∈ E x i -x j 2 2 = d 2 i j . (1)
Let n = |V| and m = |E|. The DGP is well known in the literature [START_REF] Liberti | Molecular distance geometry methods: from continuous to discrete[END_REF][START_REF] Liberti | Euclidean distance geometry and applications[END_REF], as it serves as a model for several applications (e.g. to sensor networks [START_REF] Krislock | Explicit sensor network localization using semidefinite representations and facial reductions[END_REF], molecular conformation [START_REF] Cucuringu | Eigenvector synchronization, graph ridigity and the molecule problem[END_REF][START_REF] Cassioli | An algorithm to enumerate all possible protein conformations verifying a set of distance constraints[END_REF] and more [START_REF] Billinge | Assigned and unassigned distance geometry: Applications to biological molecules and nanostructures[END_REF]). A natural extension of the DGP is the interval version (denoted by iDGP), where d : E → IR + associates intervals (instead of scalars) to edges. This variant is the one often used in applications, since intervals naturally model data uncertainty and noise [START_REF] Berger | Reconstructing a three-dimensional model with arbitrary errors[END_REF]. We adapt the naive algorithm proposed by A. Barvinok in [8, §5] to the DGP setting. Barvinok's main insight is that the solutions of a Semidefinite Programming (SDP) relaxation are "not too far" from the feasible set of the Quadratically Constrained Problem (QCP) which gives rise to the SDP relaxation: it suffices to factor the SDP solution and multiply it by a random vector having components sampled from a normal distribution. The very natural idea we propose is to use the naive algorithm first, then use this approximate solution as a starting point for a local Nonlinear Programming (NLP) solver deployed on the original QCP, hoping it will converge to a realization satisfying the feasibility of Eq. ( 1).

Barvinok's naive algorithm is a randomized algorithm based on the concentration of measure phenomenon. It applies to SDP relaxations of systems of quadratic equations, i.e. pure feasibility, equation-only QCPs having vector solutions, or, equivalently, n × 1 matrix solutions. This algorithm cannot be natively applied to the DGP, since a realization x is naturally represented by an n × K matrix x = (x ik | i ≤ n, k ≤ K) the i-th row of which is the position vector of vertex i ∈ V in R K . This matrix may in general have rank greater than one. In [START_REF] Barvinok | Measure concentration in optimization[END_REF], Barvinok gives a proof sketch which only applies to n × 1 matrices. Of course, an n × K matrix can also be represented as an n × 1 matrix where n = nK, but this would entail a QCP (and hence also an SDP relaxation) with nK × nK data matrices, which is practically prohibitive to solve. Moreover, the naive algorithm cannot be applied to the iDGP since it is constrained by quadratic inequalities rather than equations.

We make two contributions in this paper, one theoretical and the other computational. (i) We propose a generalization of Barvinok's result [START_REF] Barvinok | Measure concentration in optimization[END_REF] in two directions: the n × K case, and the case of QCP inequalities of the form d L ≤ x Qx ≤ d U (where x is n × K). (ii) We establish the practical usefulness of Barvinok's result on a set of medium and large-scale DGP instances extracted from the Protein Data Bank (PDB) [START_REF] Berman | The protein data bank[END_REF]. We also remark that our proof is detailed, and fills many gaps in Barvinok's "proof sketch" found in [START_REF] Barvinok | Measure concentration in optimization[END_REF]. The rest of this paper is organized as follows. In Sect. 2 we present Barvinok's naive algorithm. In Sect. 3 we state and prove our generalization of Barvinok's concentration of measure result. In Sect. 5 we discuss our computational results.

The naive algorithm

Consider a system of quadratic equations:

∀i ≤ m x Q i x = a i . (2)
Barvinok's naive algorithm works as follows:

1. solve the SDP relaxation [∀i ≤ mtrace(Q i X) = a i ∧ X 0] of Eq. ( 2), get a solution X * ; 2. factor X * = T T ; 3. sample each component of a vector y ∈ R n from the standard Gaussian distribution N(0, 1); 4. let x = T y.

A concentration of measure argument shows that, with probability at least 0.9,

∀i ≤ m dist(x , A i ) ≤ c X * 2 ln(n),
where

A i = {x ∈ R n | Q i x = a i }, dist(•, •) is the Euclidean distance from a point to a set, X *
2 is the largest eigenvalue of X * , and c is a constant which only depends (linearly) on log n m.

As mentioned in the introduction, applying this algorithm to the DGP requires the representation of n × K realization matrices using vectors in R nK , which, in turn, requires each Q i to be nK × nK, thereby increasing the size of matrices involved in the DGP solution. This provides our first motivation for extending the naive algorithm to work with realizations x represented naturally as n × K matrices. The second motivation is given by the desire of solving the iDGP rather than the (exact) DGP. iDGP solutions are realizations satisfying L 2 i j ≤ x ix j 2 2 ≤ U 2 i j for each edge {i, j} ∈ E, where [L i j , U i j ] is an interval weight on {i, j}.

Extension of the naive algorithm to the n × K case

Our proposed generalization of the naive algorithm, so that it can find n × K realization matrices, is very simple: we replace Step 3 with: 3. sample each component of an n × K matrix y from the standard Gaussian distribution N(0, 1/K) (where 1 K is the variance). We shall see that this has the desired effect. We call the algorithm obtained by this replacement extended naive algorithm.

Concentration of measure and the isoperimetric inequality

In [START_REF] Barvinok | Measure concentration in optimization[END_REF], Barvinok gives the following definition: Definition 1. A measure space (X, µ) has the measure concentration property if there exist positive constants c and C (called concentration constants) such that the following holds: for any Lipschitz function f :

X → R with | f (x) -f (y)| < Mdist(x, y)
for some M > 0 and all x, y ∈ X, and for any ε > 0:

µ x ∈ X | | f (x) -E( f )| > ε < c e -Cε 2 /M 2 (3)
where E( f ) is the expectation of f defined as E( f ) = X f (x)dµ.

A notable example of measure concentration is given by the Gaussian distribution in R n . Let µ be the probability measure on R n with Gaussian density defined by φ(x) = (2π) -n/2 e -x 2 /2 . It is a well known fact that there is measure concentration with c = 2 and C = 2/π 2 .

The well-known isoperimetric inequality states that, if µ has the measure concentration property, then for any closed set A that is large enough, then even a small neighborhood A(ε) of A (defined as A(ε) = {x ∈ X | dist(x, A) ≤ ε} for any ε ≥ 0) can contain almost the whole measure of the space X. This result is stated as [START_REF] Barvinok | Measure concentration in optimization[END_REF]Lemma 1.1.2]; the proof is said to "go along the lines of Sect. 7.6 in [START_REF] Milman | Asymptotic theory of finite dimensional normed spaces[END_REF]". Here we will present a direct proof for it, using elementary arguments. Proposition 2. Assume (X, µ) has measure concentration with constants c, C, and let A ⊆ X be a closed set. For any given p ∈

(0, 1), if ε ≥ √ log c-log µ(A)+ √ log c-log p 1 2 √ C then we have µ(A(ε)) > 1 -p.
Proof. Let us define the function f : X → R + by f (x) = dist(x, A) = min y∈A dist(x, y). It is well-known that f is a 1-Lipschitz function. Then the measure concentration property (3) of (X, µ) implies that, for the given ε > 0:

µ x ∈ X | f (x) -E( f )| > ε/2 < c e -Cε 2 /4 . (4)
Note that, by the choice of ε in the hypothesis, we have c e -Cε 2 /4 < max{µ(A), p}, since it is equivalent to ε > max

√ log c-log µ(A) 1 2 √ C , √ log c-log p 1 2 √ C
. Thus, (4) implies that

µ x ∈ X | f (x) -E( f )| ≤ ε/2 > 1 -µ(A). (5)
On the other hand

µ x ∈ X f (x) = 0 = µ x ∈ X | dist(x, A) = 0 = µ(A). (6)
Therefore, by ( 5) and ( 6), there must exist at least one

x ∈ A such that | f (x) -E( f )| ≤ ε 2 .
But in this case, f (x) = 0, and thus we have

E( f ) ≤ ε 2 . Now we have µ(A c ε ) = µ x ∈ X | f (x) > ε ≤ µ x ∈ X | | f (x) -E( f )| > ε 2 < c e -Cε 2 4
< p, where A c ε is the complement of the set A ε . This concludes the proof.

Analysis of the naive algorithm

In this section we prove a generalization of Barvinok's analysis of his naive algorithm which applies to n × K matrices and systems of quadratic inequalities such as:

∀i ≤ m L i ≤ x Q i x ≤ U i . (7)
Consider the SDP relaxation of Eq. ( 7):

∀i ≤ m L i ≤ trace(Q i X) ≤ U i X 0. ( 8
)
For technical reasons we assume that the number m of constraints of Eq. ( 2) is bounded by a polynomial of nK, i.e. m ≤ (nK) s for some fixed number s > 0. This is a natural assumption in many cases, e.g. we have s = 2 in the DGP. With this assumption, then log m is proportional to log(nK). We remark that the quadratic form in Eq. ( 7) can be written as ∀i ≤ m Q i (x) = x Q i x = 1≤ j, ≤n q i j x j • x , where q i j are real coefficients and z • y stands for the inner product of the vectors z and y. By the naive algorithm, we let x = T y, where T is an n × n factor of the SDP solution and y is an n × K random matrix sampled from the Gaussian distribution N(0, 1 K ). Its j-th row vector is defined by x j = T j y for 1 ≤ j ≤ n, where T j is the j-th row of T . Our objective is to prove that x is an approximate solution for the quadratic system [START_REF] Berger | Reconstructing a three-dimensional model with arbitrary errors[END_REF]. In other words, let A be the set of all feasible solution for (7); we shall then prove that dist(x , A) is small with high probability. The following lemma states that the expectation of quadratic functions Q i (x ) satisfies the system [START_REF] Berger | Reconstructing a three-dimensional model with arbitrary errors[END_REF].

Proposition 3. Let X be the set of all n×K real matrices. Given above notations, we have

X Q i (T y) µ(dy) = trace(Q i X * ) ∈ [L i , U i ]. Proof. For each 1 ≤ i ≤ m we have X Q i (T y) µ(dy) = X 1≤ j, ≤n q i j T j y T y µ(dy) = 1≤ j, ≤n q i j X T j (yy )T µ(dy).
Since entries in the matrix y are sampled independently, for each j, ≤ n and h, k ≤ K, we have:

X y jk y h µ(dy) =        1 K if j = and k = h, 0 otherwise.
For each j, ≤ n, moreover, we note that T j (yy )T is a linear combination of all terms of the form y jk y h (for h, k ≤ K). Therefore, in the integral below, we can simply cancel all the quantities y ik y jh out whenever k h or i j. This yields:

X T j (yy )T µ(dy) = X T j          y 1 2 . . . 0 . . . . . . . . . 0 . . . y n 2          T µ(dy) (9)
Since, for all p ≤ n, y p 2 = K k=1 y 2 pk , Eq. ( 9) can then be rewritten as

n p=1 X T jp T p y p 2 µ(dy) = n p=1 K k=1 X T jp T p y 2 pk µ(dy) n p=1 K k=1 T jp T p X y 2 pk µ(dy) =1/K = 1 K n p=1 K k=1 T jp T p = T j T = X * j .
Therefore, we have

X Q i (T y) µ(dy) = j, ≤n q i j X * j = trace(Q i X * ) ∈ [L i , U i ] (10)
due to the feasibility of X * for the system [START_REF] Barvinok | Measure concentration in optimization[END_REF].

Note that, in Eq. [START_REF] Milman | Asymptotic theory of finite dimensional normed spaces[END_REF], trace(Q i X * ) is also the trace of the quadratic form Q i (T y) with respect to the variable y: indeed, we have

Q i (T y) = (T y) Q i (T y) = y (T Q i T )y, whence trace(Q i (T y)) = trace(T Q i T ) = trace(Q i T T ) = trace(Q i X * ).
In order to proceed, we will need to use the following result, which was proved in Barvinok's paper (a Fact in page 51). Here we extend the result to the n × K dimension case, and provide a proof for completeness.

Lemma 4. Let Q(y) = y Qy be a quadratic form on y ∈ R n×K , A -= {y ∈ R n×K | Q(y) ≤ trace(Q)}, and A + = {y ∈ R n×K | Q(y) ≥ trace(Q)}.
Then for all sufficiently large n we have µ(A -) > 1 3nK and µ(A + ) > 1 3nK .

Proof. Denote nK by N. Without loss of generality, we can assume that trace(Q) ≥ 0. Consider a coordinate system (x 1 , . . . ,

x N ) such that Q is diagonal, i.e. Q(x) = c 1 x 2 1 + . . . + c N x 2
N . Now consider a norm-preserving cyclic transform σ : (x 1 , x 2 , . . . , x N ) → (x 2 , x 3 , . . . , x N , x 1 ). With the convention that

x k = x k-N if k > N, then we have for all vectors x ∈ R N : N k=1 Q(σ k (x)) = N k=1 N i=1 c i x 2 i+k = N i=1 c i N k=1 x 2 i+k = trace(Q) x 2 . Let B N = {x ∈ R N | x ≤ √ N}. Then , for any x ∈ B N , we have N k=1 Q(σ k (x)) ≤ N trace(Q).
It follows that there must be at least one index k ≥ 0 such that Q(σ k (x)) ≤ trace(Q), which means that σ k (x) ∈ A -. The ball B N admits a partition into equivalent classes based on σ (i.e. each class of the form {x, σ(x), . . . , σ N-1 (x)}). Therefore µ(A -) ≥ µ(B N )

N . Similarly, we can also prove that µ(A

+ ) ≥ 1-µ(B N ) N . It is known that lim N→∞ µ(B N ) = 1 2 : indeed, this limit is equal to lim N→∞ x 2 ≤N 1 (2π) N/2 e -x 2 2 dx = 1 2 = lim N→∞ Prob(X 1 + . . . + X N ≤ N),
where X 1 , . . . , X N ∼ χ 2 (1) are i.i.d. random variables. Denote by µ and σ the expectation and variance of X 1 (note that µ = 1 since it is equal to the expectation of the standard normal distribution). The above limit can now be written as

lim N→∞ Prob X 1 +...+X N -µ σ √ N ≤ 0 = Φ(0) = 1 2
, where Φ(x) is the standard normal cdf evaluated at x (due to the Central Limit Theorem). Now, since lim N→∞ µ(B N ) = 1 2 , there is N large enough such that we 1 3 ≤ µ(B N ) ≤ 2 3 . Combining with the earlier inequalities µ(A -) ≥ µ(B N )

N and µ(A

+ ) ≥ 1-µ(B N ) N
, the required inequalities follow.

We now state and prove the extension of Barvinok's concentration of measure result to the n × K case.

Theorem 5. Let X * be a solution of the SDP relaxation (8) and x be the output of the extended naive algorithm. For each i ≤ m we introduce the sets

D i = {x ∈ R n×K | Q i (x) = trace(Q i X * )}.
Then for large enough n, with probability at least 0.9, x satis-

fies dist(x , D i ) ≤ C 0 √ X * 2 ln(nK)
, where C 0 is a constant.

Proof. We define the following sets for each i ≤ m:

A - i = {y ∈ R n×K | Q i (T y) ≤ trace(Q i X * )} A + i = {y ∈ R n×K | Q i (T y) ≥ trace(Q i X * )} A i = {y ∈ R n×K | Q i (T y) = trace(Q i X * )}.
From Lemma 4, we have µ(

A - i ) ≥ 1 3nK and µ(A + i ) ≥ 1 3nK
for all sufficiently large n. Now we can apply Proposition 2 for p = 1 20m to obtain µ(

A - i (ε)) ≥ 1 -1 20m whenever ε ≥ 1 √ C ln c -ln µ(A - i ) + ln c -ln 1 20m
(the measure concentration constants can be taken to be c = 2 and C = 2/π 2 ). Since

ln c -ln µ(A - i ) ≤ ln c -ln( 1 3nK ) = ln c + ln(3nK) = O( ln(nK)) ln c -ln 1 20m = ln 20c + ln(m) ≤ ln 20 c + ln(nK) s = O( ln(nK)),
if we choose some ε = O( √ ln(nK)) with a large enough constant, then we have µ(

A - i (ε)) ≥ 1 -1 20m and similarly µ(A + i (ε)) ≥ 1 -1 20m for i ≤ m.
By the union bound we have

µ(A + i (ε) ∩ A - i (ε)) ≥ 1 -1 10m . We now claim that A + i (ε) ∩ A - i (ε) = A i (ε) for all i ≤ m. Indeed, since A i ⊂ A + i and A i ⊂ A - i , we have A i (ε) ⊆ A + i (ε) ∩ A - i (ε).
For the other inclusion, take any ŷ

∈ A + i (ε) ∩ A - i (ε). It means that there exists y + i ∈ A + i and y - i ∈ A - i such that ŷ-y + i ≤ ε and ŷ -y - i ≤ ε.
By definition it then follows that

Q i (T y - i ) ≤ trace(Q i X * ) and Q i (T y + i ) ≥ trace(Q i X * ).
Since Q i (T y) is a continuous function of y, on the interval defined by two endpoints y + i and y - i , there must exist some y such that Q i (T y ) = trace(Q i X * ), i.e. y ∈ A i . The claim follows since ŷy ≤ max{ ŷy + i , ŷy - i } ≤ ε.

From the claim above, µ(A i (ε)) ≥ 1 -1 10m for i ≤ m. Applying the union bound for all m sets A i (ε) with i ≤ m, we have

µ m i=1 A i (ε) ≥ 1 -m 1 10m = 0.9.
Hence, our randomly generated point y is in A i (ε) for each i ≤ m with probability at least 90%. Therefore

x = T y ∈ T A i (ε) ⊆ D i + ε T ⊆ D i (ε T ) ⊆ D i (ε X * 2 ),
which proves the theorem.

We remark that in Thm. 5 the constant C 0 depends on s (assumed fixed).

The DGP case: removing the dependence on K

Note that the results obtained above are for the general quadratic case. For the DGP case, the problem has a specific block structure which allows us to improve upon them. We will show in this section that the dependence on K in Lemma 4 and Theorem 5 can be removed.

The equality constraint in Eq. ( 1) corresponding to {i, j} ∈ E can be explicitly written as k x 2 ik + k x 2 jk -2 k x ik x jk = d 2 i j . This is a quadratic form on the n × K matrix variable x. If we write it in matrix form x T Q i j x = d 2 i j , the matrix Q i j has the regular structure

Q i j =                     Q n i j 0 0 . . . 0 0 Q n i j 0 . . . 0 0 0 Q n i j . . . 0 . . . . . . . . . . . . . . . 0 0 0 . . . Q n i j                    
with K diagonal blocks, where each for each i, j Q n i j is an n × n matrix with zeros everywhere except for the entries (i, i), (i, j), ( j, i), ( j, j) set to 1 -1 -1

1 . Now we can rewrite Lemma 4 for this specific block matrix, in which the dependence on K is removed:

Lemma 6. Let Q(y) = y Qy be a quadratic form on y ∈ R n×K , where Q consists of K diagonal identical symmetric blocks Q n ∈ R n×n , namely Q = diag(Q n , . . . , Q n K ). Let A -= {y ∈ R n×K | Q(y) ≤ trace(Q)} and A + = {y ∈ R n×K | Q(y) ≥ trace(Q)}.
Then for all sufficiently large n we have µ(A -) > 1 3n and µ(A + ) > 1 3n . Proof. Without loss of generality, we can assume that trace(Q) ≥ 0. By changing of the coordinate system if necessary, we can assume that each block Q n is diagonal, i.e. Q n = diag(λ 1 , . . . , λ n ). Now consider an orthogonal transform σ that sends(x 1 , x 2 , . . . , x nK ) to (x 2 , . . . , x n , x 1 ; x n+2 , . . . , x 2n , x n+1 ; x 2n+2 , . . . , x 3n , x 2n+1 ; . . . ; x nK , x nK-n+1 ).

Then σ generates a cyclic group of order n. Averaging over n shifts, we get

n k=1 Q(σ k (x)) = 1 K trace(Q) x 2 . Let N = nK and let B N = {x ∈ R N | x ≤ √ N}.
Then , for any x ∈ B N , we have

n k=1 Q(σ k (x)) ≤ N K trace(Q) = n trace(Q).
It follows that there must be at least one index k ≥ 0 such that Q(σ k (x)) ≤ trace(Q), which means that σ k (x) ∈ A -. The ball B N admits a partition into equivalent classes based on σ (i.e. each class of the form {x, σ(x), . . . , σ n-1 (x)}). Therefore µ(A -) ≥ µ(B N )

n . Similarly, we can also prove that µ(A + ) ≥ 1-µ(B N ) n . From Lemma 4, we have proved that lim N→∞ µ(B N ) = 1 2 . Therefore, when N = nK is sufficiently large, then we have

µ{y ∈ R n×K | Q(y) ≤ trace(Q)} > 1 3n and µ{y ∈ R n×K | Q(y) ≥ trace(Q)} > 1 3n as claimed.
Now we can obtain an improvement for the result in Theorem 5. We will state it without proof here because the proof is exactly the same except for the application of Lemma 6 instead of Lemma 4.

Theorem 7. Let X * be a solution of the SDP relaxation (8) for the DGP (1) and x be the output of the extended naive algorithm. For each i ≤ m we introduce the set

D i = {x ∈ R n×K | Q i (x) = trace(Q i X * )}.
Then for large enough n, with probability at least 0.9, x satisfies dist(x , D i ) ≤ C 0 √ X * 2 ln(n), where C 0 is a constant.

Note that the logarithmic term under the square root in the RHS was improved from nK to n.

A different approach to the iDGP case

In this section, we argue that with appropriate assumptions on the solution X * and the interval bounds [L i , U i ], we can obtain a different bound, which depends on ln m instead of ln nK, which might . We will use the following result, which is proved in [START_REF] Barvinok | A Course in Convexity[END_REF].

Proposition 8. Let Q : R N → R be a positive semidefinite quadratic form, i.e. Q(x) = x Qx for a positive semidefinite matrix Q. Let Q = µ 2 1 + . . . + µ 2 N be the Schatten 2-norm of Q, in which µ 1 , . . . , µ N are the eigenvalues of Q. Then 1. For any τ ≥ 0 we have

µ x ∈ R N | Q(x) < trace(Q) -τ ≤ exp - τ 2 4 Q 2 . ( 11
)
2. For any τ ≥ 0 such that τµ i ≤ Q 2 for all i ≤ N, we have

µ x ∈ R N | Q(x) > trace(Q) + τ ≤ exp - τ 2 8 Q 2 . ( 12
)
Now we apply this proposition to the naive algorithm. Since X * is a feasible solution of the SDP relaxation of the iDGP, then by [START_REF] Barvinok | Measure concentration in optimization[END_REF] we have L i ≤ trace(Q i X * ) ≤ U i for all i ≤ m. We define some parameters λ i such that λ i ≤

min U i -trace(Q i X * ) trace(Q i X * ) , trace(Q i X * )-L i trace(Q i X * )
for i ≤ m. By this definition,

we have L i ≤ trace(Q i X * ) -λ i trace(Q i X * )
. Therefore, by applying [START_REF] Barvinok | A Course in Convexity[END_REF] for τ = λ i trace(Q i X * ), we have

µ y ∈ R n×K | Q i (T y) < L i ≤ µ y ∈ R n×K | Q i (T y) < trace(Q i X * ) -λ i trace(Q i X * ) ≤ exp - λ 2 i trace(Q i X * ) 2 4 Q i X * 2 ,
which will be smaller than or equal to 1 3 if we have λ i ≥ √ 4 ln 3

Q i X * trace(Q i X * )
. Next, we do the same with the inequality in [START_REF] Barvinok | Math 710: Measure Concentration[END_REF]. By the definition of λ i , we have U i ≥ trace(Q i X * )+λ i trace(Q i X * ), therefore by applying [START_REF] Barvinok | Math 710: Measure Concentration[END_REF] for τ = λ i trace(Q i X * ), we have

µ y ∈ R n×K | Q i (T y) > U i ≤ µ y ∈ R n×K | Q i (T y) > trace(Q i X * ) + λ i trace(Q i X * ) ≤ exp - λ 2 i trace(Q i X * ) 2 8 Q i X * 2 ,
which will be smaller than or equal to 1

√ 3 if we have λ i ≥ √ 4 ln 3 Q i X *
trace(Q i X * ) (similarly to the previous case). Note that, to ensure the assumption in ( 12) is satisfied, we have to further assume that λ i trace(Q i X * )σ max (Q i X * ) ≤ Q i X * 2 as a sufficient condition, where σ max (A) stands for the maximum singular value of A. In conclusion, the above conditions are equivalent to 2

√ 3 ≤ min U i -trace(Q i X * ) Q i X * , trace(Q i X * ) -L i Q i X * , Q i X * 2 trace(Q i X * )σ max (Q i X * ) . (13)
If this condition holds, then by using the union bound we have

µ y ∈ R n×K | L i ≤ Q i (T y) ≤ U i ≥ 1 - 1 3 - 1 √ 3 = 2 - √ 3 3 . (14)
Proposition 9. Let x be the output of the extended naive algorithm and assume that Eq. (13) holds. For each

1 ≤ i ≤ m, denote D i = {x ∈ R n×K | L i ≤ Q i (x) ≤ U i }.
Then with the probability at least 90%, x satisfies the dist(x , D i ) ≤ C 0 √ ln m X * 2 , where C 0 is a positive universal constant.

Proof. For each 1 ≤ i ≤ m, we define a closed set

A i = {y ∈ R n×K | L i ≤ Q i (T y) ≤ U i }. From the estimation in (14), we have µ(A i ) ≥ 2- √ 3 3
for all i ≤ m. Then for any ε and 0

< p < 1 such that log c -log µ(A i ) + log c -log p 1 2 √ C ≤ log c -log 2- √ 3 3 + log c -log p 1 2 √ C = O( log( 1 p ) ≤ ε, we have µ(A i ( )) > 1 -p (by Proposition 2). Take p = 1 10m
and apply the union bound, then we have µ m i=1 A i (ε) > 1m p = 0.9. Hence, our randomly generated point y is in A i (ε) for each i ≤ m with probability at least 90%; thus

x = T y ∈ T A i (ε) ⊆ D i + ε T ⊆ D i (ε T ) ⊆ D i (ε √ X * 2 ). Note that we can choose ε = O( log 1 p ) = O( log m). Therefore dist(x , D i ) ≤ O( √ ln m X *
), which proves the theorem.

We remark that in Prop. 9 the constant C 0 does not depend on s as in Thm. 5.

Computational results

Testing probabilistic algorithms based on the phenomenon of concentration of measure involves the difficulty of having to determine the minimum instance size starting from which the algorithm will start producing useful outputs. Most proofs involve unknown constants (e.g. c, C in Eq. ( 3)). Valid values for these constants can be determined (see e.g. [START_REF] Barvinok | Math 710: Measure Concentration[END_REF]), but such values often make empirical applications difficult or impossible. This is why computational tests are sometimes employed in the attempt to determinate practically useful values [START_REF] Venkatasubramanian | The Johnson-Lindenstrauss transform: An empirical study[END_REF].

The objective of our computational testbed is to ascertain whether Barvinok's naive algorithm is beneficial or not. We therefore considered two heuristics, the only difference between them being the activation/deactivation of this algorithm as a sub-step. We call bvk the version of this algorithm with Barvinok's naive algorithm, and sdpnlp the version without. Both heuristics have the following structure:

1. compute the solution X * of the SDP relaxation of the DGP or iDGP instance; 2. map the n × n matrix X * into an n × K matrix x ; 3. let x * be the n × K realization obtained by applying a local NLP solver to the starting point x .

The bvk algorithm is obtained by replacing Step 2 with Barvinok's naive algorithm, while sdpnlp replaces Step 2 with Principal Component Analysis (PCA):

• factor X * as V ΛV, where Λ = diag(λ 1 , . . . , λ n ) are the eigenvalues, and V is the matrix of (column) eigenvectors;

• since X * 0, λ i ≥ 0 for all i, which means that √ Λ is real; • let Λ = diag(λ 1 , . . . , λ K , 0, . . . , 0), and let x = V √ Λ .

It is quite usual to compare DGP algorithms on solution quality and CPU time [START_REF] Liberti | Euclidean distance geometry and applications[END_REF][START_REF] D'ambrosio | New error measures and methods for realizing protein graphs from distance data[END_REF]. The two most established error measures are the mean distance error (MDE), defined as 1 |E| {i, j}∈E x ix j 2d i j , and the largest distance error (LDE), defined as max {i, j}∈E x ix j 2d i j . The equivalent measures for the iDGP replace x ix j 2d i j in the above formulae with max(0, L i jx ix j 2 ) + max(0, x ix j 2 -U i j ) [START_REF] Liberti | Molecular distance geometry methods: from continuous to discrete[END_REF][START_REF] D'ambrosio | New error measures and methods for realizing protein graphs from distance data[END_REF]. Both errors only take into account the discrepancy of the computed solution with the input data. While a low MDE can still yield a significantly different realization (if only locally), a low LDE usually gives a higher confidence of the reconstructed realization being close to the input data.

As a local NLP solver, we chose the "root" option out of the solvers in the scipy library [START_REF] Jones | SciPy: Open source scientific tools for Python[END_REF]. Limited to the DGP, we found this solver to be superior (both CPU-time and qualitywise) to the other solvers in the library, except on really tiny instances where it appears to fail often. In all cases, the total CPU time is dominated by the SDP solution and by the local NLP solution; this is natural as optimization solvers are largescale software packages implementing complicated algorithms, while both Barvinok's naive algorithm and PCA are quite simple.

Note that Step 2 is stochastic in bvk and deterministic in sdpnlp. It therefore makes sense to consider the best outcome (w.r.t. LDE) out of a certain number of iterations of Barvinok's algorithm, while PCA need only be run once. While this might appear to penalize the former as regards the CPU time, note that its worst-case complexity is dominated by the matrix product (resulting in O(Kn 2 )), while PCA requires a matrix factorization, which is O(n 3 ) in general. We empirically set the maximum number of iterations of Barvinok's algorithm to 500.

The computational results, obtained using a Python 2.7 code on an Intel i5 CPU with 8GB RAM running Windows 10, are presented for the DGP in Table 1, and for the iDGP in Table 2. For the DGP, it is quite obvious from Table 1 better solutions than sdpnlp. It is easy to see that the latter is faster than the former; it is, however, also easy to notice that bvk seems to "catch up" for instances of increasing size. the experiments on the iDGP presented in Table 2 are not as clear. We notice that three instances cause a memory overflow error on both heuristics, which is a consequence of the NLP solver. Again, the sdpnlp heuristic is superior to bvk as concerns CPU time; and the ratio between the two values again becomes smaller as the size increases. In terms of solution quality, sdpnlp is superior to bvk for all but the three largest instances, and slightly inferior for the two largest. While the discrepancies in LDE are too small to derive any meaningful conclusion, we remark that concentration of measure phenomena always "kick in" from some instance size onwards, but it is difficult to establish the threshold size. The larger instances we tested all yielded memory overflows, due to both SDP and NLP solvers. An anonymous referee suggested a direction for future work which we find very interesting. Our idea for deploying a local descent method from the starting point x found by Barvinok's naive algorithm was given by the geometric intuition that, if the distance from x to the feasible set of the given QCP is bounded with high probability, then a local descent from x might be all that is needed to find a feasible point of the QCP. The referee asks the following question: are there natural conditions under which a gradient descent from x is guaranteed (perhaps with some probability) to perform nicely?

While researchers in mathematics and theoretical computer science (TCS) produce many algorithms with interesting theoretical properties, their adaptation to applications and consequent testing is scattered at best. Some TCS researchers have been heard, during talks at conferences, making various statements to the effect that their algorithms are not for practical use, and that their value is wholly theoretical. We beg to differ: whether these algorithms have a practical value or not is simply not known (for lack of testing). While we realize that no amount of testing can ever prove anything, an indication of usefulness might remarkably increase the impact of a "wholly theoretical algorithm" (witness e.g. the Goemans-Williamson randomized approximation algorithm for max cut [START_REF] Goemans | Improved approximation algorithms for maximum cut and satisfiability problems using semidefinite programming[END_REF]). We hope this paper will contribute to increase the impact of the wonderfully elegant, simple and yet effective "naive algorithm" that A. Barvinok invented twenty years ago.

Table 1 :

 1 that bvk yields Computational comparison of Barvinok's algorithm and PCA on DGP instances.

			MDE		LDE	CPU
	instance	bvk	sdpnlp	bvk	sdpnlp	bvk	sdpnlp
	names	0.00	0.11	0.07	1.00	39.33	22.44
	pept	0.00	0.10	0.03	1.81	83.91	56.65
	C0020pdb	0.00	0.12	0.01	2.72	76.73	49.39
	1guu-1	0.03	0.00	0.26	0.08	370.73	322.66
	1guu-4000	0.03	0.12	0.73	1.15	415.66	397.87
	1guu	0.02	0.01	0.29	0.33	305.54	267.52
	res 5000	0.00	0.15	0.00	2.24	84.19	53.36
	res 2000	0.00	0.07	0.00	1.46	85.26	53.39
	res 0	0.00	0.00	0.00	0.01	93.08	62.64
	res 3000	0.00	0.01	0.00	1.08	88.51	53.43
	res 1000	0.00	0.10	0.00	3.05	87.88	52.98
	res 2kxa	0.00	0.15	0.00	2.92	764.34	713.35
	C0030pkl	0.00	0.11	0.07	2.19	1178.73	1024.86

Table 2 :

 2 This is consistent with the fact that we run a fixed number of iterations of Barvinok's naive algorithm (which is O(Kn 2 )) and only one iteration of PCA (which is O(n 3 )). The outcome of Computational comparison of Barvinok's algorithm and PCA on iDGP instances.

			MDE		LDE	CPU
	instance	bvk	sdpnlp	bvk	sdpnlp	bvk	sdpnlp
	names	0.04	0.00	2.11	0.07	53.91	36.86
	pept	0.01	0.00	0.46	0.40	133.28	99.60
	C0020pdb	0.02	0.00	1.64	0.42	112.38	79.22
	1guu-1	0.02	0.01	1.09	0.58	500.64	440.50
	1guu-4000	0.03	0.02	1.49	1.49	522.19	461.53
	1guu			memory overflow
	res 5000	0.01	0.00	0.69	0.08	30764.21	30465.16
	res 2000	0.01	0.00	1.78	0.10	33017.88	32713.91
	res 0	0.00	0.00	0.11	0.11	22897.14	22619.79
	res 3000	0.00	0.00	0.05	0.08	26095.91	25846.81
	res 1000	0.00	0.00	0.05	0.07	27790.87	27542.96
	res 2kxa			memory overflow
	C0030pkl			memory overflow

Preprint submitted to Operations Research Letters July 6, 2018

$ We are grateful to A. Barvinok for helpful suggestions and discussions, as well as to an anonymous referee for helping us improve the paper. The first author (LL) has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement n. 764759 (ITN "MINOA").