Oralie Cattan

Sahar Ghannay

Christophe Servan

Sophie Rosset

Benchmarking Transformers-based models on French Spoken Language Understanding tasks

Keywords: French Language, Language Models, Spoken Language Understanding, Benchmarking, Model costs

In the last five years, the rise of the self-attentional Transformerbased architectures led to state-of-the-art performances over many natural language tasks. Although these approaches are increasingly popular, they require large amounts of data and computational resources. There is still a substantial need for benchmarking methodologies ever upwards on under-resourced languages in data-scarce application conditions. Most pre-trained language models were massively studied using the English language and only a few of them were evaluated on French. In this paper, we propose a unified benchmark, focused on evaluating models quality and their ecological impact on two well-known French spoken language understanding tasks. Especially we benchmark thirteen well-established Transformer-based models on the two available spoken language understanding tasks for French: MEDIA and ATIS-FR. Within this framework, we show that compact models can reach comparable results to bigger ones while their ecological impact is considerably lower. However, this assumption is nuanced and depends on the considered compression method.

Introduction

The development of pre-trained language models based on Transformer architectures [START_REF] Vaswani | Attention is all you need[END_REF], such as BERT [START_REF] Devlin | BERT: Pretraining of deep bidirectional transformers for language understanding[END_REF] has recently led to significant progress in the field of natural language processing (NLP). However, the trend of training large pre-trained language models on ever larger corpora, with an ever-increasing amount of parameters has raised questions related to the usability of these approaches [START_REF] Cattan | On the Usability of Transformers-based models for a French Question-Answering task[END_REF]. Because these models require considerable computational resources, major efforts have been made to develop compact models in order to reduce the cost of using them. Compact models offer alternatives to the energyintensive models with comparable performances while reducing their computational complexity and size. These models also allow to solve some industrial problems related to online speech processing. In particular, some applications (speech recognition, speech to text, etc.) have some known problems associated to network latency, transmission path difficulties, or privacy concerns.

Spoken Language Understanding (SLU) in language dialogue systems refers to the task of producing a semantic analysis and a formalization of the user's utterance. SLU traditionally encompasses the processes of determining a broad range of information conveyed in dialogue such as identifying the domain, the intent and the concepts of the conversation.

Benchmarking models have been extensively used in the performance evaluation of NLP-based systems [START_REF] Guo | Benchmarking meaning representations in neural semantic parsing[END_REF][START_REF] Le | FlauBERT: Unsupervised language model pre-training for French[END_REF][START_REF] Farha | Benchmarking transformer-based language models for Arabic sentiment and sarcasm detection[END_REF]. With respect to our task of interest, recent research has focused on evaluating word embeddings representations. Word embeddings have proven to be effective in capturing semantic relationships between words. They are also an essential element of deep learning-based architectures.

In [START_REF] Ghannay | What is best for spoken language understanding: small but task-dependant embeddings or huge but out-of-domain embeddings?[END_REF], contextual (ELMo [START_REF] Peters | Deep contextualized word representations[END_REF]) and flat (Word2Vec [START_REF] Mikolov | Efficient estimation of word representations in vector space[END_REF], GloVe [START_REF] Pennington | GloVe: Global vectors for word representation[END_REF] and Fast-Text [START_REF] Bojanowski | Enriching word vectors with subword information[END_REF]) representations have been evaluated in order to investigate different input representations and their influence on the results obtained in SLU tasks. They highlighted the competitiveness of Word2Vec and ELMO on the French SLU corpus: MEDIA. More recently, [START_REF] Ghannay | Neural networks approaches focused on French spoken language understanding: application to the MEDIA evaluation task[END_REF] investigated the transferability of two French pre-trained BERT models [START_REF] Devlin | BERT: Pretraining of deep bidirectional transformers for language understanding[END_REF] and their integration in BiLSTM and BiLSTM+CNN-based architectures. They obtained state-of-the-art results on MEDIA using CamemBERT [START_REF] Martin | CamemBERT: a tasty french language model[END_REF] model. Lately, [START_REF] Xu | End-to-end slot alignment and recognition for cross-lingual NLU[END_REF] gave a comparison of various cross-lingual transfer approaches including the use of the multilingual BERT encoder, evaluated on MultiATIS++, a multilingual corpus for the task of language understanding. They show that mBERT brings substantial improvements on multilingual training and cross-lingual transfer tasks, yielding up to 1.4% of relative improvements on the French subcorpus.

While there are now many English corpora and models for various tasks, there are still very few resources in French. Direct comparisons between SLU models are difficult because of the lack of a unified evaluation framework and size diverse Transformer-based models in French.

Contributions1 : In this study we are investigating the use of transformer-based architectures for French language in solving an SLU problem, a concept detection task [START_REF] Bonneau-Maynard | Results of the French Evalda-Media evaluation campaign for literal understanding[END_REF]. We also assess the impact related to model compacteness on SLU performance and their ecological impact (e.g.: their CO2 cost). We focus on the benchmarking of the existing French and multilingual Transformer-based models on two French SLU corpora: MEDIA and ATIS-FR. Models considered for this benchmark are detailed in section 2, the experiments, the data, the model optimization and the experimental protocol are described in section 3. Section 4 presents the results, and section 5 proposes an analysis of some chosen results.

Transformer models considered

In order to compare the models performances on French SLU tasks, we highlight in this section and in Table 1 some characteristics of the considered models.

Among the language models that have been proposed in recent years we consider the most commonly used multilingual benchmarking models that have "reasonable" resource consumption: XLM-R [START_REF] Conneau | Unsupervised cross-lingual representation learning at scale[END_REF] and mBERT [START_REF] Devlin | BERT: Pretraining of deep bidirectional transformers for language understanding[END_REF], as well as the French models: CamemBERT [START_REF] Martin | CamemBERT: a tasty french language model[END_REF] and FlauBERT [START_REF] Le | FlauBERT: Unsupervised language model pre-training for French[END_REF] 2 since these are the only two available large models for French at the time of writing. We also evaluate compact models that are scalable Transformer-based models. We exploit distil-mBERT [START_REF] Sanh | DistilBERT, a distilled version of BERT: smaller, faster, cheaper and lighter[END_REF], a distilled version of mBERT, small-mBERT [START_REF] Abdaoui | Load what you need: Smaller versions of mutililingual BERT[END_REF] a mBERT model whose original vocabulary has been reduced to the French language and FrALBERT [START_REF] Cattan | On the Usability of Transformers-based models for a French Question-Answering task[END_REF], a recently released model for French. As mentioned in [START_REF] Kaplan | Scaling laws for neural language models[END_REF], the performance of language models is determined by multiple factors such as the pre-training objective, the layers specification or the dataset size.

Indeed, pre-training objectives vary by model and may depend on the downstream task being solved [START_REF] Joshi | SpanBERT: Improving Pre-training by Representing and Predicting Spans[END_REF]. BERT-like models adopt the masked language modeling (MLM) and next sentence prediction (NSP) objectives. MLM is a fill-in-theblank task consisting of predicting tokens of the input sequence that have been masked whereas NSP is a binary classification task to predict whether two segments are adjacent in the original text. XLM-R differs from mBERT solely in the pre-training procedure, eliminating the NSP task to address its ineffectiveness. This is also the case for CamemBERT and FlauBERT models. FrALBERT pre-training objective consists of MLM and Sentence Order Prediction (SOP) objective. SOP effectively models inter-sentence coherence by predicting whether a sentence order in a given sentence pair is swapped or not.

For monolingual and multilingual language modeling, the quality and representativeness of the large-scale datasets used are important for efficient modeling and for producing good generalization. Most models are pre-trained on the content of either Wikipedia, the web pages gathered by Common Crawl (CC) or a combination of diverse corpora. This represents a few gigabytes to several hundreds or even several terabytes of text. The other noteworthy difference is the vocabulary size with several tens of thousands of tokens for monolingual models to several hundred for multilingual models.

Models like Text-To-Text Transfer Transformer (T5) [START_REF] Raffel | Exploring the limits of transfer learning with a unified text-to-text transformer[END_REF] or GPT-3 [START_REF] Brown | Language models are few-shot learners[END_REF], have been excluded from our study, since they pose several usability problems (emphasized in the introduction). Although the GPT-3 and T-5 models are models for English, the data used for their pre-training contains many other languages, however, these languages represent only a small portion of the dataset compared to English (93% in word count for GPT-3). Furthermore, in terms of size, GPT-3 and T5 are huge com-pared to the size of established general NLP models such as BERT, and these earlier models are already quite expensive to run on GPUs. For example, GPT-3 is trained on billions of parameters, 470 times larger than the BERT model. We highlight the question of their usability in the introduction section.

In terms of availability finally, BERT is an open-source tool and easily available for users to access and fine-tune according to their needs and solve various downstream tasks. GPT-3 on the other hand is not open-sourced. It has limited access to users.

Experiments

Experiments are conducted on two well known French SLU task: MEDIA and ATIS-FR.

Datasets

The French MEDIA3 corpus, composed of 1258 transcribed dialogues, which is about hotel reservation and information [START_REF] Bonneau-Maynard | Results of the French Evalda-Media evaluation campaign for literal understanding[END_REF]. The corpus was manually annotated with semantic concepts characterized by a label and its value. There are in total 76 semantic labels. The corpus is split into three parts: a training corpus composed of 13k sentences, a development corpus composed of 1.3k sentences, and a test corpus composed of 3.5k sentences. The MEDIA task is also known as one of the most challenging slot-filling tasks, according to [START_REF] Béchet | Benchmarking benchmarks: introducing new automatic indicators for benchmarking spoken language understanding corpora[END_REF].

The French version of the Air Travel Information System (ATIS) corpus from the recent MultiATIS++ extension [START_REF] Xu | End-to-end slot alignment and recognition for cross-lingual NLU[END_REF], named ATIS-FR, concerns flight information [START_REF] Price | Evaluation of spoken language systems: The atis domain[END_REF]. The ATIS-FR corpus corresponds to the manual translation of the original ATIS sentences (in English). It is composed of 84 semantic labels, and is split into three parts: a training corpus composed of 4.5k sentences, a development corpus composed of 490 sentences, and a test corpus composed of 893 sentences.

Experimental protocol

Spoken language understanding (SLU) is considered here as a sequence labeling task that assigns a concept label (from a predefined set) to each token of a sentence. We follow BIO-tagging scheme, where each concept is associated with two labels, Blabel (for Beginning) and I-label (for Intermediate). Finally, the O-label (for Other) identifies the non-concept tokens. SLU performance is evaluated in terms of F-measure or F1 and Concept Error Rate (CER). The CER score is the official metric used in the MEDIA campaign [START_REF] Bonneau-Maynard | Results of the French Evalda-Media evaluation campaign for literal understanding[END_REF] which is estimated in the same way as the classical word error rate but applied to semantic concepts instead of words (the lower the better). The significant results are marked with a star and measured using the 95% confidence interval.

We rely in this study on the standard approach introduced in [START_REF] Devlin | BERT: Pretraining of deep bidirectional transformers for language understanding[END_REF], that corresponds to a model-based transfer learning method, used to facilitate the modeling of the target task with the knowledge learned from the language modeling task. Specifically, it consists of adding on top of the pre-trained model a token-level classifier with two hidden linear layers (with ReLU activations and dropout) that takes as input the last hidden state of the input sequence and outputs probabilities over concepts. In our experiments, we use Adam optimizer and conduct an automated hyperparameter optimization on a development dataset.

This optimization is based on the population-based learning algorithm [START_REF] Jaderberg | Population based training of neural networks[END_REF], with proven efficiency, and in which a population of models and their hyperparameters are jointly optimised. Among the hyperparameters considered are the number of training epochs from 5 to 100, the batch size in the interval of 8 and 32 or the learning rate in the range between 1 and 5.

Results

This section presents the evaluation results of the different Transformer-based pre-trained models described in section 2 on both French datasets MEDIA and ATIS-FR.

MEDIA results

Performances on the MEDIA test dataset are presented in the two first columns of the Table 2. Monolingual models scores are very close and vary from 89 to 90 of F1 while CER varies between 7.5 and 8.6. FlauBERTbase gets the worst F1 score (89.0) while CamemBERTbase, Wiki 4 Gb, pre-trained only on 4 Gb of Wikipedia, gets the best F1 score (90.0) with a CER score of 8.4. Considering the CER, as the official metric of the ME-DIA task, the best model is CamemBERTbase, CCNet 135 Gb, a base CamemBERT model pre-trained with CCNet on 135 Gb of text. It obtains the lowest CER, at 7.5 for an F1 at 89.9.

The multilingual models (mBERTbase, distill-mBERTbase, small-mBERTbase-fr, XLM-Rbase and XLM-Rlarge) obtained CER scores ranging from 10.1 to 8.0, for the worst one, distill-mBERTbase, and the best one XLM-Rlarge, respectively. Regarding the F1 scores, they are comparable to the French monolingual models, for XLM-R models with 89.5 for the base one and 89.9 for the large one. On the other side, multilingual BERT and its distilled and small versions obtained performances comparable to the FlauBERTbase, with 88.9 F1 at best.

ATIS-FR results

The performance of the monolingual models on the ATIS-FR task in terms of F1 varies between 92.5 and 94.1 and between 3.3 and 5.3 of CER (

Analysis

In deph analysis

Benchmarking Transformer-based models on the two French SLU tasks (MEDIA & ATIS-FR) allows to observe some trends.

In both tasks, the CamemBERT models perform the best, in terms of Concept Error Rate (CER). FrALBERT obtains comparable results to CamemBERTbase, Wiki 4 Gb, this may come from that the two models are trained on the same kind and amount of data (Wikipedia 4G). We also observe that FrALBERT has comparable performances to XLM-Rbase in both tasks even if the training data and structure are both different. In addition, we notice the underachievement of distill-mBERTbase compared to other BERT models and especially to FrALBERT. Finally, FlauBERTbase has comparable results to CamemBERT models in the MEDIA task, but the CER score of FlauBERT in the ATIS-FR task, is significantly worse than the CamemBERT models, which led us to go deeper in the analysis.

We propose to focus the analysis on the most representative labels of each task. For MEDIA we focus on The first trend we observed is an underperformance of the distill-mBERTbase on named entity tags on both tasks. For instance, STATE NAME, CITY NAME in ATIS or NOM-HOTEL, LOCALISATION-VILLE for MEDIA have up to 2 F1 points less than the other models. Note that these tags are ones of the most frequent in both tasks. On the other side, we could not observe such a big trend by comparing larger models and compact models, or multilingual models versus French monolingual models. Differences between models are very small, for instance the CamemBERTbase, Wiki 4 Gb model will be better than the FrALBERTbase, Wiki 4G on the MEDIA labels NOMBRE-CHAMBRE (93,14 versus 91,85 of F1 respectively), NOM-HOTEL (82,76 versus 79,60 of F1 respectively).

When diving deeper in the MEDIA task, we can detect some small trends when we look at the whole results. For in- stance, even if the OBJET label is one of the most frequent, it seems that all models have difficulties to correctly detect it. In the same way, the label NOM is highly difficult, even for one of the best model (CamemBERTbase, CCNet 135 Gb). It seems that the biggest difference between models occurs in their ability to transfer their knowledge according to the amount of data used for pretraining. In this way, the models trained with the biggest corpora are enable to handle the best, the most infrequent labels. What we observe in the MEDIA task can also be observed in the ATIS-FR task, even if the overall performance of the models is higher. This leads us to suggest that a better sample of the training data could bring interesting results. Moreover, recent works in few-shot learning for slot-filling approaches [START_REF] Cattan | On the cross-lingual transferability of multilingual prototypical models across NLU tasks[END_REF] enable models to perform at high levels with a strong subsampling of examples, which in this case could be a plus.

Ecological and computational costs

We conducted an impact study of Transformer models we considered by measuring the ecological impact of these models [START_REF] Henderson | Towards the systematic reporting of the energy and carbon footprints of machine learning[END_REF]. The table 3 presents the time to process one fine-tuning epoch, and one inference step, associated to the energy consumed and the CO2 produced, for each task. Each score is the mean of five runs for each step.

CamemBERTlarge, CCNet 135 Gb and XLM-Rlarge produce the most over all tasks and steps.

FlauBERTbase, CamemBERTbase, CCNet 135 Gb and mBERTbase produce nearly the same amount of CO2 during fine-tuning and inference steps and for the both tasks (MEDIA and ATIS-FR), which lead us to correlate the amount of parameters and their impact in terms of energy consumed and CO2 produced.

FrALBERTbase, Wiki 4 Gb is the model which uses the less energy (in both fine-tuning and inference steps), and produces the less CO2 of all models. Especially, in the inference step, the FrALBERT model consumes nearly the half of other compact models: distill-mBERTbase and small-mBERTbase-fr.

This study shows compact models like FrALBERT, which obtained comparable results to bigger models (especially with large BERT models in Table 2), have also an important lower ecological impact than big Transformer models. For instance, the impact of FrALBERT is more than 5 time lower than the impact of large models such as XML-Rlarge and CamemBERTlarge in the inference step for both tasks and in the fine-tuning step, the FrALBERT model is more than 3 times lower than the XML-Rlarge model. Surprisingly, the compression approaches of compact models may not have the same impact. We can observe that distill-mBERTbase is a little bit different from small-mBERTbase with a resource requirement and thus an impact a little bit lower but still quite comparable to the mBERTbase model, but with significantly lower performances, compared to mBERT.

Conclusion

Transformer-based architectures are currently the state-of-theart model for many NLP tasks. In this study we have proposed to benchmark the existing French and multilingual Transformer-based pre-trained models, for the purpose of comparing their performance on two French SLU corpora: ME-DIA and ATIS-FR. We have also assessed the ecological impact related to model compacteness on SLU performances and conducted an extensive side-by-side comparison of thirteen recently proposed efficient Transformer models.

The experimental results show that these tasks are very challenging even for large models. Experimental results show that for both tasks, the CamemBERT models perform the best, in term of Concept Error Rate (CER). The best CER results on MEDIA and ATIS-FR are respectively achieved by CamemBERTbase, CCNet 135 Gb and CamemBERTlarge, CCNet 135 Gb. In both tasks, the French compact model FrALBERT obtains comparable results to the large model CamemBERT (base and Wiki 4 Gb configuration). Moreover, this model achieves comparable performances to multilingual models in both tasks and outperforms the distilled version.

Then, our detailed analysis of F1 scores provides interesting model insights in each task. From those analyses we observe that the biggest difference between models occurs in their ability to transfer their knowledge according to the amount of pre-trained data.

Finally, the ecological study conducted on these models shows that compact models have significantly less ecological impact compared to big Transformer models in both fine-tuning and inference steps. This less CO2 is a plus in a context of reducing our ecological impact but it also means less energy consumed, which can reach more than 5 time less between the FrALBERT model and large BERT models.

We plan to open source our code and benchmarks to facilitate future benchmarking, research and model development.

  between 70.15 and 95.76. In the ATIS task, we focus on CITY NAME, AIRLINE NAME, DE- PART TIME.PERIOD OF DAY,DEPART DATE.DAY NAME, TOLOC.CITY NAME, FROMLOC.CITY NAME labels and their F1 varies between 59.41 and 100.

Table 1 :

 1 Model characteristics, in terms of pre-training settings, sources of data and models size.

	Model	Objectives	Data	Vocabulary size Tokenization	# parameters Model size
	FlauBERTbase	MLM	24 French subcorpora (71 Gb of text)	68,729 BPE	138 M	553 Mb
	CamemBERTbase	MLM	French OSCAR (138 Gb of text)	32,005 SentencePiece	110 M	445 Mb
	CamemBERTbase	MLM	French CCNet (135 Gb of text)	32,005 SentencePiece	110 M	445 Mb
	CamemBERTbase	MLM	French OSCAR (4 Gb of text)	32,005 SentencePiece	110 M	445 Mb
	CamemBERTbase	MLM	French CCNet (4 Gb of text)	32,005 SentencePiece	110 M	445 Mb
	CamemBERTbase	MLM	French Wikipedia (4 Gb of text)	32,005 SentencePiece	110 M	445 Mb
	CamemBERTlarge	MLM	French CCNet (135 Gb of text)	32,005 SentencePiece	335 M	1.35 Gb
	FrALBERTbase	MLM and SOP French Wikipedia (4 Gb of text)	32,005 SentencePiece	12 M	50 Mb
	XLM-Rbase	MLM	CC-100 (2.5 Tb of text)	250,002 BPE	278 M	1.12 Gb
	XLM-Rlarge	MLM	CC-100 (2.5 Tb of text)	250,002 BPE	559 M	1.24 Gb
	mBERTbase	MLM and NSP Wiki-100	119,547 WordPiece	177 M	714 Mb
	small-mBERTbase FR MLM and NSP Wiki-100	24,495 WordPiece	104 M	420 Mb
	distil-mBERTbase	MLM and NSP Wiki-100	119,547 WordPiece	134 M	542 Mb

Table 2 )

 2 . FlauBERT gets the worst F1 and CER scores while the large version of the Camem-BERT model gets the best F1 and CER scores. FrALBERT and CamemBERTbase, Wiki 4 Gb obtain similar performance with 0.3 points of F1 and 0.1 point of CER difference. The best model is CamemBERTlarge, CCNet 135 Gb, which obtains the lowest CER, at 3.3 for an F1 at 94.1. The F1 scores vary from 88.1 (CER at 6.0) for the distilled version of mBERT to 93.6 (CER at 5.0) for mBERTbase. The two versions of XLM-R obtained comparable results in terms of CER and F1 scores to CamemBERTbase, Wiki 4 Gb and FrALBERTbase, Wiki 4 Gb.

	MEDIA	ATIS-FR

* 87.5 * 10.1 * 88.1 *

Table 2 :

 2 SLU performances on the MEDIA and ATIS-FR test dataset. Results are given in terms of F-measure (F1) and Concept Error Rate (CER). Significant results are marked with a star.

Table 3 :

 3 Estimation of fine-tuning and inference costs on MEDIA and ATIS-FR corpora.

This work has been made possible thanks to the Saclay-IA computing platform

We did not include the small CamemBERT models[START_REF] Micheli | On the importance of pre-training data volume for compact language models[END_REF] and the large version of FlauBERT, the former not being available and the latter not converging in our experiments despite our efforts.

MEDIA is available for academic use: https://catalogue.elra.info/ELRA-S0272.


Antoine Le Ravallec
email: antoine.le-ravallec@grenoble-inp.fr

João Patrice Garcia

João Carlos

Azevedo Goncalves

Loïc Vincent

Jean-Marc Duchamp

Philippe Benech

Carlos Azevedo

Azevedo Gonçalves

A Compact 28-nm FD-SOI CMOS 76-81 GHz Automotive

Keywords: automotive, radar, mm-Wave, receiver, phase control, CMOS, varactor, SOI

I. INTRODUCTION

First embedded automotive radars appeared in the 1970s to reduce road traffic crashes [START_REF] Meinel | Evolving Automotive Radar -from the very beginnings into the future[END_REF]. Automotive radars are a key factor for Advanced Driver Assistance Systems (ADAS) to provide safer driving, whatever the weather and the road conditions. Therefore, robustness, reliability, long range and peripheral covering with high resolution and accuracy are required for automotive radars [START_REF] Strohm | Development of Future Short Range Radar Technology[END_REF], [START_REF] Murad | Requirements for next generation automotive radars[END_REF].

The cost reduction and the high level of integrability of nanometer CMOS technology are attractive for automotive radars due to the increasing number of sensors inside vehicles. However, due to wafer's variability in advanced CMOS SOI nodes with Process, Voltage and Temperature variations (PVT) [START_REF] Ioannidis | Statistical analysis of dynamic variability in 28nm FD-SOI MOSFETs[END_REF], amplitude and phase calibrations are mandatory.

Radar systems determine the range, the speed and the Angle of Arrival (AoA) of a target. The focus of this work is on the accuracy of the AoA determination. The azimuth, or the AoA, of the reflected signal from a target is determined by the phase difference between each receiver chain signal in phased array systems. The intrinsic phase shift difference between receiver chains is a source of inaccuracy. Thus, for an accurate AoA determination, the phase variations must be controlled in the receiver circuit. Analog solutions to calibrate the phase are usually done with phase-shifters, however active phaseshifters are power consuming and passive ones are lossy and require significant area [START_REF] Kim | A 4-Bit Passive Phase-Shifter for Automotive Radar Application in 0.13 µm CMOS[END_REF]. In addition, it is very complex to achieve an accuracy of less than one degree with phaseshifters at these frequencies [START_REF] Pepe | [END_REF]. Another method is to use an This paper presents a new accurate phase control technique directly included in the front-end LNA. A compact 76-81 GHz receiver has been achieved using a 28-nm FD-SOI CMOS technology. Section II describes the phase control technique and the receiver architecture, and Section III presents the measurement results. Finally, Section IV concludes this work.

II. RECEIVER ARCHITECTURE

Frequency Modulated Continuous Wave (FMCW) radar principle and Direct Conversion architecture are now commonly used for automotive radar applications. FMCW principle simultaneously provides the range and the speed using Doppler effect in a short computing time. Fig. 1 shows the proposed receiver architecture including a phase control solution embedded inside the LNA. The receiver is composed of a two-stage LNA, an LO driver, a passive mixer with its load also used as a low-pass filter (LPF) and a baseband amplifier (BB). These blocks are described in the following sub-sections.

A. Phase control technique

The front-end blocks are the most susceptible to cause phase variations. Monte Carlo simulations of 2000 trials with mismatch and process variations on the front-end indicate phase variations of ± 5° per amplifier stage. Thus, the proposed solution consists in correcting the phase deviations directly inside the LNA. This technique neither requires a significant silicon area nor power consumption to control the receiver phase variations. The solution includes an LNA with varactors used as variable neutralizing cross-capacitors to tune the phase, shown in Fig. 2 (a). The FD-SOI process allows threshold voltage (Vth) modulation through body biasing. Its back-gate access provides each varactor two different voltage controls, which gives both large and fine tuning of its capacitance. Therefore, FD-SOI technology confers a better tuning range than bulk technology.

B. LNA

The topology of the LNA is a differential two-stage common source with capacitive neutralization and resistive feedback, shown in Fig. 2 (a). The differential topology theoretically cancels common mode noise and common mode voltage. Two stages are used as a trade-off between gain and linearity. Both stages are identical to use the same varactors and therefore to have an identical phase tuning at each stage.

To compensate Miller effect and boost the gain, crosscoupling capacitors technique is used. The capacitor 𝐶 𝑛 in Fig. 2 (a) is a MOS varactor designed to tune the capacitance value used to neutralize the gate-drain capacitance 𝐶 𝐺𝐷 of M1-4 transistors and therefore to adjust the phase of the amplifier. The 𝐶 𝐺𝐷 capacitance has a strong impact on the phase. This technique also impacts the gain and the amplifier stability. Varactors offer a non-negligible degree of control on the LNA. Resistive feedback is used to reduce input impedance for the matching and stabilize the amplifier at lower frequencies. The costs of using varactors and resistive feedback are around 1 dB of gain reduction and around 0.8 dB NF degradation.

C. Mixer

The mixer represented in Fig. 2 The mixer is loaded by a 100 Ω differential resistor to make a trade-off between linearity and voltage conversion gain. Capacitors were added to the resistive load to filter unwanted RF signals, as a low-pass filter. A common mode voltage 𝑉 𝐶𝑀 of 450 mV is applied to supply the following baseband amplifier.

D. LO Driver

The LO driver is used as a buffer for the LO signal. It is the same amplifier topology than the LNA in Fig. 2 (a) but with only one stage and a fixed varactors capacitance value. As the LO signal is sinusoidal, there is a short time where both transistors of the same side of the mixer are ON, which degrades the NF of the mixer. Thus, LO signal must have a large voltage swing to shorten the common ON time. LO signal symmetry is important for NF, ICP1dB and conversion gain performances of the mixer. The input and output baluns of the LO driver are carefully designed to reduce asymmetries.

E. Baseband Amplifier

The baseband amplifier is a classic common source differential amplifier with a source follower stage, shown in Fig. 2 (c). This amplifier is designed for LNA-Mixer testability to compensate for the passive mixer losses and to improve NF measurements. The BB amplifier only consumes 3 mW from a 0.9 V supply voltage.

F. MOS varactors

The varactor 𝐶 𝑛 inside the LNA in Fig. 2 (a), used to tune the phase, is composed of two NMOS transistors like in [START_REF] Zhang | A wide range 60 GHz VCO using back-gate controlled varactor in 22 nm FDSOI technology[END_REF] with their drain and source nodes connected to a common node arbitrarily named 𝑉 𝑐 . Their back-gate accesses are also connected to a common node arbitrarily named 𝑉 𝑏 . Thus, there are two voltage control nodes for a more accurate capacitance tuning and therefore a more accurate phase tuning. The backgate control node allows a very fine tuning. 𝑉 𝑐 and 𝑉 𝑏 voltage control ranges are respectively 0 to 0.9 V and 0 to 1.8 V. The MOS transistors in varactor 𝐶 𝑛 are matched to the LNA MOS transistors to present the same PVT variations. 𝐶 𝑚𝑎𝑥 and 𝐶 𝑚𝑖𝑛 values are carefully chosen to ensure LNA stability and to offer the maximum phase control.

III. EXPERIMENTAL RESULTS

The proposed receiver was fabricated in 28-nm FD-SOI CMOS technology from STMicroelectronics. The active area of the receiver depicted in Fig. 3 is only 0.057 mm², including decoupling capacitors.

Measurements were realized at 25°C. The LNA was measured in standalone and the receiver path was measured in two different configurations as illustrated in Fig. 1: LNA-Mixer measurements were done at the outputs of the mixer load, named MixOUT, and LNA-Mixer-Baseband Amplifier measurements at IFOUT outputs. The receiver measurements were performed at 10 MHz intermediate frequency (IF).

Wafer level measurements in Fig. 4 on the two-stage LNA report 5.5 dB NF at 81 GHz with only 0.9 dB of variation within the 76-81 GHz frequency band. A very good correlation is observed between the simulated NF and the measured NF on four different dies, with 0.2 dB measurement accuracy. Depending on the varactors tuning, measurements of the LNA also exhibit 17.1 dB gain and an input compression point of -11.3 dBm in the 76-81 GHz frequency band. It consumes 29.4 mW from 0.9 V supply voltage and has a 1.8 V back-gate biasing.

The receiver consumes 41.7 mW, including the LO driver and the baseband amplifier. Voltage gain measurements illustrated in Fig. 5 (a) indicates that the LNA-Mixer of the receiver achieves a voltage conversion gain GC of 12.5 dB at 78.5 GHz and an ICP1dB of -12.6 dBm. A good correlation between the simulated and measured conversion gains is observed as well. The 1 dB difference is mainly explained by the measurement accuracy at these frequencies. The LO-to-RF isolation is 39.9 dB and both RF and LO input reflection coefficients are below -14.7 dB in the 76-81 GHz frequency band as illustrated in Fig. 5 (b).

Due to LO jitter and synchronization issues between the IF output signals and a frequency reference, the phase variations have been measured on the LNA in stand-alone, where the phase control is directly implemented, for an accurate phase In this band, the maximum range of phase control is 22° (at 78.5 GHz) and the minimum is 17.6° (at 81 GHz). At 78.5 GHz, shown in Fig. 6 (b), the varactor's drain/source node 𝑉 𝑐 , used for coarse tuning, has a phase control resolution of 1.5°/100 mV and its back-gate node 𝑉 𝑏 , used for fine tuning, has a phase control resolution of -0.2°/100 mV. The 𝑉 𝑐 = 0 V / 𝑉 𝑏 = 1.8 V configuration in Fig. 6 corresponds to the highest capacitance value and 𝑉 𝑐 = 0.9 V / 𝑉 𝑏 = 0 V to the lowest capacitance value of the varactor. Measurements on four different dies of the same wafer in the 76-81 GHz band indicate very few variations of the phase control resolution: ±0.2° on the coarse tuning and ±0.1° on the fine tuning; and few variations of the maximum phase control range: ±1°. Table 1 shows that the addition of the varactor offers a tuning element with minimum noise, linearity and gain degradation, for a given frequency between 76 and 81 GHz with an IF of 10 MHz.

The front-end receiver exhibits performances in line with the best-in-class state-of-the-art receivers reported in Table 2. The proposed phase control technique achieves a fine 0.2 o phase control resolution while managing a trade-off between linearity, power consumption and occupied area, which are key factors in automotive radar applications. IV. CONCLUSION A compact 28-nm FD-SOI CMOS receiver path for 76-81 GHz automotive radar applications was presented. A new compact analog phase control solution has been implemented in the receiver through the LNA using MOS varactors as variable neutrodyne capacitors. This solution exhibits a 22° phase control range and a resolution accuracy of 0.2° at 78.5 GHz, which compensates for phase variations from mismatch and process variations, simulated at ± 5° per amplifier stage. Phase measurements over four dies indicate a good reproducibility of this solution. The receiver front-end LNA achieves a minimum 5.5 dB NF at 81 GHz. NF measurements among different dies also point out the design robustness. The receiver path presents a 24.9 dB voltage conversion gain GC with -25.2 dBm ICP1dB, a power consumption of 41.7 mW and an active area of only 0.057 mm². The receiver path without the baseband amplifier (LNA-Mixer only) presents a 12.5 dB voltage conversion gain GC with -12.6 dBm ICP1dB, a power consumption of 38.7 mW and an active area of 0.046 mm². The demonstrated receiver path shows state-of-the-art performances and includes an accurate phase control solution. This contributes to the reduction of both analog area and power consumption for a precise angle of arrival determination. Indeed, a major part of the phase control is done with compact varactors in the frontend of the receiver. The proposed solution guarantees dies alignment in post-production and therefore improves the production yield of radars circuits.

Fig. 1 .

 1 Fig. 1. Receiver block diagram I/Q architecture but it increases the design complexity and power consumption. Digital solutions using Direction of Arrival (DoA) estimation algorithms such as MUSIC can be very accurate but require high Signal-to-Noise Ratio (SNR > 15 dB) [7]. Furthermore, noise and interferences degrade the accuracy of the DoA estimation: the receiver must be very linear with a low noise contribution.This paper presents a new accurate phase control technique directly included in the front-end LNA. A compact 76-81 GHz receiver has been achieved using a 28-nm FD-SOI CMOS technology. Section II describes the phase control technique and the receiver architecture, and Section III presents the measurement results. Finally, Section IV concludes this work.

Fig. 2 .

 2 Fig. 2. Schematic of the 2-stage LNA with varactors (a), the passive mixer with its load (b) and the baseband amplifier (c)

  (b) is a double-balanced passive mixer. The double-balanced topology provides lower losses and better RF-to-IF and LO-to-IF isolations than unbalanced or single-balanced topologies. Passive mixers are very linear and benefit from low on-state resistance of MOS transistors used as switches. Other advantages of passive mixers are the absence of DC current consumption and a reduced flicker noise. The mixer has a simulated ICP1dB of 5 dBm and estimated conversion losses of 4.5 dB from LNA and LNA-Mixer measurement results in Section III. Increasing M5-8 transistors width reduces their on-state resistance but increases parasitic capacitances. The back-gate biasing is used to reduce threshold voltage of the four mixing transistors M5-8, therefore reducing their on-state resistance and improving both linearity and noise performances.

Fig. 3 .

 3 Fig. 3. Die photograph

Fig. 5 .Fig. 6 .

 56 Fig. 5. Measured LNA-Mixer gain versus RF input power at 78.5 GHz (a) and measured RF and LO input reflection coefficients versus frequency (b)

Table 2 .

 2 Receiver performances comparison

		This work	RFIC 19 [9]	ISSCC 21 [10]	TMTT 13 [11] TMTT 21 [12]	JSSC 17 [13]
	Technology	28-nm FD-SOI CMOS	40-nm CMOS	45-nm CMOS	65-nm CMOS	28-nm FD-SOI CMOS	130-nm SiGe
	Frequency (GHz)		76-81	77-81	76-81	76-77	77	76-81
	Integrated RF & analog blocks	2-stage LNA, passive mixer, LO driver	2-stage LNA, passive mixer, LO driver, BB amplifier	3-stage LNA, mixer, 4-stage VGA, LDO, PLL	2-stage LNA, passive I/Q mixer, VGA	5 stage LNA, sub-harmonic mixer, passive mixer	Mixer, VGA, LPF	Mixer, BB multiplexer, IF buffer
	ICP1dB (dBm)	-12.6	-25.2	-22.3 (VGA off)	-10	-15 5	-23.5	>1
	NF_DSB (dB)	6.1 (sim.), 5.5 (LNA)	6.5 (sim.), 5.5 (LNA)	7.3 (at 77 GHz LO), 6 (at 81 GHz LO)	10.8 4	10	8.2	7-8
	GC (dB)	12.5	24.9	30.8	N.A.	16	33	15
	Phase control / resolution (°)	22 max. / 0.2 1 (at 78.5 GHz)	22 max. / 0.2 1 (at 78.5 GHz)	N.A.	N.A.	N.A.	N.A.	N.A.
	Power (mW)	38.7, 29.4 2	41.7, 29.4 2	143 3	3500	28.2 6	27 7	1800 (1TX/6RX)
	Supply (V)	0.9, 1.8 (body biasing)	1.1 (RF), 1.8 (BB)	1/1.2/1.8	1.2	1	3.3
	Active area (mm²)	0.046	0.057	0.8	22 (analog area)	2.3 (with pads)	0.96	12.7 (2TX/6RX)
	Notes	1 LNA measurements 2 Without LO driver	Packaged die, 3 LNA, mixer, BB amplifier	4RX/3TX chip, 4 Supposed DSB	5 LNA only 6 Without LO blocks	7 Without LO buffer	Wafer level measurements, 2TX/6RX chip

ACKNOWLEDGMENT

The authors would like to thank STMicroelectronics members of Crolles for their valuable comments and support.


Bruno Barracosa

Julien Bect

Héloïse Dutrieux Baraffe

Juliette Morin

Josselin Fournel

Emmanuel Vazquez
email: emmanuel.vazquez@centralesupelec.fr

Bayesian multi-objective optimization for stochastic simulators: an extension of the Pareto Active Learning method

Keywords: Multi-objective simulation optimization, Bayesian optimization, Stochastic simulator, Gaussian process modeling, Kriging

This article focuses on the multi-objective optimization of stochastic simulators with high output variance, where the input space is finite and the objective functions are expensive to evaluate. We rely on Bayesian optimization algorithms, which use probabilistic models to make predictions about the functions to be optimized. The proposed approach is an extension of the Pareto Active Learning (PAL) algorithm for the estimation of Pareto-optimal solutions that makes it suitable for the stochastic setting. We named it Pareto Active Learning for Stochastic Simulators (PALS). The performance of PALS is assessed through numerical experiments over a set of bi-dimensional, bi-objective test problems. PALS exhibits superior performance when compared to a selection of alternative approaches based on scalarization or random-search.

Introduction

We consider the problem of multi-objective optimization with stochastic evaluations, also known as multi-objective simulation optimization [START_REF] Hunter | An introduction to multiobjective simulation optimization[END_REF]. More precisely, let f 1 , . . . f q be q real-valued objective functions defined on a search domain X ⊂ R d , and assume that these objectives can only be observed with some random errors: each evaluation of the objectives at a given point x ∈ X in the search domain yields random responses Z j = f j (x) + ε j , 1 ≤ j ≤ q, where the ε j s are zero-mean random variables (the distribution of which may depend on x), mutually independent and independent from past evaluations. Our objective is to estimate the Pareto-optimal solutions of the problem:

min f 1 , . . . , f q .

(1)

This setting has many applications in the industry when stochastic simulators are used to make decisions about uncertain systems with several, often conflicting, criteria. Examples can be found in all areas of engineering and science, for instance plant breeding [START_REF] Hunter | Maximizing quantitative traits in the mating design problem via simulation-based Pareto estimation[END_REF], aircraft maintenance [START_REF] Mattila | Maintenance scheduling of a fleet of fighter aircraft through multi-objective simulation-optimization[END_REF] or electric network planning [START_REF] Dutrieux | An innovative method to assess solutions for integrating renewable generation into distribution networks over multi-year horizons[END_REF].

In this work, we choose to focus on Bayesian optimization algorithms, which use probabilistic models to make predictions about the functions to be optimized. One of the classical algorithms for Bayesian deterministic multi-objective optimization is the ParEGO algorithm [START_REF] Knowles | ParEGO: A hybrid algorithm with on-line landscape approximation for expensive multiobjective optimization problems[END_REF]. It relies on a scalarization approach to extend the very popular single-objective Efficient Global Optimization (EGO) algorithm [START_REF] Jones | Efficient global optimization of expensive black-box functions[END_REF], based on the Expected Improvement (EI) criterion. Another scalarization approach, called multi-attribute Knowledge Gradient [START_REF] Astudillo | Multi-attribute Bayesian optimization under utility uncertainty[END_REF], uses the Knowledge Gradient (KG) criterion instead of the EI criterion. Several Bayesian multi-objective optimization methods that do not rely on scalarization have been proposed as well. These are, for instance, the well-known Expected Hypervolume Improvement (EHVI) algorithms [START_REF] Emmerich | Single-and multiobjective evolutionary optimization assisted by Gaussian random field metamodels[END_REF][START_REF] Feliot | A Bayesian approach to constrained single-and multi-objective optimization[END_REF], the Stepwise Uncertainty Reduction (SUR) approach suggested in [START_REF] Picheny | Multiobjective optimization using Gaussian process emulators via stepwise uncertainty reduction[END_REF], the Predictive Entropy Search for Multi-objective Bayesian Optimization (PESMO) algorithm [START_REF] Hernández-Lobato | Predictive Entropy Search for Multi-objective Bayesian Optimization with Constraints[END_REF], the Max-value Entropy Search for Multi-objective Bayesian Optimization (MESMO) [START_REF] Belakaria | Max-value entropy search for multi-objective Bayesian optimization with constraints[END_REF], and the Pareto Active Learning (PAL) algorithm [START_REF] Zuluaga | Active learning for multi-objective optimization[END_REF].

In this article, we choose to focus on PAL because it is easy to implement and inexpensive in terms of computational time, contrarily to other Bayesian approaches based on computationally-intensive criteria for point selection. The main idea of PAL consists in using confidence intervals from the posterior distribution of a Gaussian process to balance exploration and exploitation, in a way inspired by the Gaussian Process Upper Confidence Bound rule (GP-UCB) [START_REF] Srinivas | Gaussian process optimization in the bandit setting: no regret and experimental design[END_REF] for single-objective optimization. [START_REF] Zuluaga | Active learning for multi-objective optimization[END_REF] present PAL as a multi-objective Bayesian optimization tool directed at cases where "evaluating the objective function is expensive and noisy". However, as will be discussed in this article, PAL is fundamentally an algorithm designed for the optimization of deterministic objective functions over a finite input space, and significant limitations arise when considering truly stochastic (random) evaluations.

This article brings two contributions: (1) an extension of the PAL algorithm, called PALS, making it suitable for stochastic evaluations, and (2) a numerical benchmark comparing the performance of the proposed approach with several competing approaches based on scalarization and random-search.

The article is organized as follows. Section 2 introduces in more detail the problem of multi-objective optimization for stochastic simulators. Section 3 describes the original PAL algorithm. Section 4 presents the proposed extension for stochastic simulators. Section 5 details the numerical experiments and their main results. Section 6 draws conclusions and suggests several perspectives. A separate supplementary material file is available.

Bayesian multi-objective optimization for stochastic simulators

Consider a black-box stochastic simulator whose outputs are noisy evaluations of q latent functions f 1 , . . . , f q : X → R over a search space X ⊂ R d , d ∈ N. The vector of function values (f 1 (x), . . . , f q (x)) will be denoted by f (x). Suppose that the functions f 1 , . . . , f q represent (possibly conflicting) costs that we want to minimize.

In a multi-objective optimization problem, the goal is to identify the solutions that represent the best possible trade-offs among the q objectives. This is usually defined using the Pareto domination rule: for any two points x, x ∈ R d with images z = (z 1 , . . . , z q ) = f (x), z = (z 1 , . . . , z q ) = f (x ) ∈ R q , we write z ≺ z (z dominates z in the sense of Pareto), when z j ≤ z j for all j, with at least one of the inequalities being strict. The set P of all non-dominated points is called the Pareto set:

P = {x ∈ X : x ∈ X, f (x ) ≺ f (x)} .
(2)

The Pareto front F ⊂ R q is by definition the image f (P) of the Pareto set by the objective functions. Figure 1 illustrates the Pareto front of a finite set of evaluations of two objective functions.

The goal of multi-objective optimization is to obtain a good approximation P n of P using a sequence (X 1 , . . . , X n ), n ≥ 1, of evaluation points and the corresponding evaluation results. (The use of capital letters to denote the evaluation points indicates that these points depend on past evaluation results.)

In this article, we assume a stochastic setting, where evaluations are corrupted by additive, normally distributed, homoscedastic noise. In other words, instead of

z 1 z 2 z 3 z 4 z 5 z 6 f 1 (x) f 2 (x) Figure 1: Example of a Pareto front F = {z 1 , z 2 , z 3 }, in a bi-objective example with finite search domain X = {x 1 , . . . , x 6 }.
observing the values of f 1 , . . . , f q at X 1 , . . . , X n , we observe random responses

Z i,j = f j (X i ) + ε i,j , 1 ≤ i ≤ n, 1 ≤ j ≤ q,
where the random variables ε i,j are mutually independent, with Gaussian distributions N (0, σ 2 j ) for some unknown σ 2 j > 0. Using a Bayesian framework provides a practical solution to handle this stochastic setting. The main idea is to use a prior probability distribution, often simply called the prior, to model each objective f j as a sample path of a random process ξ j . For convenience, we assume moreover that the ξ j s are independent Gaussian processes, since it makes it possible to easily obtain the posterior distributions of the models ξ j conditional on the observations Z i,j . In this case indeed, the ξ j s remain independent and Gaussian a posteriori (i.e., given the evaluation results). We denote by µ n,j and k n,j their posterior mean and covariance functions, which can be derived by solving systems of linear equations (see, e.g. [START_REF] Chiles | Geostatistics: modeling spatial uncertainty[END_REF][START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF].

Posterior distributions based on past observations can be used simultaneously to choose additional evaluation points, and to generate predictions of the Pareto front. Pareto front predictions can be generated through conditional simulations (Figure 2), making it possible not only to estimate the Pareto front, but also to quantify the remaining uncertainty [START_REF] Binois | Quantifying uncertainty on Pareto fronts with Gaussian process conditional simulations[END_REF].

As pointed out in the introduction, several Bayesian multi-objective optimization algorithms have been proposed in the literature, but few of them can actually be used under a stochastic framework. [START_REF] Astudillo | Multi-attribute Bayesian optimization under utility uncertainty[END_REF] propose an algorithm that takes as input a user preference and converges to a particular solution of the Pareto front (i.e., does not provide an estimate of the entire Pareto set and front). The PESMO algorithm [START_REF] Hernández-Lobato | Predictive Entropy Search for Multi-objective Bayesian Optimization with Constraints[END_REF] on the other hand, which has been proposed in the machine learning literature, can be used to jointly estimate the Pareto set and front, as considered in this article. However, this algorithm has a high implementation complexity and a non-negligible execution cost-it relies at each iteration on conditional simulations of the GP models, a genetic algorithm and an expectation-propagation algorithm.

Unlike PESMO, the PAL algorithm proposed by [START_REF] Zuluaga | Active learning for multi-objective optimization[END_REF] is a simple algorithm, but it is not as such suitable for the stochastic case. In the next sections, we introduce the original PAL algorithm, discuss its limitations, and propose modifications to deal with the stochastic case.

Remark 1. It is worth noting that several algorithms from the literature of budget allocation / ranking and selection, namely MOCBA [START_REF] Lee | Finding the nondominated Pareto set for multi-objective simulation models[END_REF], SCORE [START_REF] Pasupathy | Stochastically constrained ranking and selection via score[END_REF], and M-MOBA [START_REF] Branke | A new myopic sequential sampling algorithm for multi-objective problems[END_REF] algorithms, use similar ideas to those outlined above, but with one notable difference: they assume no correlation between the values of the objective functions at distinct points in the search space. Recently, an algorithm using a mix of Bayesian optimization and ranking and selection approaches was proposed by Rojas Gonzalez et al. (2020) and named SK-MOCBA.

The original PAL algorithm

In the following, we assume that X ⊂ R d is a (possibly large) finite set, corresponding for instance to a discrete grid or a large random sample from a uniform distribution on a bounded set. The PAL algorithm builds a sequence (X n ) n≥1 of evaluation points using, at each iteration, a classification of all the points x ∈ X of the search space as potentially Pareto-optimal or not, based on the posterior distributions of the Gaussian processes ξ j . More precisely, the algorithm partitions X into three sets at iteration n: a set P n of points that are deemed Pareto-optimal, a set N n of points that are deemed dominated, and the set U n = X \ (P n ∪ N n ) of points that remain unclassified. Then, the PAL algorithm improves the approximation of the Pareto set by choosing an additional evaluation point in P n ∪ U n , where uncertainty is maximal (in a sense defined below).

To build P n , U n and N n , each x of the search space is associated to a region R n (x) ⊂ R q , n ≥ 0, in the objective space, which quantifies the uncertainty about the values of the objectives at that point. For each n ≥ 0, the region R n (x) is constructed as follows. First, for each x ∈ X and each j ∈ {1, . . . , q}, compute the posterior means µ n,j (x) and posterior variances σ 2 n,j (x) = k n,j (x, x) of the Gaussian random variables ξ j (x) modeling the unknown values of the f j s at x, given observations Z 1,j , . . . , Z n,j (for n = 0, they are equal to the prior means and variances). Then, define hyper-rectangles Q n (x) ⊂ R q capturing the uncertainty of the predictions, as illustrated in Figure 3:

Q n (x) = z ∈ R q : µ n (x) -β 1/2 σ n (x) ≺ z ≺ µ n (x) + β 1/2 σ n (x) , (3)
where µ n (x) = (µ n,1 (x), . . . , µ n,q (x)), σ n (x) = (σ n,1 (x), . . . , σ n,q (x)), and β > 0 is a positive scaling factor. Finally, for each x ∈ X, the regions used for classification are defined as

R n (x) = {µ n (x)} if x ∈ {X 1 , . . . , X n } , R n-1 (x) ∩ Q n (x) otherwise, (4)
with the convention that R -1 (x) = R q . (NB: in the original article by [START_REF] Zuluaga | Active learning for multi-objective optimization[END_REF], a point x ∈ X already visited is artificially assigned a zero posterior variance; here, we prefer to state that the corresponding region R n (x) collapses to µ n (x), which is equivalent.) For each R n (x), define an optimistic outcome R min n (x) and a pessimistic outcome R max n (x) as the best possible outcome and the worst possible outcome within R n (x) according to the Pareto domination rule (see Figure 4). Then, each point x ∈ X is classified in P n , U n or N n according to the following rules illustrated in Figure 5:

P n = x ∈ X x ∈ X, R min n (x ) + ≺ R max n (x) -, (5)
where ∈ R q + is a vector of margin parameters controlling the acceptation in P n and N n . In other words, x ∈ P n when the ( -shifted) pessimistic outcome R max n (x) is not dominated by the ( -shifted) optimistic outcome R min n (x ) of any other point x . Conversely,

f 1 (x) f 2 (x) µ n (x) -β 1/2 σ n (x) µ n (x) + β 1/2 σ n (x) µ n (x) Q n (x)
N n = x ∈ X ∃x ∈ X \ {x} , R max n (x ) -≺ R min n (x) + , (6)
or in other words, x ∈ N n when the ( -shifted) optimistic outcome R min n (x) is dominated by at least one of the the ( -shifted) pessimistic outcome R max n (x ) of another point x ; and finally U n = X \ (P n ∪ N n ).

According to [START_REF] Zuluaga | Active learning for multi-objective optimization[END_REF], the margin parameter provides flexibility around the uncertainty region, and can be defined globally or individually for each objective.

The parameter β plays the role of a scaling factor for the uncertainty regions and is defined at each iteration. [START_REF] Zuluaga | Active learning for multi-objective optimization[END_REF] suggest using an increasing sequence β n defined at each iteration n as

β n = 2 log(q|X|π 2 n 2 /6δ) , (7)
with δ ∈ [0, 1]. The proposed approach allows to ensure an upper bound of the hypervolume error with probability 1 -δ, by specifying δ. Actual convergence and upper bound on the maximal number of needed iterations are presented for the case where intersecting regions are used, and a squared exponential kernel is adopted with parameters that are not re-estimated at each iteration (refer to Corollary 2 in the original paper for details). The last step of the algorithm consists in selecting a new evaluation point X n+1 among the non-visited points in the union P n ∪ U n for which uncertainty is maximal, in the sense of the Euclidean distance between the pessimistic and the f 1 (x) optimistic outcomes (i.e., the diameter of R n (x)):

f 1 (x) f 2 (x) R min n (x) R max n (x) Q n (x) R n-1 (x) R n (x)
f 2 (x) R min n (x) R max n (x) f 1 (x) f 2 (x) R max n (x) R min n (x)
X n+1 = argmax x∈(Un∪Pn)\Sn R min n (x) -R max n (x) 2 , ( 8
)
where S n denotes the set of visited points. Algorithm 1 summarizes the PAL algorithm.

Remark 2. It appears from the right-hand side of (8) that PAL expects the objectives to be suitably normalized, in order to be comparable. Indeed, the criterion at any given point involves a sum of differences of objective values:

R min n (x) -R max n (x) 2 2 = q j=1 R min n,j (x) -R max n,j (x) 2 .
We will assume for simplicity that such a normalization has been made beforehand. When this is not possible, the range of the objective functions can be learned on-the-fly as done, e.g., by [START_REF] Feliot | A Bayesian approach to constrained single-and multi-objective optimization[END_REF].

Remark 3. In practice, in Bayesian optimization, it is advisable to start the optimization procedure using a set of n 0 preliminary evaluations, which make it possible to select the parameters of the Gaussian process model before sequential selection of points. This initial set of results is called "training set" by [START_REF] Zuluaga | Active learning for multi-objective optimization[END_REF], or more traditionally initial design of experiments. In this case, the index n in Algorithm 1 corresponds to the number of simulations performed after the initial design, and the posterior means and variances at iteration n are in fact computed with respect to n + n 0 evaluation results.

PALS: modification of the PAL algorithm for stochastic simulators

The original PAL algorithm presents limitations that make it unsuitable for stochastic evaluations. To circumvent these limitations, we propose several modifications, which can in fact be seen as simplifications of the original algorithm. The simplifications, together with an extension to batch evaluations, are detailed below. As in PAL, it is assumed that X is a finite set. The resulting algorithm, summarized in Algorithm 2, is called PALS (S for Stochastic).

Modification 1: Taking into account prediction variance at visited points. By artificially setting the prediction variance at visited points to zero (equivalently, by collapsing R n (x) to µ n (x) as in (4)), [START_REF] Zuluaga | Active learning for multi-objective optimization[END_REF] prevent the PAL algorithm from visiting any point of the search space more than once. This behavior is also enforced by the explicit requirement (see ( 8)) that X n+1 should belong to the set (P n ∪ U n ) \ S n , where S n denotes the set of visited points. In fact, this approach can be seen as using a "nugget" component [START_REF] Matheron | Principles of geostatistics[END_REF] in the GP prior, which would be appropriate for a very irregular deterministic Algorithm 1: The original PAL algorithm Input : input space X; GP priors; q ; β n Output : predicted Pareto set P

P -1 = ∅, N -1 = ∅, U -1 = X, S -1 = ∅ R -1 (x) = R q for all x ∈ X n = 0, all classified = false repeat Obtain µ n (x) and σ n (x) for all x ∈ X, {µ n (X i ) = (Z i,1 , . . . , Z i,q ) and σ n (X i ) = 0, ∀X i ∈ S n } R n (x) = R n-1 (x) ∩ Q µn,σn,βn (x) for all x ∈ X P n = P n-1 , N n = N n-1 , U n = ∅ forall x ∈ U n-1 do if x ∈ X \ {x} : R min n (x ) + ≺ R max n (x) -then P n = P n ∪ {x} else if ∃x ∈ X \ {x} : R max i (x ) -≺ R min i (x) + then N n = N n ∪ {x} else U n = U n ∪ {x} if U n = ∅ then all classified = true else Choose X n+1 = argmax x∈(Un∪Pn)\Sn R min n (x) -R max n (x) 2 S n+1 = S n ∪ {X n+1 } Sample Z n+1,j = f j (X n+1 ) + ε n+1,j
for all j ∈ {1, . . . , q} n = n + 1 until all classified; return P = P n simulator, but is not for a stochastic simulator-where replicated simulations at a given point yield independent and identically distributed random responses. To circumvent this limitation, we remove the requirement that X n+1 should belong to X \ S n and use the correct posterior variance at all points, visited or not; more explicitly, we replace ( 4) and ( 8) with

R n (x) = R n-1 (x) ∩ Q n (x) . ( 9
)
and

X n+1 = argmax x∈Un∪Pn R min n (x) -R max n (x) 2 . ( 10
)
As a consequence, PALS can visit a given point several times and thus the number of iterations of the algorithm is not bounded as in PAL by the cardinality of the search domain. We introduce a user-defined maximum budget of evaluations n max , which will constitute our main stopping criterion (since it is often extremely costly, in the problems that we consider, to reach a situation where all the points are classified in P n or in N n ). We note, however, that keeping the cardinality of U n as an additional stopping criterion is the only reason why the distinction between U n and P n is still relevant.

Modification 2: Removing the intersection of uncertainty regions. [START_REF] Zuluaga | Active learning for multi-objective optimization[END_REF] define R n at each iteration as the intersection of R n-1 and Q n in (4) in order to get a non-increasing sequences of regions-a monotonicity property which plays a role in the proof of their theoretical results. However, problematic situations can arise as a result of this definition, where Q n (x) is not contained in R n-1 (x) for some x at some iteration n. When this happens, R n (x) is empty and the subsequent behavior of the algorithm is not properly defined since the pessimistic and optimistic outcomes at x, R max n (x) and R min n (x), are not defined. In the PALS algorithm we remove intersections to address this issue, and thus further simplify (9) to

R n (x) = Q n (x). ( 11
)
As a consequence, the subsets P n and N n defined by ( 5)-( 6) loose their monotonicity property as well, and thus membership of all the points in X to P n , N n or U n must be re-evaluated at each iteration.

Remark 4. As an alternative, it is possible to keep the original idea of using intersections and simply "correct" the sets R n (x) so that R n (x) always contains µ n (x). This can be achieved by taking R n (x) to be the smallest hyperrectangle that contains both R n-1 ∩ Q n (x) and {µ n (x)}. Note that doing so artificially fixes the problem of empty regions, but nonetheless destroys the monotonicity property, which was the original motivation of [START_REF] Zuluaga | Active learning for multi-objective optimization[END_REF] for introducing intersections. For the sake of completeness, this approach will also be assessed in our experimental section.

Modification 3: Batches of simulations. When the output variability is large, a single simulation provides a very noisy evaluation of the quantities of interest, which are the expected responses f j of the simulator, and therefore yields by itself little progress in the optimization procedure. A natural idea to gain more information, when parallel computing resources are available, is to perform several simulations at the same point at each iteration of the algorithm. An approach of this type has been proposed for instance by [START_REF] Binois | Replication or exploration: Sequential design for stochastic simulation experiments[END_REF].

We denote by k the size the simulation batches, i.e., the number of simulations performed at the same point at each iteration.

The number of simulation results available after n iterations with batch size k is equal to N = nk, and thus grows quickly if k is large (for instance k = 200 or more in Section 5). Such a large number of results poses apparently a challenge for Gaussian process regression, whose computational complexity is O(N 3 ). However, this complexity can be reduced, without any approximation, since Gaussian process regression with Gaussian noise only uses the empirical mean and variance of the simulation results at each point x ∈ X where at least one simulation was made (see, e.g., [START_REF] Dutrieux | The informational approach to global optimization in presence of very noisy evaluation results. Application to the optimization of renewable energy integration strategies[END_REF][START_REF] Binois | Practical heteroscedastic Gaussian process modeling for large simulation experiments[END_REF]: a careful implementation leveraging this fact has a complexity that scales cubically with the number of "visited" points, which does not grow with k.

Remark 5. Another approach, not pursued in this article, would consist in sampling simultaneously multiple locations at each iteration, as proposed by [START_REF] Habib | A multi-objective batch infill strategy for efficient global optimization[END_REF], for instance.

Predicting the Pareto front and set. PAL uses the condition U n = ∅ as a stopping criterion, and then estimates the unknown Pareto set using the set P n obtained at the final iteration. With the introduction of a budget-based stopping condition, however, PALS typically stops before all points are classified. A first approach to estimate the Pareto front and set consists in applying Pareto domination rule directly to the posterior means of the GP models, in a "plug-in" manner. This approach, however, does not take the residual uncertainty into account. An alternative approach would be to generate conditional sample paths of the Gaussian processes instead to represent possible realizations of the objective functions, and then use them to estimate for each point in the objective space the probabilities of that point being dominated (also called empirical attainment function in [START_REF] Binois | Quantifying uncertainty on Pareto fronts with Gaussian process conditional simulations[END_REF]), and for each point in the input space, the probabilities of belonging to the Pareto set (also called coverage probabilities). An illustration is provided for the Pareto front and set in Figure 6. It is worth noting that a high confidence about the Pareto front location does not imply a high confidence about the location of the Pareto set, and vice versa.

Although providing enriched information about the Pareto front and set, this approach entails a higher computational cost. Therefore, for the remaining sections, the plug-in approach is used to compute estimates of the Pareto front and the Pareto set (at the final iteration, but also at each iteration for the EI as an infill criterion which is suitable only for deterministic observations. The first version, hereafter named ParEGO-EI m , uses the EI m criterion defined in [START_REF] Vazquez | Global optimization based on noisy evaluations: An empirical study of two statistical approaches[END_REF], which considers the Expected Improvement with respect to the best prediction over the visited points (as seen used in the implementation in the PESM branch of Spearmint (Spearmint, 2016)). The second one, hereafter named ParEGO-KG, uses the Knowledge Gradient (KG) sampling criterion instead, which is suitable for a stochastic setting [START_REF] Frazier | The knowledge-gradient policy for correlated normal beliefs[END_REF].

For the PALS approach, when an increasing β is considered (cf. Equation ( 7)), the parameter δ is taken equal to 0.05 as in the original article. We consider = 0 for all experiments, given that this parameter has not shown a significant effect in our numerical experiments (details available in the supplementary material).

Numerical experiments are performed over a set of nine bi-objective and bi-dimensional test problems subject to additive homoscedastic Gaussian noise and defined over a finite 21 × 21 input space (see Appendix A). A total of 200 optimization runs are performed for each problem and optimization method.

All the methods except PRS rely on GP modeling: the GP priors on the two objectives are independent, with a constant unknown mean and a Matérn 5/2 covariance function. The unknown mean is integrated out using an improper uniform prior ("ordinary kriging"), and the parameters of the covariance function are estimated at each iteration by the Restricted Maximum Likelihood (ReML) method (see, e.g., [START_REF] Santner | The Design and Analysis of Computer Experiments[END_REF].

Parameters are initialized based on an initial design which is built by choosing among 1000 initializations of 20 randomly selected points so as to maximize the minimal distance among them. The initial observations are obtained by performing 10 evaluations at each initial point. A total budget n max of 50000 evaluations is allowed at each experiment in addition to the initial 20 × 10 = 200 evaluations. At each iteration, a batch of k evaluations is performed at the selected point. The Pareto set and front predictions are generated, for all methods, using the plug-in approach described in Section 4.

The performance of the methods is assessed using two metrics, averaged over the 200 runs: 1) the volume of the symmetric difference between the predicted dominated region and the true dominated region, and 2) the misclassification rate, meaning the percentage of points incorrectly classified in the Pareto set prediction. See Appendix B for details. These metrics provide quantitative information on the error between the prediction and the correct Pareto set or front.

Unless otherwise mentioned, the setup used for experiments is as follows: no intersection considered, a constant value for β corresponding to the 0.50 probability interval, a null value for , and a value of 200 for k.

Results

Each of the numerical experiments presented in this section has been carried out on the nine test problems. For brevity, only some representative results are presented; the remaining results are available in the supplementary material. Influence of the parameter β of PALS. We are interested in testing the usefulness of an increasing β value for the algorithm's performance. For this purpose, several constant β values were chosen to generate intervals with different coverage probabilities p. More precisely, in this context, we define the coverage probability as the probability that, for any particular objective j, at any point x ∈ X, the function value lies within the region µ j (x) ± β 1/2 σ j (x). Since the posterior is Gaussian, this definition implies that √ β = Φ -1 (0.5 + 0.5p), where Φ denotes the standard normal cumulative distribution function. These values are compared with the increasing-β approach proposed in the original PAL algorithm (cf. Equation ( 7)).

The choice of β impacts the size of the uncertainty region around the prediction. A small β allows the algorithm to focus only on the most promising regions, by quickly excluding those that seem less interesting. On the other hand, a large β favors exploration by taking longer to exclude points from being visited.

In terms of algorithm performance, an example is presented in Figure 7 for problem g 2 , for the volume of the symmetric difference error metric, comparing the increasing-β approach with fixed β for coverage values of 0.10, 0.25, 0.50, 0.75, 0.90, and 0.99.

When looking at the other problem results (as presented in the supplementary material), the choice of the ideal β seems to be problem-dependent and varies depending on whether the metric used focuses on the quality of Pareto front or Pareto set predictions. However, constant values corresponding to 0.50 and 0.75 probability regions have shown an overall interesting performance.

Additionally, experiments were performed to assess the impact of the noise level on the choice of the ideal β. Figure 8 presents the same problem as Figure 7 but with two different noise levels. In our benchmark, reduced noise levels seem to render the β choice less important given the similar performance observed among choices. Small coverage levels (10% or 25%), although promising to accelerate classification, should be avoided given the low performance when considering the Pareto front prediction quality metric. Coverage levels between 50% and 90% have shown satisfying results across a multitude of problems and noise levels for both metrics. The increasing-β approach (in the form proposed by the PAL authors) brings no significant advantage when compared to a fixed β for the classification error metric, but usually underperforms for the volume of the symmetric difference metric.

Considering the results obtained, a β with constant value corresponding to the 0.50 probability interval is used in PALS in the following experiments. About intersections. In the PALS algorithm, intersecting hyper-regions does not seem to accelerate or improve the results. Figure 9 presents a comparison between the non-intersecting and intersecting variants of PALS by showing for the different problems the average performance metrics of using (or not) intersections, normalized by the reference PRS performance. It turns out that average results in the middle and final iterations are very similar, with a slightly better performance of PALS without intersection for the misclassification rate metric. Detailed results can be consulted in supplementary material. Given these results, only the non-intersecting version of PALS will be considered for the remainder of the article.

Influence of the batch size. When using a large budget of evaluations, the results show no significant impact of the batch size k in the average metrics obtained at the final iteration. However, it is worth noting that a smaller batch size for most problems allowed a faster reduction in metric value at initial iterations. Figure 10 depicts the comparison between the average metrics for problem g 5 for k values of 200, 500, 1000 and 2000. This acceleration of the initial convergence can also be observed for the other problems on the plots available in the supplementary material. Given the results, k = 200 will be used for the remainder of the article.

Comparison with other approaches. This section compares the results of the different approaches: PRS, an alternative random search approach (we will name it "Concentrated Random Sampling" (CoRS)), two modified ParEGO approaches, and PALS. For PALS, the default parameters are used: a β with constant value corresponding to the 0.50 probability interval, an of 0, no intersections, and a k value of 200.

Figure 11 presents the evolution of the performance metrics on problem g 8 for the considered methods. Performance illustration for the remaining problems is available in the supplementary material. We can observe that PALS presents a good performance level on both metrics, allowing for better Pareto set and front estimations in most cases. A summary of the numerical results obtained when comparing the proposed PALS approach and the reference methods over the set of benchmark problems is presented in Table 1. PALS consistently outperforms the other methods in terms of the ability to better estimate the Pareto set (as seen through the misclassification rate metric). For Pareto front estimation (see volume of the symmetric difference metric), the most efficient algorithm depends on the problem considered, with PALS presenting an overall very good performance-although not being the best performer for problems g 7 and g 9 , it places second with a metric value close to the leader. Additionally, PALS performs significantly better than pure random search on both metrics for all problems.

Table 1: Average metrics (value in percentage), for the volume of the symmetric difference (V d ) and misclassification rate (M ) at final iteration, for PRS, CoRS, ParEGO-EIm, ParEGO-KG, and PALS. The best metric values are highlighted in bold with a gray background. Metrics at 10% of the best metric are highlighted with a light gray background and considered as acceptable performance.

PRS

CoRS ParEGO-EIm ParEGO-KG PALS

g V d M V d M V d M V d M V d M g 1 0.

Conclusions and future work

In this article, an extension of the Pareto Active Learning (PAL) algorithm was presented, which makes it suitable for stochastic simulators with possibly high output variability. This new version was named Pareto Active Learning for Stochastic Simulators (PALS) and showed promising results in a numerical benchmark on nine bi-objective test problems in dimension two. Performance was compared with a Pure Random Search (PRS) approach, a Concentrated Random Sampling (CoRS) approach and two modified ParEGO algorithms. For the considered benchmark, PALS is the best algorithm to estimate the Pareto set solutions of the problem. It is also the best performer for Pareto front estimation for most of the problems. The modified ParEGO algorithm named ParEGO-KG also presents promising results when considering the Pareto front estimation.

We believe that, for future work, comparing the proposed algorithm with alternative approaches such as SK-MOCBA [START_REF] Rojas Gonzalez | A multiobjective stochastic simulation optimization algorithm[END_REF] or PESMO [START_REF] Hernández-Lobato | Predictive Entropy Search for Multi-objective Bayesian Optimization with Constraints[END_REF] should be envisaged. Additionally, we could be interested in comparisons with algorithms from the Multi-objective Evolutionary Algorithms (MOEA) literature. Interesting perspectives could also include creating hybrid algorithms that could leverage the different advantages of each algorithm.

A difference in PALS with regards to PAL is that each point selected during the optimization process is evaluated several times. This is especially advantageous when multiple parallel computing units are available or when the infill criterion to choose the next point to evaluate is expensive to compute with regards to the actual simulation time. Although the numerical experiments have shown no significant difference at final iteration between different batch sizes, the intermediate performances were different. Results show that smaller batch sizes have better performance at initial stages, which seems to indicate a possible interest in having a noise-level dependent adaptive batch size [START_REF] Lyu | Adaptive batching for Gaussian process surrogates with application in noisy level set estimation[END_REF] during the optimization process.

A main component in the PAL approach is the β parameter that defines the size of the uncertainty region considered around the model's prediction. In the original work, an increasing β is proposed. Large parameter values tend to promote exploration, while smaller values promote exploitation. Numerical experiments have shown, for the problems considered, that selecting a constant value for the parameter could provide better results than the original approach. Nevertheless, the concept of using an increasing β parameter seems relevant, but further studies should be performed to determine a suitable initial value and rate of change depending on the problem structure. Therefore, the ideal β selection seems to depend on the problem and noise level, and should still be further studied.

Finally, even though PALS is an algorithm suitable for the case where the input space is finite, we believe that the main idea could be applied-with some adaptations-to continuous input spaces.

Figure 3 :

 3 Figure 3: Bi-objective example of a prediction uncertainty region Qn(x).

Figure 4 :

 4 Figure 4: Bi-objective example of the construction of a classification region Rn(x) (filled region), and the respective optimistic R min n (x) and pessimistic R max n (x) outcomes.

Figure 5 :

 5 Figure 5: Illustration of the classification of a Pareto-optimal point (left) presented in light gray, and the classification of a non-Pareto-optimal point (right) presented in dark gray, while unclassified points remain as non-filled square dashed outlines.

Figure 6 :

 6 Figure 6: Bi-objective and bi-dimensional example of a normalized heat map representing the probability of a point being dominated for the objective space (top) and the probability of belonging to the Pareto set (bottom) constructed from a probabilistic model (built from the same data and model as in Figure 2, using 40 sample paths). Level curves shown for 5%, 50%, and 95% levels.

Figure 7 :

 7 Figure 7: Comparison of average volume of the symmetric difference for problem g 2 , between the increasing-β approach and fixed β for coverage values of 0.10, 0.25, 0.50, 0.75, 0.90, and 0.99.

Figure 8 :

 8 Figure 8: Average volume of the symmetric difference for problem g 2 for different β values, with the problem noise standard deviation divided by 10 (top) and by 100 (bottom)

Figure 9 :

 9 Figure 9: Performance comparison between PALS algorithm with and without intersection, normalized by the reference PRS performance. Average performance, for each problem for the volume of the symmetric difference (top) and misclassification rate (bottom), considering 200 random initialization, obtained after a budget of 25000 evaluations (circles) and 50000 evaluations (diamonds).

Figure 10 :Figure 11 :

 1011 Figure 10: Average performance metrics volume of the symmetric difference (top) and misclassification rate (bottom) for problem g 5 , when considering k of 200, 500, 1000 and 2000.

Figure A. 1 :

 1 Figure A.1: Pareto front solutions for test problems.

Table A .

 A 1: Definition of functions f 1 to f 15 -a * x 1 . * x 1 + b * x 1 -6). 2 + c * cos(x 1 ) + 10 f 5 100(x 2 -x 2 1 ) 2 + (1 -x 1 ) 2 f 6 to f 15 c 1 + c 2 x 1 + c 3 x 2 + c 4 x 1 x 2 + c 5 x 2 1 + c 6 x 2 2 + c 7 x 2 1 x 2 + c 8 x 1 x 2 2 + c 9 x 3 1 + c 10 x 3 2Table A.2: Coefficients used for functions f 6 to f 15 This section introduces the test problems used for the numerical experiments of Section 5. All the test functions are presented in Table A.1 and are defined in [0, 1] 2 . Functions f 1 and f 2 mimick a real-life problem presented by Dutrieux et al. (2015b). Functions f 3 , f 4 (Branin) and f 5 (Rosenbrock) are well-known test functions, rescaled to [0, 1] 2 . Functions f 6 to f 15 correspond to randomly generated third degree polynomials, whose coefficients are presented in Table A.2. In our experiments, the values of these 15 test functions are scaled to [0, 1], assuming that the minimum and maximum values are known (cf. Remark 2).A total of nine problems denoted by g 1 to g 9 are proposed based on the previously defined functions (see TableA.3). For each problem, the input set X is a 21 × 21 regular grid on[0, 1]

	f		Definition
	f 1 f 2		780000 + 110000x 1 -12000x 2 -36000x 1 x 2 + 280000x 2 1 + 50000x 2 2 0.83 + 0.17x 1 -0.015x 2 -0.0038x 1 x 2 + 0.061x 2 1 + 0.0011x 2 2
	f 3		e (0.36(x1+x2)) + 0.6x 1 + 1.2x 2 2 + 3 sin(0.8πx 1 )
	f 4		(x 2 Function coefficients
		f 6	f 7	f 8	f 9	f 10	f 11	f 12	f 13	f 14	f 15
	c 1	0.36 0.68 0.094 0.61 -0.38 -0.19 0.78 -0.45 -0.45 0.75
	c 2	8.1	-9.4	-7.2	5	8.5	4.8	6	7.8	-9.3	7.4
	c 3	7.5	9.1	7	2.3	1.4	2.1	-4.7	-7.7	-3.5	-8.2
	c 4	-83	-2.9	49	-5.3	63	42	90	28	14	-98
	c 5	26	-60	68	30	81	56	-85	34	-9.7	15
	c 6	-80	72	-49	-66	96	77	-82	-31	22	-31
	c 7	-440 160	630	-170 -120	410	600	-500 -880 -450
	c 8	94	-830 -510	-99	-780	360	890	-170 -370	-62
	c 9	920 -580	860	-830 -480	150	370	-480	550	780
	c 10 930 -920 -300	430	-180	-16	-740	530	390	-260
	Appendix A. Benchmark problems

Table A .

 A 3: Definition of nine test problems using the test functions presented in Tables A.1-A.2. For problems g 5 to g 9 the input argument is shifted by x 0 : if function f k is indicated in the table, then the actual objective function is x → f k (x -x 0 ). The last two columns provide respectively the cardinality of the Pareto set and the proportion of Pareto-optimal solutions in the 21x21 input set. 10 f 11 2.1 × 10 3 3.2 × 10 2 (0.5, 0.5) 67 15% g 8 f 12 f 13 1.4 × 10 4 1.6 × 10 3 (0.3, 0.8)(0.6, 0.6) 63 14% g 9 f 14 f 15 3.7 × 10 3 2.0 × 10 4

		objectives	noise variance
	g	f g 1	f g 2	σ 2 1	σ 2 2	x 0	| P | | P | / | X |
	g 1 f 1	f 2	3.6 × 10 9 3.9 × 10 -3		136	31%
	g 2 f 4	f 3	3.1 × 10 2 4.8 × 10 3		10	2%
	g 3 f 4	f 5	3.1 × 10 2 5.7 × 10 8		12	3%
	g 4 f 3	f 5	4.8 × 10 3 5.7 × 10 8		7	2%
	g 5 f 6	f 7	7.0 × 10 2 5.6 × 10 3	(0.5, 0.5)	60	14%
	g 6 f 8	f 9	5.8 × 10 2 3.1 × 10 3	(0.5, 0.5)	22	5%
	g 7 f (0.3, 0.8)	36	8%

Algorithm 2: The PALS algorithm Input : input space X; GP priors µ 0,j , σ 0,j , k j , ∀j; j ; β i ; n max ; k. Output : predicted Pareto set P. n = 0, all classified = false repeat Obtain µ n (x) and σ n (x) for all x ∈ X R n (x) = Q µn,σn,βn (x) for all x ∈ X P n = ∅,

until all classified or n ≥ n max ; return P purpose of algorithm monitoring and, in our benchmarks, for error assessment).

Numerical experiments

Methodology

In this section, PALS is compared with a set of alternative optimization methods: 1) a pure random search approach (PRS), where evaluation points are drawn from the uniform probability distribution on X; 2) an alternative random search approach, where evaluation points are drawn from a model-based probability distribution concentrated on "promising" regions; and 3) two scalarization approaches, which are modifications of the ParEGO algorithm [START_REF] Knowles | ParEGO: A hybrid algorithm with on-line landscape approximation for expensive multiobjective optimization problems[END_REF].

More precisely, the alternative random search approach consists in drawing the next evaluation point from a distribution proportional to the probability of that point being misclassified, according to the posterior distributions of the Gaussian processes ξ j .

For the scalarization techniques, we consider two modified versions of the ParEGO algorithm suitable for a stochastic setting. ParEGO relies on the use of We introduce the performance metrics proposed to assess prediction quality. The volume of the symmetric difference measures a global error in the estimation. We use it to assess Pareto front estimates. The misclassification rate is used as a global error metric to evaluate Pareto set estimates.

Appendix B.1. Volume of the symmetric difference

Let the "dominated region" defined between a Pareto front F and a reference point R be the set of all the points that simultaneously dominate R, and are dominated by a member of F: D(F) = ∪ y ∈F {y ∈ R q : y ≺ y ≺ R} . Recalling that our objective functions are all scaled to [0, 1], we take R = (1.1, 1.1) in our (bi-objective) experiments. The volume of the dominated region is denoted by V (F), and is illustrated in Figure B.1 (left).

We can then define the volume of the symmetric difference of the hypervolumes defined between the true Pareto front and the estimation, and a reference point R. We denote it V d (F , F): the misclassification rate M (P , P) is defined as the proportion of points that differ between the prediction and true sets: