Venue: NeurIPS 2021 Conference

Reviews:

final_score: rating: 6: Marginally above the acceptance threshold
review_content: 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.
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.

The 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. 

As 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.

The 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.

In 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).

Some missing reference in the multi-output GP literature:

- “Kernels for vector-valued functions: A review”, MA Alvarez, L Rosasco, ND Lawrence, Foundations and Trends in Machine Learning, 2012.
- “Efficient Modeling of Latent Information in Supervised Learning using Gaussian Processes”, Z Dai, MA Álvarez, ND Lawrence, NeuRIPS 2017.

limitations_and_societal_impact: The authors have discussed the various limitations of the proposed method.
review_title: null

final_score: rating: 8: Top 50% of accepted NeurIPS papers, clear accept
review_content: 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.
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.

The 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.

This 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.

Minor comments:
- The "Single-Output" time cost for Matheron's rule in Table 1 seems to be missing a bracket, or some symbols.
- 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.
- You don't say what the shaded area represents, e.g. is it 1 or 2 standard deviations?
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_title: null

final_score: rating: 7: Good paper, accept
review_content: 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.
main_review: **Strengths:**
- 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.
- 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.
- 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.

**Weaknesses:**
- Some prior work has not been discussed:
    * Uhrenholt AK, Jensen BS. Efficient Bayesian optimization for target vector estimation. AISTATS 2019
    * Chowdhury SR, Gopalan A. No-regret algorithms for multi-task bayesian optimization. AISTATS 2021
    * Matsui K, Kusakawa S, Ando K, Kutsukake K, Ujihara T, Takeuchi I. Bayesian Active Learning for Structured Output Design. arXiv preprint 2019
    * 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.
    * Multi-objective BO references can be expanded: see the differentiable EHVI paper.

Minor concerns:
- P4 L158: should k(v_a, v_a') be k(v_(i_a), v_(j_a)) instead for a=1 to k?
- L228-230: Are the tables wrongly referenced?
limitations_and_societal_impact: Limitations and societal aspects have been adequately discussed in the appendix.
review_title: null

final_score: rating: 7: Good paper, accept
review_content: 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.
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.

Consequently, 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.

The limitations discussed in supplementary A should be included in the main text. 

(**) 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).




Typos :
- B.1: there are i indices in sums but that’s it.
- L729: strucutre


### Post-rebuttal comments ###
I 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.
limitations_and_societal_impact: Yes
review_title: null

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)