Bayesian Optimization with High-Dimensional Outputs
Authors: Wesley Maddox
Maximilian Balandat
Andrew Gordon Wilson
Eytan Bakshy
Venue: NeurIPS 2021 Poster (2021)
Keywords: multi-task Gaussian processes
Bayesian optimization
Gaussian processes

Abstract: Bayesian optimization is a sample-efficient black-box optimization procedure that is typically applied to a small number of independent objectives. However, in practice we often wish to optimize objectives defined over many correlated outcomes (or “tasks”). For example, scientists may want to optimize the coverage of a cell tower network across a dense grid of locations. Similarly, engineers may seek to balance the performance of a robot across dozens of different environments via constrained or robust optimization. However, the Gaussian Process (GP) models typically used as probabilistic surrogates for multi-task Bayesian optimization scale poorly with the number of outcomes, greatly limiting applicability. We devise an efficient technique for exact multi-task GP sampling that combines exploiting Kronecker structure in the covariance matrices with Matheron’s identity, allowing us to perform Bayesian optimization using exact multi-task GP models with tens of thousands of correlated outputs. In doing so, we achieve substantial improvements in sample efficiency compared to existing approaches that model solely the outcome metrics. We demonstrate how this unlocks a new class of applications for Bayesian optimization across a range of tasks in science and engineering, including optimizing interference patterns of an optical interferometer with 65,000 outputs. 

Bayesian Optimization with High-Dimensional Outputs
Wesley Maddox New York University wjm363@nyu.c edu
Maximilian Balandat
Andrew Gordon Wilson New York University andrewgw@cims nyu.
Eytan Bakshy Facebook eytanefb.com
Facebook ballndattffb.com
Abstract
Bayesian Optimization is sample-efficient black-box optimization procedure that is typically applied to problems with small number of independent objectives. However, in practice we often wish to optimize objectives defined over many correlated outcomes (or "tasks"). For example, network operators may want optimize the coverage of a cell tower network across a dense grid of locations. Similarly, engineers may seek to balance the performance of robot across dozens of different environments via constrained or robust optimization. However, the Gaussian Process (GP) models typically used as probabilistic surrogates for multi- task Bayesian Optimization scale poorly with the number of outcomes, which greatly limitis their applicability. We devise efficient technique for exact multi- task GP sampling that combines exploiting Kronecker structure in the covariance matrices with Matheron's 's identity, allowing to perform Bayesian Optimization using exact multi-task GP models with tens of thousands of correlated outputs. In doing so, we achieve substantial improvements sample efficiency compared to existing approaches that only model aggregate functions of the outcomes. We demonstrate how this unlocks new class of applications for Bayesian Optimization across range of tasks in science and engineering, including optimizing interference patterns of an optical interferometer with more than 65,000 outputs.
Introduction
Many problems science and engineering involve reasoning about multiple, correlated outputs. For example, cell towers broadcast signal across an area, and thus signal strength is spatially correlated. In randomized experiments, treatment effects on multiple outcomes are naturally correlated due to shared causal mechanisms. Without further knowledge of the internal mechanisms (i.e., in bolack-box" setting), Multi-task Gaussian processes (MTGPs) are natural model for these types of problems as they model the relationship between each output (or "task") while maintaining the gold standard predictive capability and uncertainty quantification of Gaussian processes (GPs). Many downstream analyses require more of the model than just prediction; they also involve sampling from the posterior distribution to estimate quantities interest. For instance, we may interested the performance of complex stock trading strategy that requires modeling different stock prices jointly, and want to characterize conditional value at risk (CVaR) [49), which generally requires Monte Carlo (MC) estimation strategies [1). Or, we want to use MTGPs in Bayesian Optimization (BO), method for sample-efficien optimization of black-box box functions. Many state the art BO approaches MC acquisition functions [60 4 which require sampling from the posterior distribution over
new candidate data points.
Drawing posterior samples from MTGPs means sampling over all tasks and all new data points, which typically scales multiplicatively in the number of tasks (f) and test data points (n), e.g. like On [52) 6 For problems with more than few tasks, posterior sampling thus quickly becomes
35th Conference on Neural Information Processing Systems (NcurIPS 2021).

---

intractable due to the size of the posterior covariance matrix. This is especially problematic in the case of many real-world problems that can have hundreds or thousands of correlated outputs that should be modelled jointly in order to achieve the best performance.
For instance, the cell tower signal maps in Fig- ure I each contain 2,500 outputs (pixels), In this problem. we aim jointly tune the down- tilt angle and transmission power of the anten- nas on each cell tower (locations shown in red) to optimize a global coverage quality metric, which is known function of power and interfer- ence at each location [21]. Since simulating the power and interference maps given a parameter- ization is computationally costly, traditionally one might apply BO to optimize the aggregate metric. At its core, this problem is composite BO problem 3 4. so we expect an approach that models the constituent outcomes at each pixel individually to achieve higher sample ef- ficiency. However, modelling each pixel using existing approaches used for BO is completely
Power
Interference
Figure 1: Map of radio signal power and interfer- ence for fixed locations of cell towers (red dots). Outcomes vary smoothly with respect to the tow- ers' down-tilt angle and transmission power. Our goal is to optimize statistics of these maps to maximize the overall signal coverage across an area while minimizing interference.
intractable in this setting, as would have to train and sample from MTGP with 5,000 tasks.
To remedy the poor computational scaling with the number of tasks, we exploit Matheron's rule for sampling from GP posterior distributions [13 59 We derive efficient method for MTGP sampling that exploits Kronecker structure inherent to the posterior covariance matrices, thereby reducing the complexity of sampling from the posterior to become effectively additive in the combination of tasks of data points, i.e. O(n + {s), as compared to (n Our implementation of Matheron's rule draws from the exact posterior distribution and does not require random features or inducing points, unlike decoupled sampling [59]. More specifically, our contributions are as follows:
We propose an exact sampling method for multi-task Gaussian processes that has additive time costs in the combination of tasks and data points, rather than multiplicative (Section 3).
We demonstrate empirically how large-scale sampling from MTGPs can aid in challenging multi-objeccive, constrained, and contextual Bayesian Optimization problems (Section 4
We introduce method for efficient posterior sampling for the High-Order Gaussian Process (HOGP) model [64]. allowing it to be used for Bayesian Optimization (Section (32)). This advance allows us to more efficiently perform BO on -igeeniinnsional outputs such as images including optimizing PDEs, optimizing placements of cell towers for cell coverage, and tuning the mirrors of an optical interferometer which optimizes over 65,000 tasks jointly (Section 4.3
The rest of the paper is organized as follows: First, in Section we review GPs, MTGPs, and sampling procedures from the posterior both GPs and MTGPs. In Section 3 we review Matheron's 's rule for sampling from GP posteriors and explain how to employ it for efficient sampling from MTGP models including the HOGP model. Section 4 we illustrate the utility of our method on a wide suite of problems ranging from constrained BO to the first demonstration of large scale composite BO with the HOGP, Please see Appendix A for discussion of the limitations and broader impacts of our work. Our code is fully integrated into Borch,,, see battps://bottte bo_ so_with hogp and https::///// httpl:::////tttt// ch, org/ butttrrrllla composite _ateeeeb___oot____eeeeeeeeeeote_ for tutorials.
org/tutorials composite_
Background
2.1 Bayesian Optimization
In Bayesian Optimization (BO), the goal is to minimize an preenivveto--evllatee black-box function, ie.. finding min f fz), by constructing a surrogate model to emulate that function. Gaussian processes (GPs) are often used as surrogates due to their flexibility and well-calibrated uncertainty estimates. BO optimizes acquisition function defined on the predictive distribution of the surrogate model to select the next point(s) to evaluate on the true function. These acquisition functions are often written as intractable integrals that are typically evaluated using Monte Carlo (MC) integration

---

660 3). MC acquisition functions rely on posterior samples from the surrogate model, which should support fast sampling capabilittess for efficient optimization 4. BO has been applied throughout machine learning, engineering, and the sciences, and many extensions to the setting described above exist. We focus on multi-task BO (MTBO), where (x) composed of several correlated tasks
55 15.
There are many sub-classes of MTBO problems: constrained BO uses surrogate models to optimize an objective subject to black-box constraints (29] 25] (30] contextual BO models a single function
26 multi-objective BO aims to explore 24 17 and composite BO considers the setting of a differentiable objective function defined on the outputs of a vector-valued black-box function [ 4.. In all of these problems, the setting is similar: several outputs are modelled by the surrogate, whether the output is a constraint, another objective, or separate context. As the outputs are often correlated, multi-task Gaussian processes, which model the relationships between the outputs in data-efficient manner, are a natural and common modeling choice.
that varies across different contexts or environments 38
a Pareto frontier across several objectives [36 37
23
2.2 Gaussian Processes
Single Output Gaussian Processes: We briefly review single output GPs, see Rasmussen and Williams [48] for a more detailed introduction. We assume that y = f(z) + E, f ~ GP(O, (z)x, and e ~ N(0 o²), where f is the noiseless latent function and are noisy observations of f with standard deviation a. Ka (2 x²) the kernel with hyperparameters 0 (we will drop the dependence on 0 for simplicity); we use KAB == (o (A, B) to refer to the evaluated kernel function on data points A and B, matrix which has size AA| BB) The predictive distributions over new data points, Tesa, is given by the conditional distribution of the Gaussian distribution. That is, p(f(Zee)) 0) N ( Eipp), where D = yX. is the training dataset of size n = |X] and
1 = KIXEX (Karin + +0²1)=yy.
(1)
Dj = Kreetue - Krrex (Kuain + o 0²2) Kx
xx zme
(2)
with Kiain := Kxx. For simplicity, we will drop the subscripts f in all future statements. Computing the predictive mean μ² and variance Σ* requires Onn²) time and O(n²) space when
using Cholesky decompositions for the linear solves (48). Sampling is usually performed by
f(xee)((Y = v) = n² (E*)1/22,
where 2 ~ (N,, Computing samples at Thest test points from the predictive distribution costs O(n² + sm + n2" nestt + neag), computed by adding up the cost of all of the matrix vector multipli- cations (MVMs) and matrix solves. For fixed Xlest we can incur the cubic terms only once by re-using Cholesky factorizations of Knain + 1 and E* for each sample.
To reduce the time complexity, we can replace all matrix solves with n steps of conjugate gradients (CG) and the Cholesky decompositions with rank r < n Lanczos decompositions (an approach called LOVE (45).. These change the major scaling from ne3 down to rn² and the overall
time complexity to Orn? + srMeea + TTN rm [45] 28) In general, r < n is used and is
usually accurate nearly numerical precision 28 We provide additional details in Appendix [
One straightforward way of modelling multiple outputs is to consider each output as an independent GP, modelled in batch together, either with shared or independent hyperparameters (48] 28 25). However, there are two major drawbacks to this approach: (i) the model is not able to model correlations between the outputs, and (ii) if there are many outputs then maintaining a separate model for each can result in high memory usage and slow inference times. To remedy these issues, Higdon et al. (33] propose the PCA-GP, using principal component analysis (PCA) to project the outputs to a low-dimensional subspace and then use batch GPs to model the lower-dimensional outputs.
Multi-Output Gaussian Processes:
We consider multi-task Gaussian processes (MTGPs) with the intrinsic co-regionalization model (ICM), which considers the relationship between responses as the product of the data and task features We focus on this model due to its popularity and simplicity, leaving a similar derivation of the linear model of co-regionalization to Appendix C. Given inputs and x belonging to tasks i
32 6
and j. respectively, the covariance under the ICM is k[ 1), [a', (J) =n(x, (() Given n
data points X with the n associated task indices I, the covariance is Hadamard product of the data
3

---

Table 1: Time complexities for posterior sampling in single-output, multi-task, and high-order Gaussian Process (HOGP) models. Time complexities shown in blue are our contributions that have not yet been considered by the literature. Standard sampling from MTGPs scales multiplicatively in the combination of the number of tasks, L, and the number of data points, n, while using Matheron's rule reduces the combination to effectively become additive in these components.
Model
Distributional (Standard) (Eq.
With Matheron's 's rule (Eq.
Singlc-Output O((nn + T )(a
Onn" + meese)
On + neea)
Multi-Task
O(1 + Thea) + (*)
HOGP
O(1re + (reaa) I1 d?)
O(1 + moaa 2izz df)
kernel and the task kernel, Kunin = Xxx @ KII, and the response y is still of size n. We term this implementation of multi-task GPs "Hadamard MTGPs" in our experiments. In general, there is no easily exploitable structure in this model.
If we observe each task at each data point (the so-called block design case), the covariance matrix becomes Kronecker structured, e.g. Kuain = Xxx @ Kr. where Xxx is the data covariance matrix and Kr is the / x / task covariance matrix between tasks [6), and we now have nl scalar responses. To simplify our exposition, we assume that K is full-rank (this is not required as we can use pseudo-inverses in place of inverses). Thus, the GP prior is vec(y) (O, Kxx Kr), where y is a matrix of size n x 1, and vec(y) is a vector of shape nt. The GP predictive distribution is given by p(f* [ D) N (μ" ))),, where
x¹ = =KKuai & (x) @ Kr +
Σ. -==(0 =KKmar & KT) - (Kreiix @ K @ Kr +
s & Kr).
(4)
The kernel matrix on the training data, Kxx @ Kr. is of size nt x nl, which under standard (Cholesky-baste based) approaches yields inference cubic and multiplicative in n and t, that is O(n
However, the Kronecker structure can be exploited to compute the posterior mean and variance in Ont((n + t) + n + t²), which is dominated by the individual cubic terms (50 54].
Sampling from the posterior distribution in 4 produces additional computational challenges as we must compute a root (e.g. Cholesky) decomposition of Σ", which naively costs (( 33 plus
an additional Cholesky decomposition of [xxx @ Kr + +0² which similarly costs ((())
time (6) Thus, the time complexity of drawing s samples is multiplicative in n and t, O((nt)
(Restt) + s(nn)²² + (Reea)))) )). Using CG and LOVE reduces the complexity; see Appendix B.4
High-Order Gaussian Processes: Recently, Zhe et al. (64] proposed the high-order Gaussian process (HOGP) model, a MTGP designed for matrix- and tensor-valued outputs. Given out- puts y E R4 Rd1x-xxd (c.g. a matrix or tensor), the covariance is the product of the data dimension and each output index (i and if respectively) k([z, i,, , [x', "i, :%) = ((),,')k();; r44;) ).. k(V,,, , 1f()) ) where i1, ".,.... i are the indices for the output tensor and U1,-.. Uk are latent parameters that are optimized along with the kernel hyper-parameters. Thus, Kr in the MTGP framework is replaced by chain of Kronecker products, that the GP prior is vec(y) ~ (0, K Xxx& @ K₂ 8 ... & K2). Exploiting the Kronecker structure, computation of posterior mean, posterior variance and hyper- parammta learning takes On + Dk time as dj = n. Zhe et al. [64] demonstrate that the HOGP outperforms PCA- 331. among other high-dimensional output GP methods, with respect to prediction accuracy. However, their experi- ments measure only the error of the predictive mean, rather than quantities that use the predictive variance (such as the negative log likelihood or calibration), and they do not provide way sample from the posterior distribution.
Diz de + IT d,)
2.3 Matheron's 's Rule for Single- Task Gaussian Processes
Matheron's 's rule provides an alternative way of sampling from GP posterior distributions: rather than decomposing the GP predictive covariance for sampling, one can jointly sample from the prior distribution over both train and test points and then update the training samples using the observed data. Matheron's rule is well known in the geostatistics literature where it is used for "prediction by conditional simulation"
[4] [3]. Wilson et al. [59] revitalized interest in Matheron's 's rule within the
4

---

machine learning community by developing a decoupled sampling approach that exploits Matheron's s rule to use both random Fourier features (RFFs) and inducing point approaches in the context of sampling from the posterior in single task GPs. Wilson et al. [61]] extended decoupled sampling to use combinations of RFFs, inducing points, and iterative methods. They applied these approaches to approximate Thompson sampling, simulations of dynamical systems, and deep GPs.
Matheron's rule states that if two random variables, X and Y, are jointly Gaussian, then
xI(Y= v) - X Cov(X, ()CvXy ())
where Cov(a, b) is the covariance of a and b and denotes equality in distribution 20, 59 The identity can easily be shown by computing the mean and variance of both sides. Following Wilson ct al. [59], we can use this identity to draw posterior samples
MIY + =) 2 KaEx (Kxx + 10 Y c),,
(5)
where e ~ (00 I2 f is the random variable and fo (Y + 3) is the conditional random variable we are interested in drawing samples from. To implement this procedure, we first draw samples from the joint prior (of size + Tles):
(f.Y)~ ~ N (0.
(o.(kx Kxx K X
K KmniX
)):
(6)
K KrnssTiu
For shorthand, we denote the joint covariance matrix in (6) by Kyoint- We next sample ē ~ NNO (1) o (y-Y Y--
and compute: f= f n X n so that f is our realization of the random variable f*I(Y++ +E= y). The total time requirements are then ((n + mes)) + n²),, which is slightly slower than O + n²)). Thus, sampling from the
K K x xXxx- xx +
E). The solve is against matrix of size
single-task GP posterior using (5) is slower than the standard method based on (1) and (2).
3 Matheron's Rule for Multi-Tast Sampling
While the time complexity of Matheron-based posterior sampling is inferior for single-task GPs, the opposite holds true for multi-task GPs, provided that one exploits the special structure in the covariance matrices. In the following, we demonstrate the advantages in using Matheron-based MTGP sampling in Section 3.T extend it to the HOGP (44) in Section 3.2 and identify further advantages of it for Bayesian Optimization in Section 3.3 This approach allows for further pre-computations than distributional sampling and that it maintains the same convergence results.
3.1 Extending Matheron's 's Rule for Multi-task GPs
The primary bottleneck that we wish to avoid is the multiplicative scaling in n (or nsa)) and f. Ideally, we hope to achieve posterior sampling that is additive in n and I, similar to how Kronecker matrix vector multiplication is additive in its components. Unlike in the single task case, Matheron's rule can substantially reduce the time and memory complexity sampling for MTGPs. For brevity we limit our presentation here to the core ideas, please see Appendix C for further discussion of the implementation, as well as Appendix B.] for a description of the Kronecker identities we exploit.
The covariance Kjuin in (6) is structured as the Kronecker product of the joint test-train covariance matrix and the inter-task covariance; that is, Kyemmt K[x, ,X,)) @ Kr, where K K(x, (XX,))(XX,aa (N.a appends Neat rows and columns to the joint training data covariance matrix Kxx. To sample from the prior distribution, we need to compute a root decomposition of Kyoint-
We assume that we have pre-computed R.T |xx with either a Cholesky decompostion (n33) time) or a Lanezos decomposition (0(rn?) time). Then, we update R to compute RR x ((zue.)) x). Chevalier et al. [2] used similar strategy to update samples in the context of single task kriging. Following Jiang et al. 34 Prop. 2], computing it from R costs O( + rmlo)) time, dominated by the it terms for small Teest Using Cholesky decompo- sition, this is (Ttss?)) time following the same procedure. We then have
Kjoin K (xnea)) X) 8 Kr (R L)
(7)
5

---

out of 22- Hatamand
2.5
Tire t1
-2-
Time
12
12
21
320
130
200
200
400
930
Narter of Toxt Forts
Mambbr Parts
STot Teal Pora
(a) Test points, GPU
(b) Tasks, GPU
(c) Test points, CPU
(d) Tasks, CPU
Figure 2: Timings for distributional sampling with Hadamard and Kronecker MTGPs as well as Matheron's rule sampling for MTGPs as the number of test points vary for fixed tasks (a,c) and the number of tasks vary for fixed test points (b,d) on single Tesla V 100 GPU (a,b) and single CPU (c,d). The multiplicative scaling of the number of tasks and data points creates significant timing and memory overhead, causing the Kronecker and Hadamard implementations memory very quickly for all but the smallest numbers tasks and data points, whereas sampling using Matheron's 's rule is efficient even in the many- task large-data regime. The plots show mean and two standard errors over 10 trials on the GPU, and 6 trials on the CPU.
where LL = Kr. To use the root to sample from the joint random variables, (J, Y), we need to compute only (f. (R @ L )z, where ~ (N, Ane): this computation is a Kronecker matrix vector multiplication (MVM), which costs just (nnt (n + t)) time. Thus, the overall cost of sampling to compute the joint random variables is On + + nt(n + f) + neett?²) time, reduced Orn? + rt? + rnf + r2t TRess))) if using Lanezos decompositions throughout. only needs to be computed once, so samples new test points only require re-computing i and further MVMs.
Computing the solve w: = xXxx @ Kr+ o² ) (y- Y - <) takes O(n + 13 + nt(n + f)) time using sigen-decompositions and Kronecker MVMs. The cost of the iigen-eccomposition dominates the Kronecker MVM cost so the major cost in this is O(n +t²). Finally, there is one more Kronecker MVM (K KT) w which takes (n + (haat)) time. We then only need reshape f to be of the correct size.
Therefore, the total time complexity of using Matheron's 's rule to sample from the MTGP posterior is Onn + for small Thess as the cubic terms dominate due to the Cholesky and eigen deceenddcoosiiinnnss We show the improvements from using Matheron's rule in Table where the dominating terms now additive in the combination of the tasks and the number of data points, rather than multiplicative. Memory complexities, which also much reduced, we emphasize that sampling in this manner is exact to floating point precision as the solves we use all exact by virture of being computed with eigen-decompositions.
provided in Table 2 in Appendix C Finally,
3.2 Extension to HOGP
The HOGP model can be seen as a special case of the general procedure for sampling MTGPs. We replace Kr with kernels across each tensor dimension of the response, that K K, Therefore, the time complexities for sampling go from cubic dependence, n II ds to
(n II_,,,,, d,) + (n' + > d?) dependence. The latter will usually be dominated by the addi- tive terms for k S 5, as generally
much larger than the tensor sizes, hence their product will generally be less than n². Prior to this work, sampling from the posterior was infeasible for all but the smallest HOGP models due to the time complexity. By using Matheron's rule, we can sample from HOGP posterior distributions over large tensors, unlocking the model for much wider range of applications such as BO with composite objectives computed on high-dimensional simulator outputs.
3.3 Usage in Bayesian Optimization
The primary usage of efficient multi-task posterior sampling that we consider is that of Bayesian optimization. Here, we want to use and optimize Monte Carlo acquisition functions that require many samples from the posterior of the surrogate model.

---

At a high level, Bayesian optimization loops look like the following procedure that we repeat after initializing several random points, First, we fit MTGPs by maximizing the log marginal likelihood (see Appendix [B.3) Then, we draw many posterior samples from the MTGP in the context of computing MC acquisition functions, e.g. (A.5). We use gradient-based optimization to find the points æ which optimize the acquisition x (#;; V) After choosing these points, eeandd we then query the function,
returning Meand and updating our data with these points. Finally, we continue back to the top of the loop, and re-train the MTGP.
The convergence guarantecs for Sample Average Approximation (SAA) of MC acquisition functions from Balandat et al. [4] still apply if posterior samples drawn via Matheron's rule. Consider acquisition functions of the form ax y) = E [((((x)) Y== = y] with (>=a((([[]) a()([)))) where S i..d. samples drawn from the model posterior at E Re using Matheron's 's rule. Then, under sufficient regularity conditions, the optimizer arg max, & (z; y) converges the true optimizer arg max, (zz; y) almost surely as N + oc. For a more formal statement and proof of this result see Appendix
Convergence Results:
h: 1 IR, and their MC approximation d (z;; y) N
Experiments
We first demonstrate the computational efficiencies gained by posterior sampling using Matheron's rule. While this contribution is much more broadly useful, here we focus on the benefits it provides for BO. Namely, show that accounting for correlations between outcomes in the model improves performance in multi-objective and large-scale constrained optimization problems. Finally, we perform composite BO with tens of thousands of tasks using HOGPs with Matheron- based sampling. Additional experiments on contextual policy optimization are given in Appendix D.2] All plotted lines represent the mean over repeated trials with shading represent two standard deviations of the mean across trials.
4.1 Drawing Samples from Multi-Task Models
To demonstrate the performance gains of sam- pling using Matheron's 'ss rule, we first vary the number of test points and tasks for fixed num- ber of samples and training points. Following Feng et al. (26). we consider a multi-task ver- sion of the Hartmann-6 function, generated by splitting the response surface into tasks using slices of the last dimension. In Figures and 2c we vary the number of test points for 5 and 50 tasks, demonstrating that Matheron's 's rule pling is faster and more memory efficient than distributional sampling with either Kronecker or Hadamard MTGPs. In Figures 2h and 2d] we vary the number of tasks for tixed test points, again finding that distributional sampling is only competitive for 5 tasks on the GPU. See Ap- pendix D for sampling with LOVE predictive covariances.
Random
CEVV-MTGP HHHi-Batch
aParego-PCA (EHVLPCA
@Parago-MTGP
oParego Batch
8900
6300
6000
2300
100
206
300
400
108
200
200
(a) C2DTLZZ, 2
(b) OSY, 10
Figure Constrained multi-objective Bayesian Optimization tasks. MTGPs outperform batch models in both the (a) small batch (q 2) and the (b) large batch (q 10) setting. The latter was previously computationally infeasible for MTGPs.
4.2 Multi-Objective Bayesian Optimization
We next consider constrained multi-objective BO, where the goal is to find the Pareto Frontier, i.e., the set of objective values for which no objective can be improved without deteriorating another while satisfying all constraints. To measure the quality of the Pareto Frontier, we compute the hypervolume (HV) of the non-dominated objectives [62] The optimization problem is made more difficult by the presence of black-box constraints (and hence additional outcomes) that must be modelled. We use MC batch versions of the ParEGO and EHVI acquisition functions (qParEGO and
7

---

1.0
330
250
320
06
230
+10
120
1.4
Baldd OP
00
240
Feasibility, Lunar Lander, 50
(b) Feasibility, Lunar
(c) Merit, Lunar Lander, m
(d) Vehicle mass, MOPTA08, 68
Lander, m 100
50
Figure 4: Scalable constrained Bayesian Optimization on Lunar Lander = 50, 100 (a-c) and on MOPTA0S (d), On all three problems, a multi-task GP provides better solutions with better data efficiency. The batch GP reaches feasibility on lunar lander with 50 constraints (a) and a competitive solution (c) but requires more trials, while on 100 constraints, it simply runs out memory while the MTGP succeeds. On MOPTA08, the MTGP reaches better solution faster (d).
qEHVI) To generate new candidate points, qParEGO maximizes expected improvement using random Chebyshev scalarization of the objectives, while qEHVI maximizes expected hypervolume improvement. As far as we are aware, Shah and Ghahramani [52] are the only authors to investigate the use of MTGPs in combination with multi-objective optimization, but only consider 2-4 tasks in the sequential setting (generating one point at a time). Here, we use full rank inter-task covariances with LKJ priors 40] which find to work well even in the low data regime. Our Matheron-based sampling scales to large batches and tasks, and more sample-efficient on both tasks.
We compare Matheron sampled MTGPs to batch independent MTGPs on the C2DTLZ2 [ (2 objectives, constraint for total of modelled tasks) and OSY 43 objectives. constraints for total of 8 modelled tasks) test problems. On OSY, we also compare to PCA GPs [33 due to the larger number outputs. Following Daulton et [17] we use both qParEGO and qEHVI with q = for C2DTL772 and optimize for 200 iterations. In Figure 3a we see that the MTGPs outperform their batch counterparts by coming closer to the known optimal HV. OSY, Figure (3b] we plot the maximum HV achieved for each method, using a batch size of 10, optimizing for 30 iterations, where again we see that the MTGPs significantly outperform their batch counterparts as well as the PCA GPs, which stagnate quickly.
4.3 Scalable Constrained Bayesian Optimization
We next extend scalable constrained Bayesian Optimization [SCBO, 25 state of the algorithm for constrained BO in high- digh-dimensional problems, to use MTGPs instead of independent GPs. In constrained BO, the goal minimize the objective, f. subject to black box constraints, c; e.g..
arg min f(x) s.t. (x(x) S O, E ... ....
SCBO is method based on trust regions and uses batched independent GPs model the outcome and transformed constraints. We compare to their results on their two largest problems - the 12- dimensional lunar lander and the 124-dimensional MOPTA08 problem, using the same benchmark To extend their approach, we replace the independent GPs with a single MTGP model with full rank ICM kernel over objectives and constraints.
Lunar Lander: first consider the lunar lander problem with both 50 and 100 constraints from the OpenAl Gym g. Following Eriksson and Poloczek (25] we initialize with 150 data points and use TuRBO with Thompson sampling with batches of q = 20 for a total 1000 function evaluations and repeat over 30 trials. Using multi-task GP reduces the number of iterations to achieve at least a 50% chance of feasibility by about 75 steps for the 50 constraint problem (Figure single GPU and never achieves feasibility, as indicated in Figure 4b) In Figure we show the best merit (f()) I (0) achieved where the MTGPs are able to achieve feasibility in fewer samples, but do not reach significantly higher reward. Wall clock times for the m 50 constraint problem, a table of the steps to achieve feasibility, and comparison to PCA-GPs
4a). On the 100 constraint problem, the batch GP runs out of memory after 350 steps
33 are given in Appendix [
MOPTA08: We next compare to batch GPs on the MOPTA08 benchmark problem (35) which has 68 constraints that measure the fcasibility of a vehicle's 's design, while each dimension involves gauges,

---

EI-HOGP-CF + GP
EI-HOGP-CF
Random
EI-CF
Fundin Eniksttors
110
15
29
10
(a) Pollutant concentration
(b) PDE control (c) Cell-tower coverage
(d) Optics
(t= (x=3x4,d=0) 4)
(t 2 64 x 64, d 4) (t 2 50 x 50, 30) (t = 16 x 64 x 64, 4)
Figure 5: Performance of BO with and without HOGP- based composite objectives (EI-HOGP-CF) on four scientific problems. EI-HOGP-CE outperforms a standard BO model directly on the metric itself (EI) and a random baseline (Random). Composite BO with the HOGP also outperforms composite BO with independent GPs (EI-CF). Our smooth latents for the HOGP (EI-HOGP + GP) typically outperform the HOGP itself. EI-CF is only feasible on the smallest problem.
materials, and shapes. We use Thompson sampling to acquire points with a batch size of q = 10, 130 initial points, and optimize for 2000 iterations repeating over 9 trials. The results are shown in Figure 4d where we again observe that SCBO with MTGPs significantly improves both the time to feasibility and the best overall objective found. Using MTGPs would have been computationally infeasible in this setting without our efficient posterior sampling approach.
4.4 Composite Bayesian Optimization with the HOGP
Finally, we push well beyond current scalability limits by extending BO to deal with many thousands of tasks, enabling us to perform sample-efficient optimization in the space of images. Composite BO a form of BO where the objective is of the form max, (h(x)), where h a multi-output black-box function modelled by a surrogate and g is cheap to evaluate and differentiable. Decomposing single objective into constituent objectives in this way can provide substantial improvements in sample complexity [3]. Balandat et al. 44 gave convergence guarantees for optimizing general sampled composite acquisition function objectives that extend to MTGP models under the same regularity conditions. However, both works experimentally evaluate only batch independent multi-output GPs.
We compare to three different baselines: random points (Random), expected improvement on the metric (EI), and batch GPs optimizing EI in the composite setting (EI-CF). We consider the HOGP (64 with Matheron's 's rule sampling (HOGP-CH CCF) as well as an extension of HOGPs with prior over the latent parameters that encourages smoothly varying latent dimensions (HOGP-CF + GP): see Appendix C.2 for details. More detailed descriptions of the problems are provided in Appendix
Chemical Pollutants:
Following Astudillo and Frazier (3), we start with a simple spatial problem observed on 3 4 grid originally defined in
in which environmental pollutant concentrations Bliznyuk et al. [5). The goal is to optimize a set of four parameters to achieve the true observed value by minimizing the mean squared error of the output grid to the output grid of the true parameters.
The results, over 50 trials, are shown in Figure Sa where we find that the HOGP models with these few tasks outperform both independent batch GPs (but slightly) and BO on the metric itself.
Optimizing PDEs: As a larger experimental problem, we consider optimizing the two diffusivity and two rate parameters of spatial Brusselator PDE (solved in py-pde 655 to minimize the weighted variance of the PDE output as an example of control of a dynamical system. Here, we solve the PDE on 64 64 grid, producing output solutions of size 2 x 64 64. Results over 20 trials are shown in Figure 5b where the HOGP models outperform EI fit on the metric and the random baseline.
Following Dreifuerst et al. (21). we optimize the simulated "coverage map" resulting from the transmission power and down-tilt settings of 15 cell towers (for a total of 30 parameters) based on scalarized quality metric combining signal power and inference at each location so as to maximize total coverage, while minimizing total interference. To reduce model complexity, we down-sample the simulator output to 50 x 50, initializing the optimization with 20 points. Figure (5 presents the results over 20 trials, where the HOGP models with composite objective El outperform EI, indicating that modeling the full high-dimensional output valuable.
Cell-Tower Coverage:
Optical Interferometer:
Finally, we consider the unnng of an optical interferometer by the alignment of two mirrors as in Sorokin et al. [53]. Here, the problem is to optimize the mirror coordinates

---

to align the interferometer so that it reflects light without interference. There is a sequence of 16 different interference patterns and the simulation outputs are 64 x 64 images (a tensor of shape 16 x 64 x 64). Thus, we jointly model 65,536 output dimensions. Scaling composite BO to problem of this size would impossible without the step change in scalability our method provides. Results are shown in Figure 5c over 20 trials, where we find that the HOGP-CF GP models outperform EI, with the HOGP + GP under-performing (perhaps due to high frequency variation the images).
ovee
Across all of our experiments, we consistently find that composite BO is considerably more sample efficient than BO on the metric itself, with significant computational improvements gained from using the HOGP as compared to batch GPs, which are infeasible much beyond batch sizes of 100. Furthermore, our structured prior approach for the latent parameters of the HOGP tends to outperform the random initialization strategy in the original work of Zhe et al. [64].
Conclusion
We demonstrated the utility of Matheron's rule for sampling the posterior in multi-task Gaussian processes. Combining Matheron's 's rule with scalable Kronecker algebra enables posterior sampling in One + (a) rather than the previous (nn)tS) time. This renders posterior sampling from high- order Gaussian processes [[4]] practical, for the first time unlocking Bayesian Optimization with composite objectives defined on dimensional outputs. This increase in computational efficiency dramatically reduces the time required to do multi-task Bayesian Optimization, and thus enables practitioners to achieve better automated Bayesian decision making. While we focus on the application to Bayesian Optimization in this work, our contribution is much broader and provides a step change in scalability all methods that in involve sampling from MTGP posteriors. We hope in the future to explore stronger inter-task covariance priors to make MTGP model fits even more sample efficient.
Acknowledgements
WJM. AGW are supported by Amazon Research Award, NSF I-DISRE 193471, NIH RO1 DA048764-01A1. NSF IIS-1910266, and NSF 1922658 NRT-HDD:FUTURE Foundations, Transla- tion, and Responsibility for Data Science. WJM was additionally supported by an NSF Graduate Research Fellowship under Grant No. DGE- 183932,, and performed part of this work during an internship at Facebook. We would like to thank Paul Varkey for providing code and David Eriksson, Qing Feng, and Polina Kirichenko for helpful discussions.
References
[1] Álvarez, M. A., Rosasco, L.,, Lawrence, D., et al. (2012). Kernels for vector valued functions: A review. Foundations and Trends& in Machine Learning, 195):195-266.
[2] Andrade-Pachece R., Rerolle, F., Lemoine, J., Hermandez, L.. Meité, A., Juziwelo, L., Bibaut, F., van der Laan, M. J., Arnold, B. F., and Sturrock, H. (2020). Finding hotspots: development of an adaptive spatial sampling approach. Scientific reports. 10(1):1-12.
[3] Astudillo, R. and Frazier, P. (2019). Bayesian Optimization of Composite Functions. In International Conference on Machine Learning, pages 354-363. PMLR. ISSN: 2640-3498.
[4] Balandat, M., Karrer, B., Jiang, D., Daulton, S., Letham, B., Wilson, A. G., and Bakshy, E. (2020). BoTorrh: A Framework for Efficient Monte Carlo Bayesian Optimization. In Advances in Neural Information Processing Systems, volume 33.
[5] Bliznyuk, N., Ruppert, D., Shoemaker, C., Regis, R., Wild, S., and Mugunthan, P. (2008). Bayesian calibration and uncertainty analysis for computationally expensive models using op- timization and radial basis function approximation. Journal of Computational and Graphical Statistics, 17(2):27-2244
[6] Bonilla, V., Chai, K., and Williams, C. (2007). Multi-task Gaussian Process Prediction. In Advances in Neural Information Processing Systems, volume 20, pages 153-160.
10

---

[7] Boustati, A., Damoulas, T., and Savage, R. S. (2019). Non-linear multitask learning with deep gaussian processes. arXiv preprint 1rXiv:1905.12400
[8] Brockman, G., Cheung, V., Pettersson, L., Schneider, J., Schulman, J., Tang, J., and Zaremba, W (2016). Openai gym.
[9] Bruinsma, W. Perim, E., Tebbutt, W., Hosking, S., Solin, A., and Turner, R. (2020). Scalable exact inference in multi-output gaussian processes. In International Conference on Machine Learning, pages 1190-1201. PMLR.
[10] Cakmak, S., Astudillo Marban, R., Frazier, P., and Zhou, E. (2020). Bayesian optimization of risk measures. In Advances in Neural Information Processing Systems, volume 33, pages 20130-20141. Curran Associates, Inc.
[11] Char, I., Chung, Y., Neiswanger, W., Kandasamy, K., Nelson, A. O., Boyer, M., Kolemen, E., and Schneider, J. (2019). Offline contextual bayesian optimization. In Advances in Neural Information Processing Systems, volume 32, pages 4627-4638.
[12] Chevalier, C.,, Emery, X., and Ginsbourger, D. (2015). Fast update of conditional simulation ensembles. Mathematical Geosciences, 47(7):71--79..
[13] Chiles, J.-P. and Delfiner, P. (2009). Geostatistics: modeling spatial uncertainty, volume 497.
John Wiley & Sons.
[14] Chilès, J.-P. and Lantuéjoul, C. (2005). Prediction by conditional simulation: models and algorithms. In Space, Structure and Randomness, pages 39-68. Springer.
[15] Chowdhury, S. R. and Gopalan, A. (2021). No-regret algorithms for multi-task bayesian optimization. In International Conference on Artificial Intelligence and Statistics, pages 1873- 1881. PMLR.
[16] Dai, Z., Álvarez, M. A., and Lawrence, N. D. (2017). Efficient modeling of latent information in supervised learning using gaussian processes. In Proceedings of the 31st International Conference on Neural Information Processing Systems, pages 5137-5145.
[17] Daulton, S., Balandat, M., and Bakshy, E. (2020). Differentiable Expected Hypervolume Improvement for Parallel Multi Sulti-Objective Bayesian Optimization. Advances in Neural Information Processing Systems, 33.
[18] de Fouquet, C. (1994). Reminders on the conditioning kriging. In Geostatistical simulations, pages 131-145 Springer.
[19] Deb, K. (2019). Constrained multi-objective evolutionary algorithm. In Evolutionary and Swarm Intelligence Algorithms,, pages 85-118. Springer.
[20] Doucet, A. (2010). A Note on Efficient Conditional Simulation of Gaussian Distributions. https://www. cs, ubc. caaccaa/aaaaa awwaudddoucetttee.e...
[21] Dreifuerst, R. M., Daulton, S., Qian, Y., Varkey, P., Balandat, M., Kasturia, S., Tomar, A., Yazdan, A., Ponnampalam, V., and Heath, R. W. (2020). Optimizing coverage and capacity in cellular networks using machine learning. arXiv preprint arXiv:2010.137100
[22] Emery. X. (2007). Conditioning simulations of gaussian random fields by ordinary kriging. Mathematical Geology, 39(6):607- -623.
[23] Emmerich, M. (2005). Single-and multi-objective evolutionary design optimization assisted by gaussian random field metamodels. dissertation, Universität Dortmund.
[24] Emmerich, M. Deutz, A. H., and Klinkenberg, J. W. (2011). Hypervolume-basee expected improvement: Monotonicity properties and exact computation. In 2011 IEEE Congress of Evolutionary Computation (CEC), pages 2147-2154. IEEE.
[25] Eriksson, and Poloc7ek, M. (2021). Scalable constrained bayesian optimization. In Interna- tional Conference Artificial Intelligence and Statistics. PMLR.

---

[26] Feng, Q.. Letham, B., Mao, H., and Bakshy, E. (2020). High-Dimension Contextual Policy Search with Unknown Context Rewards using Bayesian Optimization. Advances in Neural Information Processing Systems, 33.
[27] Feydy, J., Glaunès, J., Charlier, B., and Bronstein, M. (2020). Fast geometric learning with symbolic matrices. Advances in Neural Information Processing Systems, (3(4):
[28] Gardner, J., Pleiss, G., Weinberger, K. Q., Bindel, D., and Wilson, A. G. (2018). GPyTorch: Blackbox Matrix-Matrix Gaussian Process Inference with GPU Acceleration. In Advances in Neural Information Processing Systems, volume 31, pages 7576-7586.
[29] Gardner, J. R., Kusner, M. J., Xu, Z. E., Weinberger, K. Q. and Cunningham, J. P. (2014). Bayesian optimization with inequality constraints. In International Conference on Machine Learning, volume 32, pages 937-945. PMLR.
[30] Gelbart, M. A., Snoek, J., and Adams, R. P. (2014). Bayesian optimization with unknown constraints. In Uncertainly in Artificial Intelligence, volume 30, pages 250-259. AUAI Press.
[31] Golub, G. H. and Van Loan, C. F. (2013). Matrix Computations. The Johns Hopkins University Press, edition.
[32] Goovaerts, al. (1997). Geostatistics for natural resources evaluation. Oxford University Press on Demand.
[33] Higdon, D., Gattiker, J., Williams, B., and Rightley, M. (2008). Computer model calibration using high-dimensional output. Journal of the American Statistical Association, 10344)2) 583.
[34] Jiang, S., Jiang, D., Balandat, M., Karrer, B., Gardner, J., and Gamett, R. (2020). Efficient Nommyopic Bayesian Optimization via One-Shot Multi-Step Trees. In Advances in Neural Information Processing Systems, volume 33.
[35] Jones, D. R. (2008). Large-scale multi-disciplinary mass optimization in the auto industry. In MOPTA 2008 Conference (20 August 2008).
[36] Khan, N., Goldberg, D. E., and Pelikan, M. (2002). Multi-objective bayesian optimization algorithm. In Proceedings of the 4th Annual Conference on Genetic and Evolutionary Computation, pages 684-684.
[37] Knowles, J. (2006). Parego: A hybrid algorithm with on-line landscape approximation for expensive multiobjective optimization problems. IEEE Transactions on Evolutionary Computation, 101))500-66
[38] Krause, A. and Ong. C. S. (2011). Contextual gaussian process bandit optimization. In Advances in neural information processing systems, pages 2447-2455.
[39] Larocque, G., Dutilleul, P., Pelletier, B., and Fyles, J. W. (2006). Conditional gaussian co- simulation of regionalized components soil variation. Geoderma, 134(1-2):1--6.
[40] Lewandowski, D., Kurowicka, D.,, and Joe, H. (2009). Generating random correlation matrices based on vines and extended onion method. Journal of multivariate analysis, 1989922988920000
[41] Lorch, L., Kremer, H., Trouleau, W., Tsirtsis, S., Szanto, A., Schölkopf, B., and Gomez- Rodriguez, M. (2020). Quantifying the Effects of Contact Tracing, Testing, and Containment Measures in the Presence of Infection Hotspots. arXiv e-prints, page arXiv:2000.00641.
[42] Nguyen, T. V., Bonilla, E. V., et al. (2014). Collaborative multi-output gaussian processes. In Uncertainty in Artifical Intelligence, volume 30, pages 643-652. AUAI Press.
[43] Osyezka, A. and Kundu, S. (1995). A new method to solve generalized multicriteria optimization problems using the simple genetic algorithm. Structural optimization, 10(2):94-99.
[44] Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J., Chanan, G., Killeen, T., Lin, Z., Gimelshein, N., Antiga, I.,, et al. (2019). Pytorch: An imperative style, hegf-peeformance deep learning library. Advances in Neural Information Processing Systems, 32:8026-8037.
12

---

[45] Pleiss, G., Gardner, J. R., Weinberger, K. Q. and Wilson, A. G. (2018). Constant-Tine Predictive Distributions for Gaussian Processes. In Artificial Intelligence and Statistics, volume arXiv: 10::18000 883.06058 (es, stat].
[46] Pleiss, G., Jankowiak,, M., Eriksson, D., Damle, A., and Gardner, J. (2020). Fast Matrix Square Roots with Applications to Gaussian Processes and Bayesian Optimization. In Advances in Neural Information Processing Systems, volume 33.
[47] Rakitsch, B., Lippert, C., Borgwardt, K., and Stegle, O. (2013). It all in the noise: Efficient multi-task gaussian process inference with structured residuals. In Advances in Neural Information Processing Systems, volume 26, pages 4474. Curran Associates, Inc.
[48] Rasmussen, C. E. and Williams, C. K. I. (2008). Gaussian processes for machine learning. Adaptive computation and machine learning. MIT Press, Cambridge, Mass., 3. print edition.
[49] Rockafellar, R. T., Uryasev, S., al. (2000). Optimization of conditional value-at-risk, Journal of risk, 2:21-42.
[50] Rougier, J. (2008). Efficient emulators for multivariate deterministic functions. Journal of Compulational and Graphical Statistics, 17(4)8827-843.
[51] Saatchi, Y. (2011). Scalable inference for structured Gaussian process models. PhD thesis, University of Cambridge.
[52] Shah, A. and Ghahramani, Z. (2016). Pareto frontier learning with expensive correlated objectives. In Proceedings The 33rd International Conference on Machine Learning, volume 48, pages 1919-1927, New York, New York, USA. PMLR
[53] Sorokin, D., Ulanov, A., Sazhina, E., and Lvovsky, A. (2020). Interferobot: aligning an optical interferometer by a reinforcement learning agent. arXiv preprint arXiv:2006.02252.
[54] Stegle, O., Lippert, C., Mooij, J., Lawrence, N., and Borgwardt, K. (2011). Efficient inference matrix-variate gaussian models with iid observation noise. Advances in Neural Information Processing Systems, 24:1-9.
[55] Swersky, K., Snoek, J., and Adams, R. P. (2013). Multi-task Bayesian optimization. In Advances Neural Information Processing Systems, volume 26, pages 2004-2012, Red Hook, NY, USA. Curran Associates Inc.
[56] Uhrenholt, A. K. and Jensen, B. S. (2019). Efficient bayesian optimization for target vector estimation. In The 22nd International Conference Artificial Intelligence and Statistics, pages 2661-2670. PMLR.
[57] Wang, K. A., Pleiss, G., Gardner, R., Tyree, S., Weinberger, K. Q., and Wilson, A. G. (2019). Exact Gaussian Processes on a Million Data Points. In Advances in Neural Information Processing Systems. arXiv: 1903.08114.
[58] Wilson, A. G., Gilboa, E., Cunningham, J. P., and Nehorai, A. (2014). Fast kernel learning for multidimensional pattern extrapolation. In Advances in Neural Information Processing Systems, volume 27, pages 3626-3634.
[59] Wilson, J., Borovitskiy, V., Terenin, A., Mostowsky, P., and Deisenroth, M. (2020a). Efficiently sampling functions from Gaussian process posteriors. In International Conference on Machine Learning, pages 10292-10302. PMLR. ISSN: 2640-3498.
[60] Wilson, J. Hutter, F., and Deisenroth, M. P. (2018). Maximizing acquisition functions for Bayesian optimization. In Advances in Neural Information Processing Systems 31, pages 9905- 9916.
[61] Wilson, J. T., Borovitskiy, V., Terenin, A., Mostowsky, P., and Deisenroth, M. (2020b). Pathwise conditioning of gaussian processes. arXiv preprint arXiv:2011.04002
13

---

[62] Yang, K., Emmerich, M.. Deutz, A., and Fonseca, C. M. (2017). Computing 3-d expected hypervolume improvement and related integrals in asymptotically optimal time. In 9th Interna- tional Conference Evolutionary Ceitt-CCrteeionn Optimization Volume 10173, EMO 2017, page 685-700, Berlin, Heidelherg. Springer- Springer-Verlag.
[63] Zhang, H. (2007). ization models. Environmetrics: The official journal of the International Environmetrics Society, 18(2):125-139
estimation for multivariate spatial linear coregional-
[64] Zhe, S., Xing, W., and Kirby, R. M. (2019). Scalable High-Order Gaussian Process Regression. In The 22nd International Conference on Artificial Intelligence and Statistics, pages 261 1-2620. PMLR. ISSN: 2640-3498.
[65] Zwicker, D. (2020). py-pde: A python package for solving partial differential equations. Journal of Open Source Software, 5(48):2158.
14


Reviews:
[{"review_id": "oO3a44_qDhz", "reviewer": "Reviewer_MXiX", "summary": "The authors devise a novel method for exact multi-task Gaussian process (GP) sampling that dramatically improves time costs, from multiplicative to additive in the combination of tasks and data points. Their method works by exploiting Kronecker structures in covariance matrices and then applying Matheron's identity. This allows them to perform Bayesian optimisation (BO) for tens of thousands of correlated outputs, which was hardly possible previously.", "questions": "", "limitations": "", "rating": 8, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The authors set out to solve the problem of multi-task BO with an extremely large number of tasks, i.e. correlated outputs. While this idea is a known and important problem, I am not aware of any work that has attempted to solve multi-task BO with this many tasks. The authors combine two known techniques, i.e. that of exploiting Kronecker structures in covariances and using Matheron's identity, and apply them to the domain of multi-task GPs, yielding impressive results for posterior sampling times. The authors adequately cited related work in Section 1 and 2, but I believe this could have been improved by having a separate section that summarises related approaches more concisely.\n\nThe derivations appear technically sound to me, with well-supported theoretical claims, in particular with regards to sampling time costs. The experimental results were extensive and impressive. The authors tested their method on a variety of problems that each showcased different strengths of their method, e.g. sampling efficiency, applicability for multi-objective and constrained BO problems, as well as extremely large output dimensions.\n\nThis paper was clearly-written and well-motivated. The results are important, in my opinion, as they allow practitioners to perform BO in settings where there are an extremely large number of output dimensions. I believe this allows practitioners to start thinking about using BO in such settings which were previously infeasible.\n\nMinor comments:\n- The \"Single-Output\" time cost for Matheron's rule in Table 1 seems to be missing a bracket, or some symbols.\n- Although the color schemes used in the experimental section look good, I think some of the colors are hard to distinguish. E.g. the purple and magenta lines in Figure 1, the Random, qParego-Batch and qEHVI-Batch lines in Figure 3, as well as the Random and EI lines in Figure 5.\n- You don't say what the shaded area represents, e.g. is it 1 or 2 standard deviations?", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The authors extensively discussed the limitations of their work in the supplementary materials (a summary of this would have been helpful in the main text). This discussion also includes a very reasonable paragraph on societal impacts.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "4"}, "extra_text": {}}, {"review_id": "AQ09Qb4gt8N", "reviewer": "Reviewer_Dpyw", "summary": "The authors propose to scale Gaussian process modeling and sampling to high-dimensional outputs with an intrinsic coregionalization model (ICM). They do it by leveraging both a Kronecker structure on the outputs (all tasks observed everywhere) and the so-called Matheron rule to update prior samples. High-order GP models are considered as well. A variety of application examples are presented to illustrate the benefits of the method and pushing the state of the art in terms of number of outputs when MC acquisition functions are used.", "questions": "", "limitations": "", "rating": 7, "confidence": 5, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The paper reads well but the methodological part is too short. Crucial points are left in the supplementary while this material raises some concerns and needs proof-reading.\n\nConsequently, the main problem --in my opinion-- is that I don’t see how to implement the method and reproduce the results without resorting to the authors code. Actual Kronecker solves, determinant expressions are not given (just a reference to [41] in supplementary). Only statements about complexity are provided (e.g., L221).  Then several cases are described, e.g., input-varying noise (or correlation in the noise), that is, not \\sigma^2 I but generic T_l as discussed in C.1, while it is unclear which ones can really be done (**). I think this can be solved by providing a detailed (complete) pseudo-code of the approach that should be in the main text (the one for HOGPs can be in supplementary). If space is needed, some of the applications can be entirely moved to supplementary as they are direct applications of existing MC acquisition functions.\n\nThe limitations discussed in supplementary A should be included in the main text. \n\n(**) In fact, from C.1.1: some of these options are more like ideas for future research than actual results:  “We do not know of an efficient solve of sums of more than two Kronecker products, and we do not have a strong implementation of approximation methods or specialized preconditioners for solves of the form (K train + T_l ) −1 z” . It is contradicting previous statements. Either really show that it works or leave it as a perspective (the supplemental is already long enough).\n\n\n\n\nTypos :\n- B.1: there are i indices in sums but that’s it.\n- L729: strucutre\n\n\n### Post-rebuttal comments ###\nI would like to thank the authors for their detailed response, which I think should be fully integrated in the paper and its supplementary to be more self contained. In particular, I appreciate the addition of the Matheron's rule Posterior Sampling for MTGPs algorithm (requiring expressions that where not referenced or given before). Assuming that all the changes proposed on the reviews will be made (perhaps adding another round of proof reading on the inevitably numerous notations), I am happy to increase my score.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Yes", "review_text": "", "extra_scores": {"ethical_concerns": "No", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "6"}, "extra_text": {}}, {"review_id": "84h8fPQZCpg", "reviewer": "Reviewer_zeBf", "summary": "This paper presents the idea of scalable multi-task Gaussian processes using a trick known as Matheron's rule. Efficient sampling from multi-task GPs is essential in solving many black box optimization problems with high-dimensional outputs including multi-objective constrained BO and composite functions. Matheron's rule is used to reduce the dependency of the number of tasks from a multiplicative factor to an additive factor. Experimental results are presented validating the benefits of the approach wrt. the computational resource requirements and also sample complexity.", "questions": "", "limitations": "", "rating": 7, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "**Strengths:**\n- The problem studied in this paper is quite significant to the BO community. BO with high dimensional outputs is a recurring problem in many applications. Traditional ways of modeling multi-task GPs do not scale well - this paper presents a scalable approach to this problem. The paper is very well written with a clear description of the method and a sound motivation.\n- The main contribution of the paper is using Matheron's rule to reduce the dependence of the number of tasks from a multiplicative factor to an additive one, leading to a large speedup. This is also verified experimentally showing at least an order of magnitude speedup without any approximations. Although Matheron's rule has been utilized in prior work, it has not been used for optimization of high-dimensional output functions, to my knowledge. As such, this contribution seems novel.\n- The proposed method has been applied to a variety of settings - constrained multi-objective BO, and composite black box functions. Other than the computational speedup, experimental results also indicate improved sample complexity when modeling multi-output GPs as MTGPs rather than standard batched GPs. Overall the experimental results seem strong and the conclusions are interesting.\n\n**Weaknesses:**\n- Some prior work has not been discussed:\n    * Uhrenholt AK, Jensen BS. Efficient Bayesian optimization for target vector estimation. AISTATS 2019\n    * Chowdhury SR, Gopalan A. No-regret algorithms for multi-task bayesian optimization. AISTATS 2021\n    * Matsui K, Kusakawa S, Ando K, Kutsukake K, Ujihara T, Takeuchi I. Bayesian Active Learning for Structured Output Design. arXiv preprint 2019\n    * Likelihood free inference is also a closely related topic. For instance, see: Gutmann MU, Corander J. Bayesian optimization for likelihood-free inference of simulator-based statistical models. JMLR 2016.\n    * Multi-objective BO references can be expanded: see the differentiable EHVI paper.\n\nMinor concerns:\n- P4 L158: should k(v_a, v_a') be k(v_(i_a), v_(j_a)) instead for a=1 to k?\n- L228-230: Are the tables wrongly referenced?", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Limitations and societal aspects have been adequately discussed in the appendix.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "6"}, "extra_text": {}}, {"review_id": "23dwu-8YNCw", "reviewer": "Reviewer_FR3J", "summary": "This paper presents an efficient sampling method for multi-output Gaussian processes and demonstrates the benefits of the proposed sampling method on Bayesian optimization with high dimensional outputs. The proposed sampling method is based on the Matheron’s rule, which was brought into the attention of the ML community by Wilson et al. (2020), and extends the formulation into the multi-output GP scenario.", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "This paper applies the Matheron’s rule to multi-output GPs and derives an efficient sampling method for multi-output GP in the exact formulation. It demonstrates the benefits of efficient multi-output GP sampling in the context of BO with high dimensional outputs. The derivation is not surprising but is a solid piece of work that can benefit the GP and BO community.\n\nThe paper is well written and easy to follow, while containing lots of details. In experiments, the paper compared the proposed method with baselines on multi-objective BO, constrained BO and BO with a composite objective function. \n\nAs the main contribution of this paper is the efficient sampling algorithm, which contains several numerical approximations, it would be nice to compare the drawn samples to the samples drawn from the exact multi-output GP formulation.\n\nThe trick of using Eigen decomposition of Kronecker product covariance matrix has been exploited in KISS-GP (Wilson&Nickisch 2015). With the Eigen decomposition trick, one can draw efficient samples on finite locations, which is sufficient for standard acquisition functions such as EI and UCB. I wonder what is the main improvement of the proposed method compared to this simple solution.\n\nIn the proposed method, the coregionalization matrix (the covariance matrix for tasks) is estimated with maximum likelihood. Due to the large number of parameters in the coregionalization matrix for high dimensional output, it often suffers from overfitting. Previous multi-task GP modeling or multi-task BO approaches tend to perform Bayesian inference for the coregionalization matrix such as MCMC sampling or variational inference. It would be interesting to compare the proposed method with other multi-task BO approaches (although those approaches might not be as scalable).\n\nSome missing reference in the multi-output GP literature:\n\n- “Kernels for vector-valued functions: A review”, MA Alvarez, L Rosasco, ND Lawrence, Foundations and Trends in Machine Learning, 2012.\n- “Efficient Modeling of Latent Information in Supervised Learning using Gaussian Processes”, Z Dai, MA Álvarez, ND Lawrence, NeuRIPS 2017.\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The authors have discussed the various limitations of the proposed method.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "2"}, "extra_text": {}}]

Author rebuttal:

Meta-review: The reviewers have provided thoughtful and constructive comments. They have responded to the authors' feedback and have reached consensus to recommend acceptance. I hope the authors will take the reviewers' comments to heart and encourage them to incorporate their thoughts in preparing the camera-ready version of their manuscript.
Decision: Accept (Poster)