Efficient spatial designs using Hausdorff distances and Bayesian optimization

This is a self-archived version of an original article. This version may differ from the original in pagination and typographic details. Author(s): Title: Year: Version: Copyright: Rights: Rights url: Please cite the original version: CC BY-NC-ND 4.0 https://creativecommons.org/licenses/by-nc-nd/4.0/ Efficient spatial designs using Hausdorff distances and Bayesian optimization © 2021 the Authors Published version Paglia, Jacopo; Eidsvik, Jo; Karvanen, Juha Paglia, J., Eidsvik, J., & Karvanen, J. (2022). Efficient spatial designs using Hausdorff distances and Bayesian optimization. Scandinavian Journal of Statistics, 49(3), 1060-1084. https://doi.org/10.1111/sjos.12554 2022 Received: 3 February 2020 Revised: 29 March 2021 Accepted: 30 June 2021 DOI: 10.1111/sjos.12554 ORIGINAL ARTICLE Efficient spatial designs using Hausdorff distances and Bayesian optimization Jacopo Paglia1Jo Eidsvik1Juha Karvanen2 1Department of Mathematical Sciences, Norwegian University of Science and Technology, Trondheim, Norway 2Department of Mathematics and Statistics, University of Jyvaskyla, Jyväskylä, Finland Correspondence Jacopo Paglia, Department of Mathematical Sciences, Norwegian University of Science and Technology, Trondheim 7034, Norway. Email: [email protected] Funding information Jacopo Paglia’s and Jo Eidsvik’s work are supported by the KPN project 255418/E30: "Reduced uncertainty in overpressures and drilling window prediction ahead of the bit (PressureAhead)", of the Norwegian Research Council and the DrillWell Centre (AkerBP, Wintershall, ConocoPhillips and Equinor). Juha Karvanen’s work is supported by Grant number 311877 "Decision analytics utilizing causal models and multiobjective optimisation" (DEMO), of the Academy of Finland. Abstract An iterative Bayesian optimization technique is presented to find spatial designs of data that carry much information. We use the decision theoretic notion of value of information as the design criterion. Gaussian process surrogate models enable fast calculations of expected improvement for a large number of designs, while the full-scale value of information evaluations are only done for the most promising designs. The Hausdorff distance is used to model the similarity between designs in the surrogate Gaussian process covariance representation, and this allows the suggested algorithm to learn across different designs. We study properties of the Bayesian optimization design algorithm in a synthetic example and real-world examples from forest conservation and petroleum drilling operations. In the synthetic example we consider a model where the exact solution is available and we run the algorithm under different versions of this example and compare it with existing approaches such as sequential selection and an exchange algorithm. KEYWORDS Bayesian optimization, decision-making, Hausdorff distance, value of information Abbreviations: EI, expected improvement; GP surrogate, Gaussian process; PoV posterior value; PV prior value; VOI, value of information This is an open access article under the terms of the Creative Commons Attribution-NonCommercial-NoDerivs License, which permits use and distribution in any medium, provided the original work is properly cited, the use is non-commercial and no modifications or adaptations are made. © 2021 The Authors. Scandinavian Journal of Statistics published by John Wiley & Sons Ltd on behalf of The Board of the Foundation of the Scandinavian Journal of Statistics. Scand J Statist. 2021;1–25. wileyonlinelibrary.com/journal/sjos 1 2PAGLIA et al. 1INTRODUCTION This paper is inspired by challenging decision situations in the earth and environmental sciences. In these situations, data are gathered to support decisions about resource management. Data acquisition and processing is often costly, and it is then important to choose the sampling design wisely. There exist several common design or information criteria, see for example, Ryan et al. (2016) for a recent review. For decision-makers, value of information (VOI) analysis is useful in this context (Abbas & Howard, 2015; Eidsvik et al., 2015), as it is directly connected with the information gain associated with the decision situation and it provides a bound on the expected monetary amount one should be willing to pay for data to aid in resolving this decision situation. We focus on designing experiments in spatial domains. Here, Kriging interpolation (Stein, 2012) is often used to propagate the effects of observations based on spatial correlations. When choosing the criterion for the selection of the sampling design it is important to keep in mind the scope of the experiment. One could, for example, be interested in choosing a design that spread well across the domain, also called spatially balanced designs (Grafström et al., 2012; Stevens Jr & Olsen, 2004). In the study presented here, however, the focus is on finding a design that maximizes the VOI, and where the spatial balance may arise as a consequence of the spatial modelling of the variables of interest. Using VOI analysis, we aim to provide the decision-maker with efficient survey designs including the optimal number of measurement locations and their spatial configuration. We assume that the spatial domain is discretized to a grid so that there is a finite set of possible observation locations. Moreover, we limit scope to static designs (Diggle & Lophaven, 2006; Dobbie et al., 2008; Huan & Marzouk, 2013), where the experimental configuration is selected once, at the onset of data gathering. The alternative is sequential data gathering, where the design can be adapted based on the observations made in the first (batches of) measurements (Binois et al., 2019; Drovandi et al., 2013; Eidsvik et al., 2018), but this is not always possible in practical experimental planning, which must comply with project management and budgetary limitations. As pointed out by several others, this design problem is not trivial as the number of possible designs grows combinatorially fast. Royle (2002) proposed a random exchange algorithm to search for the optimal design. García-Ródenas et al. (2020) presented an interesting overview of some of the main algorithms for finding efficient designs. Weaver et al. (2016) and Overstall and Woods (2017) applied Bayesian optimization to focus the search for good designs. We use a Gaussian process (GP) surrogate model enabling fast computation of the expected improvement (EI) in Bayesian optimization. This is combined with techniques from search algorithms, to find efficient spatial designs. As was also regarded as a possibility in Ginsbourger et al. (2016) in the context of computer experiments, the current paper presents an approach for using the Hausdorff distance between various designs. The contribution of our work is using this to correlate outcomes of similar site configurations, within a realistic statistical model, and combining this in an algorithm for quickly locating valuable sampling designs. Even though our focus is on spatial decision situations and design, we believe that this approach can also be applicable to other big-data challenges (Drovandi et al., 2017) and active learning approaches (Bouneffouf, 2016; Settles, 2012), where the challenge is more related to which data to process for learning and improved classifications. In Section 2 we describe the spatial design problem in mathematical detail and define the VOI criterion which we use as a practically relevant information measure. In Section 3 we outline the Bayesian optimization approach using Hausdorff distances to borrow information among PAGLIA et al. 3 North I I I I II II II III III III III III East FIGURE 1 Illustration of a spatial domain split in 40 regional units of varying size. Three different designs are indicated (Design I, II, and III) of different cardinality and spatial allocation [Colour figure can be viewed at wileyonlinelibrary.com] similar designs. In Section 4 we study the properties of the methodology via simulations. In Section 5 we show results on forestry and petroleum examples to demonstrate possible applications of the methods. Section 6 has closing remarks on the methodological contributions presented here, including viable opportunities future work. 2SPATIAL DESIGN OF EXPERIMENTS 2.1 Spatial survey designs We consider a situation as illustrated in Figure 1, with a spatial phenomenon allocated to a two-dimensional domain divided in grid cells or sites. The distances between the sites are defined as the Euclidean distances between the centers of the cells. The approach presented in the paper can be extended to higher dimensions with minor changes. The spatial variables of interest are represented at nsites, denoted s1,…,snwith si= (northi,easti),i=1,…,n. In our applications, these sites have a particular interest to the decision-maker. For instance, in the forestry example, the governmental institute must choose at each of the nsites whether this forest unit should be harvested or left for conservation. Because there is much at stake and uncertain outcomes, the decision-maker is likely to benefit from doing surveys at (a subset of) the sites. Data can be gathered at any of the nsites in our description, and a design defines a subset of these nsites where the data collection will be conducted (other cases can be constructed similarly, see e.g. Section 5.2). The possible spatial designs then include no sites, single sites, couples, triplets, and so on, up to all nsites in the design. We denote these by 𝒟=⋃n i=0𝒟i, defined by; 𝒟0=∅,no sites in design, 𝒟1={(s1),(s2),…,(sn)},onesiteindesign, 4PAGLIA et al. 𝒟2={(s1,s2),(s1,s3),…,(sn−1,sn)},two sites in design, ⋮⋮ 𝒟n={(s1,…,sn)},all sites in design. There are npossible designs of cardinality one, (n 2)possible designs of cardinality two, etc. This means that there are 2npossible designs in 𝒟. We will further denote a general design by D∈𝒟and its cardinality by |D|. The sites in this design are then sD,1,…,sD,|D|. The number of sitessharedbydesignsCand Dis |C∩D|, while the number of sites in at least one of the designs is |C∪D|. In our setting we compare the information gain obtained by different designs, and it makes sense that similar spatial designs contain almost the same information. In Figure 1 three different designs are indicated (I, II, and III). Designs I and II appear very similar in the spatial allocation of survey sites even though they have different cardinalities (three and four). Most likely, Design I will not have much to offer over Design II, unless there is much noise in the data or large gain in capturing additional covariate information which could be important for predictive purposes. Say, in the forestry example, a biologist would spend time doing one more experiments in Design I, at an extra cost. But unless she learns substantially more about the model, there is not much additional spatial information in Design I compared with doing just the three measurements in Design II. The last survey plan, Design III, is spatially very different from the others because it allocates the measurements in the central parts of the domain. The value of this design could be very different from that of Design I and II. To find the optimal design one must evaluate the information gain and cost for all possible design sets, but in practice one can only evaluate it for a fraction of all possible designs. We suggest a statistical approach for this optimization problem, where we utilize the similarity of spatial designs to estimate the information gain. 2.2 Value of information The goal of spatial design of experiments is to choose a valuable survey plan for information gathering. This choice must balance expected information gain with the cost of data acquisition and processing. To evaluate the expected information gain associated with designs, one must formulate a value or utility function. In the applications that we consider here, it is relatively straightforward to relate the question about information gain to an underlying decision situation, meaning that data are only valuable when their outcome can materialize in different decisions. For instance, in the forestry example the underlying decision is to conserve forest units or not, and data can help the decision-maker to decide one or the other, depending on what the information reveals. Managers are further often willing to phrase these decision situation in terms of monetary units, and then the VOI which gives the expected gain in information is directly comparable to the cost of data gathering. If the VOI exceeds this cost, the experiment is worthwhile and the decision-maker should commit to gather the information, if the budget permits the cost. We next define the VOI formally through a model for the random variables of interest, the decision alternatives and the information gathered by a chosen design. The variables of interest are denoted by x=(x1,…,xn),wherexi=x(si),i=1,…,n.Inour context these are directly tied to a decision situation and connected to economic values. For PAGLIA et al. 5 instance, they can be random profits allocated to forest units or loss associated with a drilling operation. Note that other parameters will be important in the statistical modeling of the phenomenon of interest, such as regression parameters and covariance function parameters, but in our setting they are only used in the construction of a realistic statistical model for the phenomenon that is studied, and in particular for the variables of interest x. Assuming a continuous sample space for the variables of interest, we denote its probability density function by p(x),with marginal density p(xi)for each sites si. The decision alternatives are generally denoted by a∈,whereis the set of all possible alternatives. In some situations, the alternatives decouple (Eidsvik et al., 2015), involving for instance local decisions about harvesting units in our forest conservation example. In general, the prior value (PV), without any additional information, is defined as the value from doing the optimal decisions. Assuming a risk-neutral decision maker (Abbas & Howard, 2015), the PV is calculated from expected values as follows; PV =max a∈{E(𝜈(x,a))},E(𝜈(x,a)) = ∫𝜈(x,a)p(x)dx.(1) Here, 𝜈(x,a)represents the value function, which could be quite general, but in our application it is the monetary profits associated with choice a∈when the variable outcome is x.Inthe forestry example, the decision-maker will choose to conserve the sites that have high preservation value, while the others are harvested. It is difficult to make decisions under uncertainty, and one can choose to purchase information that facilitate decision-making. We here let yDdenote the data gathered by design D∈𝒟. This data is relevant to the decision situation in the sense that it will provide information about the variable of interest x. In the applications below, the model for data is given as a conditional probability density or mass function p(yD|x), and the marginal model for data is then p(yD)=∫p(yD|x)p(x)dx. When the data are available, the conditional value (CV) is CV(yD)=max a∈{E(𝜈(x,a)|yD)},(2) and the expected posterior value (PoV) before the data gathering is obtained by taking the expectation of expression (2) over the possible data outcomes: PoV(D)=EyD[CV(yD)]=EyD[max a∈{E(𝜈(x,a)|yD)}].(3) The VOI is defined as the difference between the expected PoV in (3) and the PV in (1): VOI(D)=PoV(D)−PV.(4) The goal is to choose a valuable design D. Keeping in mind that data comes with a cost, we should compare the VOI with the cost C(D)of design D. This means that the objective is to optimize D∗=argmaxD∈𝒟I(D),I(D)=VOI(D)−C(D).(5) 6PAGLIA et al. Other objectives are of course possible. For instance, a decision-maker might have a fixed budget for the design, and the goal would then be to maximize the VOI among all designs that have a cost less than the budget. With large opportunities for data gathering, it is extremely difficult to find the optimal design. First, the complexity grows extremely fast with the number of sites. Second, in common settings, the calculation of the information design criterion in (5) for a fixed design typically requires quite a bit of computational effort as is emphasized by the complexity of the integral maximum expectation expressions required in (3). In practice one must often turn to heuristic approaches to such design problems (García-Ródenas et al., 2020. We suggest to use a statistical approximation strategy that evaluates I(D)in (5) only for a few promising designs which are extracted by a fast Bayesian optimization approach building on GP surrogate models and EI. 3BAYESIAN OPTIMIZATION FOR DESIGNS We develop a Bayesian optimization approach to guide the search for the maximum of I(D) in (5). We combine computational search algorithms with the EI acquisition criterion to select which designs to evaluate in an iterative optimization workflow. In doing so, we suggest to model the information measure I(D)using a GP surrogate model. This is in line with the common approaches for Bayesian optimization (Brochu et al., 2010; Frazier, 2018). The benefits of using a GP surrogate for the information measure is that it enables: •efficient model updating based on evaluations (Section 3.1), •learning across different but similar designs (Section 3.2), •computing EI in closed form, to focus on evaluating promising designs (Section 3.3), •framing a useful algorithmic description of the overall procedure (Section 3.4). 3.1 GP surrogate The information gain I(D)is represented by a GP surrogate model. This relies on mean and variance–covariance specifications of the information gain for input designs. In the current setting with Bayesian optimization, the GP surrogate model is updated sequentially when more evaluations become available. When mdesigns D(1),…,D(m)have been evaluated, the knowledge is denoted = {(I(j),D(j)); j=1,…,m}. By standard multivariate Gaussian theory, the conditional distribution for the information measure at design Dis then Gaussian with mean and variance 𝜇(D;)=𝜇+kt D,K−1 (I()−𝜇1), 𝜎2(D;)=𝜎2(1−kt D,K−1 kD,).(6) Here, I()=(I(1),…,I(m))tis the length mvector of information gain evaluations, Kthe m×mcorrelation matrix between evaluations of designs, kD,the length mvector of correlations between the evaluations and the information gain for design D,and1is a length mvector of 1 entries. The representation requires specification of the mean 𝜇and variance 𝜎2,whichare PAGLIA et al. 7 assumed constant for all designs. It further needs a valid correlation function specification K(C,D) between two different designs Cand D(see Section 3.2). 3.2 Distance between designs The correlation function gauges the similarity between designs, as defined via kD,and Kin (6). This specification of a correlation function is a common task in spatial statistics and Bayesian optimization over a regular input space. In our setting with spatial designs, it is not obvious how to assign this correlation function, and a main contribution of this paper is to formulate a distance measure between designs which is useful in the context of Bayesian optimization. Our proposed distance measure for this task is the Hausdorff distance which is presented next, but we also outline other distance measures below to discuss this topic in a more general context. Throughout this description, we consider two general designs D=(sD,1,…,sD,|D|)and C=(sC,1,…,sC,|C|). For two sites siand sj,welet‖‖si−sj‖‖be the Euclidean distance between the two sites. The Hausdorff distance is commonly used to measure the distance between curves, images, or point sets (Huttenlocher et al., 1992). In our context it represents the maximum of the minimal distances from sites in one set to sites in the other set, and it hence measures similarity of designs: h=distH(D,C)=max {hH(D,C),hH(C,D)},(7) hH(D,C)= max i=1∶|D|{min j=1∶|C|‖‖sD,i−sC,j‖‖}.(8) Figure 2 illustrates several designs of size 1, …4. For each subplot the maximum distances from sites in one set to the other is calculated and shown. One design Dis marked as circle, the other design Cis marked as cross. The Hausdorff distance in (7) is printed in the displays, and hH(D,C)and hH(C,D)are indicated. We note that in some cases the maximum distances from one set to the other are identical (upper right display and bottom middle display), but for most of these site configurations this symmetry is not present. For instance, in the upper left display, the circle is relatively close to the southernmost site in the cross set, but the northernmost site in the cross set is quite far from the circle. Similarly, in the centre display, both sites in the circles set are close to a site in the cross set, but one site in the crossset is far from the closest site inthe circle set. In the bottom-middle display the designs are rather similar, and the distance is small (h=0.113). In all the right displays, the designs are very different, and the distances are large. Based on Hausdorff distances for sets like that displayed in Figure 2, it seems to be a useful way to measure the difference between designs. We next present some alternative distances (Fujita, 2013; Min et al., 2007) that could be useful for our purpose, and discuss their pros and cons. One alternative distance is defined by the minimum of the distances between design sites of the two sets. The main problem with this distance is that designs with sites in common are not separated because the distance in this case will be zero. In addition, this distance is not a proper metric because the triangular inequality does not hold. We will hence discard this distance—it is not suitable for our purpose. Similarly it is possible to define another distance considering the maximum of the distances between designs sites instead. Again, this is not a proper distance metric, and it is not convenient for our measure of similarity because the distance between two equal sets is greater than 0. Yet another candidate is the Jaccard distance (Levandowsky & Winter, 1971) which is defined via the relative counts of 8PAGLIA et al. FIGURE 2 Hausdorff distance between various designs, his the Hausdorff distance between the two sets, marked as circle for Dand cross for C. The solid lines represent the maximum of the minimal distances between Cand D, while the dashed line represents the maximum of the minimal distances between Dand C[Colour figure can be viewed at wileyonlinelibrary.com] sites that are not shared in the designs. We believe this could be a very sensible way to consider dissimilarity between designs for features or class covariates, but not for our type of applications because it does not account for the spatial distance between sites in the sets. Fujita (2013) describes another metric based on the average distances between the two designs: distF(D,C)= 1 |D∪C||D|∑ sD,i∈D∑ sC,j∈C∖D‖‖sD,i−sC,j‖‖+1 |D∪C||C|∑ sD,i∈D∖C∑ sC,j∈C‖‖sD,i−sC,j‖‖.(9) This metric seems to work sensibly for our purpose, and we study the possibility of using distFin Section 4. Yet another possibility is the modified Hausdorff distances. A variant is defined as the average of minimum (or alternatively the minimum squares) distances. Dubuisson and Jain (1994) show that the modified Hausdorff distance is a valid tool for object matching. The problem with this measure in our application is that it smooths the effect of outlier sites, whereas we believe that even a single outlier site could add valuable information to the design, giving knowledge of a larger area, and that should then have an important impact on the distance. In summary, we use the regular Hausdorff distance hin (7) to model design dissimilarities. When building the covariance matrix in the GP surrogate formulation (6) it is important to PAGLIA et al. 15 -6 -4 -2 0 2 4 6 8 10 x -6 -4 -2 0 2 4 6 8 y Best design Sequential selection Algorithm 1 Exchange algorithm FIGURE 4 Representation of the spatial Hausdorff distance for the best 3000 designs (grey) in a two-dimensional Euclidean space using multidimensional scaling technique. The pink star represents the best design and the green diamond the best design from the sequential selection. The red dots and lines represents the path of Algorithm 1 while the blue crosses and lines represents the results from the exchange algorithm [Colour figure can be viewed at wileyonlinelibrary.com] 4.7 4.8 4.9 5 VOI(D)-C(D) 104 0 0.5 1 1.5 2 2.5 3 Density 10-3 Algorithm 1 Exchange algorithm (a) 4.7 4.8 4.9 5 VOI(D)-C(D) 104 (b) 4.65 4.7 4.75 4.8 4.85 4.9 4.95 5 VOI(D)-C(D) 104 (c) FIGURE 5 Comparison of the performances of Algorithm 1 (red) and exchange algorithm (dashed blue) over 100 replicate restarts. The Bayesian optimization approach is able to get large values of I(D)after few iterations. When the number of evaluations grows the exchange algorithm starts to perform well. (a) 250 evaluations of I(D)(in €); (b) 500 evaluations of I(D)(in €); (c) 800 evaluations of I(D)(in €) [Colour figure can be viewed at wileyonlinelibrary.com] 16 PAGLIA et al. TABLE 1 Sensitivity analysis of different inputs on the algorithm performance. The column “Algorithm 1” represents the highest I(D)among the 10 replicates, with its rank in parenthesis, the column “%best 100” indicates the fraction of times we get a score among the best 100 information gain values over the replicates. Note that low and high 𝜎xrepresent different models with varying value of information and hence also different values of I(D) 𝝈2 xParameter re-estimation Iterations (Tmax)Algorithm1%best 100 Low off 5 €33,011 (3) 100% High off 5 €85,740 (1) 20% Low on 5 €33,012 (1) 100% High on 5 €84,486 (18) 60% Low off 15 €33,011 (3) 100% High off 15 €85,365 (3) 100% Low on 15 €33,012 (1) 100% High on 15 €85,740 (1) 100% Overall, higher prior variability appears to be more difficult for the algorithm, especially with few iterations, where only 20 and 60% of the runs resulted in the top 100 ranking designs. Still, one lucky restart run without any re-estimation ended up with the highest rank. The performance clearly increases with iterations since all restarts are in the top 100 ranking after 15 iterations, no matter prior variance high/low or re-estimation on/off. Re-estimation of parameters gives improved performance, especially when there is much prior uncertainty in the profits model. The running time is a bit more than three times larger for 15 iterations compared with 5, because of the growing matrix expressions in the GP surrogate calculations. Doing parameter re-estimation after every batch also takes some additional time, but the empirical variogram calculation and least squares matching is very fast. The basic assumption of our work is that for real-world settings, most of the computer time is spent on evaluating I(D), and the number of such evaluations is considered to be the main computational restriction. 5EXAMPLES 5.1 Forestry This example regards forest management and conservation (Eyvindson et al., 2017; Kangas et al., 2008). In this application the decision-maker must choose to conserve forest stands or not. The decision-maker is here a governmental institute that has a budget for conservation. The forest stands are owned by private owners who may harvest the timber unless the forest is conserved. In order to conserve a forest stand, the institute must pay a compensation (bi) to the forest owner. When a forest stand is conserved, the ecological benefit (ri) is proportional to a biodiversity indicator. The study is inspired by data analyzed by Eyvindson et al. (2019). The data consist of 70 forest stands (sites) of various size from the Satakunta region in southwest Finland (Figure 6a). Each stand is classified according to the age class (1: <=80, 2: 81 −95, 3: 96 −110, 4: >110 years). Figure 6b sketches forest stands, with colors identifying the different age groups. PAGLIA et al. 17 (a) East North 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 <=80 >110 (b) FIGURE 6 Study area for forest conservation. (a) The region of interest is located in the southwest of Finland. The red identify the geographic location of the Satakunta region (Wikimedia Commons, 2010). (b) Forest stands of different size and age numbered from 1 to 70. Each color identify a different age group, orange: <=80, yellow: 81–95, violet: 96–110, green: >110 years. The hatched regions correspond to the best design [Colour figure can be viewed at wileyonlinelibrary.com] The possible alternatives for the decision-maker are ={ai;i=1,…,n},whereai= {not conserve,conserve}={0,1}at stand i(Eyvindson et al., 2019). We let xibe the log-intensity of the number of wood inhabiting fungi at stand i=1,…,n. This number is a commonly used biodiversity indicator. The value function is then 𝜈(xi,ai)={0ai=0, rexi−biai=1.(15) Thelog-intensitiesarehere modeled with a multivariateGaussiandistributionoverthestands. In doing so, the age of the forest stand is treated as a covariate zin the simulation study. The mean and the covariance matrix of log intensities vector variable xare computed by double mean and variance over the regression uncertainty (Section 4). Designs are constructed to gather information that can assist the decision maker. There are age-dependent inventory costs, so it is important to plan wisely and obtain effective designs at a low overall cost. The measurements of species richness in fungi are defined with a Poisson likelihood function yi|xi∼Poisson (exi),(16) assuming conditional independence between the stands and constant area for each inventory. The VOI is here defined by VOI(D)= n ∑ i=1 EyD[max{0,E(rexi−bi|yD)}]− n ∑ i=1 max{0,E(rexi−bi)},(17) 18 PAGLIA et al. and the information gain I(D)is obtained as the difference VOI(D)−C(D)where C(D)denotes the inventory costs, obtained by accumulating costs over all design sites. The inventory cost depends on the age of the forest (1: €5100; 2: €5800; 3: €6200; 4: €5600). We use the method developed by Evangelou and Eidsvik (2017) to compute the VOI in (17). This is an approximation where the evaluation of VOI(D)relies on iterative matrix linearizations and refitting Gaussian approximations going into the Laplace approximation. Even though this evaluation can be done in reasonable time for a candidate design, it is relatively time-demanding, and as there is a total of 1.18 ⋅1021 possible designs, it is not feasible to calculate the VOI for all of them to find the optimal design. Instead, we use the suggested Bayesian optimization method to find efficient designs. We initiate the algorithm by evaluating m0=50 random designs of various cardinalities. The GP surrogate model parameters (𝜎,𝜃H)are specified for each of 10 such re-starts. Based on this, the approximate 20 percentiles of these parameters are (€76,750, 0.25) and 80 percentiles are (€136,140, 5.4). At each batch of size m=50, new evaluations are selected using the EI acquisition function. The parameters are re-estimated at each batch, and after 15 iteration the approximate 20 percentiles are (€42,430, 0.07) and 80 percentiles are (€135,500, 1.27). The results show that one gets lower SDs and correlation range by re-estimation at each batch. In the actual optimization, the re-estimation tends to give slightly faster improvements for I+in the algorithm. For the average value 𝜇of the GP surrogate model for I(D), the initial approximate 20 percentile is €325,320 and 80 percentile is €354,270. After 15 iterations the approximate 20 percentile is €301,520 and 80 percentile is €352,430. It is perhaps surprising that the 20 percentile decreases when the goal is to maximize the function. However, the exploration elements of the algorithm can also lead to the evaluation of very poor designs, and possibly outliers with very small I(D). The results from the Bayesian optimization are shown in Figure 7. Here, we plot the optimum evaluation so far, for each restart (different lines). Even though we do not know the optimal solution in this case, the results improve over batches and this indicates that the algorithm finds efficient designs within a few batch iterations. Table 2 provides a list of the top five largest values of I(D)obtained by running the algorithm, together with the associated design. The best designs have many sites in common, even though the designs have different cardinalities. In this way the algorithm spots the sites that carry more information. The highest value of I(D)corresponds to a rather large set of |D|=26. This design is illustrated using the hatched areas in Figure 6b, where we observe that the stands of the design tend to spread and cover both the geographical region and also the various age levels. Similar to what was done in the simulation study, we also run both the exchange algorithm and the sequential selection method. The exchange algorithm gets a largest replicate information gain of only I(D)= €491,760 after 800 evaluations, and is not doing so well in this case. The sequential selection algorithm gives I(D)=€552,930 with 2485 evaluations. The associated design is D=(1,2,3,4,5,8,15,16,19,24,26,32,40,41,44,47,49,57,59,60,66,69).Inthis example the sequential selection algorithm hence performs better than the iterative Bayesian optimization in Table 2, at a cost of extra VOI evaluations to find the sequential solution. With this in mind, we added the sequential solution to the evaluations of the Bayesian optimization method, and continued to run that algorithm. We then achieved slightly larger information gain for designs very similar to the one detected with the sequential search, but no significant improvement. We hence suspect that the sequential method gives a near optimal solution for this example. PAGLIA et al. 19 0 5 10 15 Batch iteration 4.2 4.3 4.4 4.5 4.6 4.7 4.8 4.9 5 5.1 VOI(D)-C(D) 105 FIGURE 7 Performance of the Bayesian optimization algorithm for the application in forest conservation. The algorithm gives improved results over batches for the replicate restarts [Colour figure can be viewed at wileyonlinelibrary.com] TABLE 2 The best five designs obtained over 10 re-start replicates of the algorithm in the forestry example, listed in descending order Design D I(D) 1, 7, 8, 11, 12, 18, 19, 24, 28, 30, 35, 38, 40, 41, 42, 43, 46, 47, 48, 49, 50, 51, 55, 59, 64, 67 €507,420 3, 11, 12, 13, 15, 16, 17, 18, 21, 22, 26, 28, 30, 32, 39, 44, 47, 48, 49, 50, 52, 56, 57, 59, 62 €500,320 1, 4, 6, 8, 9, 11, 13, 15, 16, 17, 20, 24, 25, 26, 27, 28, 31, 32, 34, 44, 45, 48, 49, 50, 51, 53, 57, 62, 63, 66, 68 €495,860 2, 5, 6, 8, 9, 13, 14, 17, 19, 20, 23, 24, 30, 32, 35, 37, 40, 42, 44, 46, 49, 51, 53, 54, 56, 57, 61, 62, 66, 69 €494,190 8, 9, 12, 13, 18, 19, 22, 23, 26, 27, 32, 33, 36, 38, 40, 42, 51, 52, 53, 56, 60, 61, 62, 65, 68 €482,710 5.2 Petroleum drilling risks This example regards decision making during drilling operations in the petroleum industry (Lothe et al., 2019; Paglia et al., 2020). We study a drilling situation in the Alvheim oil field located in the central part of the North Sea, on the Norwegian continental shelf (Figure 8a). The field is divided in 68 compartments (Figure 8b) separated by faults. The circles in Figure 8b represent wells, and Figure 8c highlights these to indicate the decision site (black star) and the potential data 20 PAGLIA et al. Alvheim North Sea Oslo (a) (b) (c) FIGURE 8 The study area is an offshore oil and gas field in the central part of the North Sea. (a) Geographical location of Alvheim. The rectangle indicates the position of the field. (b) Map view of the oil field, circles indicate the locations of wells for data gathering. Different colours are used to identify different geological compartments. (c) Three-dimensional view of the location of the measurements in red and the decision site in black star [Colour figure can be viewed at wileyonlinelibrary.com] gathering site (red dots) in a three-dimensional plot. In the following we describe this decision situation and the opportunities for data gathering to make improved decisions. Drilling operations at the Alvheim field are characterized by the risk of overpressure, which occurs when the pore pressure in the rock exceeds the hydrostatic pressure. To prevent big hazards, the drilling mud pressure must be calibrated. We study a specific layer, located at about 3700 m depth, as marked in black in Figure 8c. This layer is composed of mainly shale rocks and believed to be at drilling risk. The decision maker, which is the petroleum company in this case, must decide if it is safe enough to just keep drilling, or if they should set casing to strengthen the well because of a high risk of blowout. The alternatives are ={0,1}= {keep drilling,set casing}. To set casing is an expensive operation and it will reduce the borehole diameter, so the decision-maker is interested in trying to postpone this operation, if not necessary because of very high risk. The value function of the decision is 𝜈(x,a=0)=−c0(LB(x)−x)and 𝜈(x,a=1)=−c1(LB(x)−x),wherexis the unknown pore pressure variable and LB is the lower bound of the mud weight drilling window. Here, c0is the cost when one keeps drilling, while c1is the additional cost of casing. We note that the costs stretch the value functions so that it becomes more valuable to set casing instead of continued drilling for some values of pore pressure. Critically, the mud weight is used during the drilling of a well to exert a pressure on the borehole wall and avoid well collapse, and it is difficult to make decisions when the pore pressure is not known. Please keep in mind that there are a number of other parameters that would also affect LB, but pore pressure is an important parameter that is always taken into consideration during drilling operation (Moos et al., 2004). Figure 8c (red) shows the possible measurement sites. The design will entail any combination of these sites, and the measurements gathered at the design sites will be informative of the pore pressure where they are made and at other sites via the statistical model formulation. There are five wells where accurate measurements of pore pressure can be gathered in four different layers. The cost of data acquisition is assumed to be the same for each well. However it will be cheaper to obtain more information from the same well, since the tools for gathering the measurements have been already been placed into the well. PAGLIA et al. 21 TABLE 3 The largest five designs obtained over 10 restarts of the algorithm in the petroleum drilling risk example, listed in descending order Design (D)I(D) 5, 6, 8, 14 €9,115,100 6, 7, 8, 15 €9,115,100 5, 7, 8,11 €9,115,100 3, 5, 6, 7, 8, 19, 20 €9,114,900 3, 5, 6, 7, 8, 9, 10, 12 €9,114,600 Based on seismic data from the region along with geological simulations of pressure build-up and release (Borge, 2000; Lothe, 2004; Paglia et al., 2019) we fit an initial multivariate Gaussian model for the pore pressure at the well location and at neighboring wells. The pore pressure measurements yDin neighboring wells will then be indicative of the pore pressure in the target well via correlations. These modeling assumptions simplify the computation of the conditional expectation of pore pressure in the well of interest, given data obtained in the other wells. The main challenge of the VOI evaluation in the current setting is then to compute the expectation of the nonlinear value function and to find its integral over all possible data yD. These expectations required for the PV (expression (1)) and the PoV (expression (3)) are here calculated with numerical approximations of the integrals. This entails rather time demanding computations for the value function 𝜈over discretized levels of pore pressure, where the LB is then finally obtained with a spline interpolation. We let Δxjand ΔyDdenote the distances between two consecutive discretized levels of xand yD, respectively. The VOI approximation can then be computed from PV =max a∈{E(𝜈(x,a))} ≈max a∈{∑ j 𝜈(xj,a)p(xj)Δxj}, PoV(D)=EyD[max a∈{E(𝜈(x,a)|yD)}]≈∑ yd max a∈{∑ j 𝜈(xj,a)p(xj|yD)Δxj}p(yD)ΔyD, VOI(D)=PoV(D)−PV.(18) Figure 9 shows the performance of Algorithm 1 for this application. As in the first example we do not know the design with the largest information gain value, but we observe a convergence of the method toward larger information gain as more batches are evaluated. Table 3 shows the five largest values of I(D)in descending order, together with the associated design. Once one starts gathering data at a specific well, acquiring more data at other depth is relatively cheap. This implies that the cost effective designs suggest to explore more then one depth for a single well. We notice from the table that there are designs with the same large value of I(D), but with different sites configurations. All the best designs have four sites and have some common sites. We believe that these sites are highly informative. In this case the sequential selection method gives a solution of €9,114,900 with 210 iterations, which is rather good but not as high as the best one achieved with the Bayesian optimization method. The exchange algorithm obtains €9,114,600 with 800 iterations. That small monetary 22 PAGLIA et al. 0 5 10 15 Batch Iteration 9.1136 9.1138 9.114 9.1142 9.1144 9.1146 9.1148 9.115 9.1152 VOI(D)-C(D) 106 FIGURE 9 Performance of the Bayesian optimization algorithm for a drilling application. Over the 10 restarts, the algorithm converges to a large value of VOI-C [Colour figure can be viewed at wileyonlinelibrary.com] differences for the three methods is due to the symmetry assumption of cost being the same for sites at the same depth in different wells. 6CLOSING REMARKS The main purpose of this study is to develop an algorithm that can assist a decision-maker in choosing a spatial design for collecting information. The methodology has its applications in earth sciences, where data are often distributed over a spatial domain. We illustrated the approach by presenting one example from the forestry and one from petroleum. We have adopted the Hausdorff distance to model dissimilarities between designs, and incorporate this in the kernel of a GP surrogate model. We demonstrated its use in examples. We believe that, depending on the field where the methodology is applied, other metrics can work as well. The Jaccard distance and the metrics introduced by Fujita (2013), can be valid alternatives to the Hausdorff distance. For instance, the Jaccard distance could work in situations where we are not too interested in spatial distance between designs, but in a more machine learning oriented context where one must select the appropriate number of sets for training, and because data may come from different sources it is important to guide the active learning wisely to cover appropriate features (Settles, 2012). The developed methodology could also be applied in subset selection problems such as the choice of individuals in epidemiological follow-up studies (Reinikainen et al., 2016) or in genotyping (Karvanen et al., 2009). It would be very useful to gain insight in PAGLIA et al. 23 theoretical properties of the algorithm, at least in special situations, or potentially to connect the mix of design combinations in the algorithm to some useful theoretical properties. The value function could be extended to include a parametric form. For instance, a possibility is to focus on learning the regression parameters 𝜷in the simulation study, or other parameters. This will go beyond only expectedeconomic outputs, and rather contain additional (hybrid) terms in the prior and PoV related with the standard deviation in profits. The common space-covering geometric designs are sometimes considered robust because they do not use target specific prediction purposes. It is possible to blend our approach with other design constructions to balance multiple criteria. It is possible to extend the study considering more challenging probability distributions for data, where the computation of VOI becomes more difficult. The combination of the VOI analysis with Bayesian optimization techniques gives us an efficient way to find satisfactory data gathering scheme. With the Bayesian optimization we move the problem of evaluating VOI to that of computing EI, which requires less computational effort. The total number of evaluation of VOI is considerably reduced. In situations where computing the information gain is computationally demanding or the number of alternatives to explore is too large, the developed methodology reduces the time of computation. ORCID Jacopo Paglia https://orcid.org/0000-0003-0314-9284 Juha Karvanen https://orcid.org/0000-0001-5530-769X REFERENCES Abbas, A. E., & Howard, R. A. (2015). Foundations of decision analysis. Pearson Higher Education. Bachoc, F., Suvorikova, A., Ginsbourger, D., Loubes, J.-M., & Spokoiny, V. (2020). Gaussian processes with multidimensional distribution inputs via optimal transport and Hilbertian embedding. Electronic Journal of Statistics, 14, 2742–2772. Bhattacharjya, D., Eidsvik, J., & Mukerji, T. (2013). The value of information in portfolio problems with dependent projects. Decision Analysis,10, 341–351. Binois, M., Huang, J., Gramacy, R. B., & Ludkovski, M. (2019). Replication or exploration? Sequential design for stochastic simulation experiments. Technometrics,61, 7–23. Borg, I., Groenen, P. J., & Mair, P. (2018). Applied multidimensional scaling and unfolding. Springer. Borge, H. (2000) Fault controlled pressure modelling in sedimentary basins. [Ph.D. thesis]. Norwegian University of Science and Technology. Bouneffouf, D. (2016). Exponentiated gradient exploration for active learning. Computers,5,1. Brochu, E., Cora, V. M., & De Freitas, N. (2010) A tutorial on Bayesian optimization of expensive cost functions, with application to active user modeling and hierarchical reinforcement learning. arXiv preprint arXiv:1012.2599. Chiles, J.-P., & Delfiner, P. (2012). Geostatistics: Modeling spatial uncertainty. Wiley. Diggle, P., & Lophaven, S. (2006). Bayesian geostatistical design. Scandinavian Journal of Statistics,33, 53–64. Dobbie, M. J., Henderson, B. L., & Stevens, D. L., Jr. (2008). Sparse sampling: Spatial design for monitoring stream networks. Statistics Surveys,2, 113–153. Drovandi, C. C., Holmes, C., McGree, J. M., Mengersen, K., Richardson, S., & Ryan, E. G. (2017). Principles of experimental design for big data analysis. Statistical Science: A Review Journal of the Institute of Mathematical Statistics,32, 385. Drovandi, C. C., McGree, J. M., & Pettitt, A. N. (2013). Sequential Monte Carlo for Bayesian sequentially designed experiments for discrete data. Computational Statistics & Data Analysis,57, 320–335. Dubuisson, M.-P., & Jain, A. K. (1994) A modified Hausdorff distance for object matching. Proceedings of 12th International Conference on Pattern Recognition (Vol. 1, pp. 566–568). IEEE. 24 PAGLIA et al. Eidsvik, J., Martinelli, G., & Bhattacharjya, D. (2018). Sequential information gathering schemes for spatial risk and decision analysis applications. Stochastic Environmental Research and Risk Assessment,32, 1163–1177. Eidsvik, J., Mukerji, T., & Bhattacharjya, D. (2015). Value of information in the earth sciences: Integrating spatial modeling and decision analysis. Cambridge University Press. Evangelou, E., & Eidsvik, J. (2017). The value of information for correlated GLMs. Journal of Statistical Planning and Inference,180, 30–48. Eyvindson, K., Hakanen, J., Mönkkönen, M., Juutinen, A., & Karvanen, J. (2019). Value of information in multiple criteria decision making: An application to forest conservation. Stochastic Environmental Research and Risk Assessment,33, 2007–2018. https://doi.org/10.1007/s00477-019-01745-4 Eyvindson, K. J., Petty, A. D., & Kangas, A. S. (2017). Determining the appropriate timing of the next forest inventory: Incorporating forest owner risk preferences and the uncertainty of forest data quality. Annals of Forest Science,74,2. Frazier, P. I. (2018) A tutorial on Bayesian optimization. arXiv preprint arXiv:1807.02811. Fujita, O. (2013). Metrics based on average distance between sets. Japan Journal of Industrial and Applied Mathematics,30, 1–19. García-Ródenas, R., García-García, J. C., López-Fidalgo, J., Ángel Martín-Baos, J., & Wong, W. K. (2020). A comparison of general-purpose optimization algorithms for finding optimal approximate experimental designs. Computational Statistics & Data Analysis,144, 106844. Ginsbourger, D., Baccou, J., Chevalier, C., & Perales, F. (2016). Design of computer experiments using competing distances between set-valued inputs.InmODa 11-advances in model-oriented design and analysis (pp. 123–131). Springer. Goldberg, D. E. (1989). Genetic algorithms in search, optimization and machine learning (1st ed.). Addison-Wesley Longman Publishing Co., Inc. Grafström, A., Lundström, N. L., & Schelin, L. (2012). Spatially balanced sampling through the pivotal method. Biometrics,68, 514–520. Gramacy, R. B. (2020). Surrogates: Gaussian process modeling, design, and optimization for the applied sciences.CRC Press. Harman, R., Filová, L., & Richtárik, P. (2020). A randomized exchange algorithm for computing optimal approximate designs of experiments. Journal of the American Statistical Association,115, 348–361. Huan, X., & Marzouk, Y. M. (2013). Simulation-based optimal Bayesian experimental design for nonlinear systems. Journal of Computational Physics,232, 288–317. Huttenlocher, D. P., Rucklidge, W. J., & Klanderman, G. A. (1992) Comparing images using the Hausdorff distance under translation. Proceedings 1992 IEEE Computer Society Conference on Computer Vision and Pattern Recognition (pp. 654–656). IEEE. Kangas, A., Kangas, J., & Kurttila, M. (2008). Decision support for forest management (Vol. 16). Springer. Karvanen, J., Kulathinal, S., & Gasbarra, D. (2009). Optimal designs to select individuals for genotyping conditional on observed binary or survival outcomes and non-genetic covariates. Computational Statistics & Data Analysis, 53, 1782–1793. Krause, A., & Golovin, D. (2014). Submodular function maximization. Tractability,3, 71–104. Levandowsky, M., & Winter, D. (1971). Distance between sets. Nature,234, 34–35. Lothe, A. E. (2004) Simulations of hydraulic fracturing and leakage in sedimentary basins. [Ph.D. thesis]. University of Bergen. Lothe, A. E., Cerasi, P., & Aghito, M. (2019). Digitized uncertainty handling of pore pressure and mud-weight window ahead of bit: North Sea example. SPE Journal,25, 24. https://doi.org/10.2118/189665-PA Madan, V., Singh, M., Tantipongpipat, U., & Xie, W. (2019) Combinatorial algorithms for optimal design. Proceedings of the Conference on Learning Theory (pp. 2210–2258). PMLR. Min, D., Zhilin, L., & Xiaoyong, C. (2007). Extended Hausdorff distance for spatial objects in GIS. International Journal of Geographical Information Science,21, 459–475. Mitchell, T. J. (1974). An algorithm for the construction of "d-optimal" experimental designs. Technometrics,16, 203–210. Moos, D., Peska, P., Ward, C., and Brehm, A. (2004) Quantitative risk assessment applied to pre-drill pore pressure, sealing potential, and mud window predictions from seismic data. Proceedings of the 6th North America Rock Mechanics Symposium (NARMS) Gulf Rocks 2004. American Rock Mechanics Association.