Markos Viggiato

Ricardo Terra
email: terra@dcc.ufla.br

Henrique S C Rocha
email: henrique.rocha@inria.fr

Marco Tulio Valente

Eduardo Figueiredo
email: figueiredo@dcc.ufmg.br

Microservices

Microservices in Practice: A Survey Study

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Microservices have become largely popular in the last years together with the spread of DevOps practices and containers technologies, such as Kubernetes and Docker [Pahl 2015]. We can see a significant increase in the use of microservices architectural style since 2014 [START_REF] Klock | Workload-based clustering of coherent feature sets in microservice architectures[END_REF]], which can be verified in the service-oriented software industry where the usage of microservices has been far superior when compared to other software architecture models [START_REF] Alshuqayran | A Systematic Mapping Study in Microservice Architecture[END_REF].

Microservices are autonomous components that isolate fine-grained business capabilities. Furthermore, a microservice usually runs on its own process and communicates using standardized interfaces and lightweight protocols [START_REF] Fowler | [END_REF]Lewis 2017, Hassan et al. 2017]. In practice, microservices are widely used by large Web companies, such as Netflix, LinkedIn, and Amazon, which can be motivated by the benefits that microservices bring, e.g., the reduced time to put a new feature in operation [START_REF] Alshuqayran | A Systematic Mapping Study in Microservice Architecture[END_REF]].

There are many benefits of using microservices, such as technology diversity in a single system, better scalability, increase productivity, and ease of deployment [START_REF] Alshuqayran | A Systematic Mapping Study in Microservice Architecture[END_REF], Newman 2015]. Consequently, these benefits may improve software maintainability [START_REF] Alshuqayran | A Systematic Mapping Study in Microservice Architecture[END_REF]]. However, microservices also have their drawbacks. Usually, the services that compose the software are part of a distributed setting. Therefore, microservices could complicate some tasks such as finding a service within the network, managing the security, executing transactions, and optimizing the communication between services [START_REF] Alshuqayran | A Systematic Mapping Study in Microservice Architecture[END_REF], Yu et al. 2016].

Shedding light on microservices usage in practice is important for many reasons. We can perceive the advantages that motivate practitioners and the most important challenges faced when developing software under this architecture. This information can support decision-making about migrating systems to microservices or even start to develop an entire application under this architectural style. It can also aid software developers to understand and follow the best practices, making the microservices usage more effective. In addition, it may support the adoption of practices in software domains by practitioners and the developers' perception of the software quality [START_REF] Mori | Evaluating domain-specific metric thresholds: an empirical study[END_REF], Oliveira et al. 2018, Guimaraes et al. 2013].

Nevertheless, there are no studies that investigate how microservices are used in practice. To the best of our knowledge, existing works consist of systematic mapping studies, which summarize the progress of microservices technology so far [Alshuqayran et al. 2016, Pahl andJamshidi 2016]. By contrast, in this paper, we propose to look at microservices from a practical perspective, i.e., with the aim of understanding and reveal how practitioners are in fact using microservices. More specifically, we describe the results of a survey designed to reveal the usage of microservices in practice. First, we conducted a mapping study to identify and highlight potential advantages and challenges faced by professionals who work with microservices. Next, we surveyed developers about the findings of the mapping study, aiming at verifying whether the use of microservices in the industry is according to the best recommendations mentioned in the literature and whether the advantages and challenges found in the mapping study are in fact what practitioners face.

Microservices

Microservices are an architectural style in which the process of software development is done by using autonomous components that isolate fine-grained business functionalities and communicate one with other through standardized interfaces [START_REF] Hassan | Microservice ambients: An architectural meta-modelling approach for microservice granularity[END_REF]]. Due to an extensive use in web and cloud-based applications, we can observe a migration of some companies from the monolith architecture to microservices since the latter brings many benefits such as self-manageable (decentralized governance) and lightweight components [START_REF] Aderaldo | Benchmark Requirements for Microservices Architecture Research[END_REF]].

The purpose of microservices is to use autonomous units that are isolated one from another and coordinate them into a distributed infrastructure by a lightweight container technology, such as Docker. Usually, the adoption of this architectural model implies also in adopting agile practice, such as DevOps, which reduces the time between implementing a change in the system and transferring this change to the production environment [START_REF] Aderaldo | Benchmark Requirements for Microservices Architecture Research[END_REF].

The isolation of business functionalities is highly recommended when using microservices, and allows independent development and deployment of each microservice. Moreover, the isolation also optimizes the autonomy and the replaceability of the services. Indeed, a microservice architectural style brings many benefits for developers but is also comes with many challenges. In Section 3.1, we present a detailed description of the advantages and challenges of working with microservices.

Study Design

This study included two phases, a mapping study (Section 3.1), and a survey (Section 3.2), as described next.

Mapping Study

Initially, we performed a mapping study to collect information about microservices from blogs and articles, as well as from more traditional literature, including books and papers. Mapping studies are particularly recommended for understanding emerging fields or technologies [START_REF] Wohlin | Experimentation in Software Engineering[END_REF], which is certainly the case of microservices. In order to retrieve documents about microservices, we used three specific search strings on Google: microservices architecture, good features of microservices architecture, and bad features/parts of microservices architecture. Our intention was to gather documents regarding all aspects of microservices, such as documents proposing general terms and definitions (using the first string), and documents with the best characteristics and the main challenges faced by developers (using the second and third strings). After analyzing the retrieved documents, we identified five papers and one book from the scientific literature. Furthermore, we also considered 15 relevant documents from well-known experienced practitioners, including ten articles from websites and five documents from blogs. It is important to note that microservices have become popular in recent years, therefore there is still not a large number of relevant documents available.

The first author of this paper carefully read all the 21 relevant documents in order to extract their recurrent topics and themes. Thereupon, we classified the topics into advantages and challenges faced by developers when already using microservices architectures.

Advantages. In a system composed of multiple microservices, developers have the possibility of using many different technologies [Newman 2015]. This technology diversity is a very common characteristic in applications using microservices and it allows the use of the right tool for the right job. Moreover, the heterogeneity of technologies allows the addition of new technologies during development or maintenance in a more efficient way. For example, this is suitable for web applications where we can observe a constant and fast change in development environments and frameworks [START_REF] Ramos | Angu-larJS in the wild: A survey with 460 developers[END_REF]. Another benefit often associated to microservices is the possibility of deploying a given service independently from the others. This independent deployment could lead to a faster implementation of new features [Newman 2015]. Scalability is also mentioned as an advantage of microservices since it can be achieved on demand, scaling only the service that contains a given functionality. It is also possible to replicate specific services, instead of the entire system. Finally, maintainability is often reported as a benefit since developers can modify or replace a service without impacting the entire application.

Challenges. It is often reported that microservices demand distributed data management and hence distributed transactions, which makes their implementation much more complex. Other studies [Pahl andJamshidi 2016, Aderaldo et al. 2017] report that automated tests are extremely important in microservices, especially integration tests, which can be more complex and time-consuming. In addition, when the tests fail, it can be harder to determine which functionality has been broken [Newman 2015]. Service faults are also cited as a challenge when using microservices since the identification of a fault in a distributed setting is much harder than in a monolithic one. Finally, developers often mention that Remote Procedure Calls (RPC) are expensive and take much longer than local calls, which means that RPC may become a challenge when developing applications under microservices.

Survey Design

We designed a survey aiming at confirming (or not) the general characteristics, advantages, and challenges associated with microservices, as indicated by our mapping study. The survey has 14 questions and it is divided into three sections. The first section is about the background of the participants, including questions about experience with software development and with microservices as well as about the size of applications/number of services that they already worked with. The second section is related to definitions and trade-offs of using microservices. For instance, it includes questions about the ideal size of a microservice, and advantages and problems faced by developers when using this technology. Finally, the third section is composed of open-ended questions about the technologies used by developers when implementing microservices applications. Our intention is to identify the most popular programming languages and technologies used under microservices architecture since the literature states that many technologies can be used in microservices systems.

To find participants, we implemented an algorithm to search for microservices developers in the Stack Overflow community. We identified and retrieved nicknames from users who own questions or answers containing tag microservices. In order to collect the email, we matched the Stack Overflow nickname with the equivalent nickname at GitHub. In addition, we promoted the survey in many communities about microservices and cloudbased development, including a Google group1 , a Google Plus community2 , and a Reddit community3 . The survey remained open between June and July 2017, and we obtained 122 complete responses.

Survey Results

In this section, we describe the results obtained from the survey by presenting the participants background experience (Section 4.1), the popular languages and technologies (Section 4.2), the perceived advantages and challenges (Section 4.3), and the participants' feedback (Section 4.4).

Participants Background

As presented in Figure 1, almost 72% of the participants have at least five years of experience with software development. Furthermore, about 74% of them have more than one year of experience with microservices; more specifically, almost 67% have one to five years of experience. In addition, approximately 64% of the respondents are back-end developers while only 11.5% are DevOps.

Although not graphically illustrated, we also collected data regarding the size of the applications. 70.5% of the respondents worked with monolithic-based systems larger than 50 KLOC, and about 54% worked with microservices-based applications larger than 10 KLOC. These numbers confirm that most survey participants are not novice in the microservices field.

Considering that microservices are an emerging field, and based on the professional background of the respondents, we claim the participants have sufficient knowledge on the subject to answer the survey questions.

Most Popular Programming Languages and Technologies

Aiming at characterizing microservices applications, we asked the survey participants about the languages and technologies they usually use in their projects. We found that four programming languages are largely used: Java (33%), JavaScript through Node.js (18%), C# (12%), and PHP (8%). The answers also mention other 14 programming languages, which indicate the flexibility of microservices-based applications regarding programming languages. Regarding the most common DBMS, the results include Postgres (30%), MySQL (25%), MongoDB (20%), SQL Server (12%), and Oracle (9%). Finally, regarding the communication protocols, 62% of the participants declared they use REST over HTTP.

Advantages and Challenges of Microservices

Figure 2a presents the percentage of participants' answers about four advantages usually associated to microservices, in a scale from 1 (very important) to 4 (not important at all). For independent deployment, we can see the largest difference from score 1 to the others, which reveals a higher agreement rate for this feature when compared to the others. Yet, for the other three advantages, we can verify that more than 50% of the respondents chose scores 1 or 2, indicating these characteristics are in fact relevant when working with microservices.

From the mapping study, we also identified four main challenges faced by developers when working with microservices. In Figure 2b, we can see the responses of the survey's participants regarding these challenges. The participants agree that complex distributed transactions are a very important challenge. In fact, this challenge has the highest percentage for score 1 (32.8%). We can also observe that about 57% and 49% consider testing the whole system and service faults, respectively, as important challenges (scores 1 and 2). In contrast, challenge expensive remote calls is very important to only 13.8% of the participants. Interestingly, practitioners disagree with the literature in this last challenge since they do not find expensive remote calls a very important challenge in microservices development. In a nutshell, developers should pay special attention to distributed transactions, testing of the whole system and service faults, since these may become big problems in the system.

Feedback from Participants

We also sent another email to the developers directly contacted to answer our survey. In this follow-up message, we described the major survey results with the aim of receiving their impressions about our study and findings. We also intended to verify whether they agree or not with our results. We received nine answers; all of them with a positive feedback. In general, developers answered that our paper provides a good overview of the microservice practice. For instance, two developers commented this is really interesting and good paper. Another developer highlighted that microservices are not a "holy grail" and that monolithic can also have small and separated modules, even in distributed settings.

Threats to Validity

Some threats may affect the validity of our findings. First, the survey participants may not represent the entire population of microservices practitioners. To mitigate this risk, we put efforts in promoting the survey in many different communities to include professionals from different software ecosystems. Second, the term DevOps (which is one of the options of the survey question about the participant's role) might not be common in some contexts. For example, in small organizations, DevOps tasks such as development and delivery process automation may fall on senior developers and architects. Third, although it would be desirable for our survey the analysis of larger (w.r.t. size), stable (w.r.t. age), and specific branches of industry applications, we argue that our survey brings a broad overview since it is based on the expertise of 122 developers who work with heterogeneous microservices-based applications.

Related Work

Most of the research in microservices restrict their study to a specific domain, such as business [START_REF] Yu | A microservice based reference architecture model in the context of enterprise architecture[END_REF]. There are also researches investigating microservices by performing a systematic mapping study. For instance, a study summarized the progress of studies about microservices until 2016, and identified the gaps and requirements [START_REF] Alshuqayran | A Systematic Mapping Study in Microservice Architecture[END_REF]]. Other study taxonomically classified and compared studies of this architectural style and their application in the cloud [START_REF] Pahl | Microservices: A systematic mapping study[END_REF]. Our study follows a different route as we look at the microservices from a practical perspective. We aimed to understand and indicate how the software development industry is, in fact, using this popular architecture, and how practitioners perceive the advantages and challenges of microservices.

Conclusion

The findings of this paper indicate that practitioners usually follow the best practices for microservices reported in the literature. We also confirmed the benefits provided by microservices, such as independent deployment, ease to scale the applications, maintainability, and no commitment to a single technology stack. Last but not least, we also confirmed the challenges developers may face, such as complex distributed transactions, testing the whole system, and service faults.

However, we also found some important topics that are in disagreement. First, professionals usually work as back-end developers (64%), instead of as DevOps specialists in cross-functional teams. Second, according to the mapping study, expensive remote calls is one of the challenges that developers face when working with microservices. However, about 52% of the respondents declared that it is not an important or it is a little important issue in their systems.

As future work, we intend to conduct interviews with microservices professionals to confirm our results and to better understand if and why practitioners do not follow some best practices. We also plan to perform an industrial-scale case study with companies that adopt microservices to monitor real software developers developing microservices-based projects in order to report the problems they face and the solutions they apply.

  Figure 1. Participants' background.

Figure 2 .

 2 Figure 2. Microservices advantages and challenges

https://groups.google.com/forum/#!forum/microservices

https://plus.google.com/communities/112442985624053749478

https://www.reddit.com/r/microservices/

Acknowledgements

Our research has been supported by CAPES, FAPEMIG, and CNPq.


Victor Picheny
email: victor.picheny@inra.fr

Mickaël Binois
email: mickael.binois@chicagobooth.edu

Abderrahmane Habbal
email: habbal@unice.fr

A Bayesian optimization approach to find Nash equilibria

Keywords: Game theory, Gaussian processes, Stepwise Uncertainty Reduction

Game theory finds nowadays a broad range of applications in engineering and machine learning. However, in a derivative-free, expensive black-box context, very few algorithmic solutions are available to find game equilibria. Here, we propose a novel Gaussian-process based approach for solving games in this context. We follow a classical Bayesian optimization framework, with sequential sampling decisions based on acquisition functions. Two strategies are proposed, based either on the probability of achieving equilibrium or on the Stepwise Uncertainty Reduction paradigm. Practical and numerical aspects are discussed in order to enhance the scalability and reduce computation time. Our approach is evaluated on several synthetic game problems with varying number of players and decision space dimensions. We show that equilibria can be found reliably for a fraction of the cost (in terms of black-box evaluations) compared to classical, derivative-based algorithms. The method is available in the R package GPGame available on CRAN at https://cran.r-project.org/package=GPGame.

Introduction

Game theory arose from the need to model economic behavior, where multiple decision makers (MDM) with antagonistic goals is a natural feature. It was further extended to broader areas, where MDM had however to deal with systems governed by ordinary differential equations, the so-called differential games. See e.g., Gibbons [START_REF] Gibbons | Game Theory for Applied Economists[END_REF] for a nice introduction to the general theory and Isaacs [START_REF] Isaacs | Differential games. A mathematical theory with applications to warfare and pursuit, control and optimization[END_REF] for differential games. Recently, engineering problems with antagonistic design goals and with real or virtual MDM were formulated by some authors within a game-theoretic framework. See e.g., León et al [START_REF] León | Concurrent aerodynamic optimization of rotor blades using a nash game method[END_REF] for aerodynamics, Habbal et al [START_REF] Habbal | Multidisciplinary topology optimization solved as a Nash game[END_REF] for structural topology design, Habbal and Kallel [START_REF] Habbal | Neumann-Dirichlet Nash strategies for the solution of elliptic Cauchy problems[END_REF] for missing data recovery problems. The study of multi-agent systems or games such as poker under this setting is also quite common in the AI and machine learning communities, see e.g., Johanson and Bowling [START_REF] Johanson | Data biased robust counter strategies[END_REF], Lanctot et al [START_REF] Lanctot | No-regret learning in extensiveform games with imperfect recall[END_REF], Brown et al [START_REF] Brown | Hierarchical abstraction, distributed equilibrium computation, and post-processing, with application to a champion no-limit Texas hold'em agent[END_REF].

Solutions to games are called equilibria. Contrarily to classical optimization, the definition of an equilibrium depends on the game setting (or rules). Within the static with complete information setting, a relevant one is the so-called Nash equilibrium (NE). Shortly speaking, a NE is a fixed-point of iterated many single optimizations (see the exact definition in Section-2 below). Its computation generically carries the well known tricks and pitfalls related to computing a fixed-point, as well as those related to intensive optimizations notably when cost evaluations are expensive, which is the case for most engineering applications. There is an extensive literature related to theoretical analysis of algorithms for computing NE [START_REF] Bas ¸ar | Relaxation techniques and asynchronous algorithms for on-line computation of noncooperative equilibria[END_REF]40,[START_REF] Uryas'ev | On relaxation algorithms in computation of noncooperative equilibria[END_REF]], but very little -if any-on black-box models (i.e., non convex utilities) and expensive-to-evaluate ones; to the best of our knowledge, only home-tailored implementations are used. On the other hand, Bayesian optimization [BO,[START_REF] Mockus | Bayesian Approach to Global Optimization: Theory and Applications[END_REF] is a popular approach to tackle black-box problems. Our aim is to investigate the extension of such approach to the problem of computing game equilibria.

BO relies on Gaussian processes, which are used as emulators (or surrogates) of the black-box model outputs based on a small set of model evaluations. Posterior distributions provided by the Gaussian process are used to design acquisition functions that guide sequential search strategies that balance between exploration and exploitation. Such approaches have been applied for instance to multi-objective problems [START_REF] Wagner | On expected-improvement criteria for modelbased multi-objective optimization[END_REF], as well as transposed to frameworks other than optimization, such as uncertainty quantification [START_REF] Bect | Sequential design of computer experiments for the estimation of a probability of failure[END_REF] or optimal stopping problems in finance [START_REF] Gramacy | Sequential design for optimal stopping problems[END_REF].

In this paper, we show that the BO apparatus can be applied to the search of game equilibria, and in particular the classical Nash equilibrium (NE). To this end, we propose two complementary acquisition functions, one based on a greedy search approach and one based on the Stepwise Uncertainty Reduction paradigm [START_REF] Fleuret | Graded learning for object detection[END_REF]. The corresponding algorithms require very few model evaluations to converge to the solution. Our proposal hence broadens the scope of applicability of equilibrium-based methods, as it is designed to tackle derivative-free, non-convex and expensive models, for which a game perspective was previously out of reach.

The rest of the paper is organized as follows. Section 2 reviews the basics of game theory and presents our Gaussian process framework. Section 3 presents our main contribution, with the definition of two acquisition functions, along with com-putational aspects. Finally, Section 4 demonstrates the capability of our algorithm on three challenging problems.

Background

Games and equilibria

Nash games

We consider primarily the standard (static, under complete information) Nash equilibrium problem [NEP,[START_REF] Gibbons | Game Theory for Applied Economists[END_REF].

Definition 1 A NEP consists of p ≥ 2 decision makers (i.e., players), where each player i ∈ {1, . . . , p} tries to solve his optimization problem:

(P i ) min xi∈Xi y i (x), (1)
where y(x) = [y 1 (x), . . . , y p (x)] : X ⊂ R n → R p (with n ≥ p) denotes a vector of cost functions (a.k.a. pay-off or utility functions), y i denotes the specific cost function of player i, and the vector x consists of block components x 1 , . . . , x p (x = (x j ) 1≤j≤p ).

Each block x i denotes the variables of player i and X i its corresponding action space and X = i X i . We shall use the convention y i (x) = y i (x i , x -i ) when we need to emphasize the role of x i .

Definition 2 A Nash equilibrium x * ∈ X is a strategy such that:

(N E) ∀i, 1 ≤ i ≤ p, x * i ∈ arg min xi∈Xi y i (x i , x * -i ). (2)
In other words, when all players have chosen to play a NE, then no single player has incentive to move from his x * i . Let us however mention by now that, generically, Nash equilibria are not efficient, i.e., do not belong to the underlying set of best compromise solutions, called Pareto front, of the objective vector (y i (x)) x∈X .

Random games

We shall also deal with the case where cost functions are uncertain. Such problems belong to a family of random games that are called disturbed games by Harsanyi in [START_REF] Harsanyi | Games with randomly disturbed payoffs: A new rationale for mixed-strategy equilibrium points[END_REF]. We denote such cost functions f i (x, (ξ)), where = ( i ) : Ξ → R p is a random vector defined over a probability space (Ξ, F, P). In the following we refer to our setting as random games, but emphasize that we consider static Nash games with expectations of randomly perturbed costs. Definition 3 Assuming risk-neutrality of the players, a random Nash game consists of p ≥ 2 players, where each player i ∈ {1, . . . , p} tries to solve

(SP i ) min xi∈Xi E[f i (x, (ξ))]. (3)
Definition 4 A random Nash equilibrium x * ∈ X is a strategy such that:

(SN E) ∀i, 1 ≤ i ≤ p, x * i ∈ arg min xi∈Xi E[f i (x i , x * -i , (ξ))]. (4)
Note that since E[f i (x, (ξ))] is not random, one can equivalently set

y i = E[f i (x, (ξ))]
to formulate the solution as a (N E). The difference lies in the fact that here the cost function y i cannot be evaluated exactly.

Working hypotheses

In this work, we focus on continuous-strategy non-cooperative Nash games (i.e., with infinite sets X i ) or on large finite games (i.e., with large finite sets X i ).

Our working hypotheses are:

queries on the cost function (i.e., pointwise evaluation of the y i 's for a given x) result from an expensive process: typically, the y i 's can be the outputs of numerical models; the cost functions may have some regularity properties but are possibly strongly not convex (e.g., continuous and multimodal); the cost functions evaluations can be corrupted by noise; -X is either originally discrete, or a representative discretization of it is available (so that the equilibrium of the corresponding finite game is similar to the one of original problem).

Note that to account for the particular form of NEPs, X must realize a full-factorial design: Let us remark that we do not require for an equilibrium to exist and be unique, the case when there is no or equilibria is discussed in Section 3.2.4.

X = X 1 × . . . × X p . Given each action space X i = {x 1 i , . . . , x mi i } of size m i , X consists of all the combinations (x k i , x l j ) (1 ≤ i = j ≤ p, 1 ≤ k ≤ m i , 1 ≤ l ≤ m j ),
In the case of noisy evaluations, we consider only here an additive noise corruption:

f i (x, (ξ)) = y i (x) + i (ξ), (5)
and we assume further that has independent Gaussian centered elements: i ∼ N (0, τ 2 i ). Notice that in this case, both problems and equilibria coincide, and can be solved with the same algorithm. Hence, in the following, all calculations are given in the noisy case, while the deterministic case of the standard NEP is recovered by setting i = 0 and τ i = 0.

Furthermore, we consider solely pure-strategy Nash equilibria [as opposed to mixed-strategies equilibria, in which the optimal strategies can be chosen stochastically, see 19, Chapter 1], and as such, we avoid solving the linear programs LP or linear complementarity problems LCP generally used in the dedicated classes of algorithms à la Lemke-Howson [START_REF] Rosenmüller | On a generalization of the lemke-howson algorithm to noncooperative n-person games[END_REF].

Related work

Let us mention that in the continuous (in the sense of smooth) games setting, there is an extensive literature dedicated to the computation of NE, based on the rich theory of variational analysis, starting with the classical fixed-point algorithms to solve NEPs [START_REF] Uryas'ev | On relaxation algorithms in computation of noncooperative equilibria[END_REF][START_REF] Bas ¸ar | Relaxation techniques and asynchronous algorithms for on-line computation of noncooperative equilibria[END_REF]40]. When the players share common constraints, Nash equilibria are shown to be Fritz-John (FJ) points [11], which allows the op. cit. authors to propose a nonsmooth projection method (NPM) well adapted for the computation of FJ points. From other part, it is well known (and straightforward) that Nash equilibria are in general not classical Karush-Kuhn-Tucker (KKT) points, nevertheless, a notion of KKT condition for generalized Nash equilibria GNEP is developed in Kanzow and Steck [START_REF] Kanzow | Augmented lagrangian methods for the solution of generalized nash equilibrium problems[END_REF], which allows the authors to derive an augmented Lagrangian method to compute GNEPs.

Noncooperative stochastic games theory, starting from the seminal paper by Shapley [START_REF] Shapley | Stochastic games[END_REF], occupies nowadays most of the game theorists, and a vast literature is dedicated to stochastic differential games [START_REF] Friedman | Stochastic differential games[END_REF], robust games [45], games on random graphs, or agents learning games [START_REF] Hu | Nash q-learning for general-sum stochastic games[END_REF], among many other branches, and it is definitely out of the scope of the paper to review all aspects of the field. See also the introductory book Neyman and Sorin [44] to the basic -yet deep-concepts of the stochastic games theory.

We also do not consider games with additional specific structures, like cooperative, zero-sum stochastic or deterministic games, or repeated Nash games [START_REF] Littman | A polynomial-time nash equilibrium algorithm for repeated games[END_REF]. For these games, tailored algorithms should be used, among which are the pure exploration statistical learning with Monte Carlo Tree Search [START_REF] Garivier | Maximin action identification: A new bandit framework for games[END_REF] or multi-agent reinforcement learning MARL, see e.g., Games [START_REF] Games | Lenient learning in independent-learner stochastic cooperative games[END_REF] and references therein.

We stress here that none of the above-mentioned approaches are designed to tackle expensive black-box problems. They may even prove unusable in this context, either because they could simply not converge or require too many cost function evaluations to do so.

Bayesian optimization

Gaussian process regression

The idea of replacing an expensive function by a cheap-to-evaluate surrogate is not recent, with initial attempts based on linear regression. Gaussian process (GP) regression, or kriging, extends the versatility and efficiency of surrogate-based methods in many applications, such as in optimization or reliability analysis. Among alternative non-parametric models such as radial basis functions or random forests, see e.g., Wang and Shan [START_REF] Wang | Review of metamodeling techniques in support of engineering design optimization[END_REF], Shahriari et al [START_REF] Shahriari | Taking the human out of the loop: A review of bayesian optimization[END_REF] for a discussion, GPs are attractive in particular for their tractability, since they are simply characterized by their mean m(.) and covariance (or kernel) k(., .) functions, see e.g., Cressie [START_REF] Cressie | Statistics for spatial data: Wiley series in probability and statistics 11[END_REF], Rasmussen and Williams [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF]. In the following, we consider zero-mean processes (m = 0) for the sake of conciseness.

Briefly, for a single objective y, conditionally on n noisy observations f = (f 1 , . . . , f n ), with independent, centered, Gaussian noise, that is, f i = y(x i ) + ε i with ε i ∼ N (0, τ 2 i ), the predictive distribution of y is another GP, with mean and covariance functions given by:

µ(x) = k(x) K -1 f , (6) σ 2 (x, x ) = k(x, x ) -k(x) K -1 k(x ), (7)
where k(x) := (k(x, x 1 ), . . . , k(x, x n )) and K := (k(x i , x j ) + τ 2 i δ i=j ) 1≤i,j≤n , δ standing for the Kronecker function. Commonly, k(., .) belongs to a parametric family of covariance functions such as the Gaussian and Matérn kernels, based on hypotheses about the smoothness of y. Corresponding hyperparameters are often obtained as maximum likelihood estimates, see e.g., Rasmussen and Williams [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF] or Roustant et al [START_REF] Roustant | DiceKriging, DiceOptim: Two R packages for the analysis of computer experiments by kriging-based metamodeling and optimization[END_REF] for the corresponding details.

With several objectives, a statistical emulator for y(x) = [y 1 (x), . . . , y p (x)] is needed. While a joint modeling is possible [see e.g., 2], it is more common practice to treat the y i 's separately. Hence, conditioned on a set of vectorial observations {f 1 , . . . , f n }, our emulator is a multivariate Gaussian process Y:

Y(.) ∼ GP (µ(.), Σ (., .)) , (8)
with µ(.) = [µ 1 (.), . . . , µ p (.)], Σ = diag σ 2 1 (., .), . . . , σ 2 p (., .) , such that {µ i (.), σ 2 i (., .)} is the predictive mean and covariance, respectively, of a GP model of the objective y i . Note that the predictive distribution of an observation is:

F(x) ∼ N µ(x), Σ (x, x) + diag(τ 2 1 , . . . , τ 2 p ) . (9)
GPs are commonly limited to a few thousands of design points, due to the cubic cost needed to invert the covariance matrix. This can be overcome in several ways, by using inducing points [START_REF] Wilson | Kernel interpolation for scalable structured gaussian processes (kiss-gp)[END_REF], or local models [START_REF] Gramacy | Local gaussian process approximation for large computer experiments[END_REF][START_REF] Rullière | Nested kriging predictions for datasets with a large number of observations[END_REF]; see also [START_REF] Heaton | Methods for analyzing large spatial data: A review and comparison[END_REF] for a comparison or [START_REF] Žilinskas | Stochastic global optimization: a review on the occasion of 25 years of informatica[END_REF] for a broader discussion. In addition, here the full-factorial structure of the design space can potentially be exploited. For instance, if evaluated points also have a full-factorial structure, then the covariance matrix can be written as a Kronecker product, reducing drastically the computational cost, see e.g., [START_REF] Plumlee | Fast prediction of deterministic functions using sparse grid experimental designs[END_REF].

Sequential design

Bayesian optimization methods are usually outlined as follows: a first set of observations {X n0 , f n0 } is generated using a space-filling design to obtain a first predictive distribution of Y(.). Then, observations are performed sequentially by maximizing a so-called acquisition function (or infill criterion) J(x), that represents the potential usefulness of a given input x. That is, at step n ≥ n 0 ,

x n+1 = arg max x∈X J(x). (10)
Typically, an acquisition function offers an efficient trade-off between exploration of unsampled regions (high posterior variance) and exploitation of promising ones (low posterior mean), and has an analytical expression which makes it inexpensive to evaluate, conveniently allowing to use of-the-shelf optimization algorithms to solve Eq. [START_REF] Cressie | Statistics for spatial data: Wiley series in probability and statistics 11[END_REF]. In unconstrained, noise-free optimization, the canonical choice for J(x) is the so-called Expected Improvement [EI, 36], while in the bandit literature (noisy observations), the Upper Confidence Bound [UCB, 58] can be considered as standard. Extensions abound in the literature to tackle various optimization problems: see e.g. Wagner et al [START_REF] Wagner | On expected-improvement criteria for modelbased multi-objective optimization[END_REF] for multi-objective optimization or Hernández-Lobato et al [START_REF] Hernández-Lobato | A general framework for constrained bayesian optimization using information-based search[END_REF] for constrained problems.

Figure 1 provides an illustration of the BO principles.

q q q q q q q 0.0 0.1 0.2 0.3 0.4 0.5 0.6 J(x) q q y(x) GP model J(x) Observations New observation J maximizer
Fig. 1 One iteration of Bayesian optimization on a one-dimensional toy problem. A first GP is conditioned on a set of five observations (left), out of which an acquisition function is maximized to find the next observation. Once this observation is performed (right), the GP model and acquisition function are updated, and start pointing towards the optimum x = 0.75. Note that the acquisition is also large in unexplored regions (around x = 0.2).

Acquisition functions for NEP

We propose in the following two acquisition functions tailored to solve NEPs, respectively based on the probability of achieving equilibrium and on stepwise uncertainty reduction. Both aim at providing an efficient trade-off between exploration and exploitation.

Probability of equilibrium

Given a predictive distribution of Y(.), a first natural metric to consider is the probability of achieving the NE. From (2), using the notation x = (x i , x -i ), this probability writes:

P E (x) = P p i=1 Y i (x i , x -i ) = min x k i ∈Xi Y i (x k i , x -i ) , (11)
where {x 1 i , . . . , x mi i } denotes the m i alternatives in X i , and x -i is fixed to its value in x.

Since our GP model assumes the independence of the posterior distributions of the objectives, we have:

P E (x) = p i=1 P Y i (x) = min x k i ∈Xi Y i (x k i , x -i ) := p i=1 P i (x). (12)
Let us now introduce the notation

x i = x l i (1 ≤ l ≤ m i ).
As exploited recently by Chevalier and Ginsbourger [START_REF] Chevalier | Fast computation of the multi-points expected improvement with applications in batch selection[END_REF] in a multi-point optimization context, each P i can be expressed as

P i (x) = P   k∈mi,k =l Y i (x l i , x -i ) -Y i (x k i , x -i ) ≤ 0   . (13)
P i (x) amounts to compute the cumulative distribution function (CDF) of a Gaussian vector of size q := m i -1:

P i (x) = P (Z i ≤ 0) = Φ µ Z i ,Σ Z i (0), (14)
with

Z i = Y i (x 1 i , x -i ) -Y i (x l i , x -i ), . . . , Y i (x l-1 i , x -i ) -Y i (x l i , x -i ), Y i (x l+1 i , x -i ) -Y i (x l i , x -i ), . . . , Y i (x mi i , x -i ) -Y i (x l i , x -i ) .
The mean µ Zi and covariance Σ Zi of Z i can be expressed as:

(µ Zi ) j = µ i (x l i , x -i ) -µ i (x j i , x -i ), (Σ Zi ) jk = c ll i + c jk i -c lj i -c lk i if k, l = j and c ll i otherwise, with c jk i = σ 2 i (x j i , x -i ), (x k i , x -i ) .
Several fast implementations of the multivariate Gaussian CDF are available, for instance in the R packages mnormt (for q < 20) [START_REF] Azzalini | The R package mnormt: The multivariate normal and t distributions (version 1.5-4[END_REF] or, up to q = 1000, by Quasi-Monte-Carlo with mvtnorm [START_REF] Genz | Computation of Multivariate Normal and t Probabilities[END_REF][START_REF] Genz | mvtnorm: Multivariate Normal and t Distributions[END_REF].

Alternatively, this quantity can be computed using Monte-Carlo methods by drawing R samples Y

(1) i , . . . , Y (R) i of Y i (x 1 i , x -i ), . . . , Y i (x mi i , x -i ) , to compute Pi (x) = 1 R R r=1 1 Y (r) i (x) = min x k i ∈Xi Y (r) i (x k i , x -i ) ,
1(.) denoting the indicator function. This latter approach may be preferred when the number of alternatives m i is high (say > 20), which makes the CDF evaluation overly expensive while a coarse estimation may be sufficient. Note that in both cases, a substantial computational speed-up can be achieved by removing from the X i 's the non-critical strategies. This point is discussed in Section 3.2.3. Using J(x) = P E (x) as an acquisition function defines our first sequential sampling strategy. The strategy is rather intuitive: i.e., sampling at designs most likely to achieve NE.

Still, maximizing P E is a myopic approach (i.e., favoring an immediate reward instead of a long-term one), which are often sub-optimal [see e.g. 20, 21, and references therein]. Instead, other authors have advocated the use of an information gain from a new observation instead, see e.g., Villemonteix et al [START_REF] Villemonteix | An informational approach to the global optimization of expensive-to-evaluate functions[END_REF], Hennig and Schuler [START_REF] Hennig | Entropy search for information-efficient global optimization[END_REF], which motivated the definition of an alternative acquisition function, which we describe next.

Stepwise uncertainty reduction

Stepwise Uncertainty Reduction (SUR, also referred to as information-based approach) has recently emerged as an efficient approach to perform sequential sampling, with successful applications in optimization [START_REF] Villemonteix | An informational approach to the global optimization of expensive-to-evaluate functions[END_REF]47,[START_REF] Hernández-Lobato | Predictive entropy search for efficient global optimization of black-box functions[END_REF][START_REF] Hernández-Lobato | A general framework for constrained bayesian optimization using information-based search[END_REF] or uncertainty quantification [START_REF] Bect | Sequential design of computer experiments for the estimation of a probability of failure[END_REF][START_REF] Jala | Sequential design of computer experiments for the assessment of fetus exposure to electromagnetic fields[END_REF]. Its principle is to perform a sequence of observations in order to reduce as quickly as possible an uncertainty measure related to the quantity of interest (in the present case: the equilibrium).

Acquisition function definition

Let us first denote by Ψ (y) the application that associates a NE with a multivariate function. In the case of finite games, we have: Ψ : R N ×p → R p , for which a pseudocode is detailed in Algorithm 3. If we consider the random process Y (Eq. 8) in lieu of the deterministic objective y, the equilibrium Ψ (Y) is a random vector of R p with unknown distribution. Let Γ be a measure of variability (or residual uncertainty) of Ψ (Y); we use here the determinant of its second moment:

Γ (Y) = det [cov (Ψ (Y))] . (15)
The SUR strategy aims at reducing Γ by adding sequentially observations y(x) on which Y is conditioned. An "ideal" choice of x would be:

x n+1 = arg min x∈X Γ [Y|f = y(x)] , (16)
where Y|f = y(x) is the process conditioned on the observation f . Since we do not want to evaluate y for all x candidates, we consider the following criterion instead:

J(x) = E F (Γ [Y|F = Y(x) + ε]) , (17)
with F following the posterior distribution (conditioned on the n current observations, Eq. 9) and E F denoting the expectation over F.

In practice, computing J(x) is a complex task, as no closed-form expression is available. The next subsection is dedicated to this question.

Remark For simplicity of exposition, we assume here that the equilibrium Ψ (Y) exists, which is not guaranteed even if Ψ (y) does. To avoid this problem, one may consider an extended Ψ function equal to +∞ × I p if there is no equilibrium and to Ψ otherwise, and in Eq. 15 use the restriction of Ψ to finite values.

Approximation using conditional simulations

Let us first focus on the measure Γ when no new observation is involved. Due to the strong non-linearity of Ψ , no analytical simplification is available, so we rely on conditional simulations of Y to evaluate Γ .

Let Y 1 , . . . , Y M be independent draws of Y(X) (each Y i ∈ R N ×p ). For each draw, the corresponding NE Ψ (Y i ) can be computed by exhaustive search. We reported to Appendix C the particular algorithm we used for this step. The following empirical estimator of Γ (Y) is then available:

Γ (Y 1 , . . . , Y M ) = det [Q Y ] , with Q Y the sample covariance of Ψ (Y 1 ), . . . , Ψ (Y M ).
Now, let us assume that we evaluate the criterion for a given candidate observation point x. Let F 1 , . . . , F K be independent draws of F(x) = Y(x)+ε. For each F i , we can condition Y on the event

(F(x) = F i ) in order to generate Y 1 |F i , . . . , Y M |F i draws of Y|F i , from which we can compute the empirical estimator Γ Y 1 |F i , . . . , Y M |F i .
Then, an estimator of J(x) is obtained using the empirical mean:

Ĵ(x) = 1 K K i=1 Γ Y 1 |F i , . . . , Y M |F i .

Numerical aspects

The proposed SUR strategy has a substantial numerical cost, as the criterion requires a double loop for its computation: one over the K values of F i and another over the M sample paths. The two computational bottlenecks are the sample path generations and the searches of equilibria, and both are performed in total K × M times for a single estimation of J. Thankfully, several computational shortcuts allow us to evaluate the criterion efficiently. First, we employed the FOXY algorithm (fast update of conditional simulation ensemble) as proposed by Chevalier et al [START_REF] Chevalier | Fast update of conditional simulation ensembles[END_REF], in order to obtain draws of Y|F i based on a set of draws Y 1 , . . . , Y M . In short, a unique set of draws is generated prior to the search of x n+1 , which is updated quickly when necessary depending on the pair (x, F i ). The expression used are given in Appendix A, and we refer to [START_REF] Chevalier | Fast update of conditional simulation ensembles[END_REF] for the detailed algebra and complexity analysis.

Second, we discard points in X that are unlikely to provide information regarding the equilibrium prior to the search of x n+1 , as we detail below. By doing so, we reduce substantially the sample paths size and the dimension of each finite game, which drastically reduces the cost. We call X sim the retained subset.

Finally, Ĵ is evaluated only on a small, promising subset of X sim . We call X cand this set.

To select the subsets, we rely on a fast-to-evaluate score function C, which can be seen as a proxy to the more expensive acquisition function. The subset of X is then chosen by sampling randomly with probabilities proportional to the scores C(X), while ensuring that the subset retains a factorial form. We propose three scores, of increasing complexity and cost, which can be interleaved:

-C target : the simplest score is the posterior density at a target T E in the objective space, for instance the NE of the posterior mean (hence, it requires one NE search). C target reflects a proximity to an estimate of the NE. We use this scheme for the first iteration to select X sim ⊂ X. -C box : once conditional simulations have been performed, the above scheme can be replaced by the probability for a given strategy to fall into the box defined by the extremal values of the simulated NE (i.e., Ψ (Y 1 ), . . . , Ψ (Y M )). We use this scheme to select X sim ⊂ X for all the other iterations. -C P : since P E is faster (in particular in its Monte Carlo setting with small R) than Ĵ(x), it can be used to select X cand ⊂ X sim .

The detailed expressions of C target and C box are given in Appendix B. Note that in our experiments, C target and C box are also used with the P E acquisition function.

Last but not least, this framework enjoys savings from parallelization in several ways. In particular, the searches of NE for each sample Y can be readily distributed.

An overview of the full SUR approach is given in pseudo-code in Algorithm 1. Note that by construction, SUR does not necessarily sample at the NE, even when it is well-identified. Hence, as a post-processing step, the returned NE estimator is the design that maximizes the probability of achieving equilibrium.

Algorithm 1 Pseudo-code for the SUR approach

Require: n 0 , nmax, N sim , N cand 1: Construct initial design of experiments Xn 0 2: Evaluate yn 0 = F(Xn 0 ) 3: while n ≤ nmax do 4:

Train the p GP models on the current design of experiments {Xn, yn} 5:

if n = n 0 then 6:

estimate T E = Ψ (µ(X), the NE on the posterior mean; select X sim ⊂ X using Ctarget 7: else 8:

Select X sim ⊂ X using C box 9:

end if 10:

Generate M draws (Y 1 , . . . , Y M ) on X sim 11:

Compute Ψ (Y 1 ), . . . , Ψ (Y M ) (for C box ) 12:

Select X cand ⊂ X sim using C P 13:

Find x n+1 = arg min x∈X cand Ĵ(x) 14:

Evaluate y n+1 = F(x n+1 ) and add {x n+1 , y n+1 } to the current design of experiments 15: end while Ensure: x * = arg min x∈X cand P E (x)

Stopping criterion

We consider in Algorithm 1 a fixed budget, assuming that simulation cost is limited. However, other natural stopping criteria are available. For P E , one may stop if there is a strategy for which the probability is high, i.e., max P E (x) ≥ 1 -(with close to zero). For SUR, Ĵ is an indicator of the remaining uncertainty, so the algorithm may stop if this uncertainty is below a threshold, i.e., min J(x) ≤ .

While there always exists a Nash equilibrium for mixed strategy, in the setup we entertain there may actually be no pure strategy, or several. For P E , this is completely transparent, with a probability zero for all strategies (no NE) or several strategies having probability one (multiple NEs). For SUR, with more than one equilibrium, no change is needed, even though SUR may benefit from using a clustering method of the simulated NEs and defining local variability instead of the global Γ , Eq. [START_REF] Games | Lenient learning in independent-learner stochastic cooperative games[END_REF]. The absence of NE on the GP draws (Ψ (Y i )) can be used to detect the absence of NE for the problem at hand.

NUMERICAL EXPERIMENTS

These experiments have been performed in R [50] using the GPGame package [START_REF] Picheny | GPGame: Solving Complex Game problems using Gaussian processes[END_REF], which relies on the DiceKriging package [START_REF] Roustant | DiceKriging, DiceOptim: Two R packages for the analysis of computer experiments by kriging-based metamodeling and optimization[END_REF] for the Gaussian process regression part.

We used a fixed-point method as a competitive alternative to compute Nash equilibria, in order to assess the efficiency of our approach. It is a popular method among the audience who is familiar with gradient-descent optimization algorithms. The algorithm pseudo-code is given in Algorithm 2, and has been implemented in Scilab [START_REF]Scilab: Free and Open Source software for numerical computation[END_REF].

Algorithm 2 Pseudo-code for the fixed-point approach [START_REF] Uryas'ev | On relaxation algorithms in computation of noncooperative equilibria[END_REF] Require: P : number of players, 0 < α < 1 : relaxation factor, kmax : max iterations 1: Construct initial strategy x (0) 2: while k ≤ kmax do 3: Compute in parallel : ∀i, 1 ≤ i ≤ P, z

(k+1) i = arg min x i ∈X i J i (x (k)
-i , x i ) 4:

Update : x (k+1) = αz (k+1) + (1 -α)x (k) 5: 6:

if x (k+1) -x (k) small enough then exit 7:

end if 8: end while Ensure: For all i = 1...P , x * i = arg min

x i ∈X i J i (x * -i , x i )
In the following, performance is assessed in terms of numbers of calls to the objective function, hence assuming that the cost of running the black-box model largely exceeds the cost of choosing the points. For moderately expensive problems, the choice of algorithm may depend on the budget, as, intuitively, the time used to search for a new point may not exceed the time to simulate it. We report in Appendix D the computational times required for our approach on the three following test problems.

A classical multi-objective problem

We first consider a classical optimization toy problem (P1) as given in [START_REF] Parr | Improvement criteria for constraint handling and multiobjective optimization[END_REF], with two variables and two cost functions, defined as:

y 1 = (x 2 -5.1(x 1 /(2π)) 2 + 5 π x 1 -6) 2 + 10((1 - 1 8π
) cos(x 1 ) + 1) [START_REF] Games | Lenient learning in independent-learner stochastic cooperative games[END_REF]. Both functions are non-convex. We set x 1 = x 1 and x 2 = x 2 . The actual NE is attained at x 1 = -3.786 and x 2 = 15.

y 2 = -(10.5 -x 1 )(x 1 + 5.5)(x 2 + 0.5) - (x 2 -5.1(x 1 /(2π)) 2 -6) 2 30
- (1 -1/(8π)) cos(x 1 ) + 1 3 with x 1 ∈ [-5, 10] and x 2 ∈ [0,
Our strategies are parameterized as follow. X is first discretized over a 31 × 31 regular grid. With such size, there is no need to resort to subsets, so we use X cand = X sim = X. We set the number of draws to K = M = 20, which was empirically found as a good trade-off between accuracy and speed. We use n 0 = 6 initial observations from a Latin hypercube design [LHD,[START_REF] Mckay | Comparison of three methods for selecting values of input variables in the analysis of output from a computer code[END_REF], and observations are added sequentially using both acquisition functions. As a comparison, we ran a standard fixed-point algorithm [START_REF] Bas ¸ar | Relaxation techniques and asynchronous algorithms for on-line computation of noncooperative equilibria[END_REF] based on finite differences. This experiment is replicated five times with different initial sets of observations for the GP-based approaches and different initial points for the fixed-point algorithm. The results are reported in Table 1.

In addition, Figure 2 provides some illustration for a single run. In the initial state (top left), the simulated NEs form a cloud in the region of the actual NE. A first point is added, although far from the NE, that impacts the form of the cloud (top, middle). The infill criterion surface (Figure 2, top right) then targets the top left corner, which offers the best compromise between exploration and exploitation. After adding 7 points (bottom left) the cloud is smaller and centered on the actual NE. As the NE is quite well-identified, the infill criterion surface (Figure 2, bottom right) is now relatively flat and mostly indicates regions that have not been explored yet. After 14 additions (bottom, middle), all simulated NE but two concentrate on the actual NE. The observed values have been added around the actual NE, but not only, which indicates that some exploration steps (such as iteration 8) have been performed. Both GP approaches consistently find the NE very quickly: the worst run required 14 cost evaluations (that is, 8 infill points). In contrast, the classical fixed-point algorithm required hundreds of evaluations, which is only marginally better than an exhaustive search on the 961-points grid. Besides, two out of five runs converged to the stationary point (1, 1), which is not a NE. On this example, P E performed slightly better than SUR. Iteration 8 f1 f2 q q q q qq q q q q q q q q q Actual NE Simulated NE Current obs New obs Iteration 15 f1 f2 q q q q qq q q q q q q q q q q q q q q q q Actual NE Simulated NE Current obs New obs 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 q q q q q q q q q q q q q q Iteration 8 x1 x2 Fig. 2 Four iterations of SUR for (P1). Left and middle: observations and simulated NEs in the objective space. The hatched area represents the image of X by y. Right: level sets of the SUR criterion value in the variable space (green is better).

An open loop differential game

Differential games model a huge variety of competitive interactions, in social behavior, economics, biology among many others [predator-prey, pursuit-evasion games and so on, see 33]. As a toy-model, let us consider p players who compete to drive a dynamic process to their preferred final state. Let denote T > 0 the final time, and t ∈ [0, T ] the time variable. Then we consider the following simple process, which state z(t) ∈ R 2 obeys the first order dynamics:

ż(t) = v 0 + 1≤i≤p α i (t)x i (t) z(0) = z 0 (18)
The parameter α i (t) = e -θit models the lifetime range of each player's influence on the process (θ i is a measure of player (i)'s preference for the present). The timedependent action x i (t) ∈ R 2 is the player (i)'s strategy, which we restrict to the (finite-dimensional) space of spline functions of a given degree κ -1, that is:

x i (t) = 1≤k≤κ a (i) k × B k (t/T ) 1≤k≤κ b (i) k × B k (t/T )
, where (B k ) 1≤k≤κ is the spline basis. Now, the decision variables for player (i) is the array (a

(i) 1 , . . . , a (i) κ , b (i) 1 , . . . , b (i) κ )
(the spline coefficients), and x ∈ R κ×2×p . We used κ = 1 (constant splines, two decision variables per player) and κ = 2 (linear splines, four decision variables per player).

All players have their own preferred final states z i T ∈ R 2 , and a limited energy to dispense in playing with their own action. The cost y i of player (i) is then the following :

y i (x i , x -i ) = 1 2 z(T ) -z i T 2 R 2 + 1 2 x i 2 L 2 (0,T ) . (19)
The game considered here is an open loop differential game (the decisions do not depend on the state), and belongs to the larger class of dynamic Nash games.

-1.0 -0.5 0.0 0.5 1.0 -1.0 -0.5 0.0 0.5

1.0 z 1 z 2 q z T 2 q z T 1 q z T 3 q z T 4 z 0 initial state z(T) final state
Fig. 3 The differential game setting: The process starts at the initial state z 0 and travels during a time T to reach the state z(T ). Players targets are z

(i) T .
When the α i 's do not depend on i, and the points z i T are located on some -anycircle, then the Nash solution of the game puts the final state z(T ) on the center of the circle (provided T is large enough to let the state reach this center starting from z 0 ). In our setup, we have p = 4 players, with targets positioned on the four corners of the [-1, 1] 2 square, as illustrated by Figure 3. The state equation ( 18) is solved by means of an explicit Euler numerical scheme, with T set to 4, and 40 time steps. As for problem parameters, we chose v 0 = 0, z 0 = (0, 0.5) (non-central initial position) and θ = (θ i ) = (0.25, 0, 0.5, 0) (heterogeneous time preferences), which makes the search of NE non-trivial. The results labeled "Fixed point" presented in Table 2 are obtained by means of a popular fixed-point algorithm, as described in Algorithm-2.

The initial design space is set as [-6, 6] d (with d = 8 or d = 16) and discretized as follow. For each set of design variables belonging to a player, we first generate The data to be recovered, or missing data, are u |Γi and λ∇u.ν |Γi , which are determined as soon as one knows u in the whole Ω. The parameters λ, ϕ and Φ are given functions, ν is the unit outward normal vector on the boundary. The Dirichlet data ϕ and the Neumann data Φ are the so-called Cauchy data, which are known on the accessible part Γ c of the boundary ∂Ω and the unknown field u is the Cauchy solution. Completion/Cauchy problems are known to be severely ill-posed (Hadamard's), and computationally challenging.

Let us assume that (Φ, ϕ)

∈ H -1 2 (Γ c ) × H 1 2 (Γ c ) where H -1 2 (Γ c ) resp. H 1
(SP 1)    ∇.(λ∇u 1 ) = 0 in Ω u 1 = ϕ on Γ c λ∇u 1 .ν = η on Γ i (SP 2)    ∇.(λ∇u 2 ) = 0 in Ω u 2 = ζ on Γ i λ∇u 2 .ν = Φ on Γ c (21)
Let us define the following two costs : for any η ∈ H -1 2 (Γ i ) and ζ ∈ H

1 2 (Γ i ), J 1 (η, ζ) = 1 2 λ∇u 1 .ν -Φ 2 H -1 2 (Γc) + α 2 u 1 -u 2 2 H 1 2 (Γi) , (22)
J 2 (η, ζ) = 1 2 u 2 -ϕ 2 L 2 (Γc) + α 2 u 1 -u 2 2 H 1 2 (Γi) , ( 23
)
where α is a given positive parameter (e.g. α = 1).

Let us remark that, for the problem is severely ill posed, the classical descent algorithms for the minimization of the cost J 1 + J 2 with respect to the couple of variables (η, ζ) do not converge. Besides, the solution is extremely sensitive to small variations of φ and Φ. This originally motivated the above-described game formulation.

The fields u 1 = u 1 (η) and u 2 = u 2 (ζ) are aiming at the fulfillment of a possibly antagonistic goals, namely minimizing the Neumann gap λ∇u 1 .ν -Φ

H -1 2 (Γc)
and the Dirichlet gap u 2 -ϕ L 2 (Γc) . This antagonism is intimately related to Hadamard's ill-posedness character of the Cauchy problem, and rises as soon as one requires that u 1 and u 2 coincide, which is exactly what the coupling term u 1 -u 2 L 2 (Γi) is for. Thus, one may think of an iterative process which minimizes in a smart fashion the three terms, namely Neumann-Dirichlet-Coupling terms. From a game theory perspective, one may define two players, one associated with the strategy variable η and cost y 1 = J 1 and the second with the variable ζ and cost y 2 = J 2 , each trying to minimize its cost in a non-cooperative fashion. The fact that each player controls only his own strategy, while there is a strong dependence of each player's cost on the joint strategies (η, ζ) justifies the use of the game theory framework (and terminology), a natural setting which may be used to formulate the negotiation between these two costs.

Theorem 1 [START_REF] Habbal | Neumann-Dirichlet Nash strategies for the solution of elliptic Cauchy problems[END_REF] There always exists a unique Nash equilibrium

(η * , ζ * ) ∈ H -1 2 (Γ i )× H 1 2 (Γ i ),

Noisy data and random Nash equilibrium

In a realistic situation, the fields ϕ and Φ may be "polluted" by noise. Habbal and Kallel [START_REF] Habbal | Neumann-Dirichlet Nash strategies for the solution of elliptic Cauchy problems[END_REF] showed that the Nash equilibrium is stable with respect to small perturbations, and that the perturbed equilibrium converges strongly to the unperturbed one when the noise tends to zero. However, they did not consider directly the random Nash game:

min η E[J 1 (η, ζ, ϕ , Φ )] min ζ E[J 2 (η, ζ, ϕ , Φ )], (24)
where ϕ and Φ are perturbed values of ϕ and Φ.

Addressing this problem with a classical fixed-point algorithm à la [START_REF] Bas ¸ar | Relaxation techniques and asynchronous algorithms for on-line computation of noncooperative equilibria[END_REF] would require, for each trial pair (η, ζ), repeated evaluations of J 1 and J 2 for many different values of ϕ and Φ , which would prove extremely intensive computationally.

Implementation and experimental setup

The physical domain Ω is taken as a 2D annular structure. The accessible boundary Γ c is the outer circle, with a ray R c = 1 and the inaccessible boundary Γ i is the inner circle with a ray R i = 0.5, see Figure 4. The conductivity coefficient is λ = 1, the flux Φ = 0 and the heat field ϕ is built from an exact known solution (an academic u(x, y) = exp(x) cos(y) used for validation issues).

We use FreemFem++ [START_REF] Hecht | Freefem++ v. 2.11. User?s Manual University of Paris 6[END_REF] to develop our finite element (FE) solvers. The FE computations are performed with a P 1-triangular mesh yielding 1, 088 degrees of freedom, the outer and inner boundaries being discretized each with 60 finite element nodes. From now on, we use the same notations as above for functions to refer to their finite element approximations (values at FE nodes).

As in Section 4.2, to reduce the dimensionality of the problem, the η and ζ being originally vectors of size 60 (the number of nodes at the inner boundary), we interpolate the underlying functions with -natural-splines with κ = 8 coefficients for each quantity, resulting with a decision space of size d = 16 (instead of the original 120). Each coefficient is bounded between -3 and 3, allowing a large variety of spline shapes.

Since the Neumann condition involving Φ is known to be the most sensitive to noise, we perturb only this term with an additive noise:

Φ = Φ 0 + (ξ),
with (ξ) a white noise uniformly distributed between -0.25 and 0.25, thus resulting in a randomized vector of dimension 60 (the number of FE nodes on the outer boundary). The resulting variability on J 1 (η, ζ, Φ ) and J 2 (η, ζ, Φ ) strongly depends on the value of (η, ζ) and can be very large (with coefficients of variation, y i /τ i , from 5% to over 100%).

The original continuous space of the spline coefficients is discretized by generating first two 1000 × 8 LHD and taking all the combinations between the two designs, ending with space of size N = 10 6 potential strategies. The Bayesian optimization strategy is run with n 0 = 50 initial points (chosen by space-filling design) and 150 infill points. Since the noise is very heteroskedastic, we use five repetitions for each trial pair (η, ζ) in order to obtain a rough estimate of the noise variance τ i , and we take the mean over those five repetitions as our observations (f i 's). This results with a total budget of (50+150)×5 = 1, 000 model runs. Let us notice that computing a deterministic NE with fixed-point methods requires about five times this budget on this problem. We followed Algorithm 1 with the SUR criterion, K = M = 20, selecting Card(X sim ) = N sim = 1, 296 simulation points and Card(X cand ) = N cand = 256 candidates.

Fig. 4 The domain is an annular structure, the outer circle is the accessible boundary while the inner circle is the inaccessible one. The plot shows the isovalues of the solution u 1 at convergence.

Results

The algorithm is run five times with different initial points to assess its robustness. The results on the test problem are presented in Figure 5. In the absence of reference (ground truth to be compared to), we show the evolution of the two objectives functions during optimization. All five runs identify very quickly (after 20 iterations, which corresponds to 350 FE evaluations) a similar estimation of the equilibrium, close to the equilibrium of the noiseless problem (which is actually known to be (0, 0), Habbal and Kallel [START_REF] Habbal | Neumann-Dirichlet Nash strategies for the solution of elliptic Cauchy problems[END_REF]). However, it fails at converging finely to a single, stable value, in particular for y 2 . This is not surprising and can largely be attributed to the difficulty of the problem (with respect to the noise level and dimensionality). Hybrid approaches (using fixed-point strategies for a final local convergence) or sophisticated sampling strategies (increasing the number of repeated simulations gradually with the iterations) may be used to solve this issue (although at the price of a considerably higher computational budget).

CONCLUDING COMMENTS

We have proposed here a novel approach to solve random or deterministic Nash games with drastically limited budgets of evaluations based on GP regression, taking the form of a Bayesian optimization algorithm. Experiments on challenging synthetic problems demonstrate the potential of this approach compared to classical, derivative-based algorithms.

On the test problems, the two acquisition functions performed similarly well. P E has the benefit of not relying on conditional simulation paths, which makes it simpler to implement and less computationally intensive in most cases. Still, the SUR approach has several decisive advantages; in particular, it does not actually require the new observations to belong to the grid X sim , such that it could be optimized continuously. Moreover, it lays the groundwork for many extensions that may be pursued in future work.

First, SUR strategies are well-suited to allow selecting batches of points instead of only one, a key feature in distributed computer experiments [START_REF] Chevalier | Fast computation of the multi-points expected improvement with applications in batch selection[END_REF]. Second, other games and equilibria may be considered: the versatility of the SUR approach may allow its transposition to other frameworks. In particular, mixed equilibria, which existence are always guaranteed on discrete problems, could be addressed by using Ψ functions that return discrete probability measures. Generalized NEPs [START_REF] Facchinei | Generalized nash equilibrium problems[END_REF] could be tackled by building on existing works on Bayesian optimization with constraints [see e.g., 31, and references therein].

Finally, this work could also be extended to continuous domains, for which related convergence properties may be investigated in light of recent advances on SUR approaches [START_REF] Bect | A supermartingale approach to gaussian process based sequential design of experiments[END_REF].

A Handling conditional simulations

We detail here how we generate the draws of Y|F i to compute Ĵ(x) in practice. We employ the FOXY (fast update of conditional simulation ensemble) algorithm proposed by Chevalier et al [START_REF] Chevalier | Fast update of conditional simulation ensembles[END_REF], as detailed below.

Let Y 1 , . . . , Y M be independent draws of Y (X) (each Y i ∈ R N ×p ), generated using the posterior Gaussian distribution of Eq. ( 8), and F 1 , . . . , F K independent (of each other and of the Y i 's) draws of Y(x) + ε from the posterior Gaussian distribution of Eq. ( 9). As shown in Chevalier et al [START_REF] Chevalier | Fast update of conditional simulation ensembles[END_REF], draws of Y|F i can be obtained efficiently from Y 1 , . . . , Y M using:

Y (i) j |F (i) k = Y (i) j + λ (i) (x) F (i) k -Y (i) j (x) , (25)
with 1 ≤ i ≤ p, 1 ≤ j ≤ M , 1 ≤ k ≤ K and

λ (i) (x) = k (i) n (x, X) k (i) n (x, x)
.

Notice that λ (i) (x) may only be computed once for all Y (i) j (x).

B C(x) formulae

For a given target T E ∈ R p and x ∈ X: Then, the probability to belong to the box is:

Ctarget(x) = p i=1 φ T Ei -µ i (x) σ i (x) , (26)
C box (x) = p i=1 Φ T U i -µ i (x) σ i (x) -Φ µ i (x) -T Li σ i (x) . (27)

C Solving NEP on GP draws

We detail here a simple algorithm to extract Nash equilibria from GP draws.

D Computational time

We report here the computational time required to perform a single iteration of our algorithm for each of the three examples (not including the time required to run the simulation itself). Experiments were run on an Intel R Core TM i7-5600U CPU at 2.60GHz with 4 × 8GB of RAM.

  and we have N := Card(X) = p i=1 m i .

2 ( 1 2

 21 Γ c ) are the Sobolev spaces of Neumann resp. Dirichlet traces of functions in the Sobolev space H 1 (Ω) [see e.g. 1]. For given η ∈ H -1 2 (Γ i ) and ζ ∈ H (Γ i ), let us define u 1 (η) and u 2 (ζ) as the unique solutions in H 1 (Ω) of the following elliptic boundary value problems :

  and when the Cauchy problem has a solution u, then u 1 (η * ) = u 2 (ζ * ) = u, and (η * , ζ * ) are the missing data, namely η * = λ∇u.ν |Γi and ζ * = u |Γi .

Fig. 5

 5 Fig. 5 Convergence of the estimated Nash equilibrium for five algorithm runs. Each iteration corresponds to five calls of the FreeFem++ model.

  with φ the probability density function of the standard Gaussian variable.Let T L ∈ R p and T U ∈ R p such that ∀1 ≤ i ≤ p, T Li < T U i define a box in the objective space. Defining Ψ = [Ψ (Y 1 ), . . . , Ψ (Y M )] the p × M matrix of simulated NE, we use:∀1 ≤ i ≤ p T Li = min Ψ i,1...M and T U i = max Ψ i,1...M .

Table 1

 1 P1 convergence results.

	Strategy	Evaluations required Success rate
	P E	9-10	5/5
	SUR	8-14	5/5
	Fixed point	200-1000	3/5

Table 3

 3 Algorithm 3 Pseudo-code for Nash equilibria extractionRequire: P : number of players, Y: draw of Y(X) of size N × P , I = {1, . . . , N } 1:for 1 ≤ i ≤ N do 2: if i ∈ I then 3:for 1 ≤ j ≤ P do Average CPU times required for one iteration of the GP-based algorithm on the different test problems.

	4:	Find K = 1 ≤ k ≤ N, x k -j = x i -j
	5:	Find l * = min 1≤l≤|K|	Y(x l ) j
	6:	for 1 ≤ l ≤ |K| do
	7:	if Y(x l ) j > Y(x l * ) j then
	8:		I ← I \ l
	9:	end if
	10:	end for
	11:	end for
	12:	end if
	13: end for
	Ensure: I
	Case / Criterion	Pnash SUR
		Case 4.1	4s	20s
	Case 4.2 (κ = 1)	13s	40s
	Case 4.2 (κ = 2)	28s	82s
		Case 4.3	9s	12s

Acknowledgements The authors acknowledge inspiration from Lorentz Center Workshop "SAMCO -Surrogate Model Assisted Multicriteria Optimization", at Leiden University Feb 29 -March 4, 2016. Mickaël Binois is grateful for support from National Science Foundation grant DMS-1521702.

a 17-point LHD. For κ = 1, this means four LHD in the spaces (x 1 , x 2 ), (x 3 , x 4 ) and so on. Then, we take all the combinations of strategies, ending with a total of Card(X) = N = 17 4 = 83, 521 possible strategies.

In both cases, we followed Algorithm 1 with n 0 = 80 or 160 initial design points and a maximum budget of n max = 160 or 320 evaluations (depending on the dimension). We chose Card(X sim ) = N sim = 1296 simulation points and Card(X cand ) = N cand = 256 candidates. For the number of draws to compute the SUR criterion, we chose K = M = 20.

A fixed-point algorithm based on finite differences is also ran, and the experiment is replicated five times with different algorithm initializations. Note that, here, the actual NE is unknown beforehand. However, since all runs (including the fixed-point algorithm ones) converged to the same point, we assume that it is the actual NE.

We show the results in Table 2. For the lower dimension problem (d = 8), all runs found the solution with less than 100 evaluations (that is, 80 initial points plus 20 infills). This represents about 1% of the total number of possible strategies. In this case, SUR appears as slightly more efficient than P E . For the higher dimension problem (d = 16), more evaluations are needed, yet all runs converge with less than 240 evaluations. As a comparison, on both cases the fixed-point algorithm is more than 60 (resp. 20) times more expensive.

PDE-constrained example: data completion

Problem description

We address here the class of problems known as data completion or data recovery problems. Let be Ω a bounded open domain in R d (d = 2, 3) with a sufficiently smooth boundary ∂Ω composed of two connected disjoint components Γ c and Γ i , with the latter being inaccessible to boundary measurements. For details, see [START_REF] Habbal | Neumann-Dirichlet Nash strategies for the solution of elliptic Cauchy problems[END_REF] whence the present example is excerpt.

Let us focus, for illustration, on the particular case of steady state heat equation. The problem is formulated in terms of the following elliptic Cauchy problem : [START_REF] Ginsbourger | Towards Gaussian process-based optimization with finite time horizon[END_REF]