François Gonard

Marc Schoenauer

Michèle Sebag

Algorithm Selector and Prescheduler in the ICON challenge

Algorithm portfolios are known to offer robust performances, efficiently overcoming the weakness of every single algorithm on some particular problem instances. Two complementary approaches to get the best out of an algorithm portfolio is to achieve algorithm selection (AS), and to define a scheduler, sequentially launching a few algorithms on a limited computational budget each. The presented Algorithm Selector And Prescheduler system relies on the joint optimization of a pre-scheduler and a per instance AS, selecting an algorithm well-suited to the problem instance at hand. ASAP has been thoroughly evaluated against the state-of-the-art during the ICON challenge for algorithm selection, receiving an honourable mention. Its evaluation on several combinatorial optimization benchmarks exposes surprisingly good results of the simple heuristics used; some extensions thereof are presented and discussed in the paper.

Introduction

In quite a few domains related to combinatorial optimization, such as satisfiability, constraint solving or operations research, it has been acknowledged for some decades that there exists no universal algorithm, dominating all other algorithms on all problem instances. This result, referred to as No Free Lunch theorem [START_REF] Wolpert | No free lunch theorems for optimization[END_REF], has prompted the scientific community to design algorithm portfolios addressing the various types of difficulties involved in the problem instances, i.e., such that at least one algorithm in the portfolio can efficiently handle any problem instance [START_REF] Huberman | An economics approach to hard computational problems[END_REF][START_REF] Gomes | Algorithm portfolios[END_REF]. Algorithm portfolios thus raise a new issue, that of selecting a priori an algorithm well suited to the application domain [START_REF] Leyton-Brown | A portfolio approach to algorithm selection[END_REF]. This issue, referred to as Algorithm Selection (AS) and first formalized by Rice [START_REF] Rice | The algorithm selection problem[END_REF], is key to the successful transfer of algorithms outside of research labs. It has been tackled by a number of authors in the last years (more in section 2).

Algorithm selection comes in different flavors, depending on whether the goal is to yield an optimal performance in expectation with respect to a given distribution of problem instances (global AS), or an optimal performance on a particular problem instance (per instance AS). Note that the measure of performance depends on the domain (e.g., time-to-solution in satisfiability, or time to reach the optimal solution up to a given precision in optimization 4 ). This paper focuses on the per-instance setting, aimed at achieving peak performance on every problem instance.

In some domains, it is often the case that some problems can be solved in no time by some algorithms. It thus makes sense to allocate a part of the computational budget to a pre-scheduler, sequentially launching a few algorithms with a small computational budget each. The pre-scheduler is expected to solve "easy" instances in a first stage; in a second stage, AS is only launched on problem instances which have not been solved in the pre-scheduler phase. Note that the prescheduler enables to extract some additional information characterizing the problem at hand, which can be used together with the initial information about the problem instance, to support the AS phase.

This paper presents the Algorithm Selector And Prescheduler system (ASAP), aimed at algorithm selection in the domain of combinatorial optimization (section 3). The main contribution lies in the joint optimization of both a pre-scheduler and a per-instance algorithm selector. The extensive empirical validation of ASAP is conducted on the ICON challenge on algorithm selection [START_REF] Kotthoff | ICON challenge on algorithm selection[END_REF]. This challenge leverages the Algorithm Selection library [START_REF] Bischl | ASlib: A Benchmark Library for Algorithm Selection[END_REF], aimed at the fair, comprehensive and reproducible benchmarking of AS approaches on 13 domains ranging from satisfiability to operations research (section 4).

The comparative empirical validation of ASAP demonstrates its good performances comparatively to state-of-art pre-schedulers and AS approaches (section 5), and its complementarity with respect to the prominent Zilla algorithms [START_REF] Xu | Satzilla: portfolio-based algorithm selection for sat[END_REF]. The paper concludes with a discussion of the limitations of the ASAP approach, and some perspectives for further research.

Related work

Algorithm selectors

The algorithm selection issue, aimed at selecting the algorithm best suited to the problem at hand, was first formalized by Rice [START_REF] Rice | The algorithm selection problem[END_REF] as follows. Given a problem space mapping each problem instance onto a description x thereof (usually x in IR d ) and the set A of algorithms in the portfolio, let us denote G(x, a) a performance model, mapping each (x, a) pair onto the performance of algorithm a onto problem instance x. AS most naturally follows from such a performance model by selecting for each problem instance x the algorithm a with optimal G(x, a).

AS(x) = arg max a∈A {G(x, a)} (1)
The performance model is most usually built by applying machine learning approaches onto a dataset reporting the algorithm performances on a comprehensive set of benchmark problem instances (with the exception of [START_REF] Gagliolo | Algorithm portfolio selection as a bandit problem with unbounded losses[END_REF], using a multi-armed bandit approach). Such machine learning approaches range from k-nearest neighbors [START_REF] O´mahony | Using case-based reasoning in an algorithm portfolio for constraint solving[END_REF] to ridge regression [START_REF] Xu | Satzilla: portfolio-based algorithm selection for sat[END_REF], random forests [START_REF] Xu | Satzilla2012: improved algorithm selection based on cost-sensitive classification models[END_REF], collaborative filtering [START_REF] Stern | Collaborative expert portfolio management[END_REF][START_REF] Mısır | Algorithm selection as a collaborative filtering problem[END_REF], or learning to rank approaches [START_REF] Oentaryo | Algorithm selection via ranking[END_REF].

As expected, the efficiency of the machine learning approaches critically depends on the quality of the training data: the representativity of the problem instances used to train the performance model, and even more importantly, the description of the problem instances. Considerable care has been devoted to the definition of descriptive features in the SAT and Constraint domains [START_REF] Xu | Features for SAT[END_REF].

Schedulers

Besides AS, an algorithm portfolio can also take advantage of parallel computer architectures, by launching several algorithms working independently or in cooperation on the considered problem instance (see e.g. [START_REF] Yun | Learning algorithm portfolios for parallel execution[END_REF], [START_REF] Kadioglu | Algorithm selection and scheduling[END_REF]). Schedulers embed the parallel solving strategies in a sequential computer architecture, by considering a sequence of κ (algorithm a i , time-out τ i ) pairs, where the problem instance is successively tackled by algorithm a i with a computational budget τ i , until being solved. Notably, the famed restart strategy, launching a same algorithm with different random seeds or different initial conditions can be viewed as a particular case of scheduling strategy [START_REF] Gomes | Algorithm portfolios[END_REF]. Likewise, AS can be viewed as a particular case of scheduler with κ = 1 and τ 1 set to the overall computational budget.

As shown by [START_REF] Xu | Satzilla2012: improved algorithm selection based on cost-sensitive classification models[END_REF], schedulers and AS can be combined together along a multi-stage process, where a scheduler solves easy instances in a first stage, and remaining instances are handled by the AS and tackled by the selected algorithm in the next stage. [START_REF] Kadioglu | Algorithm selection and scheduling[END_REF] build per-instance schedules where the AS is one of the component algorithms incorporated.

Overview of ASAP

This section first discusses the rationale for the ASAP approach, before detailing the pre-scheduler and AS modules in ASAP.V1. Extensions thereof, forming ASAP.V2, are presented thereafter.

Analysis

A benchmark suite most generally involves easy and hard problem instances. The difficulty is that the hardness of a problem instance depends on the considered algorithm. As shown on Fig. 1 in the case of the SAT11-HAND dataset (section 4), while several algorithms might solve 20% of the problem instances within seconds, the oracle (selecting the best one out of these algorithms for each problem instance) solves about 40% of the problem instances within seconds.

Accordingly, one might want to launch each one of these algorithms for a few seconds each on each problem instance: after this stage, referred to as pre-scheduler stage, circa 40% of the overall problem instances would be solved.

Definition 1 (Pre-scheduler). Let A be a set of algorithms. A κ-component pre-scheduler, defined as a sequence of κ (algorithm a i , time-out τ i ) pairs,

((a i , τ i ) κ i=1 ) with (a i , τ i ) ∈ A × R + , ∀i ∈ 1, . . . , κ
sequentially launches algorithm a j on any problem instance x until either a j solves x, or time τ j is reached, or a j stops without solving x. If x has been solved, the execution stops. Otherwise, j is incremented while j ≤ κ.

Note that a pre-scheduler contributes to reduce the impact of the AS failures. An AS failure is manifested as a problem instance x for which the AS selects an inappropriate algorithm (requiring much computational resources to solve x or even failing to solve it), although there exists another algorithm which could have solved x in no time. Since a pre-scheduler increases the chance for each problem instance to be solved in no time, everything else being equal, it therefore mitigates the chances and impact of AS failures.

After this discussion, the ASAP system involves two modules, a pre-scheduler and an AS. The pre-scheduler is meant to solve as many problem instances as possible in a first stage, and the AS takes care of the remaining problem instances. A primary decision concerns the division of labor between the two modules: how to split the available runtime between the two, and how many algorithms are involved in the pre-scheduler (parameter κ). It is clear that the number of problem instances solved by a module will increase with its computational budget, everything else being equal; the pre-scheduler and the AS modules are interdependent. For simplicity and tractability however, the maximal runtime allocated to the pre-scheduler is fixed to T max ps (10% of the overall computational budget in the experiments, section 5), and the number κ of algorithms in the prescheduler is set to 3. [START_REF] Kadioglu | Algorithm selection and scheduling[END_REF] and [START_REF] Malitsky | Algorithm portfolios based on cost-sensitive hierarchical clustering[END_REF] use a most close setup for their fixed-split selection schedules, except they do not constrain the AS component to take the last part of the schedule.

Given T max ps and κ, ASAP tackles the optimization of the pre-scheduler and the AS modules. It is clear that both optimization problems remain inter-dependent: the AS should mostly focus on the problem instances which are not solved by the pre-scheduler, while the pre-scheduler should symmetrically focus on the problem instances which are most uncertain or badly identified by the AS. Formally, this interdependence is handled as follows:

-A performance model G(x, a) is built for each algorithm over all training problem instances, defining AS init (Eq. 1); -A pre-scheduler is built to optimize the joint performance (pre-scheduler, AS init ) over all training problem instances;

-Another performance model G2(x, a) is built over all training problem instances, using an additional boolean feature that indicates for each problem instance whether it was solved by the above pre-scheduler; let AS post denote the AS based on performance model G2(x, a).

ASAP finally is composed of the pre-scheduler followed by AS post .

ASAP.V1 pre-scheduler

Let (a i , τ i ) κ i=1 denote a pre-scheduler, with overall computational budget T ps = κ i=1 τ i , and let F ((a i , τ i ) κ i=1 ) denote the associated domain-dependent performance. ASAP.V1 considers for simplicity equal time-outs (a i = T ps κ , i = 1 . . . κ). The pre-scheduler is thus obtained by solving the following optimization problem: max T ps ≤T max ps ,a1,...aκ

F (a i , T ps κ ) κ i=1 (2)
This mixed optimization problem is tackled in a hierarchical way, determining for each value of T ps the optimal κ-uple of algorithms a 1 . . . a κ . Thanks to both small κ values (κ = 3 in the experiments) and small number of algorithms (≤ 31 in the ICON challenge, section 4), the optimal κ-uple is determined by exhaustive search conditionally to the T ps value.

The ASAP.V1 pre-scheduler finally relies on the 1-dimensional optimization of the overall computational budget T ps allocated to the pre-scheduler. In all generality, the optimization of T ps is a multi-objective optimization problem, e.g. balancing the overall number of problems solved and the overall computational budget. Multi-objective optimization commonly proceeds by determining the so-called Pareto front, made of non-dominated solutions. In our case, the Pareto front depicts how the performance varies with the overall computational budget, as illustrated on Fig. 2, where the performance is set to the number of solved instances.

In multi-objective decision making [START_REF] Das | On characterizing the "knee" of the pareto curve based on normal-boundary intersection[END_REF][START_REF] Branke | Finding knees in multi-objective optimization[END_REF], the choice of a solution on the Pareto front is tackled using post-optimal techniques [START_REF] Deb | Multi-objective evolutionary algorithms: Introducing bias among pareto-optimal solutions[END_REF], including: i) compromise programming, where one wants to find the point the closest to an ideal target in the objective space; ii) aggregation of the objectives into a single one, e.g., using linear combination; or iii) marginal rate of return. The last heuristics consists of identifying the so-called "knees", that is, the points where any small improvement on a given criterion is obtained at the expense of a large decrease on another criterion, defining the so-called marginal rate of return. The vanilla marginal rate of return is however sensitive to strong local discontinuities; for instance, it would select point A in Fig. 2. Therefore, a variant taking into account the global shape of the curve, and measuring the marginal rate of improvement w.r.t. the extreme solutions on the Pareto front is used (e.g., selecting point K instead of point A in Fig. 2).

ASAP.V1 algorithm selector

As detailed in section 3.1, the AS relies on the performance model learned from the training problem instances. Two learning algorithms are considered in this paper: random forests and k-nearest neighbors. One hyper-parameter was adapted for each ML approach (all other hyperparameters being set to their default value, using the Python scikit-learn library [START_REF] Pedregosa | Scikit-learn: Machine learning in Python[END_REF]), based on a few preliminary experiments: 35 trees are used for the RandomForest algorithm and the number of neighbors is set to k = 3 for the k-nearest neighbors. In the latter case, the predicted value associated to problem instance x is set to the weighted sum of the performance of its nearest neighbors, weighted by their relative distance to x:

F(x, a) = i ||x -x i ||F(a, x i ) i ||x -x i ||
where x i ranges over the 3 nearest neighbors of x. Features were normalized (zero mean, unit variance) before selecting the neighbors.

A main difficulty comes from the descriptive features forming the representation of problem instances. Typically, the feature values are missing for some groups of features, for quite a few problem instances, due to diverse causes (computation exceeded time limit, exceeded memory, presolved the instance, crashed, other, unknown). The lack of feature value is handled by i) replacing the missing value by the feature average value; ii) adding to the set of descriptive features 7 additional boolean features per group of initial features, indicating whether the feature group values are available or the reason why they are missing otherwise5 .

ASAP.V2

Several extensions of ASAP.V1 have been considered after the closing of the ICON challenge, aimed at exploring a richer pre-scheduler-AS search space while preventing the risk of overfitting induced by a larger search space.

We investigated the use of different time-outs for each algorithm in the pre-scheduler, while keeping the set of algorithms (a 1 , . . . , a κ ) and the overall computational budget T ps . The sequential optimization strategy (section 3.2), deterministically selecting T ps as the solution with maximal average return rate, exhaustively determining the κ-uple of algorithms conditionally to T ps , is thus extended to optimize the (τ 1 , . . . τ κ-1 ) vector conditionally to κ-1 i=1 ≤ T ps , using a prominent continuous black-box optimizer, specifically the Covariance-Matrix Adaptation-Evolution Strategy (CMA-ES) [START_REF] Hansen | Reducing the time complexity of the derandomized Evolution Strategy with Covariance Matrix Adaptation[END_REF].

This extended search space is first investigated by considering the raw optimization criterion F raw ((a i , τ i ) κ i=1 ) defined in section 3.2, that is, the cumulative performance of ASAP over all training problem instances. However a richer search space entails some risk of overfitting, where the higher performance on data used to optimize ASAP (training data) is obtained at the expense of a lower performance on test data. Generally speaking, the datasets used to train an AS are small ones.

A penalized optimization criterion is thus considered:

F L2 ((a i , τ i ) κ i=1 ) = F ((a i , τ i ) κ i=1 ) + w κ i=1 τ i - T ps κ 2
which penalizes the L 2 distance between the (τ i ) vector and the uniform time outs (τ i = T ps κ ). The rationale for this penalization is to prevent brittle improvements on the training set due to opportunistic adjustments of the τ i s, at the expense of stable performances on further instances. The penalization weight w is adjusted using a nested CV process.

A randomized optimization criterion is also considered. By construction, the fitness function aggregates the performances of all training problem instances. As the training problem instances sample the problem domain, this fitness defines a noisy optimization problem. Sophisticated approaches have been proposed to address this noisy optimization issue in non-convex optimizationbased machine learning settings (see e.g. [START_REF] Heidrich-Meisner | Hoeffding and Bernstein races for selecting policies in evolutionary direct policy search[END_REF]). Another approach is proposed here, based on the bootstrap principle: in each CMA-ES generation, the set of n problem instances used to compute the performance is uniformly drawn with replacement from the n-size training set. In this manner, each optimization generation considers a slightly different optimization objective noted F rand , thereby discarding the insignificant improvements.

Finally, a probabilistic optimization criterion is considered, handling the ASAP performance on a single problem instance as a random variable with a triangle-shape distribution (Fig. 3) centered on the actual performance p(x), with support in [p(x) -θ, p(x) + θ], and taking the expectation thereof. The merit of this triangular probability distribution function is to allow for an analytical computation of the overall fitness expectation, noted F df p . 4 Experimental setting: The ICON challenge

ASlib data format

Due to the difficulty of comparing the many algorithm selection systems and the high entry ticket to the AS field, a joint effort was undertaken to build the Algorithm Selection Library (ASlib), providing comprehensive resources to facilitate the design, sharing and comparison of AS systems [START_REF] Bischl | ASlib: A Benchmark Library for Algorithm Selection[END_REF]. ASlib (version 1.0.1) involves 13 datasets, also called scenarios (Table 1), gathered from recent challenges and surveys in the operations research, artificial intelligence and optimization fields. The interested reader is referred to [START_REF] Bischl | ASlib: A Benchmark Library for Algorithm Selection[END_REF] for a more comprehensive presentation. Each dataset includes i) the performance and computation status of each algorithm on each problem instance; ii) the description of each problem instance, as a vector of the expert-designed feature values (as said, this description considerably facilitates the comparison of the AS systems); iii) the computational status of each such feature (e.g. indicating whether the feature could be computed, or if it failed due to insufficient computational or memory resources). Last but not least, each dataset is equi-partitioned into 10 subsets, to enforce the reproducibility of the 10 fold CV assessment of every AS algorithm.

The ICON Challenge on Algorithm Selection

The ICON Challenge on Algorithm Selection, within the ASlib framework, was carried on between February and July 2015 to evaluate AS systems in a fair, comprehensive and reproducible manner6 . Each submitted system was assessed on the 13 ASlib datasets [START_REF] Bischl | ASlib: A Benchmark Library for Algorithm Selection[END_REF] with respect to three measures: i) number of problem instances solved; ii) extra runtime compared with the virtual best solver (VBS, also called oracle); and iii) Penalized Average Time-10 (PAR10) which is the cumulative runtime needed to solve all problem instances (set to ten times the overall computational budget whenever the problem instance is unsolved).

As the whole datasets were available to the community from the start, the evaluation was based on hidden splits between training and test set. Each submitted system provides a datasetdependent, instance-dependent schedule of algorithms, optionally preceded by a dataset-dependent presolver (single algorithm running on all instances during a given runtime before the per-instance schedule runs). Each system can also, in a dataset-dependent manner, specify the groups of features to be used (in order to save the time needed to compute useless features).

Two baselines are considered: the oracle, selecting the best algorithm for each problem instance; and the single best (SB) algorithm, with best average performance over all problem instances in the dataset. The baselines are used to normalize every system performance over all datasets, associating performance 0 to the oracle (respectively performance 1 to the single best), supporting the aggregation of the system results over all datasets.

Experimental validation

Comparative results

Table 2 reports the results of all submitted systems on all datasets (the statistical significance tests are reported in Fig. 5). The general trend is that zilla algorithms dominate all other algorithms on the SAT datasets, as expected since they have consistently dominated the SAT contests in the last decade. On non-SAT problems however, zilla algorithms are dominated by ASAP RF.V1.

The robustness of the ASAP approach is demonstrated as they never rank last; they however perform slightly worse than the single best on some datasets. The rescaled performances of ASAP RF.V1 is compared to zilla and autofolio (Fig. 4, on the left), demonstrating that ASAP RF.V1 offers a balanced performance, significantly lower than for zilla and autofolio on the SAT problems, but significantly higher on the other datasets; in this respect it can be viewed as a low-risk system.

Sensitivity analysis

The sensitivity analysis conducted after the closing of the challenge compares ASAP.V2 (with different time-outs in the pre-scheduler) and ASAP.V1, and examines the impact of the different optimization criteria, aimed at avoiding overfitting: the raw fitness, the L2-penalized fitness, the randomized fitness and the probabilistic fitness (section 3.4).

The impact of the hyper-parameters used in the AS (number of trees set to 35, 100, 200, 300 and 500 trees in the Random Forest) is also investigated.

Table 3 summarizes the experimental results that each ASAP.V2 configuration would have obtained in the ICON challenge framework, together with the actual submissions results, including On the left: per-dataset performances of ASAP RF.V1 (balls, dotted line), zilla (no marker, dashed line) and autofolio (triangles, solid line) scaled to the range of performance of all submitted systems. As a comparison, the per-dataset best submitted system (small balls, solid line) and ASAP RF scores before rescaling are depicted on the right. systems that were not competing in the challenge: llama-regr and llama-regrPairs from the organizers, and autofolio-48 which is identical to autofolio but with 48h time for training (12h was the time limit authorized in the challenge) [START_REF] Kotthoff | ICON challenge on algorithm selection[END_REF].

The significance analysis, using a Wilcoxon signed-rank test, is reported in Fig. 5. A first result is that all ASAP.V2 variants improve on ASAP.V1 with significance level 1%. A second result is that ASAP.V2 with the probabilistic optimization criterion is not statistically significantly different from zilla, autofolio and zillafolio.

A third and most surprising result is that the difference between the challenge-winner zilla and most of ASAP.V2 variants is not statistically significant. A new hybrid algorithm selection approach, the ASAP system was presented in this paper, combining a pre-scheduler and a per-instance algorithm selector. ASAP.V1 introduced a selector learned conditionally to a predetermined schedule so that it focuses on instances that were not solved by the pre-scheduler . ASAP.V2 completes the loop as it re-adapts the schedule to the new AS. The main message is that the scheduler and the AS must be optimized jointly to reflect the division of labor achieved by these two components ASAP.V1, thoroughly evaluated in the ICON challenge on algorithm selection (ranked 4th) received an honourable mention, due to its novelty and good performance comparatively to the famed and long-known Zilla algorithms.

The ASAP.V2 extension achieved significantly better results along the same challenge setting. It must be emphasized that these results must be comforted by additional experiments on fresh data. A main lesson learned is the importance of the regularization, as the amount of available data does not permit to consider richer AS search spaces without incurring a high risk of overfitting. The probabilistic performance criterion successfully contributed to a more stable optimization problem. Further research work will be devoted to extending this criterion.

Fig. 1 .

 1 Fig. 1. Percentage of solved instances vs. runtime on the SAT11-HAND dataset, for 5 algorithms and the oracle (selecting the best algorithm out of 5 for each problem instance).

Fig. 2 .

 2 Fig.2. Among a set of Pareto-optimal solutions, solution A has the best marginal rate of return; solution K, which maximizes the average rate of return w.r.t. the extreme solutions of the Pareto front (maximizing angle γ), is the knee selected in ASAP.

Fig. 3 .

 3 Fig.3. Schedule execution difference between punctual and triangular pdf. On the left (punctual pdf), schedule stops during step 1. On the right (triangular pdf, part of the instance is not solved during step 1 and the schedule executes until step 2 solves the rest.

  Fig.[START_REF] Deb | Multi-objective evolutionary algorithms: Introducing bias among pareto-optimal solutions[END_REF]. On the left: per-dataset performances of ASAP RF.V1 (balls, dotted line), zilla (no marker, dashed line) and autofolio (triangles, solid line) scaled to the range of performance of all submitted systems. As a comparison, the per-dataset best submitted system (small balls, solid line) and ASAP RF scores before rescaling are depicted on the right.

Fig. 6

 6 Fig. 6 details per dataset the performance improvement between ASAP.V2 (500 trees, f L 2 version) and ASAP.V2 (500 trees, f df p version) and on the other hand ASAP.V1 RF (35 trees). Note that ASAP.V2 outperforms the per-dataset best submission to the challenge for 3 datasets.

Fig. 6 .

 6 Fig.6. Per-dataset performances of ASAP RF.V2 with 500 trees, optimized using f L 2 (on the left with triangles) and using f df p (on the right with squares) scaled to the range of performance of the challenge submitted systems. As a comparison, ASAP RF.V1 appears with the balls and dashed line.

Table 1 .

 1 ASlib

		datasets (V1.0.1)
	dataset	# instances # algorithms # features
	ASP-POTASSCO	1294	11	138
	CSP-2010	2024	2	86
	MAXSAT12-PMS	876	6	37
	PREMARSHALLING-ASTAR-2013	527	4	16
	PROTEUS-2014	4021	22	198
	QBF-2011	1368	5	46
	SAT11-HAND	296	15	115
	SAT11-INDU	300	18	115
	SAT11-RAND	600	9	115
	SAT12-ALL	1614	31	115
	SAT12-HAND	767	31	115
	SAT12-INDU	1167	31	115
	SAT12-RAND	1362	31	115

Table 2 .

 2 Normalized performances of submitted systems, aggregated across all folds and all measures (the lower, the better). Ranks of zilla (challenge winner) and ASAP RF.V1 (honourable mention) are given in parenthesis. Numbers were computed from the challenge outputs.

		ASAP RF.V1	ASAP kNN.V1	autofolio	flexfolio	sunny	sunny-presolv	zilla	zillafolio
	ASP-POTASSCO	0.294 (2) 0.359 0.299 0.314 0.37 0.336 0.319 (5) 0.283
	CSP-2010	0.146 (1) 0.247 0.288 0.223 0.263 0.406	0.2 (3)	0.157
	MAXSAT12-PMS	0.168 (4) 0.159 0.45 0.149 0.166 0.224 0.201 (5) 0.233
	PREMARSHALLING-ASTAR-2013 0.349 (4) 0.369 0.359 0.307 0.325 0.296 0.374 (7) 0.385
	PROTEUS-2014	0.16 (4) 0.177 0.222 0.056 0.134 0.103 0.245 (8) 0.223
	QBF-2011	0.097 (2) 0.091 0.169 0.096 0.142 0.162 0.191 (7) 0.194
	SAT11-HAND	0.341 (4) 0.318 0.342 0.342 0.466 0.464 0.328 (3) 0.302
	SAT11-INDU	1.036 (5) 0.957 0.875 1.144 1.13 1.236 0.905 (2) 0.966
	SAT11-RAND	0.104 (6) 0.09 0.046 0.226 0.116 0.088 0.053 (2) 0.067
	SAT12-ALL	0.392 (5) 0.383 0.306 0.502 0.509 0.532 0.273 (1) 0.322
	SAT12-HAND	0.334 (5) 0.31 0.256 0.434 0.45 0.467 0.272 (2) 0.296
	SAT12-INDU	0.955 (6) 0.919 0.604 0.884 1.074 1.018 0.618 (3) 0.594
	SAT12-RAND	1.032 (5) 1.122 0.862 1.073 1.126 0.97 0.779 (1) 0.79

Table 3 .

 3 Optimized pre-scheduler performances8 aggregated across all datasets, all splits and all measures (the lower, the better). The hyperparameters for f L 2 and f df p were chosen after preliminary experiments using the cross validation provided with ASlib. For each configuration of the selector, the best-evaluated fitness function appears in bold fitness function (if relevant) f L 2 f df p f rand fraw Wilcoxon signed-rank test p-value between ASAP.V2 variants with 500 trees and every other systems including the single best algorithm (SB) and the virtual best solver (VBS). In particular, the different fitness functions do not differ significantly from each other in most cases.

	none

This research work has been funded by the French Program "Investissements d'Avenir".[START_REF] Deb | Multi-objective evolutionary algorithms: Introducing bias among pareto-optimal solutions[END_REF] One often considers the joint problems of selecting an algorithm and the optimal hyper-parameters thereof, referred to as Algorithm Configuration (AC), as the choice of the hyper-parameter values governs the algorithm performance. AC is outside the scope of the paper and will not be further considered.

The increase in the overall number of features is handled by an embedded feature selection mechanism, removing all features with negligible importance criterion (< 10 -5 in the experiments) in a separately learned 10-trees random forest regression model.

The codes of all submitted systems and the results are publicly available, http://challenge.iconfet.eu/challengeas

As the CSP-2010 dataset gives the choice between only two algorithms, the pre-scheduler consists of a single algorithm running for the whole time devoted to the pre-scheduler. For that particular case, all ASAP RF.V2 variants with the same selector hyperparameter are identical.


Keywords:

In-line process control in microelectronic manufacturing requires real-time and non-invasive monitoring techniques. Among the different metrology techniques, the scatterometry, based on the ellipsometric signatures (i.e stokes coefficients vs. wavelengths) of light scattered by a patterned structures, seems to be well adapted to address this kind of problem. For instance, the shape evolution of a grating in real time during etch processing can be monitored. Traditionally, the direct problem of defining the shape and computing the signature is dealt with modal methods. However, the inverse problem can not be solved as easily. Different classes of algorithms have been introduced (evolutionary, simplex), but the method of library searching is more commonly used in industry. The main limit of this method is the acquisition time of data for different wavelengths. Then, a lack of data can leads to the method failure and several database patterns can be matched to the experimental data. In this article, a technique for real time reconstruction of grating shape variation using dynamic scatterometry is presented. The different tools to realize this reconstruction, such as Modal Method by Fourier Expansion, regularization technique and mostly the specific soft and hardware architecture are also developed. Results of a dynamic experiments will illustrated these presentations.

INTRODUCTION

A part of future advances in microelectronic is conditioned by the capability for industry to monitor the manufacturing process of their chips. Since the dimensions of transistors become smaller and smaller, an important effort has to be made in the field of metrology. CD metrology techniques must be able to characterize very short dimensions with a good accuracy. AFM and XSEM techniques are suitable to reach these accuracies, but these are invasive or long time consuming and actually, not well adapted to the industry requirements. Indeed, the next step in the metrology improvement is the development of a technique able to follow the evolution of a pattern during a dynamic process. It means that the metrology tool must be placed in situ, the time of acquisition of data must be very fast and this technique must be non invasive. This objective is not a simple task but first experimental tests can be carried out for simplest processing (as described below). Based on the ellipsometric measurement of light diffracted by a grating, the scatterometry technique appears to be promising for in line process monitoring. The scatterometry is an optical and a non-invasive technique of metrology. It measures the change of polarization subsequent to the scattering of a ray of light. The used ellipsometer gives signatures for a single angle but for several wavelengths. The signatures resulting of the measurement are compared with simulated signatures coming from modelling for the best match. The simulated signatures are achieved by the Modal Method by Fourier Expansion starting from a given topography. This is the direct problem. On the opposite, the inverse problem consists in finding the shape of the grating thanks to the signatures and this is not a simple task. Different classes of algorithms have been developped (evolutionary, simplex), 1 but the method of library searching is practically the more commonly used in industry. Moreover, this method is real-time complient as the solving time is perfectly predictible. The work presented in this article will deal with improvements made on this technique for the case of dynamic scatterometry. A signature library is built-up using different combinations of profile parameters (CD, height, sidewall angle, etc) and the experimental signature is compared with the library for the best match. Although generating the library requires a lot of computation time, solving the inverse problem can be computed quickly because pattern matching is reduced to a simple numerical database search. The constraint is the limit in wavelength resolution. Indeed, in order to follow the real time process evolution, the acquisition of data is limited. This poor resolution in wavelength can leads to matching problems (i.e. several database can matched to the experimental data). In this paper, we present a technique which does not need to solve the full inverse problem at each time step but which keep a good accuracy on results. The main idea is to achieve a first and complete signature measurement for the whole scale of wavelengths at t=0 before processing, and partial signatures at next steps to intelligently refine the library search using Tikhonov time regularization. This technique assumes that only a smooth variation of the shape over time is realistic. This approach is implemented as a general C++ architecture dedicated to scatterometry. The system integrates multiprocessing (shared memory) capability for high-speed signature generation and the use of the Graphic Processin Unit (GPU, video card processor) as a numerical database accelerator for the library-based inverse problem. Due to the inherent vector architecture of the GPU, the entire library search operation can be made in only a few clock cycles. In a short first part, a brief description of the direct problem is presented. The second part deals with the real time regularisation description and the third part describes the software and the hardware tools. Finally, the first results are shown with a specific example. This application deals with the time evolution of the shape of a pattern resist when heated.

DIRECT PROBLEM SOLVING

As mentioned in the introduction part, the direct problem is used to establish simulated signatures from a given topography. It allows building the data library with parameter variations. Some modelling tools are suitable to address the direct problem. In this work, the Modal Method by Fourier Expansion is chosen. A full description of this method can be found elsewhere and will not be detailed here. Only a brief description is given. The Modal Method by Fourier Expansion is well adapted for the topography of the scatterometry problem. [START_REF] Fuard | Measurement of residual thickness using scatterometry[END_REF] MMFE is based on a simple approach where Maxwell equations are solved thanks to a matrix system. TE or TM polarization are treated using the same formalism and an improved version allows to obtain good results for TM polarization in terms of convergence speed. In our study only the TE polarization has been addressed. This method implicitly considers periodic structures. The solution of Maxwell equations is reduced to the solution of an algebraic eigenvalue problem in a discrete Fourier space. The derivation of the matrix operator involves two steps: the electromagnetic field is expanded into Floquet-Fourier series, the permittivity of the materials which appears in the Maxwell equation is here a periodic function and is expanded into Fourier series. The topography of the structure which is studied is cut in different layers. These different layers are linked by a S matrix system which accounts for the boundary condition providing a good convergence even for thick or absorbing layers. In this rigorous method, the only approximation is due to the number of modes retained for the truncation of the modal expansion. In EUV lithography, due to the short wavelength, the ratio λ/d is very small. Consequently the number of mode is important according to the following equation:

sin(θ m ) = sin(θ) -m λ d
where θ 0 is the incident angle, m is the number of mode and θ m is the angle of the mode m.

This large number of mode leads to a compromise between computation time and results accuracy. In practice, to address complex topographies and for a retained number of mode of 51, the precision is good and the computation of the reflected field is performed in around 20 seconds with an Intel Xeon 3.4GHz processor. For the scatterometry problems, the number of mode must not be important in order to generate, as fast as possible, the database. Typical library computation time, depending on m. The database contains 65.636 signatures, computed for 28 wavelengths with a gratting shape divided in 32 layers. The library size in memory is 40MBytes.

m = 0 m = 1 m = 2 m = 3 m = 10 1 h 3h 5h30 10h 5d
Depending on the number of layer : (for m = 3)

n = 32 n = 22 n = 12 10 h 6h30 3h30
For a fixed m and n, the computation time and library size are linear, with respect to the number of wavelength for a signature and the number of signature contained in the library.

It's clear that computation time has to be dramatically reduced. We are working to make the code run on multiproccessor systems. Since a library is nothing but a collection of signatures that can be computed independently, the parallel efficiency is close to 100%. We have veryfied that on a quad-Intel Xeon computer, using an OpenMP library.

REAL-TIME REGULARIZATION

Dynamic scatterometry can not be reduced to successive static problems : a high frequency acquisition and a short acquisition time (relatively to process characteristic time) are required to follow the shape variation. Thus, standard ellipsometers that usually gives full defined signature in seconds are not adapted. However, we can use multi-way ellipsometers : ellipsometric signals can be acquired quickly but only for a selected set of wavelength. The number of wavelengths is low (so that a mere static scatterometry resolution is inefficient), but the grating shape can be reconstructed if we take into account the whole set of partial signatures and the regularity (the smoothness) of the shape parameters (height, CD, side wall angle, etc.). The principle of regularization can indeed be summarized as a balance for the choice of the best parameters. This balance is made between, on an hand, the proximity of each acquired signature to the chosen signature in the library and on an other hand, the regularity (smoothness) of the parameters variation over time.

Practically, the way we proceed is :

• at t = 0, i.e. before the process, the grating shape is supposed to be perfectly known. For that we can use, for example, static scatterometry, where a full defined signature is acquired. There is no problem of computing time saving here.

• then, at each time step, we acquire a low-resolution and low-precision signature. This implies bad accuracy in database look up. The technique of regularization begin with this : we don't pick up only one Nearest Neighbor (NN) to the acquired signatures in the database but the 4 nearest neighbors. So, even if the first NN is far to be the good solution (it can be an artifact), we assume that a better solution is included in the 4 NN. The algorithm of regularization whom results are presented here selects at each time step (but relatively to the previous time steps) the best set of parameters. A standard Tikhonov-based regularization 3 is then computed :

Xmin   X -X exp 2 + β window X ′2 dt   where X = (H, CD, sidewall angle...)
The goal is to have, for each parameter describing the grating shape, a smooth and realistic variation.

Notice that, for the full regularization (integrating the "smart" library search), only two parameters have to be optimized.

HIGH PERFORMANCE COMPUTING

We've shown that the library search method is the only way to solve the inverse problem taking into account the "real-time" constraint. But the library used for this method is quite large : about half a GByte is a common size for industrial production. Large quantity of data implies two challenges :

-The computing time of this library should be improved at maximum; on a single (but powerfull : Xeon 3GHz) core, such a computation can last several weeks.

-Once computed, dealing with this data for quick library look-up is also a challenge : finding the Nearest Neighbor of an acquired signature inside a huge signature library should be made fast enough comparing to acquisition frequency of experimental data.

For the first point, we are working to make the code that computes libraries run on multiprocessor systems. Since a library is nothing but a collection of signatures that can be computed independently, the parallel efficiency is close to 100%. We have veryfied that on a dual dual-core Intel Xeon computer, using an OpenMP library.

For the second point, we use as a co-processor the GPU (Graphics processing Unit) of our graphic card.

Computer game industry has a growing need to 3D graphics capability. Nowadays, video game picture rendering is made by a special stream processor, very powerfull but very different from standard CPUs. Indeed, GPUs are said to be stream or S.I.M.D. (Single Instruction, Multiple Data)processor because of their ability to compute a whole set of data (like a picture) with a single instruction. Contrary to CPUs, who, in their simplest form compute one data by instruction. For comprehensive explanation of what can be done efficiently with a GPU, read the paper. [START_REF] Buck | Gpgpu: General-purpose computation on graphics hardware[END_REF] Using this vector processor is a bit tricky but can be very efficient for certain class of problem, the vector architecture however make it perfect for library look up. Finally, if we enumerate the advantages we have :

-For the problem of large library look-up, GPU are far faster than CPU.

-Main CPU in discharged from such computation and can be used for the second step of regularization.

-The power of GPU is increasing far fastest than GPU (see bellow)

-A graphic card with the most powerfull GPU of the moment only cost 700$...

We have said that using GPU is a bit tricky because these processors are made only for graphics. If we want to compute data with, we need to think the signature database as a collection of big pictures in video memory. And computing in database (for instance, the calculation of the distance between the acquired signature and the pictures in the database) can be view as a visual effect to be applied to the plan where the picture is mapped. The visual effect is made by a lot of sub-processors inside the GPU core called shaders. In recent GPUs, shaders are programmable and we use this to develop a technique for multi-NN search. This technique is based on the work of Benjamin Bustos et al. [START_REF] Benjamin Bustos | A graphics hardware accelerated algorithm for nearest neighbor search[END_REF] Results Here we present some results for our problem, obtained with the NVidia GeForce 8800GTX GPU. 6

The interest of using GPUs over CPU is clear on this chart. For large signature libraries (400 -700 MB), GPU make the work about 10 times faster and even more than CPUs. An important point is that using GPUs allow us to reconstruct grating shape in time ; fast enough regarding to device acquisition frequency. Using GPUs is a way to solve 10 or 20 signature frame per second, while regularizing...

RESULTS AND COMMENTS

Here we show two examples of application for dynamic scatterometry where real-time regularization is used. The first one is for demonstration purpose: we solve the inverse problem on simulated data ( i.e.. solution of the direct problem.). The second example uses a real set of scatterometry signature coming from a 32 way ellipsometer. The process we want to monitor is a reflow of patterned resist. Both of these experiment will allow us to make remarks and comments.

Results for simulated data

In this simulated experiment, we want to monitor a grating shape of which geometry is modeled by three parameters : the height H, the CD (critical dimension the width) and α the side wall angle. These parameter variation are arbitrary but includes successively a smooth part, a corner and then a flat part.

For the signature computation using the MMFE engine, each grating is modeled as a stack of 10 resist layers (neb22), lying on a silicon substrate. The angle of incidence of the light to be diffracted is set to 70. The truncate order for Fourier expansion is set to 10. These settings will be the same for the library computation.

Once the movie is generated, we add a noise to each frame (each signature) S ′ 1 = (1 + θ 1 )S 1 , S ′ 2 = (1 + θ 2 )S 2 where θ1, θ2 ∈ [0; 0.1] (5%) in order to take into account all the uncertainties due to modelization, measure and computation.

Library We choose to generate a signature library corresponding to the three-dimensional parameter space. Values taken by the parameters are :

-CD=160.32 to CD=200, by steps of 1.28 -H=160.32 to H=200, by steps of 1.28 -α=9.5 to α=20 by steps of 0.7

The choice of parameter range and density is the most tricky part in library computation : the library size in memory is proportional to the number of parameters combination. So, a good modeling work in necessary, first, to maintain the dimension of parameter space as low as possible, and then to restrict the number of parameters values. Here we have 32 × 32 × 16 = 16384 signatures in the library. Each signatures represents 16 wavelengths. Library size in memory, given by :

S Bytes = N signature × 2N λ × sizeof(float)
is here of 4MBytes.

Dynamic reconstruction

The last step of this experiment is to feed the program at each time step with simulated frames from the movie. You can see bellow the geometrical parameters set of the grating shape as it has been found using step-by-step static scatterometry using a mere library search. Geometrical parameters appears to be the most part of the time irregulars and sometimes false. The best match in library search is not necessarily the real solution : when parameters density is poor (for instance, because the modeling is difficult and thus, parameter space dimension is important); the solution vector X can be close to real shape only by some of its component. For example, the raw library search method could privileges a solution with good H and CD and insignificant α than a more significant solution where the three parameters are closer to the reality.

Thus, the reconstruction process begins with a deeper library search where not only the first nearest neighbor is picked up, but also the three next. With four nearest neighbors at each time step, we can add another criterion for selection : the criterion of proximity to the data is kept, but it is now pondered by the regularity of the parameters. It's a way to throw away "bad" solutions coming from the library. This is the pre-regularizing step. Un peu trop confidentiel non? Picking-up the 4 nearest neighbors of each acquired signature (each frame) expands the field where realistic solution can be found: Notice that the plus symbols are sometimes far from the solution, that is why a process of regularization is needed to reconstruct realistic shape variation.

The principle of regularization poses a problem when the first derivative is supposed to be discontinuouscontinuous : corners are softened. This can't be avoided but can be limited either by taking β lower as possible, or increasing frequency of acquisition.

We can see here the fundamental difference between these two regularizations : the first one does not make regular a parameter in itself (it can even be less regular than before), however, it regularizes the whole set of parameter

Resist pattern reflow

This experiment, originally carried out for resist viscosity study is here a way to begin experimental validation for our algorithms.

A nano-imprinted grating, made of resist on silicon is heated while scatterometry signatures are acquired. The hot-plate and the 32 wavelengths ellipsometer are from Horiba/Jobin Yvon. We made the temperature increase until reflow. (∼ 90 deg C) and we try to watch pattern shape variation while reflowing.

We first chose to medel pattern shape using a mere sine wave : the amplitude is used as the only parameter. When we compare reconstructed profile shape with AFM measurements mades in the same condition and at different time step, we get :

Figure 1 .

 1 Figure 1. CPUs. vs. GPUs (S.I.M.D.)

Figure 2 .

 2 Figure 2. GPUs power, in here, GFLOPS increases far faster than CPUs power.

α

  Simulated movie Using our MMFE engine, 2 we simulate an "ellipsometric movie" of 100 frames (the process control lasts 10s and we use 10 frames per second). During this movie, the values of the three parameters are supposed to vary as such :

  Refining again with a standard Tikhonov 3 time-based regularization gives : Then a comparison with "original" simulated parameters gives :

We can clearly see that there are two parts :

-A part, beginning at 40s where our simple model seems to work.

-And an other part, in the beginning of the reflow, where actually another model is required. This is obviously because the grating shape is rectangular during the first part of the process, as we can see bellow : This is typically a problem of scatterometry : when a good modeling has been found, standard scatterometry is efficient and dynamic shape reconstruction can work.

SUMMARY AND CONCLUSION