Black-box Optimization: Bayesian Optimization and Its Multi-task Extension

The previous post derived Gaussian Process Regression from linear regression. This post uses GPR as a surrogate model to assemble the full Bayesian optimization loop, and surveys six acquisition functions (EI / EIC / NEIC / UCB / ES / PI). We then extend to the multi-task setting, covering ICM, SLFM, LMC, and PC formulations of multi-output GPs, plus the acquisition functions used in multi-task Bayesian optimization.

1. The Black-box Optimization Problem

Suppose there exists an unknown objective function $f: \mathcal{X} \rightarrow \mathbf{R}$ with the following properties:

1. Objective function characteristics

   • Black-box: the analytic form, gradients, and higher-order derivatives of $f$ are unavailable; we can only obtain observations by querying with input $\mathbf{x} \in \mathcal{X}$.
   • Noisy: the observation is $y = f(\mathbf{x}) + \epsilon$, where $\epsilon \sim \mathcal{N}(0, \tau^2)$ is additive noise (the noise is linearly superimposed on the true signal).
   • Costly to evaluate:
     • High computational cost: each evaluation of $f$ consumes substantial compute or time.
     • Limited evaluation budget: in practice the number of evaluations is strictly bounded.
   • May involve multi-objective optimization — minimising several objectives simultaneously.

2. Domain and constraints

   • The input space is $\mathcal{X} \subseteq \mathbb{R}^p$, where $p$ is the number of parameters; mixed parameter types are possible (continuous, integer, categorical).

     $$\mathcal{X} = \left\{\mathbf{x} \in \mathbb{R}^d \mid \text{parameters satisfy} \left\{\begin{array}{ll} \textbf{type constraint} & x_i \in \{\text{real, integer, categorical, ...}\} \\ \textbf{linear constraint} & \mathbf{w}^T \mathbf{x} \leqslant b \ (\text{e.g. } x_1+x_2 \leqslant 100, x_1 \leqslant x_2) \\ \textbf{nonlinear constraint} & c_j(\mathbf{x}) \leq 0, \quad j=1, \ldots, J \end{array}\right.\right\}$$

   • The constraint functions $c_j$ share the same characteristics as the objective $f$ — black-box, noisy, costly to evaluate — and can be evaluated alongside the objective.

3. Optimization objective: within a finite evaluation budget, find an approximate optimum $\mathbf{x}^* \in \mathcal{X}$ such that:

   $$\mathbf{x}^* = \arg \min_{\mathbf{x} \in \mathcal{X}} f(\mathbf{x}) \text{ subject to } \mathbf{w}^T \mathbf{x} \leqslant b \text{ and } c_j(\mathbf{x}) \leq 0, \quad j=1, \ldots, J.$$

2. Bayesian Optimization Loop

Bayesian optimization tackles the problem of an unknown and expensive-to-evaluate $f(\mathbf{x})$ by introducing a surrogate model. The surrogate $\hat{f}$ must satisfy two properties: (1) evaluating the surrogate is cheap; (2) $\hat{f}$ approximates the true function $f$ well.

To reach the optimum within as few evaluations as possible, we additionally need a function that evaluates the value of querying an unevaluated point, so that we pick the highest-value point to evaluate next. This function is called the acquisition function, and it uses the surrogate $\hat{f}$ to compute the value of each candidate point.

The basic Bayesian optimization loop:

1. Step 1: Randomly select some initial points from the parameter space and evaluate them, producing the dataset $D = (\mathbf{X}, f(\mathbf{x}))$.
2. Step 2: Fit the surrogate to the parameter space using $D$, obtaining a posterior distribution at every point.
3. Step 3: Use the acquisition function to score every point’s value based on its posterior; pick the point with the highest value (subject to constraints, if any) to evaluate next.
4. Step 4: Add the newly evaluated point to $D$. To continue, return to Step 2; to stop, go to Step 5.
5. Step 5: Pick the best already-evaluated point. Specifically:
   • If observations are noiseless (no variance), simply take the $\mathbf{x}$ with the best $f(\mathbf{x})$ from the evaluated dataset.
   • If observations are noisy, compute the predictive posterior at the evaluated points and pick the optimal $\mathbf{x}$ based on some statistic (posterior mean, median, etc.). Under constraints, the chosen point must also satisfy the constraints.

3. Common Surrogate Models

3.1 Requirements on the Surrogate Model

• The surrogate must produce a posterior distribution $p(f(\mathbf{x}) \mid \mathcal{D})$ over candidate points $\mathbf{x}$, with both a mean (prediction) and a variance (uncertainty). The predicted variance must accurately reflect the true uncertainty (especially in low-data regions).
• It must not require a pre-specified functional form (e.g. linear, polynomial); it should adaptively fit complex non-linear relationships.
• Bayesian optimization frequently updates the surrogate, so its computational complexity must be manageable.
• Real observations often contain noise $y = f(\mathbf{x}) + \epsilon$; the surrogate must be able to separate signal from noise.

3.2 Gaussian Processes

The Gaussian process is the surrogate model of choice — Gaussian process regression is an efficient way to approximate a black-box function. First, GPs are non-parametric: their model complexity adapts to the data, avoiding the structural limitations of parametric models (polynomial regression, neural networks):

• Kernel flexibility: by choosing the kernel function (e.g. RBF kernel, Matern kernel), GPs implicitly map into infinite-dimensional feature spaces, capturing non-linear relationships of arbitrary smoothness. For example:
  • The RBF kernel models infinitely differentiable functions.
  • The Matern kernel allows control over smoothness (e.g. permits some “jaggedness”).
• Data-driven fitting: GPs do not pre-specify a functional form (e.g. linear, quadratic); they learn the function shape from data via the kernel and hyperparameter optimization (e.g. maximising the marginal likelihood). This flexibility lets them fit highly non-convex, multi-modal black-box functions.

GPR also outputs a posterior distribution over the search space (mean and variance) rather than a single prediction, which pairs perfectly with acquisition functions. Emphasising the mean is conservative (Exploitation — preserve current gains); emphasising the variance is aggressive (Exploration).

The mechanics are covered in earlier sections. For kernel functions, see the kernel function chapter; for Gaussian processes and GPR, see the GPR derivation; for hyperparameter tuning of GPR, see hyperparameter optimization.

4. Common Acquisition Functions

4.1 Requirements on the Acquisition Function

The acquisition function must balance Exploitation and Exploration.

• Exploitation: sample more in regions where the current model predicts well (i.e. high-mean regions), to refine the local optimum.
  • Use case 1: the objective has a clear optimum region (single-modal or with a dominant peak).
  • Use case 2: evaluation cost is extremely high; we need fast convergence to an acceptable solution.
  • Risk: over-exploitation can trap us in a local optimum, ignoring globally promising regions.
• Exploration: sample in regions of high uncertainty (i.e. high-variance regions), to discover new potentially optimal points.
  • Use case 1: the objective may have multiple local optima (multi-modal).
  • Use case 2: the data is sparse or initial sampling is insufficient; a global search is needed.
  • Risk: over-exploration wastes resources and slows convergence.

4.2 Expected Improvement (EI)

Suppose the goal is to maximise $f(\mathbf{x})$, in a noiseless setting, and the so-far-iterated data is $D = (\mathbf{X}, \mathbf{f(X)})$. GPR gives the posterior at the candidate point $\mathbf{x}_*$ as $f(\mathbf{x}_*) \sim \mathcal{N}(\mu_{\mathbf{x}_*}, \sigma_{\mathbf{x}_*})$. Let $f^{*} = \max \limits_{i \in \{1,...,n\}} f(\mathbf{x}_i)$ be the best-so-far observation. The EI of candidate $\mathbf{x}_*$ is defined as:

$$\begin{aligned} \alpha_{\text{EI}}(\mathbf{x}_* | f^{*}) & = \mathbb{E} [\max(0, f(\mathbf{x}_*) - f^{*}) | f(\mathbf{x}_*) \sim \mathcal{N}(\mu_{\mathbf{x}_*}, \sigma_{\mathbf{x}_*})] \\ & = \mathbb{E} [\max(0, \mu_{\mathbf{x}_*} + \sigma_{\mathbf{x}_*} t - f^{*}) | t \sim \mathcal{N}(0, 1)] \\ & = \int_{-\infty}^{+\infty} \max(0, \mu_{\mathbf{x}_*} + \sigma_{\mathbf{x}_*} t - f^{*}) \phi(t) dt \\ & = \int_{\frac{f^{*} - \mu_{\mathbf{x}_*}}{\sigma_{\mathbf{x}_*}}}^{+\infty} (\mu_{\mathbf{x}_*} + \sigma_{\mathbf{x}_*} t - f^{*}) \phi(t) dt \\ & = (\mu_{\mathbf{x}_*} - f^{*}) \left[ 1 - \Phi\left(\frac{f^{*} - \mu_{\mathbf{x}_*}}{\sigma_{\mathbf{x}_*}}\right) \right] + \int_{\frac{f^{*} - \mu_{\mathbf{x}_*}}{\sigma_{\mathbf{x}_*}}}^{+\infty} \sigma_{\mathbf{x}_*} t \phi(t) dt \\ & = (\mu_{\mathbf{x}_*} - f^{*}) \left[ 1 - \Phi\left(\frac{f^{*} - \mu_{\mathbf{x}_*}}{\sigma_{\mathbf{x}_*}}\right) \right] - \frac{\sigma_{\mathbf{x}_*}}{\sqrt{2\pi}} \int_{\frac{f^{*} - \mu_{\mathbf{x}_*}}{\sigma_{\mathbf{x}_*}}}^{+\infty} e^{-\frac{t^2}{2}} d\left(-\frac{t^2}{2}\right) \\ & = (\mu_{\mathbf{x}_*} - f^{*}) \left[ 1 - \Phi\left(\frac{f^{*} - \mu_{\mathbf{x}_*}}{\sigma_{\mathbf{x}_*}}\right) \right] + \frac{\sigma_{\mathbf{x}_*}}{\sqrt{2\pi}} e^{-\frac{(f^{*} - \mu_{\mathbf{x}})^2}{2\sigma_{\mathbf{x}_*}^2}} \\ & \qquad \text{let:} \ z = \frac{f^{*} - \mu_{\mathbf{x}_*}}{\sigma_{\mathbf{x}_*}} \\ & = -\sigma_{\mathbf{x}_*} z \Phi(-z) + \sigma_{\mathbf{x}_*} \phi(z) \end{aligned}$$

• Has an analytic form; easy to compute.
• The EI of an already-observed point is 0, so we never re-sample the same point.
• Requires a definite “best observation so far”, so it only applies to the noiseless setting.

4.3 EI with Constraints (EIC)

If the optimization has $J$ non-linear constraints, the problem becomes:

$$\max_{\mathbf{x}} f(\mathbf{x}) \text{ subject to } c_j(\mathbf{x}) \leq 0, \quad j=1, \ldots, J,$$

The constraints share the same characteristics as the objective (black-box, costly, etc.); during optimization we use GPR for each constraint to obtain its posterior mean and variance. The best-so-far observation (a feasible point that satisfies the constraints) is:

$$f_c^{*} = \max \limits_{\substack{i \in \{1,\dots,n\} \\ c_j(\mathbf{x}_i) \leq 0, j = 1,\dots,J}} f(\mathbf{x}_i)$$

Assuming $\{f, c_1, ..., c_J\}$ are mutually independent:

$$\alpha_{\text{EIC}}\left(\mathbf{x}_* \mid f_c^*\right) = \alpha_{\mathrm{EI}}\left(\mathbf{x}_* \mid f_c^*\right) \prod_{j=1}^J \mathbb{P}\left(c_j(\mathbf{x}_*) \leq 0\right)$$

• In the noiseless setting, a constraint is either satisfied or not, so $\mathbb{P}\left(c_j(\mathbf{x}_*) \leq 0\right)$ is either 0 or 1.

4.4 Noisy Expected Improvement With Constraints (NEIC)

In real-world settings there are usually both constraints and noise. The objective’s observed data is $D_f = (\mathbf{X}, \mathbf{y})$ with $y_i = f(\mathbf{x}_i) + \epsilon_i$. To simplify GPR, we assume $\epsilon_{i} \in \mathcal{N}(0, \tau_i^2)$. The difficulty is that we no longer know which point in $D_f$ is the best observation. One approach is to substitute the posterior mean of each observed point — but this converges too slowly (Vazquez & Bect, 2008), and the suggested points cluster too tightly to explore the space effectively (Letham et al. 2019, arXiv:1706.07094).

Instead, we can take a weighted average over all possible outcomes based on the GPR posterior. The posterior at the observed points is:

$$\begin{aligned} \mathbf{f}^n \mid D_f & \sim \mathcal{N}\left(\boldsymbol{\mu}_f, \Sigma_f\right) \\ \mathbf{c}_j^n \mid D_{c_j} & \sim \mathcal{N}\left(\boldsymbol{\mu}_{c_j}, \Sigma_{c_j}\right), \quad j=1, \ldots, J. \end{aligned}$$

Define $D = \{D_f, D_{c_1}, ..., D_{c_J}\}$. The NEI of an unevaluated point $\mathbf{x}_*$ is:

$$\alpha_{\text{NEI}}(\mathbf{x}_* \mid D) = \int_{\mathbf{f}^n} \int_{\mathbf{c}^n} \alpha_{\text{EI}}\left(\mathbf{x}_* \mid \mathbf{f}^n, \mathbf{c}^n\right) p\left(\mathbf{f}^n \mid D_f\right) \prod_{j=1}^J p\left(\mathbf{c}_j^n \mid D_{c_j}\right) d\mathbf{c}^n d\mathbf{f}^n$$

• Note that the posterior at $\mathbf{x}_*$ itself is deterministic; only when computing EI do we need to consider the posterior at the observed points.
• No analytic form, but value and gradient can both be computed efficiently, so it is practical.
• For batch sampling (producing multiple samples in one step): produce one sample, add its posterior to $\mathbf{f}^n$ and $\mathbf{c}_n$, then produce the next sample, and so on.
• Like EI, the NEI of an already-observed point is 0, so we never re-sample.

4.5 Upper Confidence Bound (UCB)

For maximisation, UCB is defined as the upper end of the posterior interval:

$$\alpha_{\text{UCB}}\left(\mathbf{x}_*\right) = \mu_{\mathbf{x}_*} + \beta \sigma_{\mathbf{x}_*}$$

For minimisation, UCB is the lower end of the posterior interval:

$$\alpha_{\text{UCB}}\left(\mathbf{x}_*\right) = \mu_{\mathbf{x}_*} - \beta \sigma_{\mathbf{x}_*}$$

where:

• $\mu_{\mathbf{x}_*}$: posterior mean at the candidate $\mathbf{x}_*$.
• $\sigma_{\mathbf{x}_*}$: posterior standard deviation at the candidate.
• $\beta \geq 0$: balance parameter. A larger $\beta$ increases exploration; a smaller $\beta$ emphasises exploitation.

• UCB is computationally efficient — only the posterior mean and variance are needed.
• $\beta$ must be chosen carefully. A fixed $\beta$ can imbalance the double-E; an adaptive strategy is generally preferred.

Theoretical Guarantee for UCB

Under certain conditions, the UCB algorithm guarantees that the cumulative regret $R(T)$ grows sub-linearly, in which case the algorithm is guaranteed to converge to the optimum:

$$R(T) = \sum_{t=1}^T f\left(\mathbf{x}^*\right) - f\left(\mathbf{x}_t\right) \leq \tilde{O}(\sqrt{T})$$

where:

• $\mathbf{x}^*$ is the global optimum.
• $\mathbf{x}_t$ is the point chosen at step $t$.
• $\tilde{O}$ denotes an asymptotic upper bound that ignores logarithmic factors.

4.6 Entropy Search (ES)

ES aims to minimise the uncertainty about the objective’s argmax $\mathbf{x}_{\text{max}}$ (measured by information entropy). Given the observation set $D$, define:

• Current entropy:

  $$H\left(\mathbf{x}_{\max} \mid D\right) = -\int p\left(\mathbf{x}_{\max} \mid D\right) \log p\left(\mathbf{x}_{\max} \mid D\right) d\mathbf{x}_{\max}$$

• Expected entropy (the conditional entropy after sampling at $\mathbf{x}_*$):

  $$\mathbb{E}_{f(\mathbf{x}_*)} \left[ H(\mathbf{x}_{\max} \mid D \cup \{(\mathbf{x}_*, f(\mathbf{x}_*))\}) \right]$$

The ES acquisition function is the expected reduction in entropy:

$$\alpha_{\mathrm{ES}}\left(\mathbf{x}_*\right) = H\left(\mathbf{x}_{\max} \mid D\right) - \mathbb{E}_{f\left(\mathbf{x}_*\right)}\left[H\left(\mathbf{x}_{\max} \mid D \cup \left\{\left(\mathbf{x}_*, f\left(\mathbf{x}_*\right)\right)\right\}\right)\right]$$

ES is computationally expensive because $p\left(\mathbf{x}_{\max} \mid D\right)$ has no analytic form and must be approximated (Monte Carlo sampling, etc.), so it is suitable for low-dimensional problems (within roughly 10 dimensions).

4.7 Probability of Improvement (PI)

PI measures the probability that a candidate $\mathbf{x}_*$ improves over the best-so-far observation $f^*$. Based on the posterior $f\left(\mathbf{x}_*\right) \sim \mathcal{N}\left(\mu_{\mathbf{x}_*}, \sigma_{\mathbf{x}_*}\right)$:

$$\alpha_{\mathrm{PI}}\left(\mathbf{x}_* \mid f^*\right) = \mathrm{P}\left(f\left(\mathbf{x}_*\right) \geq f^*\right) = \Phi\left(\frac{\mu_{\mathbf{x}_*} - f^*}{\sigma_{\mathbf{x}_*}}\right)$$

where $\Phi$ is the standard normal cumulative distribution function (CDF). An exploration parameter $\xi \geq 0$ can be added to avoid local optima:

$$\alpha_{\mathrm{PI}}\left(\mathbf{x}_* \mid f^*, \xi\right) = \Phi\left(\frac{\mu_{\mathbf{x}_*} - \left(f^* + \xi\right)}{\sigma_{\mathbf{x}_*}}\right)$$

• Has an analytic form, depending only on the posterior mean and variance, so it is computationally efficient.
• Without $\xi$, PI tends to over-exploit (prefers high-mean, low-variance regions).
• PI ignores the magnitude of the improvement, whereas EI balances probability and magnitude — so EI is more commonly used.

5. Multi-task Gaussian Process Regression

Machine learning has a paradigm called multi-task learning, whose core idea is to share information across different tasks — relative to learning each task in isolation, this can boost performance. In an optimization context, we can use multi-task optimization to jointly optimise multiple correlated tasks for a higher payoff; or, for an expensive online optimization, we can construct biased-but-related offline optimization tasks and use multi-task learning to enhance the online optimization. In the Gaussian process regression setting this is also called Multi-Output Learning or Vector-Valued Learning.

This section follows the work of Álvarez et al. (Álvarez 2017 slides, Álvarez, Rosasco & Lawrence 2012 (arXiv:1106.6251)).

5.1 Overall Framework

The earlier single-task GPR section covers single-task Gaussian process regression. For the multi-task setting, we define a vector-valued function $\mathbf{f} = \left[f_1, \ldots, f_O\right]^T$ that follows a Gaussian process:

$$\mathbf{f} \sim \mathcal{G} \mathcal{P}(\mathbf{m}, \mathbf{K})$$

where:

• $\mathbf{m}$ is the mean function vector — an $O$-dimensional vector of functions: $[m_1(\mathbf{x}), ..., m_O(\mathbf{x})]^T$.
• $\mathbf{K}$ is the covariance function matrix — an $O \times O$ positive matrix-valued function. Each entry is a covariance function as in single-task GPR. The element at row $o$ column $o^\prime$ of $\mathbf{K}(\mathbf{x}, \mathbf{x}^{\prime})$, written $(\mathbf{K}(\mathbf{x}, \mathbf{x}^{\prime}))_{o, o^{\prime}}$, represents the correlation between tasks $f_o(\mathbf{x})$ and $f_{o^\prime}(\mathbf{x}^{\prime})$.

For the input dataset $\mathbf{X}$, the prior is:

$$\mathbf{f}(\mathbf{X}) \sim \mathcal{N}(\mathbf{m}(\mathbf{X}), \mathbf{K}(\mathbf{X}, \mathbf{X}))$$

with covariance function matrix:

$$\mathbf{K}(\mathbf{X}, \mathbf{X}) = \left[\begin{array}{ccc} \left(\mathbf{K}\left(\mathbf{X}_1, \mathbf{X}_1\right)\right)_{1,1} & \cdots & \left(\mathbf{K}\left(\mathbf{X}_1, \mathbf{X}_O\right)\right)_{1, O} \\ \left(\mathbf{K}\left(\mathbf{X}_2, \mathbf{X}_1\right)\right)_{2,1} & \cdots & \left(\mathbf{K}\left(\mathbf{X}_2, \mathbf{X}_O\right)\right)_{2, O} \\ \vdots & \cdots & \vdots \\ \left(\mathbf{K}\left(\mathbf{X}_O, \mathbf{X}_1\right)\right)_{O, 1} & \cdots & \left(\mathbf{K}\left(\mathbf{X}_O, \mathbf{X}_O\right)\right)_{O, O} \end{array}\right]$$

If the $O$ functions are mutually independent, then $\mathbf{K}(\mathbf{X}, \mathbf{X})$ is a block-diagonal matrix, with each block being the covariance matrix of one function:

$$\left[\begin{array}{cccc} \left(\mathbf{K}\left(\mathbf{X}_1, \mathbf{X}_1\right)\right)_{1,1} & 0 & \cdots & 0 \\ 0 & \left(\mathbf{K}\left(\mathbf{X}_2, \mathbf{X}_2\right)\right)_{2,2} & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & \left(\mathbf{K}\left(\mathbf{X}_O, \mathbf{X}_O\right)\right)_{O, O} \end{array}\right]$$

Likelihood:

$$p(\mathbf{y} \mid \mathbf{f}, \mathbf{x}, \Sigma) = \mathcal{N}(\mathbf{f}(\mathbf{x}), \Sigma)$$

where $\Sigma$ is a diagonal matrix with diagonal entries $\{\sigma_o^2\}_{o=1}^{O}$.

Predictive posterior:

$$p\left(\mathbf{f}\left(\mathbf{x}_*\right) \mid D, \mathbf{f}, \mathbf{x}_*, \boldsymbol{\theta}\right) = \mathcal{N}\left(\mathbf{f}_*\left(\mathbf{x}_*\right), \mathbf{K}_*\left(\mathbf{x}_*, \mathbf{x}_*\right)\right),$$

where:

$$\begin{aligned} \mathbf{f}_*\left(\mathbf{x}_*\right) & = \mathbf{K}_{\mathbf{x}_*}^{\top}(\mathbf{K}(\mathbf{X}, \mathbf{X}) + \boldsymbol{\Sigma})^{-1} \overline{\mathbf{y}}, \\ \mathbf{K}_*\left(\mathbf{x}_*, \mathbf{x}_*\right) & = \mathbf{K}\left(\mathbf{x}_*, \mathbf{x}_*\right) - \mathbf{K}_{\mathbf{x}_*}(\mathbf{K}(\mathbf{X}, \mathbf{X}) + \boldsymbol{\Sigma})^{-1} \mathbf{K}_{\mathbf{x}_*}^{\top}, \end{aligned}$$

with:

• $\boldsymbol{\Sigma} = \Sigma \otimes \mathbf{I}_N$.
• $\mathbf{K}_{\mathbf{x}_*} \in \mathbf{R}^{O \times N O}$ has entries $\left(\mathbf{K}\left(\mathbf{x}_*, \mathbf{x}_j\right)\right)_{o, o^{\prime}}$ for $j=1, \ldots, N$ and $o, o^{\prime}=1, \ldots, O$.
• $\boldsymbol{\theta}$: the set of all kernel hyperparameters.

In practice, proposing a model usually starts from an assumption — for the multi-task setting, we first assume something about the data distribution and the inter-task relationships, then derive the corresponding model. The next sections cover models derived from different assumptions, including:

• Separable Kernels and Sum of Separable Kernels (SoS kernels): A special class of multi-output kernel functions that can be formulated as a sum of products between a kernel function for the input space alone, and a kernel function that encodes the interactions among the outputs. Includes ICM (Intrinsic Coregionalization Model) and SLFM (Semiparametric Latent Factor Model) (Linear Model of Coregionalization). LMC is a generalisation of ICM and SLFM.
• Beyond Separable Kernels: more sophisticated models, such as PC (Process Convolutions). Process convolutions further generalise LMC.

Glossary

• cokriging: probabilistic models and Gaussian processes for multivariate learning were originally pioneered and developed in geostatistics, where prediction over vector-valued output data is called cokriging (multivariate co-kriging).
• autokrigeability: a term describing the special condition under which multivariate cokriging can be reduced to univariate kriging. When the data satisfies certain conditions, the covariance matrix of the multivariate cokriging model becomes diagonal (block-diagonal, indicating independence among variables), and the problem automatically decouples into multiple independent univariate kriging problems.
• Isotopic / Heterotopic Data: Isotopic means the input datasets across different tasks are exactly the same; heterotopic means the input datasets across tasks differ.

5.2 Intrinsic Coregionalization Model (ICM)

ICM assumes that multiple outputs are obtained as different linear transformations of multiple samples drawn from the same Gaussian process regression.

Two-output Setting

We derive ICM from the two-output setting. Consider two outputs $f_1(\mathbf{x})$ and $f_2(\mathbf{x})$ that share the same input space $\mathbf{x} \in \mathbf{R}^p$. The corresponding vector-valued function is:

$$\mathbf{f}(\mathbf{x}) = \left[\begin{array}{l} f_1(\mathbf{x}) \\ f_2(\mathbf{x}) \end{array}\right]$$

Covariance:

$$\begin{aligned} & \operatorname{cov}\left(\mathbf{f}(\mathbf{x}), \mathbf{f}\left(\mathbf{x}^{\prime}\right)\right) \\ & = \mathbb{E}\left\{\mathbf{f}(\mathbf{x})\left[\mathbf{f}\left(\mathbf{x}^{\prime}\right)\right]^{\top}\right\} - \mathbb{E}\{\mathbf{f}(\mathbf{x})\}\left[\mathbb{E}\left\{\mathbf{f}\left(\mathbf{x}^{\prime}\right)\right\}\right]^{\top} \\ & = \mathbb{E}\left\{\left[\begin{array}{l} f_1(\mathbf{x}) \\ f_2(\mathbf{x}) \end{array}\right]\left[\begin{array}{ll} f_1\left(\mathbf{x}^{\prime}\right) & f_2\left(\mathbf{x}^{\prime}\right) \end{array}\right]\right\} - \mathbb{E}\left\{\left[\begin{array}{l} f_1(\mathbf{x}) \\ f_2(\mathbf{x}) \end{array}\right]\right\} \left[ \mathbb{E}\left\{\left[\begin{array}{l} f_1(\mathbf{x}^\prime) \\ f_2(\mathbf{x}^\prime) \end{array}\right]\right\} \right]^T \\ & = \mathbb{E} \left\{\left[\begin{array}{ll} f_1(\mathbf{x}) f_1\left(\mathbf{x}^{\prime}\right) & f_1(\mathbf{x}) f_2\left(\mathbf{x}^{\prime}\right) \\ f_2(\mathbf{x}) f_1\left(\mathbf{x}^{\prime}\right) & f_2(\mathbf{x}) f_2\left(\mathbf{x}^{\prime}\right) \end{array}\right] \right\} - \left[\begin{array}{l} \mathbb{E}\left\{f_1(\mathbf{x})\right\} \\ \mathbb{E}\left\{f_2(\mathbf{x})\right\} \end{array}\right] \left[\begin{array}{l} \mathbb{E}\left\{f_1(\mathbf{x}^\prime)\right\} \\ \mathbb{E}\left\{f_2(\mathbf{x}^\prime)\right\} \end{array}\right]^T \\ & = \left[\begin{array}{ll} \mathbb{E}\left\{f_1(\mathbf{x}) f_1\left(\mathbf{x}^{\prime}\right)\right\} & \mathbb{E}\left\{f_1(\mathbf{x}) f_2\left(\mathbf{x}^{\prime}\right)\right\} \\ \mathbb{E}\left\{f_2(\mathbf{x}) f_1\left(\mathbf{x}^{\prime}\right)\right\} & \mathbb{E}\left\{f_2(\mathbf{x}) f_2\left(\mathbf{x}^{\prime}\right)\right\} \end{array}\right] - \left[\begin{array}{ll} \mathbb{E}\left\{f_1(\mathbf{x}) \right\} \mathbb{E}\left\{ f_1\left(\mathbf{x}^{\prime}\right)\right\} & \mathbb{E}\left\{f_1(\mathbf{x}) \right\} \mathbb{E}\left\{ f_2\left(\mathbf{x}^{\prime}\right)\right\} \\ \mathbb{E}\left\{f_2(\mathbf{x}) \right\} \mathbb{E}\left\{ f_1\left(\mathbf{x}^{\prime}\right)\right\} & \mathbb{E}\left\{f_2(\mathbf{x}) \right\} \mathbb{E}\left\{ f_2\left(\mathbf{x}^{\prime}\right)\right\} \end{array}\right] \end{aligned}$$

Based on one sampled function

1. Sample one function $\mu^1(\mathbf{x})$ from the Gaussian process $\mu(\mathbf{x}) \sim \mathcal{GP}(\mathbf{0}, k(\mathbf{x}, \mathbf{x}^{\prime}))$. The sampled function $\mu^1(\mathbf{x})$ is also called a latent function.

2. Assume $f_1(\mathbf{x})$ and $f_2(\mathbf{x})$ are obtained from the following linear transformation of $\mu^1(\mathbf{x})$:

   $$\begin{aligned} & f_1(\mathbf{x}) = a_1^1 \mu^1(\mathbf{x}) \\ & f_2(\mathbf{x}) = a_2^1 \mu^1(\mathbf{x}) \end{aligned}$$

The vector-valued function is:

$$\begin{aligned} \mathbf{f}(\mathbf{x}) & = \left[\begin{array}{l} f_1(\mathbf{x}) \\ f_2(\mathbf{x}) \end{array}\right] = \left[\begin{array}{l} a_1^1 \mu^1(\mathbf{x}) \\ a_2^1 \mu^1(\mathbf{x}) \end{array}\right] = \left[\begin{array}{l} a_1^1 \\ a_2^1 \end{array}\right] \mu^1(\mathbf{x}) \\ & \quad \, \, \text{define } \mathbf{a} = \left[\begin{array}{ll} a_1^1 & a_2^1 \end{array}\right]^{T} \\ & = \mathbf{a} \mu^1(\mathbf{x}) \end{aligned}$$

The covariance function is:

$$\begin{aligned} & \operatorname{cov}\left(\mathbf{f}(\mathbf{x}), \mathbf{f}\left(\mathbf{x}^{\prime}\right)\right) \\ & = \left[\begin{array}{ll} \mathbb{E}\left\{f_1(\mathbf{x}) f_1\left(\mathbf{x}^{\prime}\right)\right\} & \mathbb{E}\left\{f_1(\mathbf{x}) f_2\left(\mathbf{x}^{\prime}\right)\right\} \\ \mathbb{E}\left\{f_2(\mathbf{x}) f_1\left(\mathbf{x}^{\prime}\right)\right\} & \mathbb{E}\left\{f_2(\mathbf{x}) f_2\left(\mathbf{x}^{\prime}\right)\right\} \end{array}\right] - \left[\begin{array}{ll} \mathbb{E}\left\{f_1(\mathbf{x}) \right\} \mathbb{E}\left\{ f_1\left(\mathbf{x}^{\prime}\right)\right\} & \mathbb{E}\left\{f_1(\mathbf{x}) \right\} \mathbb{E}\left\{ f_2\left(\mathbf{x}^{\prime}\right)\right\} \\ \mathbb{E}\left\{f_2(\mathbf{x}) \right\} \mathbb{E}\left\{ f_1\left(\mathbf{x}^{\prime}\right)\right\} & \mathbb{E}\left\{f_2(\mathbf{x}) \right\} \mathbb{E}\left\{ f_2\left(\mathbf{x}^{\prime}\right)\right\} \end{array}\right] \\ & = \left[\begin{array}{ll} \left(a_1^1\right)^2 \mathbb{E}\left\{u^1(\mathbf{x}) u^1\left(\mathbf{x}^{\prime}\right)\right\} & a_1^1 a_2^1 \mathbb{E}\left\{u^1(\mathbf{x}) u^1\left(\mathbf{x}^{\prime}\right)\right\} \\ a^1 a^2 \mathbb{E}\left\{u^1(\mathbf{x}) u^1\left(\mathbf{x}^{\prime}\right)\right\} & \left(a_2^1\right)^2 \mathbb{E}\left\{u^1(\mathbf{x}) u^1\left(\mathbf{x}^{\prime}\right)\right\} \end{array}\right] \\ & \quad \, \, - \left[\begin{array}{ll} \left(a_1^1\right)^2 \mathbb{E}\left\{u^1(\mathbf{x}) \right\} \mathbb{E}\left\{ u^1\left(\mathbf{x}^{\prime}\right)\right\} & a_1^1 a_2^1 \mathbb{E}\left\{u^1(\mathbf{x}) \right\} \mathbb{E}\left\{ u^1\left(\mathbf{x}^{\prime}\right)\right\} \\ a^1 a^2 \mathbb{E}\left\{u^1(\mathbf{x}) \right\} \mathbb{E}\left\{ u^1\left(\mathbf{x}^{\prime}\right)\right\} & \left(a_2^1\right)^2 \mathbb{E}\left\{u^1(\mathbf{x}) \right\} \mathbb{E}\left\{ u^1\left(\mathbf{x}^{\prime}\right)\right\} \end{array}\right] \\ & = \left[\begin{array}{ll} \left(a_1^1\right)^2 & a_1^1 a_2^1 \\ a^1 a^2 & \left(a_2^1\right)^2 \end{array}\right] \left[ \mathbb{E}\left\{u^1(\mathbf{x}) u^1\left(\mathbf{x}^{\prime}\right)\right\} - \mathbb{E}\left\{u^1(\mathbf{x})\right\} \mathbb{E}\left\{u^1\left(\mathbf{x}^{\prime}\right)\right\} \right] \\ & = \left[\begin{array}{l} a_1^1 \\ a_2^1 \end{array}\right]\left[\begin{array}{ll} a_1^1 & a_2^1 \end{array}\right] \underbrace{\left[\mathbb{E}\left\{u^1(\mathbf{x}) u^1\left(\mathbf{x}^{\prime}\right)\right\} - \mathbb{E}\left\{u^1(\mathbf{x})\right\} \mathbb{E}\left\{u^1\left(\mathbf{x}^{\prime}\right)\right\}\right]}_{k\left(\mathbf{x}, \mathbf{x}^{\prime}\right)} \\ & \quad \, \, \text{define } \mathbf{a} = \left[\begin{array}{ll} a_1^1 & a_2^1 \end{array}\right]^{\top} \\ & = \mathbf{a}\mathbf{a}^{\top} k\left(\mathbf{x}, \mathbf{x}^{\prime}\right) \\ & \quad \, \, \text{define } \mathbf{B} = \mathbf{a}\mathbf{a}^{\top}, \mathbf{B} \text{ has rank one} \\ & = \mathbf{B} k\left(\mathbf{x}, \mathbf{x}^{\prime}\right)\\ & = \left[\begin{array}{ll} b_{11} & b_{12} \\ b_{21} & b_{22} \end{array}\right] k\left(\mathbf{x}, \mathbf{x}^{\prime}\right) \end{aligned}$$

Or equivalently:

$$\begin{aligned} \operatorname{cov}\left(\mathbf{f}(\mathbf{x}), \mathbf{f}\left(\mathbf{x}^{\prime}\right)\right) = \operatorname{cov}(\mathbf{a}\mu^1(\mathbf{x}), \mathbf{a}\mu^1(\mathbf{x}^\prime)) = \mathbf{a}\mathbf{a}^{\top} \operatorname{cov}(\mu^1(\mathbf{x}), \mu^1(\mathbf{x})) = \mathbf{a}\mathbf{a}^{\top} k\left(\mathbf{x}, \mathbf{x}^{\prime}\right) \end{aligned}$$

Based on two sampled functions

1. Sample two functions $\mu^1(\mathbf{x})$ and $\mu^2(\mathbf{x})$ from the Gaussian process $\mu(\mathbf{x}) \sim \mathcal{GP}(\mathbf{0}, k(\mathbf{x}, \mathbf{x}^{\prime}))$. The two are mutually independent and share the same covariance function $k(\mathbf{x}, \mathbf{x}^{\prime})$.

2. Assume $f_1(\mathbf{x})$ and $f_2(\mathbf{x})$ are obtained from the following linear transformations of $\mu^1(\mathbf{x})$ and $\mu^2(\mathbf{x})$:

   $$\begin{aligned} & f_1(\mathbf{x}) = a_1^1 \mu^1(\mathbf{x}) + a_1^2 \mu^2(\mathbf{x}) \\ & f_2(\mathbf{x}) = a_2^1 \mu^1(\mathbf{x}) + a_2^2 \mu^2(\mathbf{x}) \end{aligned}$$

The vector-valued function is:

$$f(\mathbf{x}) = u^1(\mathbf{x}) + \mathbf{a}^2 u^2(\mathbf{x})$$ $$\begin{aligned} \mathbf{f}(\mathbf{x}) & = \left[\begin{array}{l} f_1(\mathbf{x}) \\ f_2(\mathbf{x}) \end{array}\right] = \left[\begin{array}{l} a_1^1 \\ a_2^1 \end{array}\right] \mu^1(\mathbf{x}) + \left[\begin{array}{l} a_1^2 \\ a_2^2 \end{array}\right] \mu^2(\mathbf{x}) \\ & \quad \text{define } \mathbf{a}^1 = \left[\begin{array}{ll} a_1^1 & a_2^1 \end{array}\right]^T, \mathbf{a}^2 = \left[\begin{array}{ll} a_1^2 & a_2^2 \end{array}\right]^T \\ & = \mathbf{a}^1 u^1(\mathbf{x}) + \mathbf{a}^2 u^2(\mathbf{x}) \end{aligned}$$

The covariance function:

$$\begin{aligned} \operatorname{cov}\left(f(\mathbf{x}), f\left(\mathbf{x}^{\prime}\right)\right) & = \mathbf{a}^1\left(\mathbf{a}^1\right)^{\top} \operatorname{cov}\left(u^1(\mathbf{x}), u^1\left(\mathbf{x}^{\prime}\right)\right) + \mathbf{a}^2\left(\mathbf{a}^2\right)^{\top} \operatorname{cov}\left(u^2(\mathbf{x}), u^2\left(\mathbf{x}^{\prime}\right)\right) \\ & = \mathbf{a}^1\left(\mathbf{a}^1\right)^{\top} k\left(\mathbf{x}, \mathbf{x}^{\prime}\right) + \mathbf{a}^2\left(\mathbf{a}^2\right)^{\top} k\left(\mathbf{x}, \mathbf{x}^{\prime}\right) \\ & = \left[\mathbf{a}^1\left(\mathbf{a}^1\right)^{\top} + \mathbf{a}^2\left(\mathbf{a}^2\right)^{\top}\right] k\left(\mathbf{x}, \mathbf{x}^{\prime}\right) \\ & \quad \text{define } \mathbf{B} = \mathbf{a}^1\left(\mathbf{a}^1\right)^{\top} + \mathbf{a}^2\left(\mathbf{a}^2\right)^{\top} \text{ which has rank two.} \\ & = \mathbf{B} k\left(\mathbf{x}, \mathbf{x}^{\prime}\right) = \left[\begin{array}{ll} b_{11} & b_{12} \\ b_{21} & b_{22} \end{array}\right] k\left(\mathbf{x}, \mathbf{x}^{\prime}\right) \end{aligned}$$

Generalisation to $R$ samples and $O$ outputs

In ICM we assume the $O$ outputs $\{f_o(\mathbf{x})\}_{o = 1}^O$ are all obtained from $R$ independent samples of the same Gaussian process. These samples are mutually independent and share a single covariance function:

$$f_o(\mathbf{x}) = \sum_{r=1}^R a_o^r u^r(\mathbf{x})$$

Define:

$$\begin{aligned} \mathbf{a}^r & = \left[ \begin{array}{llll} a_1^r & a_2^r & ... & a_O^r \end{array} \right]^T \\ \mathbf{A} & = \left[ \begin{array}{llll} \mathbf{a}^1 & \mathbf{a}^2 & ... & \mathbf{a}^R \end{array} \right] \\ \boldsymbol{\mu}(\mathbf{x}) & = \left[ \begin{array}{llll} \mu^1(\mathbf{x}) & \mu^2(\mathbf{x}) & ... & \mu^R(\mathbf{x}) \end{array} \right]^T \end{aligned}$$

The vector-valued function is:

$$\mathbf{f}(\mathbf{x}) = \left[f_1(\mathbf{x}) \cdots f_O(\mathbf{x})\right]^{\top} = \mathbf{A} \boldsymbol{\mu}(\mathbf{x})$$

The covariance function is:

$$\begin{aligned} \operatorname{cov}\left[\mathbf{f}(\mathbf{x}), \mathbf{f}\left(\mathbf{x}^{\prime}\right)\right] & = \mathbf{A} \mathbf{A}^{\top} k\left(\mathbf{x}, \mathbf{x}^{\prime}\right) \\ & \quad \text{define } \mathbf{B} = \mathbf{A}\mathbf{A}^{\top}, \mathbf{B} \in \mathbf{R}^{O \times O}, \text{ B has rank } R \\ & = \mathbf{B} k\left(\mathbf{x}, \mathbf{x}^{\prime}\right) \end{aligned}$$

ICM properties:

• If the outputs are noiseless and the data is isotopic, then ICM yields the same result as running single-task GPR on each output independently. This automatic reduction from cokriging to kriging is called autokrigeability.

5.3 Semiparametric Latent Factor Model (SLFM)

Like ICM, SLFM assumes that multiple outputs are different linear transformations of multiple sampled functions — but here the samples come from different, mutually independent GPRs (one function sampled from each GPR). With $Q$ GPRs $\mu_q(\mathbf{x}) \sim \mathcal{GP}(\mathbf{0}, k_q(\mathbf{x}, \mathbf{x}^\prime)), q \in \{1, 2, ..., Q\}$, we get $Q$ sampled functions $\{\mu_q(\mathbf{x})\}$. $Q = 1$ corresponds to the $R = 1$ case of ICM.

As in the generalised ICM, for $O$ outputs:

$$f_o(\mathbf{x}) = \sum_{q=1}^Q a_{o, q} \mu_q(\mathbf{x})$$

Define:

$$\mathbf{a}_q = \left[\begin{array}{llll} a_{1,q} & a_{2,q} & ... & a_{O,q} \end{array}\right]^T$$

The vector-valued function is:

$$\mathbf{f}(\mathbf{x}) = \left[f_1(\mathbf{x}) \cdots f_O(\mathbf{x})\right]^{\top} = \sum_{q = 1}^{Q} \mathbf{a}_q \mu_q(\mathbf{x})$$

The covariance function is:

$$\begin{aligned} \operatorname{cov}\left[\mathbf{f}(\mathbf{x}), \mathbf{f}\left(\mathbf{x}^{\prime}\right)\right] & = \sum_{q=1}^Q \mathbf{a}_q \mathbf{a}_{q}^{\top} k_q\left(\mathbf{x}, \mathbf{x}^{\prime}\right) \\ & \quad \text{define } \mathbf{B}_q = \mathbf{a}_q \mathbf{a}_{q}^{\top}, \text{ The rank of each } \mathbf{B}_q \in \mathbf{R}^{O \times O} \text{ is one.} \\ & = \sum_{q=1}^Q \mathbf{B}_q k_q\left(\mathbf{x}, \mathbf{x}^{\prime}\right) \end{aligned}$$

5.4 Linear Model of Coregionalization (LMC)

LMC is a further generalisation of ICM and SLFM. The base assumption is the same — multiple outputs are different linear combinations of sampled functions — but ICM assumes the samples come from a single Gaussian process; SLFM assumes they come from completely different Gaussian processes; LMC assumes some samples come from the same Gaussian process and others from different ones.

For $O$ outputs:

$$f_o(\mathbf{x}) = \sum_{q=1}^Q \sum_{r=1}^{R_q} a_{o, q}^r u_q^r(\mathbf{x})$$

There are $Q$ groups of samples, each containing $R_q$ sampled functions. The $R_q$ functions in a group are mutually independent and share the same Gaussian process with covariance $k_q(\mathbf{x}, \mathbf{x}^\prime)$; the $Q$ groups are also mutually independent. This is equivalent to summing $Q$ ICMs.

Define:

$$\begin{aligned} \mathbf{a}_q^r & = \left[ \begin{array}{llll} a_{1, q}^r & a_{2, q}^r & ... & a_{O, q}^r \end{array} \right]^T \\ \mathbf{A}_q & = \left[ \begin{array}{llll} \mathbf{a}_q^1 & \mathbf{a}_q^2 & ... & \mathbf{a}_q^R \end{array} \right] \\ \boldsymbol{\mu}_q(\mathbf{x}) & = \left[ \begin{array}{llll} \mu_q^1(\mathbf{x}) & \mu_q^2(\mathbf{x}) & ... & \mu_q^R(\mathbf{x}) \end{array} \right]^T \end{aligned}$$

The vector-valued function is:

$$\mathbf{f}(\mathbf{x}) = \left[f_1(\mathbf{x}) \cdots f_O(\mathbf{x})\right]^{\top} = \sum_{q = 1}^{Q} \mathbf{A}_q \boldsymbol{\mu}_q(\mathbf{x})$$

The covariance function:

$$\operatorname{cov}\left[\mathbf{f}(\mathbf{x}), \mathbf{f}\left(\mathbf{x}^{\prime}\right)\right] = \sum_{q=1}^Q \mathbf{A}_q \mathbf{A}_q^{\top} k_q\left(\mathbf{x}, \mathbf{x}^{\prime}\right) = \sum_{q=1}^Q \mathbf{B}_q k_q\left(\mathbf{x}, \mathbf{x}^{\prime}\right)$$

where $\mathbf{B}_q$ is usually called the coregionalization matrices, with rank $R_q$.

5.5 Process Convolutions (PC)

In LMC, the value of the vector-valued function $\mathbf{f}(\mathbf{x})$ at point $\mathbf{x}$ depends only on the values of the sampled functions at point $\mathbf{x}$, not on their values elsewhere. This simplification gives LMC the form of a separable kernel. PC convolves a smoothing kernel with the sampled functions, making each output depend on the values of the sampled functions at all points — a further generalisation of LMC. Because the sampled functions come from GPR, the convolution result is still a Gaussian process, so the posterior remains analytic.

About the smoothing kernel:

Here the smoothing kernel is a class of functions, but it does not refer to a kernel function. The convolution process uses the smoothing kernel to take a weighted average over a neighbourhood, producing smoother results — this is sometimes also called a moving average function.

As in LMC, assume there are $Q$ groups of samples, each with $R_q$ sampled functions. The $R_q$ sampled functions $\mu_q^r(\mathbf{x}), i \in \{1, ..., R_q\}$ are mutually independent and share the same Gaussian process with covariance $k_q(\mathbf{x}, \mathbf{x}^\prime)$; the $Q$ groups are also mutually independent. For $O$ outputs $\left\{f_o(\mathbf{x})\right\}_{o=1}^O$:

$$f_o(\mathbf{x}) = \sum_{q=1}^Q \sum_{r=1}^{R_q} \int_{\mathcal{X}} G_{o, q}^r(\mathbf{x} - \mathbf{z}) u_q^r(\mathbf{z}) \mathrm{d} \mathbf{z} + w_d(\mathbf{x}) = \sum_{q=1}^Q f_o^q(\mathbf{x}) + w_o(\mathbf{x}),$$

where:

$$f_o^q(\mathbf{x}) = \sum_{r=1}^{R_q} \int_{\mathcal{X}} G_{o, q}^r(\mathbf{x} - \mathbf{z}) u_q^r(\mathbf{z}) \mathrm{d} \mathbf{z}$$

and:

• $\left\{w_o(\mathbf{x})\right\}_{o=1}^O$ are mutually independent Gaussian processes with mean $\mathbf{0}$ and covariance $k_{w_o}(\mathbf{x}, \mathbf{x}^\prime)$.
• For the integral above to converge, $G_{o, q}^r(\mathbf{x})$ must be a continuous function with compact support (locally non-zero, so the integral does not blow up) or be square-integrable.

The covariance function of $f_o(\mathbf{x})$ and $f_{o^\prime}(\mathbf{x}^\prime)$:

$$\left(\mathbf{K}\left(\mathbf{x}, \mathbf{x}^{\prime}\right)\right)_{o, o^{\prime}} = \sum_{q=1}^Q k_{f_o^q, f_{o^{\prime}}^q}\left(\mathbf{x}, \mathbf{x}^{\prime}\right) + k_{w_o}\left(\mathbf{x}, \mathbf{x}^{\prime}\right) \delta_{o, o^{\prime}}$$

where:

$$k_{f_o^q, f_{o^{\prime}}^q}\left(\mathbf{x}, \mathbf{x}^{\prime}\right) = \sum_{i=1}^{R_q} \int_{\mathcal{X}} G_{o, q}^i(\mathbf{x} - \mathbf{z}) \int_{\mathcal{X}} G_{o^{\prime}, q}^i\left(\mathbf{x}^{\prime} - \mathbf{z}^{\prime}\right) k_q\left(\mathbf{z}, \mathbf{z}^{\prime}\right) \mathrm{d} \mathbf{z}^{\prime} \mathrm{d} \mathbf{z}$$

PC vs LMC

When the smoothing kernel is the Dirac $\delta$ function $G_{o, q}^i(\mathbf{x} - \mathbf{z}) = a_{o, q}^i \delta(\mathbf{x} - \mathbf{z})$, PC degenerates to LMC. So:

• Degenerates to ICM: $G_{o, q}^i(\mathbf{x} - \mathbf{z}) = a_{o, q}^i \delta(\mathbf{x} - \mathbf{z}), \quad Q = 1, \quad R_q > 1$
• Degenerates to SLFM: $G_{o, q}^i(\mathbf{x} - \mathbf{z}) = a_{o, q}^i \delta(\mathbf{x} - \mathbf{z}), \quad Q > 1, \quad R_q = 1$
• Degenerates to LMC: $G_{o, q}^i(\mathbf{x} - \mathbf{z}) = a_{o, q}^i \delta(\mathbf{x} - \mathbf{z}), \quad Q > 1, \quad R_q > 1$
• Degenerates to Dependent GPs: $k_1\left(\mathbf{z}, \mathbf{z}^{\prime}\right) = \sigma^2 \delta\left(\mathbf{z}, \mathbf{z}^{\prime}\right), \quad Q = 1, \quad R_q = 1$

6. Multi-task Bayesian Optimization

Once we have fitted multiple tasks with multi-task GPR, multi-task Bayesian optimization mainly addresses how to choose the next probe point — i.e. it defines an acquisition function for the multi-task setting. The paper Swersky, Snoek & Adams 2013 — Multi-Task Bayesian Optimization extends Entropy Search to the multi-task setting. The paper Letham et al. (2019) extends Expected Improvement to the multi-task setting; its implementation is open-sourced in facebook/ax.