Charlotte Truchet
email: charlotte.truchet@univ-nantes.fr

Alejandro Arbelaez
email: alejandro.arbelaez@insight-centre.org

Florian Richoux
email: florian.richoux@univ-nantes.fr

Philippe Codognet
email: codognet@is.s.u-tokyo.ac.jp

Estimating parallel runtimes for randomized algorithms in constraint solving

come

Introduction

In the last years, parallel algorithms for solving hard combinatorial problems, such as Constraint Satisfaction Problems (CSP), have been of increasing interest in the scientific community. The combinatorial nature of the problem makes it difficult to parallelize existing solvers, without costly communication schemes. Several parallel schemes have been proposed for incomplete or complete solvers, one of the most popular (as observed in the latest SAT competitions www.satcompetition.org) being to run several competing instances of the algorithm on different processes with different initial conditions or parameters, and let the fastest process win over others. The resulting parallel algorithm thus terminates with the minimal runtime among the launched processes. The framework of independent multi-walk parallelism, seems to be a promising way to deal with large-scale parallelism. Cooperative algorithms might perform well on shared-memory machines with a few tens of processors, but are difficult to extend efficiently on distributed hardware. This leads to so-called independent multi-walk algorithms in the CSP community [START_REF] Verhoeven | Parallel Local Search[END_REF] and portfolio algorithms in the SAT community (satisfiability of Boolean formulae) [START_REF] Gomes | Algorithm Portfolios[END_REF].

However, although it is easy to obtain good Speed-up on a small-scale parallel machine (viz. with a few tens of processes), it is not easy to know how a parallel variant of a given algorithm would perform on a massively parallel machine (viz. with thousands of processes). Parallel performance models are thus particularly important for parallel constraint solvers, and any indication on how a given algorithm (or, more precisely, a pair formed by the algorithm and the problem instance) would scale on massively parallel hardware is valuable. If it becomes possible to estimate the maximum number of processes until which parallelization is efficient, then the actual parallel computing power needed to solve a problem could be deduced. This piece of information might be quite relevant, since supercomputers or systems such as Google Cloud and Amazon EC2 can be rented by processor-hour with a limit on the maximum number of processors to be used. In this context, modelling tools for the behavior of parallel algorithms are expected to be very valuable in the future.

The goal of this paper is to study the parallel performance of randomized constraint solving algorithms under the independent multi-walk scheme, and to model the performance of the parallel execution from the runtime distribution of sequential runs of a given algorithm. Randomized constraint solvers considered in this paper include Local Search algorithms for Constraint Satisfaction Problems, Local Search techniques for SAT, and complete algorithms with random components e.g., a propagation-based backtrack search with random heuristics. An important application of this work relates to the increasing computational power being available in cloud systems (e.g., Amazon Cloud EC2, Google Cloud and Microsoft Azure), a good estimate on how the algorithm scales might allow users to rent just the right number of cores. Most papers on the performance of stochastic Local Search algorithms focus on the average runtime in order to measure the performance of both sequential and parallel executions. However, a more detailed analysis of the runtime behavior could be done by looking at the runtime of the algorithm (e.g., CPU-time or number of iterations) as a random variable and performing a statistical analysis of its probability distribution. More precisely, we first approximate the empirical sequential runtime distribution by a well-known statistical distribution (e.g., exponential or lognormal) and then derive the runtime distribution of the parallel version of the solver. Our model is related to order statistics, a rather recent domain of statistics [START_REF] David | Order Statistics. Wiley series in probability and mathematical statistics[END_REF], which is the statistics of sorted random draws. Our method encompasses any algorithm whose solving time is random and makes it possible to formally determine the average parallel runtime of such algorithms for any number of processors.

For Local Search, we will consider algorithms in the framework of Las Vegas algorithms [START_REF] Babai | Monte-Carlo Algorithms in Graph Isomorphism Testing[END_REF], a class of algorithms related to Monte-Carlo algorithms introduced a few decades ago, whose runtime may vary from one execution to another, even on the same input. The classical parallelization scheme of multi-walks for Local Search methods can easily be generalized to any Las Vegas algorithm. We will study two different sets of algorithms and problems: first, a Constraint-Based Local Search solver on CSP instances, and, secondly, two SAT Local Search solvers on random and crafted instances. Interestingly, this general framework encompasses other types of Las Vegas algorithms, and we will also apply it to a propagation-based constraint solver with a random labeling procedure on CSP instances.

We will confront the performance predicted by the statistical model with actual speed-ups obtained for parallel implementations of the above-mentioned algorithms and show that the prediction can be quite accurate, matching the actual speed-up up to a large number of processors. More interestingly, we can also model both the initial and the asymptotic behavior of the parallel algorithm.

This paper extends [START_REF] Truchet | Prediction of Parallel Speed-ups for Las Vegas Algorithms[END_REF] and [START_REF] Arbelaez | Using Sequential Runtime Distributions for the Parallel Speedup Prediction of SAT Local Search[END_REF] by giving a detailed presentation of the runtime estimation model, based on order statistics, and by validating the model on randomized propagation-based constraint solvers, extensive experimental results for stochastic local search algorithms on well-known CSP instances from CSPLib and SAT instances obtained from the international SAT competition. Additionally, we provide a more detailed theoretical analysis of the reference distributions used for predicting the parallel performance.

The paper is organized as follows. Section 2 presents the existing approaches in parallel constraint solving, and formulates the question we address in the following. Section 3 details our probabilistic model for the class of parallel algorithms we tackle in this article, based on Las Vegas algorithms. Several such algorithms can be used in constraint solving. Each of the three last sections is dedicated to a specific family of Las Vegas algorithm: Constraint-Based Local Search in Section 4, SAT Local Search in Section 5 and propagationbased methods with randomization in Section 6. For each method, we apply the model to compute the parallel speed-ups, run the parallel algorithm in practice and compare and analyze the results. Section 7 presents a comparison of our probabilistic method against a method that does not involve statistical tools. A conclusion, in Section 8, ends the paper.

Parallel Constraint Solving

This section presents existing approaches for parallel constraint solving. The base methods are Constraint-Based Local Search, Local Search for SAT, and complete solvers based on propagation techniques and backtrack search. We then present the state of the art on estimating parallel speed-up for randomized solvers.

Parallel Local Search

Parallel implementation of Local Search metaheuristics [START_REF]Handbook of Approximation Algorithms and Metaheuristics[END_REF][START_REF]Metaheuristics: Progress as Real Problem Solvers[END_REF] has been studied since the early 1990s, when parallel machines started to become widely available [START_REF] Pardalos | Parallel Search for Combinatorial Optimization: Genetic Algorithms, Simulated Annealing, Tabu Search and GRASP[END_REF][START_REF] Verhoeven | Parallel Local Search[END_REF]. With the increasing availability of PC clusters in the early 2000s, this domain became active again [START_REF] Alba | Special Issue on New Advances on Parallel Meta-Heuristics for Complex Problems[END_REF][START_REF] Crainic | Special Issue on Parallel Meta-Heuristics[END_REF]. [START_REF] Moisan | Parallel Discrepancy-Based Search[END_REF] recently proposed a scalable parallel algorithm for backtracking algorithms based on load-balancing to distribute the work across multiple processes with low communication overhead. Apart from domain-decomposition methods and population-based method (such as genetic algorithms), [START_REF] Verhoeven | Parallel Local Search[END_REF] distinguishes between single-walk and multi-walk methods for Local Search. Single-walk methods consist in using parallelism inside a single search process, e.g., for parallelizing the exploration of the neighborhood (see for instance [START_REF] Van Luong | Local Search Algorithms on Graphics Processing units[END_REF] for such a method making use of GPUs for the parallel phase). Multi-walk methods, i.e., parallel execution of multi-start methods, consist in developing concurrent explorations of the search space, either independently or cooperatively with some communication between concurrent processes. Sophisticated cooperative strategies for multi-walk methods can be devised by using solution pools [START_REF] Crainic | Cooperative Parallel Variable Neighborhood Search for the -Median[END_REF], but require shared-memory or emulation of central memory in distributed clusters, thus impacting on performance. A key point is that a multi-walk scheme is easier to implement on parallel computers without shared memory and can lead, in theory at least, to linear speed-ups [START_REF] Verhoeven | Parallel Local Search[END_REF]. However, this is only true under certain assumptions and we will see that we need to develop a more realistic model in order to cope with the performance actually observed in parallel executions.

Parallel Local Search for SAT

Like for CSP, it is now currently admitted that an easy and effective manner to parallelize Local Search solvers consists in executing in parallel multiple copies of the solver with or without cooperation. We remark that nearly all parallel solvers in the parallel SAT competitions are based on the multi-walk framework 1 . The non-cooperative approach has been used in the past to solve SAT and MaxSAT instances. gNovelty+ [START_REF] Pham | gNovelty+. In: Solver description, SAT competition[END_REF] executes multiple copies of gNovelty without cooperation until a solution is obtained or a given timeout is reached; 1 The SAT community usually refer to the multi-walk framework as portfolio algorithms.

and [START_REF] Pardalos | A Parallel Grasp for MAX-SAT Problems[END_REF] executes multiple copies of GRASP until an assignment which satisfies a given number of clauses is obtained. Strategies to exploit cooperation between parallel SAT Local Search solvers have been studied in [START_REF] Arbelaez | Improving Parallel Local Search for SAT[END_REF] in the context of multi-core architectures with shared memory and in [START_REF] Arbelaez | Massivelly Parallel Local Search for SAT[END_REF] in massively parallel systems with distributed memory.

Parallel Complete Solvers

Early experiments on the parallelization of propagation-based complete solvers date back to the beginning of Constraint Logic Programming, cf. [START_REF] Van Hentenryck | Parallel Constraint Satisfaction in Logic Programming: Preliminary Results of CHIP within PEPSys[END_REF], and used the search parallelism of the host logic language [START_REF] De Kergommeaux | Parallel Logic Programming Systems[END_REF]. Most of the proposed implementations have been based on some kind of OR-parallelism, splitting the search space between different processors and relying on the sharedmemory multicore architecture as the different processors work on shared data-structures representing a global environment in which the subcomputations take place. However only very few implementations of efficient constraint solvers on such machines have ever been reported, for instance [START_REF] Perron | Search Procedures and Parallelism in Constraint Programming[END_REF] or [START_REF] Chu | Confidence-Based Work Stealing in Parallel Constraint Programming[END_REF] for a shared-memory architecture with 8 CPU processors. The Comet system [START_REF] Van Hentenryck | Constraint-Based Local Search[END_REF] has been parallelized for small clusters of PCs, both for its Local Search solver [START_REF] Michel | Distributed Constraint-Based Local Search[END_REF] and its propagation-based constraint solver [START_REF] Michel | Parallelizing Constraint Programs Transparently[END_REF]. More recent experiments have been done up to 12 processors [START_REF] Michel | Parallel and Distributed Local Search in Comet[END_REF].

Search-space splitting techniques such as domain decomposition have also been implemented in the domain of Constraint Programming, but initial experiments [START_REF] Bordeaux | Experiments with Massively Parallel Constraint Solving[END_REF] show that the speed-up goes to flatten after a few tens of processors, thus away from linear speed-up. A recent approach based on a smaller granularity domain decomposition [START_REF] Régin | Embarrassingly Parallel Search[END_REF] shows better performance. The results for all-solution search on classical CSPLib benchmarks are quite encouraging and show an average speed-up of 14 to 20 up with 40 processors w.r.t. the base sequential solver.

In the domain of combinatorial optimization, the most popular complete method that has been parallelized at a large scale is the classical branch and bound method [START_REF] Gendron | Parallel Branch-And-Bound Algorithms: Survey and Synthesis[END_REF], because it does not require much information to be communicated between parallel processes: basically only the current upper (or lower) bound of the solution. It has thus been a method of choice for experimenting the solving of optimization problems using grid computing, see for instance [START_REF] Aida | A Case Study in Running a Parallel Branch and Bound Application on the Grid[END_REF] and also [START_REF] Caromel | Grid'BnB: A Parallel Branch and Bound Framework for Grids[END_REF], which use several hundreds of nodes of the Grid'5000 platform. Good speed-ups are achieved up to a few hundreds of processors but interestingly, their conclusion is that runtimes tend to stabilize afterward. A simple solving method for project scheduling problems has also been implemented on an IBM Bluegene/P supercomputer [START_REF] Xie | Massively Parallel Constraint Programming for Supercomputers: Challenges and Initial Results[END_REF] up to 1,024 processors, but with mixed results since they reach linear speed-ups until 512 processors only, and then no improvements beyond this limit.

How to Estimate Parallel Speed-Up ?

The multi-walk parallel scheme is rather simple, yet it provides an interesting test-case to study how Las Vegas algorithms can scale-up in parallel. Indeed, runtime will vary among the processes launched in parallel and the overall runtime will be equal to the minimal runtime (i.e., "long" runs are killed by "shorter" ones). The question is thus to quantify the relative notion of short and long runs and their probability distribution. This might give us a key to quantify the expected parallel speed-up. This can be deduced from the sequential behavior of the algorithm, and more precisely from the proportion of long and short runs in the sequential runtime distribution.

In the following, we propose a probabilistic model to quantify the expected speed-up of multi-walk Las Vegas algorithms. This makes it possible to give a general formula for the speed-up, depending on the sequential behavior of the algorithm. Our model is related to order statistics, which is the statistics of sorted random draws, a rather new domain of statistics [START_REF] David | Order Statistics. Wiley series in probability and mathematical statistics[END_REF]. Indeed, explicit formulas have been given for several well-known distributions. Relying on an approximation of the sequential distribution, we compute the average speedup for the multi-walk extension. Experiments show that the prediction is quite sharp and opens the way for defining more accurate models and applying them to larger classes of algorithms.

Previous works [START_REF] Verhoeven | Parallel Local Search[END_REF] studied the case of a particular distribution for the sequential algorithm: the exponential distribution. This case is ideal and it yields a linear speed-up. Our model enables us to approximate Las Vegas algorithms by other types of distribution, such as a shifted exponential distribution or a lognormal distribution. In the last two cases the speed-up is no longer linear, but admits a finite limit when the number of processors goes toward infinity. We will see that these distributions fit experimental data for some problems.

The literature provides other probabilistic models for parallel Las Vegas algorithms. For parallelization schemes based on restart strategies, Luby et al. [START_REF] Luby | Optimal Speedup of Las Vegas Algorithms[END_REF] proposed an optimal universal strategy, which achieves the best speedups amongst the restart-based universal strategies. Although our work, which is not based on restarts, does not fit within this framework, it is related to it since it uses similar probabilistic ideas. In our case, a probabilistic model is used to model the efficiency of a parallel scheme, without modifying the algorithm. In their case, probabilistic tools are used to find an optimal restart scheme for a base algorithm. With the same idea of restart-based parallelization of Las Vegas algorithms, [START_REF] Shylo | Restart Strategies in Optimization: Parallel and Serial Cases[END_REF] goes a step further and provides a more detailed model. It clarifies in particular the cases of super-linear speed-ups, which are due, in this framework, to inefficient (or not well-tuned) sequential algorithms. Interestingly, this article also experiments with log-normally distributed sequential algorithms, which confirms our hypothesis that not only exponential distributions have to be investigated (as assumed by [START_REF] Aiex | Probability distribution of solution time in grasp: An experimental investigation[END_REF][START_REF] Aiex | TTT Plots: A Perl Program to Create Time-To-Target Plots[END_REF]).

Probabilistic Model

Randomized algorithms are stochastic processes. Their behavior (output, running time...) is non-deterministic and varies according to a probability distribution, which may or may not be known. For instance, Local Search algorithms include several random components: choice of an initial configuration, choice of a move among several candidates, plateau mechanism, random restart, etc. A complete algorithm may also have a stochastic behavior when they include random components, such as random heuristics, restarts with randomization, etc.

In the following, we first define the class of algorithms that our model encompasses. We then present our probabilistic model, considering the computation time of an algorithm (whatever it is) as a random variable, and using elements of probability theory to study its multi-walk parallel version.

Parallel Las Vegas Algorithms

The notion of Las Vegas algorithm encompasses a wide range of combinatorial solvers. We borrow the following definition from [START_REF] Hoos | Stochastic Local Search: Foundations and Applications[END_REF], Chapter 4.

Definition 1 (Las Vegas Algorithm) An algorithm A for a problem class Π is a (generalized) Las Vegas algorithm if and only if it has the following properties:

1. If for a given problem instance π ∈ Π, algorithm A terminates returning a solution s, s is guaranteed to be a correct solution of π. 2. For any given instance π ∈ Π, the runtime of A applied to π is a random variable.

This definition includes algorithms which are not guaranteed to return a solution. However in practice, we will only consider terminating Las Vegas algorithms, such as Local Search algorithms which always terminate if run for an unbounded time.

Let us now formally define a parallel multi-walk Las Vegas algorithm. It is often more convenient to consider distributions with values in R because it makes calculations easier. For the same reason, although f Y is defined in N, we will use its natural extension to R. This step is merely technical, and the probability distribution in N can be retrieved from the distribution in R.

The expectation of the computation is then defined by the standard formula for real-valued distributions:

E [Y ] = ∞ 0 tf Y (t)dt.
This formula is the extension to R of the classical expectation formula in the case of integer distributions ( ∞ 0 tf Y (t)). Assume that the base algorithm is concurrently run in parallel on n processors. In other words, over each processor the running process is a copy of the algorithm with different initial random seed. The first process that finds a solution then kills all others and the algorithm terminates. The i-th process corresponds to a draw of a random variable X i , following distribution f Y . The variables X i are thus independently and identically distributed (i.i.d.). The computation time of the whole parallel process is also a random variable, let's call it Z (n) , with a distribution f Z (n) that depends on both n and f Y . Since all the X i are i.i.d., the cumulative distribution F Z (n) can be computed as follows:

F Z (n) = P[Z (n) ≤ x] by definition = P[∃i ∈ {1...n}, X i ≤ x] because of the multiwalk rule = 1 -P[∀i ∈ {1...n}, X i > x] probability formula for the negation = 1 - n i=1 P[X i > x] because the random variables are i.i.d. = 1 -(1 -F Y (x))
n by definition which leads to:

f Z (n) = (1 -(1 -F Y ) n ) = nf Y (1 -F Y ) n-1
Thus, knowing the distribution for the base algorithm Y , one can calculate the distribution for Z (n) . In the general case, the formula shows that the parallel algorithm favors short runs, by killing the slower processes. Thus, we can expect that the distribution of Z (n) moves toward the origin, and is more peaked. As an example, Fig. 1 shows this phenomenon when the base algorithm admits a Gaussian distribution.

Expectation and Speed-up

The model described above gives the probability distribution of a parallelized version of any randomized algorithm. We can now calculate the expectation for the parallel process with the following relation:

E Z (n) = ∞ 0 tf Z (n) (t) dt = n ∞ 0 tf Y (t)(1 -F Y (t)) n-1 dt
Unfortunately, this does not lead to a general formula for E Z (n) . In the following, we will study it for different specific distributions.

To measure the gain obtained by parallelizing the algorithm on n processors, we will study the speed-up G n defined as:

G n = E [Y ] /E Z (n)
Again, no general formula can be computed and the expression of the speed-up depends on the distribution of Y . However, it is worth noting that our computation of the speed-up is related to order statistics, see [START_REF] David | Order Statistics. Wiley series in probability and mathematical statistics[END_REF] for a detailed presentation. Order statistics are the statistics of sorted random draws. For instance, the first order statistics of a distribution is its minimal value, and the k th order statistic is its k th -smallest value. For predicting the speed-up of a multi-walk Las Vegas algorithm on n processors, we are indeed interested in computing the expectation of the distribution of the minimum among n draws. As the above formula suggests, this may lead to heavy calculations, but recent studies such as [START_REF] Nadarajah | Explicit Expressions for Moments of Order Statistics[END_REF] give explicit formulas for this quantity for several classical probability distributions. Except in the case of the exponential distribution, detailed below, the formulas given for the minimum order statistics are rather complicated. In some cases, a symbolic computation may not even succeed in a reasonable amount of time. When this happens, we will perform a first step of symbolic computation, then another step of numeric integration to obtain the numerical value for the speed-up.

Case of an Exponential Distribution

Assume that Y has a shifted exponential distribution, as it has been suggested by [START_REF] Aiex | Probability distribution of solution time in grasp: An experimental investigation[END_REF][START_REF] Aiex | TTT Plots: A Perl Program to Create Time-To-Target Plots[END_REF].

f Y (t) = 0 if t ≤ x 0 λe -λ(t-x0) if t > x 0 F Y (t) = 0 if t ≤ x 0 1 -e -λ(t-x0) if t > x 0 E [Y ] = x 0 + 1/λ
Then the formula of Section 3.2 can be symbolically computed by hand:

f Z (n) (t) = 0 if t ≤ x 0 nλe -nλ(t-x0) if t > x 0 F Z (n) (t) = 0 if t ≤ x 0 1 -e -nλ(t-x0) if t > x 0
The intuitive observation of Section 3.2 is easily seen on the expression of the parallel distribution, which has an initial value multiplied by n but an exponential factor decreasing n-times faster, as shown on the curves of Fig. 2.

And in this case, one can symbolically compute both the expectation and speed-up for Z (n) : Fig. 3 shows the evolution of the speed-up when the number of processors increases. With such a rather simple formula for the speed-up, it is worth studying what happens when the number of processors n goes to infinity. Depending on the chosen algorithm, If x 0 = 0, then the expectation goes to 0 and the speed-up is equal to n. This case has already been studied by [START_REF] Verhoeven | Parallel Local Search[END_REF]. For x 0 > 0, the speed-up admits a finite limit which is

E Z (n) = nλ ∞ x0 te -nλ(t-x0) dt = x 0 + 1 nλ G n = x 0 + 1 λ x 0 + 1
x0+ 1 λ x0 = 1 + 1
x0λ . Yet, this limit may be reached slowly, and depends on the values of x 0 and λ. From the previous formula we observe that the closer x 0 is to zero, the better is the speedup. Another interesting value is the coefficient of the tangent at the origin, which approximates the speed-up for a small number of processors. In case of an exponential, it is (x 0 * λ + 1). The higher x 0 and λ, the bigger is the speed-up at the beginning. In the following, we will see that, depending on the combinations of x 0 and λ, different behaviors can be observed.

Case of a Lognormal Distribution

Other distributions can be considered, depending on the behavior of the base algorithm. We will study the case of a lognormal distribution, which is the log of a Gaussian distribution. The lognormal distribution appears in several of our experiments in Section 4 and 5. The lognormal distribution has two parameters, the mean µ and the standard deviation σ. In the same way as the shifted exponential distribution, we shift the distribution so that it starts at a given parameter x 0 . Formally, a (shifted) lognormal distribution is defined as:

f Y (t) = 0 if t < x 0 Φ(log(t-x0)) t-x0 if t > x 0
where

Φ(t) = 1 √ 2 * π e -t 2 /2
. The mean, variance and median are known and equal to e µ+ σ 2

2 , e 2µ+σ 2 (e σ 2 -1) and e µ respectively. Fig. 4 depicts lognormal distributions of Z (n) , for several n. The computations for the distribution of Z (n) and the theoretical speed-up are the same as given in Section 3.3. The computation of the expectation for the moments of order statistics of a lognormal distribution can be found in [START_REF] Nadarajah | Explicit Expressions for Moments of Order Statistics[END_REF], which gives an explicit formula with only a numerical integration step. We only recall from [START_REF] Nadarajah | Explicit Expressions for Moments of Order Statistics[END_REF] this formula for the first moment (expectation) of the first order statistics (minimum distribution), shifted by x 0 , which we are interested in.

We will need Lauricella functions defined by:

F (n) A (a; b 1 , ..., b n ; c 1 , ..., c n ; x 1 , ..., x n )) = ∞ m1=0 ... ∞ mn=0 (a) m1+...+mn (b 1 ) m1 ...(b n ) mn (c 1 ) m1 ...(c n ) mn x m1 1 ...x mn n m 1 !...m n !
Then one has:

E Z (n) = x 0 + ne 1/2 n-1 l=0 n -1 l - 1 2 l r p=0 l p - 2 √ π p * E N √ 2 -1 2 √ 2 p F (p) A ( 1 2 , ..., 1 2 ; 3 2 , ..., 3
2 ; - (N √ 2 -1) 2 8 , ..., - (N √ 2 -1) 2 8 )
where N is a standard normal random variable. In addition, [START_REF] Nadarajah | Explicit Expressions for Moments of Order Statistics[END_REF] provides pointers to routines in Mathematica to compute the terms of this formula with only one step of numerical integration. In practice, we will get our numerical results with an earlier step of numerical integration in Mathematica, with a good accuracy. This allows us to draw the general shape of the speed-up, an example being given on Fig. 5. Due to the numerical integration step, which requires numerical values for the number of processors n, we restrict the computation to integer values of n.

Methodology

In this section we detail how to perform the prediction of the parallel runtime and speed-up with respect to sequential execution from the analysis of the sequential runtime distribution. On each problem, the sequential benchmark gives observations of the distribution of the algorithm runtime f Y . Yet, the exact distribution is still unknown. It can be any real distribution, not even a classical one. In the following, we will rely on the assumption that Y is distributed with a known parametric distribution. We perform a statistical test, called Kolmogorov-Smirnov test, on the hypothesis H 0 that the collected observations correspond to a theoretical distribution. Assuming H 0 , the test first computes the probability that the distance between the collected data and the theoretical distribution does not significantly differ from its theoretical value. This probability is called the p-value.

Then, the p-value is compared to a fixed threshold (usually 0.05). If it is smaller, one rejects H 0 . For us, it means that the observations do not correspond to the theoretical distribution. If the p-value is high, we will consider that the distribution of Y is approximated by the theoretical one. Note that the Kolmogorov-Smirnov test is a statistical test, which in no way proves that Y follows the distribution. However, it measures how well the observations fit a theoretical curve and, as it will be seen in the following, it is accurate enough for our purpose.

Our benchmarks appear to fit with two distributions: the exponential distribution, as suggested by [START_REF]Statistical Methods in Experimental Physics[END_REF], and the lognormal distribution. We have also performed the Kolmogorov-Smirnov test on other distributions (e.g., Gaussian and Lévy), but obtained negative results w.r.t. the experimental benchmarks, thus we do not include them in the sequel. For each problem, we need to estimate the value of the parameters of the distribution, which is done on a case by case basis. Once we have an estimated distribution for the runtimes of Y , it becomes possible to compute the expectation of the parallel runtimes and the speed-up using the formulas of Section 3.3.

In the following, all the analyses are done on the number of iterations, because they are more likely to be unbiased, and all the mathematical computations are done with Mathematica [START_REF] Wolfram | The Mathematica Book, 5th Edition[END_REF].

Application to Constraint-based Local Search

Recently, the application of Local Search to solve combinatorial problems has attracted the attention of researcher [START_REF] Codognet | Yet Another Local Search Method for Constraint Solving[END_REF][START_REF] Van Hentenryck | Constraint-Based Local Search[END_REF], as it can tackle Constraint Satisfaction Problems (CSP) instances far beyond the reach of classical propagationbased constraint solvers. A generic, domain-independent Constraint-Based Local Search method, named Adaptive Search, has been proposed by [START_REF] Codognet | Yet Another Local Search Method for Constraint Solving[END_REF][START_REF] Codognet | An Efficient Library for Solving CSP with Local Search[END_REF]. This meta-heuristic takes advantage of the structure of the problem in terms of constraints and variables and can guide the search more precisely than a single global cost function to optimize, such as for instance the number of violated constraints. The algorithm also uses a short-term adaptive memory in the spirit of Tabu Search in order to prevent stagnation in local minima and loops. The main ideas of adaptive search can be summarized as follows:

to consider for each constraint a heuristic function that is able to compute an approximated degree of satisfaction of the goals (the current error on the constraint); to aggregate constraints on each variable and project the error on variables thus trying to repair the worst variable with the most promising value; to keep a short-term memory of bad configurations to avoid looping (i.e., some sort of tabu list) together with a reset mechanism.

We used for our experiments the reference implementation of Adaptive Search (AS) which has been developed as a framework library in C is and available as a freeware at the URL: http://cri-dist.univ-paris1.fr/diaz/adaptive/

Experimental Protocol and Benchmarks

We detail here the performance and speed-ups obtained with both sequential and parallel multi-walk Adaptive Search implementations. We have chosen to test this method on a hard combinatorial problem abstracted from radar and sonar applications (Costas Array) and two problems from the CSPLib benchmark library2 : All-Interval Series (prob007 in CSPLib), and the Magic-Square problem (prob019 in CSPLib). We refer the reader to [START_REF] Truchet | Prediction of Parallel Speed-ups for Las Vegas Algorithms[END_REF] for a complete description of the benchmarks and the local search framework.

We run our benchmarks in a sequential manner a certain number of times to evaluate the quality of the proposed model. Sequential experiments, as well as parallel experiments, have been done on the Griffon cluster of the Grid'5000 platform [START_REF] Bolze | Grid'5000: A Large Scale And Highly Reconfigurable Experimental Grid Testbed[END_REF], the French national grid for research, which contains 8,596 processors deployed on 11 sites distributed over France. In our experiments, we used the Griffon cluster at Nancy, composed of 184 Intel Xeon L5420 (Quadcore, 2.5GHz, 12MB of L2-cache, bus frequency at 1333MHz), thus with a total of 736 processes available giving a peak performance of 7.36TFlops. The following table ( One can see from Table 1 that runtimes and numbers of iterations are spread over a large interval for each benchmark, illustrating the stochasticity of the algorithm. Depending on the benchmark, there is a ratio of a few thousands times between the minimum and the maximum runtimes.

Table 2 presents the speed-up for the runtime and the number of iterations up to 256 processors for the executions of large benchmarks: Magic-Square (instance of size 200×200), All-Interval (instance of size 700) and Costas Array (instance of size 21). The same code has been ported and executed, timings are given in seconds and are the average of 50 runs. The speed-up for a given parallel execution is calculated against the mean performance of its sequential version as follows: Speed-Up = Mean(Solver on 1 processor) Mean(Solver on N processors) One can notice there is no significant difference between speed-ups in CPU-time and in number of iterations, therefore we will prefer as a time measure the number of iterations, which has the good property of not being machine-dependent. Similar speed-ups have been achieved on other parallel machines [START_REF] Caniou | Performance Analysis of Parallel Constraint-Based Local Search[END_REF]. For Magic-Square and All-Interval one can observe the stabilization point is not yet obtained for 256 processor, even if speed-ups do not increase as fast as the number of processors, i.e., are getting further away from linear speed-up. For the Costas Array Problem, our algorithm reaches linear or even super-linear speed-ups for up to 256 processors. Actually, it scales linearly far beyond this point, i.e., at least up to 8,192 processors, as reported in [24]. Speed-ups of the average runtime for Magic-Square and All-Interval look similar, but their actual runtime behaviors are different, as will be seen in the next section.

The All-Interval Series Problem

The analysis is done on 720 runs of the Adaptive Search algorithm on the instance of All-Interval series for 700 notes. The sequence of observations is written AI 700 in the following.

We test the hypothesis that the observations admit a shifted exponential distribution as introduced in Section 3.4. The first step consists in estimating the parameters of the distribution, which for a shifted exponential are the value of the shift x 0 and λ. 3 We take for x 0 the minimum observed value, x 0 = 1217. The exponential parameter is estimated thanks to the following relation: for a non-shifted exponential distribution, the expectation is 1/λ. Thus we take λ = 1/(mean(AI 700) -x 0 ), which gives λ = 9.15956 • 10 -6 .

We then run the Kolmogorov-Smirnov test on the shifted exponential distribution with these values of x 0 and λ, which answers positively (computed p-value: 0.77435). We thus admit the hypothesis that AI 700 fits this shifted exponential distribution. As an illustration, Fig. 6(a) shows the normalized histogram of the observed runtimes and the theoretical distribution.

It is then possible to symbolically compute the speed-up that can be expected with the multi-walk parallel scheme. We use the formulas of Section 3.4 with the estimated parameters and obtain a theoretical expression for the 200 000 400 000 600 000 800 000 runtime (b) Predicted speedup with its limit yellow (dashed), and the ideal speedup in pink (dashed) Fig. 6 Analysis of AI 700 speed-up. This allows us to calculate its value for different numbers of processors.

The results are given in Fig. 6(b). With this approximated distribution, the limit of the speed-up when the number of processors goes to infinity is 90.7087. One can see that, with 256 processors, the curve has not reached its limit, but comes close. Thus, the speed-up for this instance of All-Interval appears significantly less than linear (i.e., less than the number of processors).

The Magic-Square Series Problem

For the Magic-Square problem with N = 200, the observations are the number of iterations on 662 runs, with a minimum of x 0 = 6210. The Kolmogorov-Smirnov test on a shifted exponential distribution fails, but we obtain a positive result with a lognormal distribution, with µ = 12.0275 and σ = 1.3398, shifted to x 0 . These parameters have been estimated with the use of the Mathematica software. As an illustration, Fig. 7(a) shows the observations and the theoretical estimated distribution. The speed-up can be computed by integrating the minimum distribution with numerical integration techniques. The results are presented in Fig. 7(b). We can observe that the speed-up grows very fast at the origin, which can be explained by the high peak of the lognormal distribution with these parameters. Again, the speed-up is computed with a numerical integration step, and we only draw the curve for integer values of n. In this case again, the speed-up is significantly less than linear from 50 processors onwards, and the limit of the speed-up when the number of processors goes to infinity is about 71.5.

The Costas Array Problem

The same analysis is done for the runs of the AS algorithm on the Costas Array problem with N = 21. The observations are taken from the benchmark with 638 runs. The sequence of observations is written Costas 21.

This benchmark has an interesting property: the observed minimum, 3.2 • 10 5 is negligible compared to its mean (1.8 • 10 8 ). Thus, we estimate x 0 = 0 and perform a Kolmogorov-Smirnov test for a (non-shifted) exponential distribution, with λ = 1/mean(Costas 21) = 5.4 • 10 -9 . The test is positive for this exponential distribution, with a p-value of 0.751915. Fig. 8(a) shows the estimated distribution compared to the observations. The computation of the theoretical speed-up is then done in the same way as for AI 700. Yet, in this case, the observed minimum for x 0 is so small that we can approximate the observations with a non-shifted distribution, thus the predicted speed-up is strictly linear, as shown in Section 3.4. The results are given on Fig. 8(b). This explains that one may observe linear speed-up when parallelizing Costas Array.

Analysis and Discussion

Table 3 presents the comparison between the predicted and the experimental speed-ups. We can see that the accuracy of the prediction is very good up to 64 parallel processors and then the divergence is limited even for 256 parallel processors.

For the MS 200 problem, the experimental speed-up and the predicted one are almost identical up to 128 processors and diverging by 10% for 256 processors. For the AI 700 problem, the experimental speed-up is below the predicted one by a maximum of 30% for 128 and 256 processors. For the Costas 21 problem, the experimental speed-up is above the predicted one by 15% for 128 and 256 processors. It is worth noting that our model approximates the behaviors of experimental results very closely, as shown by the predicted speed-ups matching closely the real ones. Moreover we can see that on the three benchmark programs, we needed to use three different types of distribution (exponential, shifted exponential and lognormal), in order to approximate the experimental data most closely. This shows that our model is quite general and can accommodate different types of parallel behaviors.

A quite interesting behavior is exhibited by the Costas 21 problem. Our model predicts a linear speed-up, up to 10,000 processors and beyond, and the experimental data gathered for this paper confirms this linear speed-up up to 256 processors. Would it scale up with a larger number of processors? Indeed we did such an experiment up to 8,192 processors on the JUGENE IBM Bluegene/P at the Jülich Supercomputing Center in Germany (with a total 294,912 processors), and reported it in [24], the speed-up is linear up to 8,192 processors, thus showing an excellent fit between the prediction model and real data

Application to SAT Local Search

Let us now look at a different problem domain (SAT -the Satisfiability Problem for Boolean formulas) and different local search solvers.

In this section we focus our attention on two well-known problem families: random and crafted instances 4 . Moreover, we consider the two best Local Search solvers from the 2012 SAT competition5 : CCASAT [START_REF] Cai | CCASAT: Solver Description[END_REF] and Sparrow [START_REF] Balint | Improving Stochastic Local Search for SAT with a New Probability Distribution[END_REF]. Both solvers were used with their default parameters and with a timeout of 3 hours for each experiment. All the experiments were performed on the Grid'5000 platform, the French national grid for research. We used a 44-node cluster with 24 cores (2 AMD Opteron 6164 HE processors at 1.7 Ghz) and 44 GB of RAM per node.

We performed experiments with two well-known problem families of instances coming from the SAT'11 competition: random and crafted. In particular, we used 10 random instances (6 around the phase transition and 4 outside the phase transition) and 10 crafted instances. Complete details of the selected instances for the evaluation is available at [START_REF] Arbelaez | Using Sequential Runtime Distributions for the Parallel Speedup Prediction of SAT Local Search[END_REF], hereafter we denote random instances as Rand-[1 to 10] and crafted instances as crafted-[1 to 9].

In order to obtain the empirical data for the theoretical distribution (predicted by our model from the sequential runtime distribution), we performed 500 runs of the sequential algorithm. The Mathematica software [START_REF] Wolfram | The Mathematica Book, 5th Edition[END_REF], version 8.0, was used to estimate the parameters of the theoretical distributions and to integrate numerically the formulas of the lognormal distribution. In order to evaluate the accuracy of the learned statistical model, we performed 50 runs of the multi-walk parallel algorithms.

Analysis of the Parallel Speed-ups

Let us now present the empirical and predicted results for random and crafted instances. We will not perform here again a detailed analysis as in Section 3.6, but rather report directly the computed p-value when comparing the empirical runtime distribution to either a shifted exponential distribution and a lognormal distribution.

Random instances

We start our analysis with Table 4, which presents initial statistics for the sequential version of Sparrow and CCASAT. We present the minimum, maximum, and mean runtime values, as well as the outcome of the Kolmogorov-Smirnov (KS) test for two types of distributions: shifted exponential and lognormal. In the following tables, bold numbers indicate the distribution chosen to predict the performance of a given solver. As it can be observed for a few instances, the best distribution does not pass the KS test, but we still use the distribution with the best KS value for predicting the runtime of the parallel solver. . The results presented in this table are consistent with the results of the 2012 SAT competition (random category) where CCASAT greatly outperformed Sparrow. For this set of instances, we choose the shifted exponential distribution in lieu of the exponential distribution as the Min runtime value for the reference solvers is not negligible compared to its mean value across 500 executions (about 100 times smaller in the best case).

As can be seen from the table, both solvers report a tendency which indicates that the empirical data for instances around the phase transition are better approximated by a lognormal distribution; all these instances pass the KS test with a confidence level (p-value) above 0.05, except for Sparrow on rand-4. Alternatively when the KS test is not enough to decide the most appropriate distribution, we use the Q-Q plot to decide whether a given distribution fits the empirical data. The closer the empirical data is to the reference line (theoretical distribution), the better the fit. Fig. 9 shows the Q-Q plot for Sparrow to solve rand-7, unlike Fig. 9(b) where a clear deviation between the empirical data and the lognormal distribution is observed, Fig. 9(a) shows that the empirical distribution is very close to the shifted exponential distribution. We have performed the same analysis for the remaining random instances and observed a similar behaviour between the empirical data and the reference distributions It is well recognized in the SAT community that k-SAT instances around the phase transition (i.e., No of clauses/No of Variables = 4.2) are known to be difficult [START_REF] Gent | The SAT Phase Transition[END_REF]. Taking this into consideration, for instances outside the phase transition, Sparrow reports enough statistical evidence to infer that the shifted exponential distribution fits better the empirical data. For CCASAT, 3 out of 4 instances outside the phase transition are better characterized with a lognormal distribution and the remaining instance pass the KS test for the shifted exponential distribution. We now look at the parallel performance of the solvers. Table 5 shows the empirical and predicted runtime for both Sparrow and CCASAT on all instances using 48, 96, 192, and 384 processors. Detailed tables with the speedup for Sparrow and CCASAT are available in [START_REF] Arbelaez | Using Sequential Runtime Distributions for the Parallel Speedup Prediction of SAT Local Search[END_REF]. Summing up both solver report the same tendency as the empirical data. Furthermore, the speed-up factor of the references far from linear (ideal), a phenomenon well described by the predicted model. Fig. 10 depicts observations and the theoretical estimated distribution, and Fig. 11 shows a performance summary of the reference solvers to tackle an instance on the phase transition (rand-4) and another instance outside the phase transition (rand-7). The y-axis gives the probability (P[Y ≤ x]) of finding a solution in a time less or equal to x and the x-axis gives the runtime in seconds. From now on, in all figures 'Emp' stands for Empirical distribution, 'LN' stands for lognormal distribution, and 'SExp' stands for shifted exponential distribution. As expected CCASAT dominates the performance on one processor. For example to solve rand-4, CCASAT reports P[Y ≤ 16-mins] ≈ 1.0, while Sparrow reports P[Y ≤ 16-mins] ≈ 0.75. Interestingly, we observe that for CCASAT increasing the number of processors does not significantly improve the solving time. Consequently, Sparrow becomes more effective for a large number of processors. Therefore, Fig. 11(b) and 11(d) show that Sparrow is better than CCASAT when using 384 processors. Interestingly, the same pattern is observed for other random instances (see Table 5). Finally, it can also be observed that random instances around the phase transition exhibit a lower speed-up factor than the remaining random instances. For instance, the best empirical speed-up factor obtained for instances in the phase transition is 7.0 for Sparrow and 3.4 for CCASAT; and the best speed-up factor obtained for instances outside the phase transition is 114.2 for Sparrow and 50.5 for CCASAT.

Crafted instances

Let us switch our attention now to crafted instances, for which we have to treat CCASAT and Sparrow differently. For CCASAT, we were unable to find a theoretical distribution which fits the empirical data. It should be also noticed that CCASAT has mainly been designed and tuned to handle random instances. Let us look for instance at Fig. 13(a), which depicts the cumulative runtime distribution of CCASAT to solve Crafted-1 using the two reference distributions detailed in this paper (lognormal and exponential) and two extra distributions (Weibull and beta-prime). None of the theoretical distributions seems to be a good approximation of the empirical data. More precisely, the KS test reported a p-value of 2.7•10 -7 (lognormal); 7.0•10 -24 (exponential); 2.4•10 -6 (Weibull); and 6.9•10 -15 (beta-prime). Therefore, none of the theoretical distributions pass the KS test with a high-enough p-value. We also experimented with other instances and observed a similar behavior. Fig. 12(b) and 12(b) display the Q-Q plots between the empirical distribution and two theoretical distribution: exponential and log-normal. Both figures show that the data is neither exponentially nor log-normally distributed.

For Sparrow on all crafted instances, the KS test shows a much better pvalue for the exponential distribution than for the lognormal one, see Table 6. The confidence level is quite high for the instances Crafted-2,-3,-4,-5,-8,-9, with p-value up to 0.97, while the p-value is between 0.01 and 0.02 for Crafted-1,-6,-7. Also, as the minimum runtime is much smaller than the mean (at least Table 6 Sequential performance of Sparrow on crafted instances 300 times smaller), we can approximate the empirical data by a non-shifted exponential distribution [START_REF] Truchet | Prediction of Parallel Speed-ups for Las Vegas Algorithms[END_REF].

As can be seen in Table 7 the multi-walk parallel approach scales well for Sparrow on crafted instances as the number of processors increases. Indeed a nearly linear speed-up is obtained for nearly all the instances. As expected, the speed-up predicted by our model is optimal, and this result is consistent with those obtained in [33].

Analysis and Discussion

Several works have been devoted to the experimental study of parallel multiwalk extensions of Local Search algorithms [START_REF] Arbelaez | Massivelly Parallel Local Search for SAT[END_REF][START_REF] Arbelaez | From Sequential to Parallel Local Search for SAT[END_REF][START_REF] Hoos | Stochastic Local Search: Foundations and Applications[END_REF], but we presented here the first approach (to our knowledge) which applies order statistics in order to predict the parallel performance of Local Search algorithms for SAT. As in Section 4.2, the lognormal or shifted-exponential distributions need to be considered for the runtime distribution, and not only the exponential distribution as suggested by most of the literature.

Interestingly, the phase transition point of SAT instances also seems to have important consequences in the parallel performance of Local Search al- gorithms. For Sparrow at least, which is the solver with an overall best speedup factor, instances around the phase transition region are lognormally distributed, while instances outside the phase transition are shifted-exponentially distributed. As for Constraint-Based Local Search, the probability of returning a solution in no iterations is not zero in theory for SAT Local Search. How-ever, in practice a minimum number of steps is in general required to reach a solution cf. [33,[START_REF] Ribeiro | 08/17) Exploiting Run Time Distributions to Compare Sequential and Parallel Stochastic Local Search Algorithms[END_REF] for the sequential case. Therefore, experimental data may be better approximated by a shifted distribution with x 0 > 0: this is the case for the random instances. This leads to a non-linear speed-up with a finite limit, even in the case of an exponential distribution. Indeed, the experimental speed-up for both CCASAT and Sparrow on random instances is far from linear. On the contrary, Sparrow on crafted instances has a linear speed-up which could be explained by the fact that the minimal runtime is negligible w.r.t. the mean time (i.e., 1/λ for an exponential distribution). Therefore, the statistical test succeeds for x 0 0. This suggests that, in general, the comparison between the minimal time and the mean time is a key element for the study of the parallel behavior.

We do not discard that other parameters for the reference solvers would lead to other theoretical distributions (e.g., exponential distribution for random instances). In [START_REF] Kroc | An Empirical Study of Optimal Noise and Runtime Distributions in Local Search[END_REF] the authors showed that a well-tuned version of Walk-SAT is exponentially distributed for instances in the phase transition region. However, we experimented by increasing the ps (smoothing probability) parameter of Sparrow and still obtained the same theoretical distribution. In addition, when ps is too high the solver was unable to solve the instances within the 3 hour time limit. Unfortunately, CCASAT is only available in binary form, and it is not possible to experiment with other parameters for the solver.

We expect this work to have significant implications in the area of automatic parameter tuning to devise scalable Local Search algorithms. Currently, most parameter tuning tools (e.g., [START_REF] Hutter | ParamILS: An Automatic Algorithm Configuration Framework[END_REF][START_REF] Ansótegui | A Gender-Based Genetic Algorithm for the Automatic Configuration of Algorithms[END_REF]) are designed to improve the expected mean (or median) runtime, however as observed in this paper, unless the algorithms exhibit a non-shifted exponential distribution, their parallel performance is far from linear and varies from algorithm to algorithm.

Application to Propagation-based Constraint Solving

In the previous two sections we have been studying the parallel performance of incomplete algorithms based on Local Search. In this section, we shift our focus to complete algorithms to tackle CSPs. Broadly speaking, complete algorithms employ a combination of domain reduction and tree-based search method, in which at the root node all variables are associated with a given domain. At each step of the search a value is assigned to some variable. Each assignment is preceded by a look-ahead process called constraint propagation which reduces the domains of the variables. The search process is stopped when a solution has been obtained or when the complete tree has been explored. In particular, we consider wdeg [START_REF] Boussemart | Boosting Systematic Search by Weighting Constraints[END_REF] and min-dom [START_REF] Haralick | Increasing Tree Search Efficiency for Constraint Satisfaction Problems[END_REF] state-of-the-art variable selection heuristics that focus on difficult parts of the search.

Experimental Protocol and Benchmarks

In order to perform the following experimentation we equipped Gecode6 with the multi-walk framework. We also used a randomized version of the solver. At each node, we use wdeg and min-dom to select the most appropriate variable. The former selects the variable with the largest weighted degree7 value is selected (breaking ties at random), the latter selects the variables with the smallest valid domain at a given state of the search. The value for the chosen variable is selected uniformly at random. It is highly recognized that the runtime distribution for tree-based search methods for a large number of problems is heavy-tailed [START_REF] Gomes | Heavy-Tailed Phenomena in Satisfiability and Constraint Satisfaction Problems[END_REF] and adding restarts (i.e., stopping the search after a given number of backtracks, and re-initializing it) helps to alleviate the heavy-tailed phenomenon. In this paper, we use a fixed restart strategy where the restart cutoff is 250 backtracks. For each problem instance we compute the sequential performance using wdeg and min-dom, and use the heuristic with better performance (in sequential settings) to evaluate the model described above to estimate the performance of the solver in parallel.

In these experiments, we use the same computer settings as for the previous section and considered the same baseline problems as for the adaptive search framework in Section 4, that is, Magic-Square, All-Interval, and Costas Array. For Magic-Square and All-Interval we used the model provided in the Gecode distribution (version 4.2.1), and the Gecode version of Costas Array is available at http://www.hakank.org/gecode/costas_array.cpp. It is worth noting that the All-Interval model for Gecode allows the trivial solution (0, n-1, n-2, n-3, . . .), and such a solution is obtained without backtracking. Therefore, we removed the trivial solution in the model by only allowing positive numbers for the first variable in the list, i.e., adding the constraint X 0 > 0.

Analysis of the Parallel Speed-ups

Table 8 shows the runtime of the sequential version. We present the minimum, maximum, mean, and the p-value for the reference line distributions of the sequential runs for Magic-Square 15, All-Interval 18, Costas Array 17. Highlighted numbers indicate the best heuristic w.r.t. Mean value, this heuristic will be used later on for the parallel evaluation of the solver. The Kolmogorov-Smirnov test suggest that the sequential execution for Magic-Square can be characterized using a shifted exponential distribution with a p-value of 0.21 (Magic-Square 15 with wdeg), All-Interval can be characterized using an exponential distribution with a p-value of 0.75 (All-Interval 18 with min-dom) 8 , and Costas Array can be characterized using a log-normal distribution with a p-value 0.04 (Costas Array 17 with min-dom). It is worth pointing out that neither Magic-Square nor All-Interval passed the KS test for the log-normal distribution, and Costas Array did not pass KS for the (shifted/and not shifted) exponential distribution. Once again, we use the formulas of Section 3.3 with λ=1/mean and x 0 is set to the minimal value of the empirical data, and used Mathematica to estimate the parameters of the log-normal distribution for Costas Array. Hereafter we use the indicated distributions for each algorithm to build the model and estimate the parallel performance in solver. Tables 9 reports the estimated runtime (top) and speed-up (bottom) of the multi-walk version of the propagation-based constraint solving algorithm. Interestingly, this simple parallel scheme exhibits very good speedup factors for the three benchmarks: super-linear for Magic-Square with 48, 96, cores, linear for nearly all experiments for All-Interval, and a factor 122 (w.r.t. the sequential solver) for Costas Array with 384 cores.

Analysis and Discussion

We recall that the goal of this section is not to compare incomplete methods based on Local Search and complete methods based on tree search. Instead, we intend to show the application of the proposed framework to a large variety of randomized techniques to solve combinatorial problems. Interestingly, even though the search strategy used in these experiments largely differs from the experiments presented in the two previous sections, we have observed that the empirical data can be characterised according to a theoretical distribution, ie shifted exponential (Magic-Square and All-Interval) and log-normal (Costas Array) distributions, and accurately estimate the runtime of the parallel algorithm.

The runtime prediction is very close to the empirical data for the three problems with gap between the prediction and the actual execution time (Table 8) of up 13.7% for Magic-Square with 384 cores, 14% for All-Interval with 92 cores, and 8% for Costas Array with 384 cores. Furthermore, the predictions of the speedup of the solvers are very encouraging, for the Magic-Square problem it can be observed that the gap between the actual and estimated speedup is up to 13% with 24 cores, for the All-Interval problem is up to 12% with 96 cores, and for Costas Array is up to 4% with 192 cores. Finally, we would like to point out that this methodology requires independently and identically distributed executions of the sequential algorithm, so that, parallel algorithms based on search-splitting techniques [START_REF] Bordeaux | Experiments with Massively Parallel Constraint Solving[END_REF][START_REF] Chu | Confidence-Based Work Stealing in Parallel Constraint Programming[END_REF][START_REF] Moisan | Parallel Discrepancy-Based Search[END_REF] and cooperative techniques [START_REF] Arbelaez | Massivelly Parallel Local Search for SAT[END_REF] do not fit within the proposed framework.

Finally, we would like to point out that we also experimented with the losing technique in sequential settings (see Table 8) and also found out very accurate predictions of up to 5% for Magic-Square with min-dom, up to 4% for All-Interval with wdeg, and up to 2% for Costas Array with wdeg.

Empirical Estimation

In this section we study a method that only uses the empirical sequential distribution to estimate the parallel runtime distribution. This method is much simpler than our model, as it does not involve many statistical tools. Our goal is to compare this empirical method to our model.

Let {k 1 , ..., k m } ∈ K be the set of m sequential executions of a given randomised algorithm to solve a given instance. We assume K to be sorted,

P[Z (n) ≤ t i ] = 1 -(1 -F Y (t i )) n Min Distribution, Section 3.2 = 1 -1 -i m n by definition
Figure 15 shows the assigned probabilities of each element in the sorted set of sequential observations, the y-axis shows the probability of selecting each k i ∈ K (x-axis) for the sequential algorithm with 1 processor,and the derived probabilities for the empirical parallel runtime distribution with 2, 4, 8, 16, and 32 processors. As expected, increasing the number of processors also increases the chances of selecting short runs.

Hereafter, we refer to prediction method for our model based on a theoretical distribution (e.g., lognormal distribution) as shown in Sections 4, 5, and 6, and we use the term estimation method to refer to the simpler method described here. Certainly, a key factor for the accuracy of both methods is the number of sequential observations. We are interested in the robustness of both methods w.r.t. the number of observations (m) and the number of cores (n). We thus show experiments where we vary the number of sequential observations to build the probability distribution model. Table 10 shows the average discrepancy across all instances with 48, 96, 192, and 384 cores 9 . We use a baseline of 500 instances for Sparrow (random and crafted), CCASAT and the Adaptive Search, and 300 instances for experiments with the Gecode solver and gradually decrease the number of observations m by considering 50%, 25%, and 12% of the sequential data. As it can be observed in the table, the empirical evaluation is a robust alternative when a large enough number of observations is available. In particular, we observe that the empirical method slightly outperforms our statistical model in 3 out of 4 scenarios when considering 100% of the sequential data. On the other hand, our statistical model reports better results in 11 out of the 12 scenarios when considering 50%, 25% and 12% of the data. Interestingly, the difference between the two prediction models is more significant as we increase the number of cores and reduce the number of observations. For instance, when using 12% of the data and 384 cores we observe that the discrepancy w.r.t. the actual data is 67.05% for empirical estimations and only 29.29% for the prediction model. In addition, summing up all the experiments (last row in Table 10), we observe that the quality of the prediction model is more stable and better than the estimation method. Fig. 16 shows a comparison using 100% of the sequential execution for the prediction, empirical 10 , and empirical estimation methods for Sparrow to solve crafted-1 with 48 and 192 processors 11 . As it can be observed in this figure the empirical estimation greatly matches the performance with 48 processors (Fig. 16(a)). However, the quality considerably degrades when using 192 processors (Fig. 16(b)), for instance, the empirical estimation suggests that the probability of finding a solution in 50 seconds for Sparrow with 192 processors is P[t ≤ 50] ≈ 60% vs. P[t ≤ 50] ≈ 90% for the actual performance in parallel settings, that is, about 30% discrepancy.

Finally, it can be observed that our prediction model, based in order statistics, closely matches the performance for an important number of processors. It is worth noting that the performance deterioration of the empirical estimation method is even higher when using a small set of sequential observations, as observed in Table 10. Furthermore, the prediction method also allows us to predict the limit of the speedup when the number of processors tends to infinity as shown in Section 3.4 for the exponential distribution. An example of the maximal speedup is depicted for the Magic-Square problem in Section 4.1.1.

The prediciton method thus ensures a more robust accuracy. The intuitive explanation is the following: as explained above, compared to the sequential distribution, the minimum distribution is concentrated close to the origin (or x 0 ). Thus, when it is approached with the empirical method, only the samples next to the origin count, no matter what the rest of the observations are: only a few samples really count in the model. On the contrary, our method finds a curve which best approximates the whole sequential distribution, leveraging the information at hand. Hence, it is more robust when the number of observations decreases. Notice also that the empirical method does not allow us to compute the theoretical limit for the speed-up. 10 For empirical we mean the empirical runtime distribution computed with the actual empirical results of the experiments with 48 and 192 processors. 11 We have performed the same analysis with other instances and observed a similar tendency, i.e., increasing the number of processors degrades the performance of the empirical estimated probability distribution

Conclusion and Future Work

We have proposed a theoretical model for predicting and analyzing the speedups of Las Vegas algorithms and applied it for Local Search methods in two different domains, Constraint-Based Local Search and SAT, and also for propagation-based constraint solving with random labeling heuristics. Interestingly, we have observed that, for the different algorithms and the variety of instances considered in this study, the runtime distribution can be characterized using two types of distributions: exponential (shifted and non-shifted) and lognormal.

It is worth noting that our model mimics the behaviors of the experimental results very closely, as shown by the predicted speed-ups matching closely the real ones. We showed that the parallel speed-ups predicted by our statistical model are accurate, matching the actual speed-ups very well up to several hundreds of processors.

However, a limitation of our approach is that, in practice, we need to be able to approximate the sequential distribution. In addition, this distribution must be one of the distributions for which the first order statistics is known, symbolically (as the exponential) or numerically (as the lognormal). Nevertheless, recent results in the field of order statistics give explicit formulas for a number of useful distributions: Gaussian, lognormal, gamma, beta. This provides a wide range of tools to analyze different behaviors.

Another interesting extension of this work would be to devise a method for predicting the speed-up from scratch, that is, without any knowledge on the algorithm distribution. Our observations suggest that the sequential runtime of both LS and complete search are well approximated by a small number of distributions. This could be intensively tested on a wider range of problems. In particular, it is important to know wether the sequential distribution of different instances of a given problem belong to the same family. Then we can devise a method for estimating the sequential distribution based on a limited number of observations, possibly on small instances, and then estimate the parallel speed-up for larger instances. This would allow us to predict if a simple multi-walk parallelization scheme, which does not imply to modify the algorithms, is likely to be efficient, or not. In the same spirit, further research includes investigating machine learning techniques to infer the distribution from a set of experiments as small as possible.

Aknowledgments

We acknowledge that some results in this paper have been achieved using the Grid'5000 experimental testbed, being developed under the INRIA ALADDIN development action with support from CNRS, RENATER and several universities as well as other funding bodies. We acknowledge that some other results in this paper have been achieved using the PRACE Research Infrastructure resource JUGENE based in Germany at Jülich Supercomputing Centre. The authors would like to thank the anonymous reviewers for their comments and suggestions that helped to improve the paper. We would like to thank the anonymous reviewers for suggesting to compare our method vs. the empirical method depicted in Section 7.

Fig. 1

 1 Fig.1Distribution of Z(n) , in the case where Y admits a Gaussian distribution (cut on R - and renormalized). The blue curve is Y . The distributions of Z (n) are in pink for n = 10, in yellow for n = 100 and in green for n = 1000.

Fig. 2 Fig. 3

 23 Fig.2For an exponential distribution, here in blue with x 0 = 100 and λ = 1/1000, simulations of the distribution of Z (n) for n = 2 (pink), n = 4 (yellow) and n = 8 (green).

Fig. 4 Fig. 5

 45 Fig.4For a lognormal distribution (in blue), with x 0 = 0, µ = 5 and σ = 1, simulations of the distribution of Z (n) for n = 2 (pink), n = 4 (yellow) and n = 8 (green).

Fig. 7

 7 Fig. 7 Analysis of MS 200

Fig. 8

 8 Fig. 8 Analysis of Costas 21

  Seq. vs. Log-normal dist.

Fig. 9 Q

 9 Fig. 9 Q-Q Plot Analysis of Sparrow on rand-7

Fig. 10

 10 Fig. 10 Histogram for the observed runtimes for CCASAT and Sparrow to solve rand-4 and rand-7

  Fig. 11 Performance summary to solve rand-4 and rand-7

  Seq. vs. Log-normal dist.

Fig

  Fig. 12 Q-Q Plot Analysis of CCASAT on crafted-1

Fig. 13

 13 Fig. 13 Performance summary on crafted-1

Fig. 14

 14 Fig. 14 Histogram for the observed runtimes for complete search

Fig. 16

 16 Fig. 16 Predicted, Empirical, and Empirical Estimated probability distributions of Sparrow to solve crafted-1

Table 1 )

 1 shows the minimum, mean, median, maximum and standard deviation among the runtimes (in seconds) and the number of iterations for our three benchmarks.

			MS 200	AI 700	Costas 21
		time	iterations	time iterations	time	iterations
	Min	5.5	6, 210	23.3	1, 217	6.6	321, 361
	Mean	382.0	443, 970	1, 354.0	110, 393	3, 744.4 183, 428, 617
	Median	126.3	164, 042	945.4	76, 242	2, 457.4 119, 667, 588
	Max	7, 441.6 7, 895, 872 10, 243.4	826, 871 19, 972.0 977, 709, 115
	Std Dev	873.0	933, 766	1, 363.4	111, 352	3, 655.5 179, 049, 696

Table 1

 1 Sequential executions in seconds and number of iterations

Table 2

 2 Parallel Speed-ups time and number of iterations.

Table 3

 3 Comparison: experimental and predicted speed-ups

	Instance		16	speed-up on k processors 32 64 128	256
	MS200	experimental predicted	16.6 15.94 22.04 28.28 34.26 22.2 29.9 34.3	45.0 39.7
	AI700	experimental predicted	12.8 13.7	20.2 23.8	29.3 37.8	37.3 53.3	48.0 67.2
	Costas21	experimental predicted	15.8 16.0	26.4 32.0	60.0 159.2 290.5 64.0 128.0 256.0

Table 4

 4 Performance of sequential algorithms on random instances

	Instance	Alg	Min	Max	Mean	p-value Shifted exp. dist. Lognormal dist.
	Rand-1	Sparrow CCASAT	98.0 103.6	4860.0 1340.1	793.9 458.5	5.6•10 -19 8.4•10 -55	0.76 0.96
	Rand-2	Sparrow CCASAT	91.8 108.8	5447.0 1007.5 1652.9 497.9	1.8•10 -20 5.1•10 -58	0.91 0.71
	Rand-3	Sparrow CCASAT 126.48 104.8	3693.7 1125.4	797.0 359.6	2.9•10 -26 2.3•10 -108	0.64 0.90
	Rand-4	Sparrow CCASAT	162.4 132.5	3037.8 980.9	781.5 382.0	1.6•10 -38 1.5•10 -117	0.03 0.92
	Rand-5	Sparrow CCASAT	164.0 158.6	7946.3 1177.9	952.4 403.1	1.8•10 -31 5.1•10 -134	0.16 0.20
	Rand-6	Sparrow CCASAT	142.0 142.9	4955.8 890.9	763.5 354.1	1.6•10 -31 9.6•10 -137	0.64 0.46
	Rand-7	Sparrow CCASAT	35.5 10637.4 3464.2 61.6 6419.1 1801.0	0.01 1.0•10 -5	1.5•10 -4 0.13
	Rand-8	Sparrow CCASAT	23.2 10738.0 3412.9 35.9 10443.7 2007.6	0.05 0.50	4.7•10 -4 0.03
	Rand-9	Sparrow CCASAT	6.8 18.1	5935.8 1028.2 2830.4 476.8	0.23 7.0•10 -3	3.7•10 -3 0.03
	Rand-10	Sparrow CCASAT	19.0 10800.0 1726.3 19.8 4854.5 758.4	0.65 7.6•10 -10	0.15 0.18

Table 5

 5 Runtimes for random instances up to 384 proc. (processors)

	Instance		Sparrow -Runtime on k proc. 48 96 192 384	CCASAT -Runtime on k proc. 48 96 192 384
	Rand-1	Actual Predicted 133.8 110.5 163.8 140.4 125.2 113.7 160.0 143.0 122.8 92.7 78.8 137.7 120.6 106.7	112.0 95.3
	Rand-2	Actual Predicted 183.5 152.8 129.2 110.6 153.4 134.7 119.6 213.2 191.4 166.2 142.5 186.8 169.3 159.3	142.8 107.1
	Rand-3	Actual Predicted 183.5 152.8 129.2 110.6 153.4 134.7 119.6 175.9 151.2 135.8 123.5 166.7 155.6 143.5	132.2 107.1
	Rand-4	Actual Predicted 175.7 149.5 128.9 112.4 170.6 155.9 143.5 202.3 179.2 159.5 141.8 193.1 176.0 169.4	158.7 132.8
	Rand-5	Actual Predicted 185.0 155.3 132.3 114.0 179.8 164.3 151.2 219.6 201.0 182.5 161.9 212.2 191.3 176.8	165.8 140.0
	Rand-6	Actual Predicted 158.3 133.6 114.4 185.5 167.1 150.3 137.5 190.9 179.3 168.4 99.1 160.6 147.0 135.4	153.4 125.6
	Rand-7	Actual Predicted 195.8 143.0 107.3 151.2 102.7 63.8	51.1 82.3 182.8 142.6 113.7 22.9 33.7 54.3	67.8 92.2
	Rand-8	Actual Predicted	126.6 93.8	81.9 58.5	51.1 40.8	30.9 131.8 32.0 76.9	83.9 56.4	64.8 46.1	39.7 41.0
	Rand-9	Actual Predicted	33.9 28.1	18.4 17.4	13.1 12.1	9.0 9.4	45.0 38.5	31.0 29.4	22.7 23.0	16.3 18.3
	Rand-10	Actual Predicted	63.4 54.6	48.9 36.8	40.7 27.9	30.9 113.8 23.4 105.6	94.7 85.3	72.9 70.2	54.2 58.6

  . 12 Q-Q Plot Analysis of CCASAT on crafted-1

	Instance	Alg	Min	Max	Mean	p-value Exp. dist. Lognormal dist.
	Crafted-1 Sparrow	9.9 10800.0	3440.3	0.02	1.2•10 -4
	Crafted-2 Sparrow	1.0 10800.0	2711.2	0.57	1.4•10 -4
	Crafted-3 Sparrow	8.7 10800.0	3432.7	0.14	1.1•10 -3
	Crafted-4 Sparrow	2.2 10800.0	2701.6	0.11	9.6•10 -3
	Crafted-5 Sparrow	4.1 10800.0	1564.1	0.95	9.2•10 -4
	Crafted-6 Sparrow	2.9 10800.0	3599.6	0.01	1.0•10 -5
	Crafted-7 Sparrow	4.4 10800.0	3598.7	0.01	7.8•10 -6
	Crafted-8 Sparrow	3.5	5456.0 972.046	0.67	0.17
	Crafted-9 Sparrow	1.9	7876.5 1298.24	0.97	7.0•10 -3

Table 7

 7 Parallel performance of Sparrow on crafted instances (proc. stands for processors)

	Instance		Runtime on k proc. 48 96 192	384	48	Speed-Up ok k proc. 96 192	384
	Crafted-1	Actual Predicted	97.7 43.7 19.1 71.6 35.8 17.9	9.8 35.1 8.9 48.0	78.6 179.6 349.8 96.0 192.0 384.0
	Crafted-2	Actual Predicted	67.8 36.4 17.5 56.4 28.2 14.1	7.2 39.9 7.0 48.0	74.4 154.7 375.2 96.0 192.0 384.0
	Crafted-3	Actual Predicted	94.8 49.3 23.2 11.9 36.1 71.5 35.7 17.8 8.9 48.0	69.6 147.6 286.1 96.0 192.0 384.0
	Crafted-4	Actual Predicted	87.5 42.0 17.3 56.2 28.1 14.0	9.7 30.8 7.0 48.0	64.2 155.4 277.8 96.0 192.0 384.0
	Crafted-5	Actual Predicted	33.7 15.1 32.5 16.2	7.6 8.1	4.2 46.3 103.2 204.1 371.6 4.0 48.0 96.0 192.0 384.0
	Crafted-6	Actual Predicted	130.0 69.8 25.6 12.8 27.6 74.9 37.4 18.7 9.3 48.0	51.5 140.5 279.5 96.0 192.0 384.0
	Crafted-7	Actual Predicted	95.0 51.3 28.4 11.6 37.8 74.9 37.4 18.7 9.3 48.0	70.0 126.3 308.0 96.0 192.0 384.0
	Crafted-8	Actual Predicted	17.2 10.8 20.2 10.1	5.3 5.0	2.6 56.4 2.5 48.0	89.6 181.1 363.6 96.0 192.0 384.0
	Crafted-9	Actual Predicted	27.2 12.1 27.0 13.5	5.9 6.7	3.6 47.5 106.6 217.3 358.0 3.3 48.0 96.0 192.0 384.0
		1.0
		0.8
	Prob Sol	0.4 0.6			Weibull Exponential Empirical
		0.2			Log Normal
					Beta Prime
		0.0	0 2000 4000 6000 8000 10 000
				Time s

Table 8

 8 Sequential executions in seconds for complete search. * Indicates that for AI-18 with min-dom we use non-shifted distributions because the min value is considerably smaller than the mean

	probability					probability
	0.025					0.014
						0.012
	0.020					0.010
	0.015					0.008
	0.010					0.006
						0.004
	0.005					0.002
	50	100	150	200	runtime		100 200 300 400 500 600	runtime
		(a) MS 15					(b) AI 18
			probability
			0.008
			0.006
			0.004
			0.002
				500	1000	1500	2000	2500 runtime
					(c) Costas 17

Table 9

 9 Parallel runtimes in seconds (proc. stands for processors). Each cell in the performance indicates the runtime (top) and speedup (bottom) for MS with wdeg, AI and Costas with min-dom thus k i denotes the i th sequential runtime. The associated empirical sequential distribution is the probability distribution defined by assigning a probability 1/m to each sequential execution. Therefore the probability of each observation can be defined as P[Y ≤ t i ]= i/m. Considering the min distribution rule of Section 3.2 we compute the empirical parallel runtime distribution F Z n as follows:

	Problem	1 proc.		performance (time/speed-up) on k proc. 24 48 96 192 384
	MS 15	38.20	Actual Predicted	2.06 18.546 1.78 21.46	0.66 57.87 0.99 38.58	0.47 81.27 0.59 64.74	0.31 123.22 0.39 97.94	0.25 152.8 0.29 131.72
	AI 18	58.08	Actual Predicted	3.54 16.41 2.42 24.00	1.56 37.23 1.21 48.00	0.69 84.17 0.60 96.00	0.42 138.29 0.30 192.00	0.28 207.43 0.15 384.00
	Costas 17	40.39	Actual Predicted	3.41 11.83 1.93 20.92	1.64 24.62 1.09 37.05	0.93 43.43 0.68 59.39	0.49 82.42 0.47 85.93	0.33 122.39 0.36 112.191

Table 10

 10 Probability of selecting each element k i ∈ K, for |K|=100 Average Discrepancy across all instances w.r.t. the actual solution for the empirical evaluation (EE) and the predicted evaluation (PE) varying the sequential executions and number of cores

		1.0	32 proc.
		0.8	16 proc.
	P Sol	0.4 0.6	2 proc. 4 proc. 8 proc.
		0.2		1 proc.
		0.0	0	20 40 60 80 100
			Sorted seq. executions
	Fig. 15 Cores Method	100%	Discrepancy 50% 25%	12%
		48	EE PE		17.92 19.04 25.66 21.80	33.71 23.67 27.79 48.25
		96	EE PE		18.36 25.45	24.40 20.11 25.17 31.50 47.58 51.62
		192	EE PE		22.64 23.90	33.33 22.76 27.06 31.89 50.91 65.04
		384	EE PE		28.22 21.89 24.86 28.25 29.29 40.03 52.89 67.05
	ALL	EE PE		21.79 24.23	29.20 22.38 26.04 30.12 46.27 57.99

http://www.csplib.org

We also experimented with industrial instances but the tested LS solvers performed poorly and a very limited number of instances could be solved with a reasonable time limit.

http://baldur.iti.kit.edu/SAT-Challenge-2012/results.html

www.gecode.org

We use AFC, the Gecode implementation of wdeg

We use the exponential distribution in lieu of the shifted exponential because the min value is negligible w.r.t. to the mean 0.008 vs. 58.08.

For experiments with the Adaptive Search we use report values with32, 64, 128, and