Ky Vu

Pierre-Louis Poirion

Leo Liberti

Gaussian random projections for Euclidean membership problems

Keywords: Johnson-Lindenstrauss lemma, Machine Learning, Euclidean Distance Geometry, clustering

published or not. The documents may come

Gaussian random projections for Euclidean membership problems

Ky Vu, Pierre-Louis Poirion, Leo

Introduction

Random projections are useful dimension reduction techniques which are widely used in many areas such as Machine Learning [START_REF] Blum | Random projection, margins, kernels, and feature-selection[END_REF][START_REF] Fradkin | Experiments with random projections for machine learning[END_REF], Computer Science [START_REF] Indyk | Algorithmic applications of low-distortion geometric embeddings[END_REF][START_REF] Matousek | Lectures on Discrete Geometry[END_REF], Numerical Linear Algebra [START_REF] Mahoney | Randomized algorithms for matrices and data[END_REF][START_REF] Woodruff | Sketching as a tool for numerical linear algebra[END_REF] and so on. In standard algorithmic settings, assume we have an algorithm A acting on a data set X consisting of n vectors in R m , where m is large, and assume that the complexity of A depends on m and n in a way that makes it impossible to run A sufficiently fast. A random projection exploits the statistical properties of some random distribution to construct a mapping which embeds X into a lower dimensional space R k (for some appropriately chosen k) while preserving distances, angles, or other quantities used by A. Their simplicity notwithstanding, the performance of random projections is comparable to more complicated methods such as SVD and PCA [START_REF] Dasgupta | Experiments with random projection[END_REF]. One striking example of random projections is the famous Johnson-Lindenstrauss lemma [START_REF] Johnson | Extensions of Lipschitz mappings into a Hilbert space[END_REF]:

Email addresses: vukhacky@gmail.com (Ky Vu), pierre.louis.poirion@huawei.com (Pierre-Louis Poirion), liberti@lix.polytechnique.fr (Leo Liberti)

1.1 Theorem (Johnson-Lindenstrauss Lemma (JLL)) Let X be a set of n points in R m and ε > 0. Then there is a map F : R m → R d where d is O( log n ε 2 ), such that for any x, y ∈ X, we have

(1 -ε) x -y 2 2 ≤ F (x) -F (y) 2 2 ≤ (1 + ε) x -y 2 2 . (1)
Intuitively, this lemma claims that X can be projected into a much lower dimensional space whilst keeping Euclidean distances approximately the same. The main idea to prove Theorem 1.1 is to construct a random linear mapping T (called JL random mapping onwards), sampled from certain distribution families, so that for each x ∈ R m , the event that

(1 -ε) x 2 2 ≤ T (x) 2 2 ≤ (1 + ε) x 2 2 (2)
occurs with high probability. By Eq. ( 2) and the union bound, it is possible to show the existence of a map F with the stated properties [START_REF] Achlioptas | Database-friendly random projections: Johnson-Lindenstrauss with binary coins[END_REF][START_REF] Dasgupta | An elementary proof of a theorem of Johnson and Lindenstrauss[END_REF]. The JLL is a consequence of a general property of sub-gaussian random mappings T = 1 √ d U [START_REF] Matoušek | On variants of the Johnson-Lindenstrauss lemma[END_REF]. Some of the most popular choices for U are: Choices of random projection 1. U = √ mP where P is the orthogonal projections on a random k-dimensional linear subspace of R m [START_REF] Johnson | Extensions of Lipschitz mappings into a Hilbert space[END_REF]; 2. random d × m matrices with each entry independently drawn from the standard normal distribution N (0, 1) [13]; 3. random d × m matrices with each entry independently taking values +1 and -1, each with probability 1 2 [1]; 4. random d × m matrices with entries independently taking values +1, 0, -1, respectively with probability 1 6 , 2 3 , 1 6 [START_REF] Achlioptas | Database-friendly random projections: Johnson-Lindenstrauss with binary coins[END_REF]. In this paper we employ random projections to study the following metaproblem (parametrized on S):

Euclidean Set Membership Problem (ESMP). Given b ∈ R m and S ⊆ R m , decide whether b ∈ S.
This problem is important and extensively studied in optimization as feasibility problems. For example, any linear and integer programming problem can be transformed easily (by bisection arguments) into feasibility problems:

Integer and Linear Feasibility Problems. Decide whether the set {x ∈ X | Ax = b} is empty or not, where X is either R n + or Z n + . This is exactly the membership problem where S is the restricted linear span {Ax | x ∈ X}.

In this paper, we consider the general case where the data set S has no specific structure. We will use a random projection T = 1 √ d U (we assume in the whole paper than d ≥ 3), where U is a Gaussian random projection, i.e., a random d × m matrix with each entry independently drawn from the standard normal distribution N (0, 1), in our arguments to embed both b and S to a lower dimensional space, and solve its associated projected version:

Projected ESMP (PESMP). Given b ∈ R m and S ⊆ R m , and let T : R m → R d be a random projection. Decide whether T (b) ∈ T (S).

Our objective is to investigate the relationships between ESMP and PESMP. It is easy to argue that b ∈ S implies T (b) ∈ T (S). In this paper, we are only interested in the problematic case when b / ∈ S. Since T is a random projection, we want to estimate Prob(T (b) / ∈ T (S)), given that b / ∈ S. Notice that the probability does not change if we multiply T by some constant. Hence we could simply consider a Gaussian random matrix U instead of T = 1 √ d U , however for the sake of clarity in the proof, we will work with normalized Gaussian random matrix, i.e., T , in the whole paper.

Our contributions

We concentrate on the use of Gaussian random projections as a tool for dimension reduction in the general Euclidean membership problem. We consider two special cases: S is at most countable (i.e. finite or countable) and S is of low doubling dimension.

In the first case, using a straightforward argument, we prove that these two problems (ESMP and PESMP) are equivalent almost surely regardless of the projected dimension. However, this result is only of theoretical interest, possibly due to round-off errors in floating point operations which make its practical application difficult. We address this issue by introducing a threshold τ > 0 and study the corresponding "threshold" problem:

Threshold ESMP (TESMP): Given b, S, T as above. Let τ > 0. Decide whether T (b) -T (S) ≥ τ .

In the case when S may also be uncountable, we prove that ESMP and PESMP are also equivalent if the projected dimension d is proportional to some intrinsic dimension of the set S. In particular, we employ the notion of doubling dimension (defined later) to prove that, if b / ∈ S, then T (b) / ∈ T (S) almost surely as long as the projected dimension d ≥ C ddim(S), where ddim(S) is the doubling dimension of S and C is some universal constant. We extend this result to the threshold case, and obtain a more useful bound for d. It turns out that, as a consequence of that result, we are able to improve a bound of Indyk-Naor on the Nearest Neigbour preserving embeddings by a factor of log(1/δ) ε .

Applicability

Although this paper only makes theoretical contributions, the applicability of random projections is beyond doubt. They are routinely used in large-scale problems arising in clustering, for example for images and text. We are ourselves pursuing random projections in the context of linear and integer programming [START_REF] Vu | Random projections for linear programming[END_REF].

Finite and countable sets

In this section, we assume that S is either finite or countable. Let U ∈ R d×m be a random matrix drawn from Gaussian distribution, i.e. each entry of U is independently sampled from N (0, 1). It is well known that, for an arbitrary unit vector a ∈ S m-1 , U a 2 is a random variable with a Chi-squared distribution χ 2 d with d degrees of freedom [START_REF] Mood | Introduction to the Theory of Statistics[END_REF]. Its corresponding density function is

F χ 2 d (δ) ≤ (ze 1-z ) d/2 < (ze) d/2 = eδ d d/2 . (3)
Thus, we have

Prob( U a ≤ δ) = F χ 2 d (δ 2 ) < ( 3 d δ 2 ) d/2 (4)
or, more simply, Prob( U a ≤ δ) < δ d when d ≥ 3.

In the case of a normalized Gaussian random projection

T = 1 √ d U , we have Prob( T a ≤ δ) = Prob( U a ≤ √ dδ) < (3δ 2 ) d 2 (5)
Using this estimation, we immediately obtain the following result.

Proposition

Given b ∈ R m and S ⊆ R m , at most countable, such that b / ∈ S. Then, for any d ≥ 1 and for a normalized Gaussian random projection

T : R m → R d (T = 1 √ d U ) , we have T (b) / ∈ T (S) almost surely, i.e. Prob T (b) / ∈ T (S) = 1.
Proof. First, note that for any a = 0, T a = 0 holds almost certainly. Indeed, without loss of generality we can assume that a = 1. Then for any 0 < δ < 1:

Prob T (a) = 0 ≤ Prob T a ≤ δ (5)
< (3δ 2 ) d/2 → 0 as δ → 0. ( 6
)
It means that for any y = b, then T (y) = T (b) almost surely. Since the event T (b) / ∈ T (S) can be written as the intersection of at most countably many "almost sure" events T (b) = T (y) (for y ∈ S), it follows that Prob T (b) / ∈ T (S) = 1, as claimed.

2

Proposition 2.1 is simple, but it looks interesting because it suggests that we only need to project the data points to a line (i.e. d = 1) and study an equivalent membership problem on a line. This idea, stated somewhat differently in terms of random aggregations of linear constraints, appears to be part of the folklore of linear and integer programming. It turns out that this result remains true for a large class of random projections.

Proposition

Let ν be a probability distribution on R m with bounded Lebesgue density f . Let S ⊆ R m be an at most countable set such that 0 / ∈ S. Then, for a random projection T : R m → R 1 sampled from ν, we have 0 / ∈ T (S) almost surely, i.e. Prob 0 / ∈ T (S) = 1.

Proof. For any 0 = s ∈ S, consider the set

E s = {T : R m → R 1 | T (s) = 0}. If we regard each T : R m → R 1 as a vector t ∈ R m , then E s is a hyperplane {t ∈ R m | s • t = 0} and we have Prob(T (s) = 0) = ν(E s ) = Es f dµ ≤ f ∞ Es dµ = 0
where µ denotes the Lebesgue measure on R m . The proof then follows by the countability of S, similarly to Proposition 2.1. 2

Appealing as the idea may be, it does not seem to work in practice. We tested it by considering the ESMP given by the Integer Feasibility Problem (IPF) defined on the set {x ∈

Z n + ∩[L, U ] | Ax = b}, where L, U ∈ Q n and, for all i = 1, • • • , n, L i ≤ U i . Numerical experiments indicate that the corresponding PESMP {x ∈ Z n + ∩ [L, U ] | T (A)x = T (b)}, with
T consisting of a one-row Gaussian projection matrix, is always feasible despite the infeasibility of the original IPF. Since Prop. 2.1 assumes that the components of T are real numbers, the reason behind this failure is possibly due to the round-off error associated to the floating point representation used in computers. Specifically, when T (A)x is too close (but not equal) to T (b), floating point operations will consider them as a single point. In order to address this issue, we force the projected problems to obey stricter requirements. In particular, instead of only requiring that T (b) / ∈ T (S), we ensure that

dist(T (b), T (S)) = min x∈S T (b) -T (x) > τ, (7)
where dist denotes the Euclidean distance, and τ > 0 is a (small) given constant.

With this restriction, we obtain the following result. Proof. For any x ∈ S, by the linearity of T , we have:

Proposition

Given τ > 0, 0 < δ < 1 and b / ∈ S ⊆ R m ,
Prob T (b -x) ≤ τ = Prob T b -x b -x ≤ τ b -x ≤ Prob T b -x b -x ≤ τ R < ( √ 3τ ) d R d ,
due to [START_REF] Dasgupta | Experiments with random projection[END_REF]. Therefore, by the union bound,

Prob min x∈S T (b) -T (x) > τ = 1 -Prob min x∈S T (b) -T (x) ≤ τ ≥ 1 - x∈S Prob T (b) -T (x) ≤ τ > 1 -|S| ( √ 3τ ) d R d .
The RHS is greater than or equal to 1 -δ if and only if

R d ( √ 3τ ) d ≥ |S| δ , which is equivalent to d ≥ log( |S| δ ) log( R √ 3τ ) , as claimed.

2

Note that R is often unknown and can be arbitrarily small. However, if both b, S are integral and τ < 1, then R ≥ 1 and we can select d >

log |S| δ log 1 √ 3τ
in the above proposition.

Sets with low doubling dimensions

In many real-world applications, the data set S is not finite or countable, but it lies in some intrinsically low-dimensional space. There are many examples of such sets, including human motion records, speed signals, image and text data and more [START_REF] Dasgupta | Random projection trees and low dimensional manifolds[END_REF][START_REF] Bingham | Random projection in dimensionality reduction: Applications to image and text data[END_REF]. Random projections can provide a tool to extract the full information of the set, in spite of the (high) dimension of the ambient space that it is embedded in.

In this section, we will use doubling dimension as a measure for the intrinsic dimension of a set. Let denote by B(x, r) = {y ∈ S : y -x ≤ r}, i.e. the closed ball centered at x with radius r > 0 (w.r.t S). We use the following definition:

Definition

The doubling constant of a set S ⊆ R m is the smallest number λ S such that any closed ball in S can be covered by at most λ S closed balls of half the radius. A set S ⊆ R m is called a doubling set if it has a finite doubling constant. The number log 2 (λ S ) is then called the doubling dimension of S and denoted by ddim(S).

One popular example of doubling spaces is a Euclidean space R m . In this case, it is well-known that its doubling dimension is O(m) [START_REF] Verger-Gaugry | Covering a ball with smaller equal balls in R n[END_REF][START_REF] Gupta | Bounded geometries, fractals, and low-distortion embeddings[END_REF]. However, there are cases where the set S ⊆ R m is of much lower doubling dimension. It is also easy to see that ddim(S) does not depend on the dimension m. The doubling dimension is therefore a powerful tool for reducing dimensions in several classes of problems such as nearest neighbor [START_REF] Krauthgamer | Navigating nets: Simple algorithms for proximity search[END_REF][START_REF] Indyk | Nearest neighbor preserving embeddings[END_REF], low distortion embeddings [START_REF] Agarwal | Embeddings on surfaces, curves, and moving points in Euclidean space[END_REF], and clustering [START_REF] Magen | Dimensionality reductions in 2 that preserve volumes and distance to affine spaces[END_REF].

In this section, we assume that ddim(S) is relatively small compared to m. Note that, computing the doubling dimension of S is generally NP-hard [START_REF] Gottlieb | Proximity algorithms for nearly doubling spaces[END_REF], although it can be approximated within a constant factor [START_REF] Gottlieb | Proximity algorithms for nearly doubling spaces[END_REF]. For simplicity, we assume that λ S is a power of 2, i.e. the doubling dimension of S is a positive integer number.

We shall make use of the two following lemmas.

Lemma

For any b ∈ S and ε, r > 0, there is a set S 0 ⊆ S of size at most λ

log 2 ( r ε ) S such that B(b, r) ⊆ s∈S0 B(s, ε). (8)
Proof. By definition of the doubling dimension, B(b, r) is covered by at most λ S closed balls of radius r 2 . Each of these balls in turn is covered by λ S balls of radius r 4 , and so on: iteratively, for each

k ≥ 1, B(b, r) is covered by λ k S balls of radius r 2 k . If we select k = log 2 ( r ε ) then k ≥ log 2 ( r ε ), i.e. r 2 k ≤ ε. This means B(b, r) is covered by λ log 2 ( r ε ) S balls of radius ε. 2

Lemma

Let S ⊆ {s ∈ R m | s ≤ 1} be a subset of the m-dimensional Euclidean unit ball. Let T : R m → R d be a normalized Gaussian random projection. Then there exist universal constants c, C > 0 such that for d ≥ C log λ S and δ ≥ 6, the following holds:

Prob(∃s ∈ S s.t. T s > δ) < e -cdδ 2 . ( 9
)
This lemma is proved in [START_REF] Indyk | Nearest neighbor preserving embeddings[END_REF] using concentration estimations for sum of squared gaussian variables (Chi-squared distribution). In particular, we recall an important inequality, proved in [START_REF] Indyk | Nearest neighbor preserving embeddings[END_REF] : for all a ∈ R m of unit norm, all δ > 0 and a mapping T as above,

Prob (| T a -1| > δ) ≤ e -dδ 2 /8 and (10)
For the sake of completeness, we will present the original proof [START_REF] Indyk | Nearest neighbor preserving embeddings[END_REF] of the above Lemma. Here, we use an additional requirement that δ ≥ 6 instead of δ > 1, however the main argument is unchanged. Proof. Choose b = 0 ∈ S, r = 1 and ε k = 1 2 k in Lemma 3.2. By earlier convention that B(x, r) = {y ∈ S : y -x ≤ r}, obviously we have S ⊆ B(0, 1). Then there is a set S k ⊆ S of size at most λ k S such that

S ⊆ s∈S k B(s, 1 2 k ). (11)
Therefore, for any x ∈ S, there is a sequence (x n ) n∈N that converges to x, with ∀n ∈ N, x n ∈ S n and x n -x n+1 ≤ 1 2 n (this holds because each ball B(x n , 1 2 n ) is itself covered by balls of radius 1 2 n+1 with centers in S n+1 ). We claim that, if T x ≥ δ, then there must be some k ∈ N such that

T (x k -x k+1 ) ≥ δ 3 ( 3 2 ) -k
Indeed, if no such k exists, then

δ ≤ T x ≤ ∞ k=0 T (x k -x k+1 ) < δ 3 ∞ k=0 ( 3 2
) -k = δ, a contradiction. Now, if we want to neglect x, we can treat x k and x k+1 (found above) as two points u, v, in which u ∈ S k and v ∈ B(u, 1 2 k ) ∩ S k+1 . Then we have

T u -T v ≥ δ 3 ( 3 2 ) -k ≥ u -v 2 -k δ 3 ( 3 2 ) -k = δ 3 ( 4 3 ) k u -v .
Therefore, we have

Prob ∃x ∈ S s.t T (x) > δ ≤ Prob ∃x ∈ S, ∃k ≥ 0 s.t T (x k -x k+1 ) > δ 3 ( 3 2 ) -k ≤ ∞ k=0 Prob ∃u ∈ S k , ∃v ∈ B(u, 1 2 k ) ∩ S k+1 s.t T u -T v > δ 3 ( 4 3 ) k u -v ≤ ∞ k=0 Prob ∃u ∈ S k , ∃v ∈ B(u, 1 2 k ) ∩ S k+1 s.t T ( u -v u -v ) > δ 3 ( 4 3 ) k ≤ ∞ k=0 λ 2k+1 S Prob T z > δ 3 ( 4 3
) k for any unit z (by the union bound).

Since δ ≥ 6, for any k ≥ 0 we have δ 3 ( 4 3 ) k ≥ 1 + δ 6 ( 4 3 ) k . Therefore, from the inequality [START_REF] Gottlieb | Proximity algorithms for nearly doubling spaces[END_REF], the above expression is less than or equal to

∞ k=0 λ 2k+1 S Prob | T z -1| ≥ δ 6 ( 4 3 ) k ≤ ∞ k=0 λ 2k+1 S e -d δ 2 8 * 6 2 ( 4 3 ) 2k ≤ e -cdδ 2 (12)
as long as d ≥ C log(λ S ) for some universal constants c, C. 2

We now state one of the main results in this paper. It says that we can maintain the equivalence between ESMP and its projected version by using Gaussian random projections with d proportional to the doubling dimension of S.

Theorem Given b /

∈ S ⊆ R m where S is a closed doubling set. Let T : R m → R d be a normalized Gaussian random projection. Then

Prob(T (b) / ∈ T (S)) = 1 (13) if d ≥ C log 2 (λ S )
, for some universal constant C.

Proof. Let d ≥ C log 2 (λ S ) for some universal constant C (large). The proof is divided in three steps: first we build a sequence of annuli, (X k ) k∈N that cover S; then on each annulus, X k , we estimate the probability that there exists an element x ∈ X k such that T (x) = T (b). Finally, we apply the union bound on the whole set S and prove that the probability that there exists x ∈ S such that T (x) = T (b) can be as small as needed. Assume that R > 0 is the distance between b and the set S. Let ε i , ∆ i , i ∈ N and R = r 0 < r 1 < r 2 < . . . be positive scalars (their values will be defined later). For each k ∈ N * we define an annulus

X k = S ∩ B(b, r k ) B(b, r k-1 ). Since X k ⊆ S ∩ B(b, r k ), by Lemma 3.2 we can find a point set S k ⊆ S of size |S k | ≤ λ log 2 ( r k ε k ) S such that X k ⊆ s∈S k B(s, ε k ).
Hence, for any x ∈ X k , there is s ∈ S k such that x -s < ε k . Moreover, by the triangle inequality, any such

s satisfies r k-1 -ε k < s -b < r k + ε k (since x is inside the annulus X k ).
So without loss of generality we can assume that

S k ⊆ B(b, r k + ε k ) B(b, r k-1 -ε k ).
Using the union bound, we have:

Prob ∃x ∈ S s.t. T (x) = T (b) = Prob ∃x ∈ ∞ k=1 X k s.t. T (x) = T (b) ≤ ∞ k=1 Prob ∃x ∈ X k s.t. T (x) = T (b) .
Now, we will try to estimate the individual probabilities inside this sum. More precisely we will consider two cases according to the following event: for each k ≥ 1, we denote by E k the event that:

∃s ∈ S k , ∃x ∈ X k ∩ B(s, ε k ) s.t. T s -T x > ∆ k (with ∆ k ≥ 6ε k ).
Then we have

Prob ∃x ∈ X k s.t. T (x) = T (b) ≤ Prob (∃x ∈ X k s.t. T (x) = T (b)) ∧ Ēk + Prob(E k ) (14)
For the second term in [START_REF] Johnson | Extensions of Lipschitz mappings into a Hilbert space[END_REF], by the union bound, we have

Prob(E k ) ≤ s∈S k Prob ∃x ∈ X k ∩ B(s, ε k ) s.t. T s -T x > ∆ k ≤ s∈S k e -c1d( ∆ k ε k ) 2
(for c 1 a univ. constant as in Lemma 3.3)

≤ λ log 2 ( r k ε k ) S e -c1d( ∆ k ε k ) 2
.

(Note that here we must choose ∆ k ≥ 6ε k in order to apply Lemma 3.3).

For the first term in ( 14), we have

Prob (∃x ∈ X k s.t. T (x) = T (b)) ∧ Ēk ≤ Prob ∃x ∈ X k , s ∈ S k ∩ B(x, ε k ) s.t. T (x) = T (b) ∧ T (s) -T (x) ≤ ∆ k ≤ Prob ∃s ∈ S k s.t. T (s) -T (b) < ∆ k ≤ λ log 2 ( r k ε k ) S Prob T (z) < ∆ k r k-1 -ε k for some unit vector z ≤ λ log 2 ( r k ε k ) S ( √ 3∆ k r k-1 -ε k ) d (by inequality 5).
Putting all the estimations, we have obtained, together, we have:

Prob ∃x ∈ S s.t T (x) = T (b) ≤ ∞ k=1 λ log 2 ( r k ε k ) S e -c1d( ∆ k ε k ) 2 + ( √ 3∆ k r k-1 -ε k ) d . (15)
Next, we will show that there are choices of ε k , ∆ k , r k such that the RHS of ( 15) can be as small as needed.

Choices of ε k , ∆ k , r k : For some N > 1 large, we will choose ε k , ∆ k , r k as follows:

1.

ε k = ε, for some ε > 0. 2. ∆ k = N kε √ 3 3. r k = (N 2 (k + 1) 2 + 1)ε
Now the RHS of ( 15) can be rewritten as follows

∞ k=1 λ log 2 (N 2 (k+1) 2 +1) S e -c1d(N k) 2 + ( 1 N k ) d ≤ ∞ k=1 λ 3 log 2 (N k) S e -c1d(N k) 2 + ( 1 N k ) d ≤ ∞ k=1 2 3 ddim(S) log 2 (N k) e -c1d(N k) 2 + ( 1 N k ) d . (16)
Note that 2 3 ddim(S) log 2 (N k) does not increase fast enough compared to the decreasing speeds of both e -c1d(N k) 2 and ( 1N k ) d when d ≥ Cddim(S) with C large enough (and also independent of N ). Therefore, there are universal constants c 2 , c 3 > 0 such that the value of ( 16) is less than or equal to

∞ k=1 e -c2(N k) 2 + ∞ k=1 ( 1 N k ) c3d (17)
Both the two infinite sums tend to 0 when N tends to ∞. This means that Prob ∃x ∈ S s.t. T (x) = T (b) = 0, which proves our theorem. 2

Our next result in this section is an extension of Thm. 3.4 to the threshold case.

Theorem Let b /

∈ S where S ⊆ R m is a closed doubling set, T : R m → R d be a normalized Gaussian random projection, and r = min x∈S b -x . Let κ < 1 be some fixed constant. Then for all 0 < δ < 1 and 0 < τ < κr, we have

Prob(dist(T (b), T (S)) > τ ) > 1 -δ if the projected dimension is d = Ω( log( λ S δ ) 1-κ ).
Proof. We follow the same scheme, as in the proof of Theorem 3.4. Let

d ≥ C( log( λ S δ )
1-κ ) for some universal constant C (large). As before, let ε i , ∆ i , i ∈ N and r = r 0 < r 1 < r 2 < . . . be positive scalars whose values will be decided later. The annuli X k and point sets S k are also defined similarly. Using the union bound, now we have:

Prob ∃x ∈ S s.t. T (x) -T (b) < τ = Prob ∃x ∈ ∞ k=1 X k s.t. T (x) -T (b) < τ ≤ ∞ k=1 Prob ∃x ∈ X k s.t. T (x) -T (b) < τ .
Now, we will try to estimate the individual probabilities inside this sum.

For each k ≥ 1, we denote by E k the event that:

∃s ∈ S k , ∃x ∈ X k ∩ B(s, ε k ) s.t. T s -T x > ∆ k (with ∆ k ≥ 6ε k ). ( 18
)
Then we have

Prob ∃x ∈ X k s.t. T (x) -T (b) < τ ≤ Prob (∃x ∈ X k s.t. T (x) -T (b) < τ ) ∧ Ēk + Prob(E k ) (19)
For the second term in [START_REF] Matoušek | On variants of the Johnson-Lindenstrauss lemma[END_REF], from the previous proof, we already had:

Prob(E k ) ≤ λ log 2 ( r k ε k ) S e -c1d( ∆ k ε k ) 2 . ( 20
)
(Note that here we must choose ∆ k ≥ 6ε k in order to apply Lemma 3.3). Now, for the first term in [START_REF] Matoušek | On variants of the Johnson-Lindenstrauss lemma[END_REF], we have

Prob (∃x ∈ X k s.t. T (x) -T (b) < τ ) ∧ Ēk ≤ Prob ∃x ∈ X k , s ∈ S k ∩ B(x, ε k ) s.t. T (x) -T (b) < τ ∧ T (s) -T (x) ≤ ∆ k ≤ Prob ∃s ∈ S k s.t. T (s) -T (b) < ∆ k + τ (by triangle inequality) ≤    λ log 2 ( r k ε k ) S Prob T (z) < ∆ k +τ r k-1 -ε k for some unit vector z if k ≥ 2 λ log 2 ( r 1 ε 1 ) S Prob T (z) < ∆1+τ r for some unit vector z if k = 1 ≤    λ log 2 ( r k ε k ) S ( ∆ k +τ r k-1 -ε k ) d if k ≥ 2 λ log 2 ( r k ε k ) S ( ∆1+τ r ) d if k = 1.
Putting all the estimations we have obtained together, we have:

Prob ∃x ∈ S s.t. T (x) -T (b) < τ ≤ ∞ k=1 λ log 2 ( r k ε k ) S e -c 1 d( ∆ k ε k ) 2 + ∞ k=2 λ log 2 ( r k ε k ) S ( ∆ k + τ r k-1 -ε k ) d + λ log 2 ( r 1 ε 1 ) S ( ∆1 + τ r ) d . (21)
Here we separate one term out, and we will prove that the remaining expression can be made as small as wanted by certain choices of parameters. Choices of ε k , ∆ k , r k : Let N > 0 be the number such that ( 7 N + 1) = r τ . From the assumptions, we have N = 7 τ r-τ and 0 < τ < κr. Hence N < 7κr r-τ < 7κ 1-τ r < 7κ 1-κ < 7 1-κ . We will choose ε k , ∆ k , r k as follows:

1. ε k = ε = τ N , 2. ∆ k = 6 √ kε, 3. r k = (6k + 7)ε + √ k + 1 τ .
(Our purpose is to choose the parameters so that

∆ k +τ r k-1 -ε k = 1 √ k and ∆ k ≥ 6ε k ).
From this choice, it is obvious r 0 = r. Now the RHS of ( 21) can be rewritten as follows:

∞ k=1 λ log 2 (6k+7+N √ k+1) S e -c1d(36k) + ∞ k=2 λ log 2 (6k+7+N √ k+1) S ( 1 √ k ) d + + λ log 2 (13+N √ 2) S (( 6 N + 1) τ r ) d ≤ ≤ ∞ k=1 λ c3 log 2 (N (k+1)) S e -c1d(36k) + ∞ k=2 λ c3 log 2 (N (k+1)) S ( 1 √ k ) d + + λ c2 S (( 6 N + 1) τ r ) d (22)
for some universal constants c 1 , c 2 , c 3 .

It is easy to show that the expression in the big bracket is bounded above by e -c4d as long as d ≥ C log 2 (λ S ) log(

However, log(( N 6+N ) r τ ) = log( 7r 6r+τ ) ≥ log( 7 6+κ ) ≥ log(1 + 1-κ 7 ) ≥ 1-κ 7 , from the Taylor series of the logarithm function. Therefore, (24) holds if we select

d ≥ C log( λ S δ ) 1 -κ (25)
for some universal constants C. The proof follows immediately from ( 23) and ( 25) by an application of the union bound. 2

One of the interesting consequence of Theorem 3.5 is the following application to the Approximate Nearest Neighbour problem.

Corollary

For T : R m → R d be a normalized Gaussian random projection, X ⊆ R m , ε ∈ (0, 1/2) and δ ∈ (0, 1/2), we can choose

d = max O( log( 1 δ ) ε 2 ), O( log( λ S δ ) ε )
such that for every x 0 ∈ X, with probability at least 1 -δ, 1. dist (T x 0 , T (X \ {x 0 })) ≤ (1 + ε)dist (x 0 , X \ {x 0 }), 2. Every x ∈ X with x 0 -x ≥ (1 + 2ε)dist (x 0 , X \ {x 0 }) satisfies T x 0 -T x > (1 + ε)dist (x 0 , X \ {x 0 }).

Note that, this result improves the bound provided by Indyk-Naor in [START_REF] Indyk | Nearest neighbor preserving embeddings[END_REF]. In that paper, the authors give the bound for the projected dimension to be with some large constant C:

Prob [dist (T x 0 , T (X \ {x 0 })) ≤ (1 + ε)dist (x 0 , X \ {x 0 })] < δ 2 . (26)
Now, assume that dist (x 0 , X \ {x 0 }) = 1. Set b = x 0 and S = {x ∈ X :

x -x 0 ≥ 1 + 2ε} and τ = 1 + ε as in Theorem 3.5. We then have r := min x∈S x -b ≥ 1 + 2ε, which implies τ r ≤ 1+ε 1+2ε = 1 -ε 1+2ε < 1 -ε 2 . Therefore, we can choose κ = 1 -ε/2. Applying Theorem 3.5, we have: ). From (26) and (27), we conclude that the two required conditions hold for some

d = max O( log( 1 δ ) ε 2 ), O( log( λ S δ ) ε ) ,
as claimed. 2

Conclusion

In this paper we discussed the possibility of applying random projections of the type used in Johnson-Lindenstrauss lemma to a very general class of set membership problems in Euclidean spaces. Our results are directly applicable, as it suffices to pre-multiply the input vectors by a randomly generated matrix, and yield the correct answer with high probability.

  where S is a finite set. Let R = min x∈S b -x > 0 and T : R m → R d be a normalized Gaussian random projection with d ≥ log(|S|/δ) log(R/( √ 3τ )) . Then: Prob min x∈S T (b) -T (x) > τ > 1 -δ.

  λ S ) , which is significantly larger than our bound. Proof. Similar as the proof of Theorem 4.1 in[START_REF] Indyk | Nearest neighbor preserving embeddings[END_REF], we have:for d ≥ C log( 1 δ ) ε 2

  Prob(dist(T (b), T (S)) ≤ τ )

Acknowledgements

The authors are grateful to the Siebel Energy Institute for a seed grant.


Esteban Salgado

Claudio Gentile
email: gentile@iasi.cnr.it

Leo Liberti
email: liberti@lix.polytechnique.fr

Perspective cuts for the ACOPF with generators

The alternating current optimal power flow problem is a fundamental problem in the management of smart grids. In this paper we consider a variant which includes activation/deactivation of generators at some of the grid sites. We formulate the problem as a mathematical program, prove its NP-hardness w.r.t. activation/deactivation, and derive two perspective reformulations.

Introduction

The Alternating Current Optimal Power Flow (ACOPF) problem is one of the most important problems arising in the energy industry. It models the propagation of power flows in electrical grids. It is often used as second-level subproblem in bilevel problems modelling the decision of electricity prices subject to production and demands [14]. Multilevel problems with ACOPF at different time-scales are also considered [START_REF] Bacher | Power system models, objectives and constraints in optimal power flow calculations[END_REF]. The ACOPF received a lot of attention over the years, and specifically after smart grids were introduced [START_REF] Bienstock | Electrical Transmission System Cascades and Vulnerability: an Operations Research Viewpoint[END_REF].

The ACOPF asks for the best power flow over an electrical network modelled by a digraph D = (N , L ), where N is the set of buses and L the set of lines. It is well known that the natural formulation can be simplified using only voltage variables [START_REF] Josz | Application of polynomial optimization to electricity transmission networks[END_REF]. The ACOPF is usually cast as a Mathematical Programming (MP) problem over the complex numbers (which make their appearance due to the cyclic nature of alternating currents). The standard ACOPF can be reformulated as a (larger) MP over the reals, by separating real and complex parts [START_REF] Salgado | Fast relaxations for alternating current optimal power flow[END_REF].

While the standard version of the ACOPF only has continuous variables, more realistic variants include binary variables which activate/deactivate various electrical components. In this paper we consider the possibility of activating/deactivating electrical generation at some of the buses. This defines an ACOPF variant which we call ACOPF with Generators (ACOPFG) [START_REF] Ruiz | A progressive method to solve large-scale AC optimal power flow with discrete variables and control of the feasibility[END_REF].

Note that the ACOPF is NP-hard even without binary variables, as shown in [START_REF] Lehmann | AC-Feasibility on tree networks is NP-hard[END_REF]. Experimentally, however, it was found that many standard benchmarks, as well as randomly generated instances, can be solved efficiently. It is shown in [START_REF] Lavaei | Zero duality gap in optimal power flow problem[END_REF] that this happens whenever the duality gap is zero. One might then question whether the ACOPFG is NP-hard simply because of the addition of the binary activation variables. The first contribution of this paper is to prove that this is indeed the case.

While ACOPF objective functions vary in the literature, it is common to consider quadratic objectives with respect to voltage. In this paper, however, we focus on a more general objective function, quadratic with respect to active power and quartic (without cubic terms) w.r.t. voltage [START_REF] Josz | Application of polynomial optimization to electricity transmission networks[END_REF]. The second contribution of this paper is the application of two perspective reformulations (PR) to the ACOPFG with the more general (quartic) objective [START_REF] Frangioni | Perspective cuts for a class of convex 0-1 mixed integer programs[END_REF][START_REF] Frangioni | Approximated perspective relaxations: a project&lift approach[END_REF].

MP formulation

We consider the network digraph D mentioned in Sect. 1. Let n = |N | and = |L |. We identify a subset G of generator buses, and let n = |G |. We note that, in modern "smart grids", generators may produce and consume electricity. Because we are dealing with alternating currents, power is represented by a complex number. The real part is the active power while the complex part is reactive.

Notationwise, we use [α, α] to denote lower/upper bounds to a quantity, and α * to denote complex conjugate.

The parameters of our problem are as follows:

• ∀b ∈ N S b ∈ C is the power demand at bus b;

• ∀g ∈ G S g = [S g , S g ] is the (complex) interval where g can generate power if active;

• ∀b ∈ N v b = [v b , v b ]
is the (real) interval where the voltage magnitude at bus b can range; • ∀(a, b) ∈ L ı ab is the maximum current which can flow through the line (a, b); • Y is a complex n × n bus admittance matrix (it plays a role analoguous to the reciprocal of resistance in Ohm's law); • Y 0 ,Y 1 are complex × n line admittance matrices (they "encode" some electrical properties of the lines).

The decision variables are:

• ∀g ∈ G s g ∈ C is the power generated at g; • ∀g ∈ G z g ∈ {0, 1} denotes the deactivation (0) or activation (1) of generator g;

• ∀b ∈ N v b ∈ C is the voltage at bus b; • ∀(a, b) ∈ L i ab ∈ C is the current on the line (a, b).
At each generator g ∈ G , the injected complex power s g -S g = v g ∑ (g,a)∈L i * ga , and at each non-generator bus b ∈ N G , we have -S b = v b ∑ (b,a)∈L i * ba . Kirchoff's law and a generalized form of Ohm's law allow us to derive i = Y v, which implies that the RHS of the above equations can be reformulated to

v b (Y * v * ) b = ∑ (a,b)∈L v b v * a Y * ab for each b ∈ N [17]
. This allows us to express current in function of voltage and power. We obtain the following constraints:

∀g ∈ G ∑ (g,a)∈L Y * ga v g v * a = s g z g -S g (1)
∀b ∈ N G ∑ (b,a)∈L Y * ba v b v * a = -S b (2)
∀(a, b) ∈ L , ω ∈ {0, 1} ∑ h =k∈N (Y ω abh ) * (Y ω abk ) * v * h v k ≤ ı ab (3) ∀g ∈ G s g ∈ S g (4) ∀b ∈ N |v b | ∈ v b (5) ∀g ∈ G z g ∈ {0, 1}. (6)
We remark that complex power variables s only appear in Eq. ( 1) and ( 4). We can eliminate them by replacing Eq. ( 1) and ( 4) with the following inequalities:

∀g ∈ G S g z g ≤ ∑ (g,a)∈L Y * ga v g v * a + S g ≤ S g z g . (7)
Moreover, if we define z over all of N and fix z b = 0 for all b ∈ G , Eq. ( 7) quantified on N can also replace Eq. ( 2).

In the ACOPF literature [START_REF] Lavaei | Zero duality gap in optimal power flow problem[END_REF][START_REF] Kuang | Alternative LP and SOCP hierarchies for ACOPF problems[END_REF][START_REF] Chen | Bound tightening for the alternating current optimal power flow problem[END_REF][START_REF] Ruiz | A progressive method to solve large-scale AC optimal power flow with discrete variables and control of the feasibility[END_REF] we consider the following generation cost function, to be minimized:

f (s, z) = ∑ g∈G z g (c g2 (Re(s g )) 2 + c g1 Re(s g ) + c g0 ). (8)
Again we can replace s by v using Eq. ( 1) and removing constant terms in order to express Eq. ( 8) as a function of voltage: essentially, we obtain f (v, z) from Eq. ( 8) by replacing s g with ∑ (g,a)∈L Y * ga v g v * a + S g . Let F be the feasible subset of C n defined by Eq. ( 2)-( 6) and [START_REF] Frangioni | Improving the approximated projected perspective reformulation by dual information[END_REF]. We call

ACOPFG C the formulation min (v,z)∈F f (v, z).
Finally, we can obtain a real formulation as follows:

1. replace each quadratic constraint v H Mv ♦ α + jβ (where ♦ ∈ {=, ≤, ≥} and j = √ -1) by the pair of constraints

v H M + v ♦ α + jβ ∧ v H M -v ♦ α + jβ , where M + = 1 2 (M + M H ) and M -= 1 2 (M -M H ); 2.

replace each complex matrix M by the real matrix

Re(M) -Im(M) Im(M) Re(M) ;

3. replace each complex vector v by the real vector (Re(v) Im(v)) .

We call this reformulation ACOPFG R .

Complexity

Assume c g2 = 0 for all g ∈ G in Eq. [START_REF] Frangioni | Perspective cuts for a class of convex 0-1 mixed integer programs[END_REF]. By ignoring activation variables we obtain the ACOPF, which is a Quadratically Constrained Quadratic Program (QCQP).

Since the ACOPF is NP-hard [START_REF] Bienstock | Strong NP-hardness of AC power flows feasibility[END_REF], it follows by inclusion that the ACOPFG is also NP-hard. On the other hand, it was shown in [START_REF] Lavaei | Zero duality gap in optimal power flow problem[END_REF] that many practical ACOPF instances turn out to be easy rather than hard. We remark that "easy", in this setting, does not necessarily mean "in P", since the decision version of the QCQP is not known to be in NP (unless there are no quadratic constraints, in which case the problem class is known to be in NP [START_REF] Vavasis | Quadratic programming is in NP[END_REF]). The meaning of "easy" in this context is that global optima can be obtained by means of a local, rather than global, optimization procedure.

The question we answer in this section is whether the addition of the binary activation variables constitute an actual additional difficulty. To show that this need not necessarily be the case, we consider a linear system

Ax ≤ b x ≤ 1 x ∈ R n + ,    ( 9
)
where A is totally unimodular. Finding a feasible solution can obviously be done in polynomial time by the, say, interior point algorithm (irrespective of total unimodularity), and so this formulation is in P. If we add n additional binary activation variables y 1 , . . . , y n ∈ {0, 1} and n additional activation/deactivation constraints

∀ j ≤ n x j ≤ y j , (10)
then the new system has a constraint matrix:

  A 0 I n 0 I n -I n   ,
which is easily seen to also be totally unimodular [START_REF] Wolsey | Integer Programming[END_REF]. Therefore this new Mixed-Integer Linear Programming (MILP) formulation is in P. This provides an example where adding boolean activation variables does not make the underlying problem more difficult.

Having established that the question makes sense, we present a reduction of the weakly NP-complete SUBSET-SUM problem to a subclass of ACOPFG. Given an instance (σ 1 , . . . , σ n , S 0 ) of SUBSET-SUM, we must decide whether there is a subset I ⊆ {1, . . . , n} such that

∑ i∈I σ i = S 0 . (11)
This is equivalent to asking whether the following linear diophantine equation has a solution x ∈ {0, 1} n :

∑ i≤n σ i x i = S 0 . (12)
We now show that we can naturally express Eq. ( 12) using the ACOPFG formulation of Sect. 2. We consider a simple network D with G = {1, . . . , n} generators with demand S g = 0 for g ≤ n and a single non-generator bus (indexed by 0) with demand S 0 (so that N = {0, . . . , n}). The set L of lines is {(g, 0) | 1 ≤ g ≤ n}, namely each generator is linked to the only non-generator bus. Each generator g ∈ G has generation interval S G = [σ g , σ g ], i.e. each generator can either be inactive, or else, if active, must produce exactly σ g . Then Eq. ( 7) becomes:

∀g ≤ n Y * g0 v g v * 0 = σ g z g . (13)
Since we know σ g is real and positive, we arrange Y * so that the complex part of the LHS of Eq. ( 13) is zero; in particular, we arrange Y * g0 v * 0 to yield a j 2 = -1 coefficient (this can be easily done when we derive ACOPFG R ). So we get:

∀g ≤ n Re(Y * g0 v g v * 0 ) = -σ g z g . (14)
Furthermore, Eq. ( 2) is:

∑ (g0)∈L Y * g0 v 0 v * g = -S 0 ,
whence, by Eq. ( 14), we have:

∑ g≤n (-σ g z g ) = -S 0 ,
which is exactly Eq. ( 12).

Perspective reformulation

The objective function Eq. ( 8) can be restated using additional variables

p g = Re(v g ∑ (g,a)∈L v * a Y * ga + S g ). In practice use the convex constraints ∀g ∈ G s.t. c g2 > 0 p g ≥ Re(v g ∑ (g,a)∈L v * a Y * ga + S g ), (15)
which are justified by the objective function direction. We now reformulate Eq. ( 8) using these new variables:

f (p, z) = ∑ g∈G (c g2 p 2 g + c g1 p g + c g0 z g ). ( 16
)
The reformulations proposed below can all be carried out on a per-generator basis.

In the rest of this paper, we assume they are only applied to generators g ∈ G for which c g2 > 0.

The power p g is subject to the following activation constraints:

p g ≤ P g z g ∧ p g ≥ P g z g ( 17
)
where P g = Re(S g ) and P g = Re(S g ). The PR reformulation [START_REF] Frangioni | Perspective cuts for a class of convex 0-1 mixed integer programs[END_REF] can be applied to ( 16) as follows:

f (p, z) = ∑ g∈G c g2 p 2 g z g + c g1 p g + c g0 z g . ( 18
)
The function ( 18) can be optimized using the perspective cuts (PC) method [START_REF] Frangioni | Perspective cuts for a class of convex 0-1 mixed integer programs[END_REF], which works as follows: (i) first we add new variables t g representing the nonlinear part of the cost in ( 18) by considering the following constraints

t g ≥ c g2 p 2 g z g , (19)
and replacing [START_REF] Salgado | Alternating current optimal power flow with generator selection[END_REF] with f (t, p, z) = ∑ g∈G (t g +c g1 p g +c g0 z g ); (ii) then constraints [START_REF] Vavasis | Quadratic programming is in NP[END_REF] can be replaced by PCs:

t g ≥ c g2 (2 pg p g -p2 g z g ), (20)
where pg are fixed values of the real power p g varying in the feasible interval P g ≤ pg ≤ Pg when z g = 1. The addition of PCs does not add further difficulties in the problem formulation except for the condition that they should be generated iteratively as their number is not finite. We can alternatively apply the AP2R technique [START_REF] Frangioni | Approximated perspective relaxations: a project&lift approach[END_REF][START_REF] Frangioni | Improving the approximated projected perspective reformulation by dual information[END_REF], which works in two phases. The first phase is a projection where the optimal value of z g for the continuous relaxation of Eq. ( 18) subject to Eq. ( 17) is found depending on p g . The second phase is a lifting where the variables z g are lifted back. The resulting problem can be solved using an off-the-shelf MIP solver. This is equivalent to replacing [START_REF] Salgado | Alternating current optimal power flow with generator selection[END_REF] and ( 17) with:

min ∑ g∈G z g (p int g ) 2 + fg (π g + p int g ) -fg (p int g ) + c g1 p g + c g0 z g ∀g ∈ G (P g -p int g )z g ≤ π g ≤ (P g -p int g )z g ∀g ∈ G π g = p g -p int g z g ,      (21)
where fg (x) = c g2 x 2 and p int g is p int g = max(P g , min( c g0 /c g2 , P g )).

(

The final AP2R reformulation consists in Eq. ( 21), Eq. ( 2)-( 6), the complex part of Eq. ( 7), and Eq. ( 15).

Computational results

We tested PRs with 4 cuts and AP2R (implemented using AMPL [START_REF] Fourer | AMPL: A Modeling Language for Mathematical Programming[END_REF]): both on the ACOPFG formulation ACOPFG R in Sect. 2 (Table 1) and on the dual Diagonally Dominant Programming (DDP) outer-approximation proposed in [START_REF] Salgado | Alternating current optimal power flow with generator selection[END_REF] (Table 2) solved using CPLEX [START_REF]ILOG CPLEX 12.6 User's Manual[END_REF]. We compared these results with local optima of the ACOPF (all active generators) obtained by MatPower [START_REF] Zimmermann | MATPOWER: Steady-state operations, planning, and analysis tools for power systems research and education[END_REF] and by solving the ACOPFG R using Baron [START_REF] Sahinidis | [END_REF] to global optimality (within a limited CPU time of 1h 1 Results on ACOPFG R ('x': solution not found within time limit).

In the "Perspective reformulation" columns we show: number of iterations, CPU seconds (limited to 1h), PR objective value obtained on 1st iteration and final value, original objective function value at optimum, percentage of active generators at optimum. In the "AP2R" and "Mi-Quartic" columns we show CPU time, objective function value and percentage of active generators. In "Solution's distances" we report a scaled distance of the optima found by PR/AP2R w.r.t. Mi-Quartic, namely p ω -p MI-Quartic Table 2 Results on dual DDP [START_REF] Salgado | Alternating current optimal power flow with generator selection[END_REF] ('x': solution not found within time limit).

While it is clear that the tests with ACOPFG R are inconclusive, those on the dual DDP approximation give very tight bounds in relatively little time.

  ). The test set includes small to medium scale instances taken from MatPower; results on one larger-scale instance are reported in Table2. All results were obtained on an Intel i7 dual-core CPU at 2.1GHz with 16GB RAM.

	Perspective reformulation (ACOPFGR)		AP2R (ACOPFGR)	MATPOWER		Mi-Quartic	Solution's distances
	Instance it time first value last value % active real value time %active value real value	value	time %active best value Persp/Mi-Quartic AP2R/Mi-Quartic
	WB2 2 24 878.182 878.182	100	878.182	13	100	878.18 878.182	877.78	13	100	878.182	0	0
	WB3 2 169 417.244 417.244	100	417.244 109 100 417.244 417.244	417.25	108 100	417.244	0	0
	WB5 2 3600 947.056 947.056	100	947.056 2269 100 947.056 947.056	1082.33	2454 100	947.056	0	0
	6ww 1 3600 2913.58 2913.58	x	2881.28 3600 100 10948.7 3135.18	3134.35	3600 100	3135.18	0.1584	0
	case9 2 3600 2062.65 5115.72	100	5430.38 3600 66.7 10948.74 7335.42	5296.69	3600 100	5296.69	0.2561	0.6234
	case14 1 3600 5250.22 5250.22	x	5375.94 3600 80	6589.72 5287.72	8081.53	3600 60	5476.90	0.1370	0.8304
	case30 1 3600 430.906 430.906	x	536.307 3600 66.7 503.508 503.508	576.89	3600 83.3 515.807	0.4234	0.3503
	Table

Acknowledgements. This paper has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement n. 764759. The second author is partially supported by MIUR PRIN2015 project no. 2015B5F27W. The last author (LL) is grateful for financial support by CNR STM Program prot. AMMCNT-CNR n. 80058 dated 05/12/2017.