Pith. sign in

REVIEW 3 major objections 5 minor 39 references

Derivation of Output Correlation Inferences for Multi-Output (aka Multi-Task) Gaussian Process

T0 review · 3 major / 5 minor · reviewed 2026-08-10 · deepseek-v4-flash

Challenge this review Re-run · record.json Download PDF Read on arXiv ↗

Pith's one-line read A full step-by-step derivation shows how to build and differentiate the multi-task Gaussian process, giving practitioners a way to check their implementations.

desk verdict A well-organized tutorial re-derivation of MTGP, but the central gradient formula in Eq. (16) has a sign error and the Monte Carlo EM objective in Step 2 is misspecified; both need fixing before this can be a reliable reference. read the letter →

arxiv 2501.07964 v4 pith:HISQE6LE submitted 2025-01-14 cs.LG cs.AIstat.ML

Shuhei Watanabe This is my paper · ORCID
classification cs.LGcs.AIstat.ML
keywords multi-taskGaussianprocessmulti-outputKroneckerstructureEMalgorithmmarginallikelihoodgradientsBayesianoptimizationcoregionalization
verification ladder T0 review T1 audit T2 compute T3 formal

The pith

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The reading

This paper aims to make the multi-task Gaussian process (MTGP) derivable from scratch. Its central claim is that the EM updates and gradient-based hyperparameter optimization used in the standard MTGP model can be written down explicitly, step by step, with all intermediate algebra shown. If correct, practitioners can verify their MTGP code against these formulas instead of trusting an opaque implementation. The value is practical: multi-task GP is a common tool in Bayesian optimization, and the paper's derivations remove a barrier to using it correctly.

What carries the argument

The central object is the Kronecker-structured covariance $K_{xf}=K_x\otimes K_f$ combined with the noise term $S=I_N\otimes\Sigma$. The argument is carried by standard linear-algebra identities: $(A\otimes B)^{-1}=A^{-1}\otimes B^{-1}$, $|A\otimes B|=|A|^M|B|^N$, and the trace identity $\operatorname{vec}(C^\top)^\top(A\otimes B)\operatorname{vec}(C^\top)=\operatorname{tr}(A C B^\top C^\top)$. These convert high-dimensional vectorized expressions into matrix products of size $N\times N$ and $M\times M$, which is what makes both the EM updates and the gradient formulas derivable and computationally efficient.

What would settle it

Take a small synthetic two-output, two-input MTGP with chosen hyperparameters, compute the marginal log-likelihood of Eq. (13) directly, and compare its finite-difference gradient to the closed-form gradient formulas in Eqs. (16)–(17); any mismatch beyond numerical tolerance would disprove the derivation.

Watch

Extended reading notes

Core claim

On Bonilla et al.'s MTGP model, a zero-mean GP prior over latent functions with covariance $K_x \otimes K_f$ and observation noise $I_N\otimes\Sigma$, the complete-data log-likelihood can be transformed using Kronecker identities and trace identities into a tractable objective. The paper derives the EM updates $\hat{\sigma}_i^2 = \frac{1}{N}\sum_n (y_{n,i}-f_{n,i})^2$, $\hat{K}_f = F^\top K_x^{-1}F/N$, and the reduced objective for $\boldsymbol{\theta}$ as $M\log|K_x| + N\log|F^\top K_x^{-1}F|$. It likewise derives the marginal log-likelihood gradients for the gradient approach, including a Cholesky parameterization of $K_f$ and a projection matrix $P$ that extends the same formulas to settings where some outputs are missing at some inputs, preserving $O(N^3+M^3)$ complexity with iterative solvers.

Load-bearing premise

The paper's formulas hold for the specific MTGP model where the prior covariance between any two outputs is separable as a product of an input kernel and a fixed output kernel; if output correlations change across the input space, or are otherwise not captured by this product structure, the derived updates and gradients no longer apply.

Editorial extensions

If this is right

  • The derived formulas give practitioners a concrete reference for checking MTGP implementations, particularly the EM updates and the marginal-likelihood gradients.
  • Because the gradient approach can be written with explicit derivatives, automatic differentiation packages can be used directly, and the closed-form derivatives document what the autodiff is computing.
  • The partial-observation extension via the projection matrix $P$ covers realistic datasets where some outputs are missing at some inputs, without changing the asymptotic complexity when combined with iterative linear solvers.
  • The diagonalization of the Kronecker structure recovers the $O(N^3+M^3)$ training cost for the full-observation case, the regime on which practical MTGP implementations rely.

Reading between the lines

Editorial extensions of the paper, not claims the author makes directly.

  • The derivations suggest that any implementation of MTGP that fails to reproduce these closed-form updates on a small synthetic dataset very likely has a coding error, so checking against the formulas is a fast smoke test.
  • The same trace-and-Kronecker machinery presumably extends to other separable multi-output models, such as linear coregionalization models, and the projection-matrix trick may carry over to those settings as a way to handle missing outputs.
  • The paper's reliance on the separable prior $K_x\otimes K_f$ implies that for data where output correlations change with input location, the formulas should not be expected to hold, and a non-separable prior would need its own derivation.
Share X Bluesky LinkedIn Reddit HN

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, and a circularity audit.

Referee Report

3 major / 5 minor

Summary. The paper presents detailed derivations of the multi-output (multi-task) Gaussian process model of Bonilla et al. (2007), covering the EM algorithm for the output-correlation kernel and the gradient-based maximization of the marginal log-likelihood. It also discusses Kronecker-structured efficient computations, partial-observation extensions, and practical considerations for Bayesian optimization. The stated purpose is to provide friendly, verifiable derivations of the MTGP formulations and their gradients.

Significance. If the derivations were correct, the paper would be a useful reference for practitioners implementing MTGP, particularly because it fills derivation gaps in the original treatment and connects to modern automatic differentiation and scalable linear algebra. I credit the paper for including explicit proofs of the linear algebra identities and for discussing partial-observation extensions. However, the value of the paper is entirely contingent on mathematical accuracy, and the presence of several incorrect key formulas currently undermines that value.

major comments (3)
  1. [Section 4.2.1, Eq. (16)] The derivative of the marginal log-likelihood is stated as ∂L_mll/∂z = -(1/2)(tr(QG) + y^T Q G Q y). Re-deriving from Eq. (13) using d(A^{-1})/dz = -A^{-1} A' A^{-1} gives ∂L_mll/∂z = -(1/2) tr(QG) + (1/2) y^T Q G Q y. The sign of the data-fit term is wrong. Because Eq. (16) is the main deliverable of the gradient derivation in Section 4.2, a practitioner following the paper would move the data-fit gradient in the wrong direction during hyperparameter optimization.
  2. [Section 4.2.1, Eq. (17)] The listed derivatives of K_xf + S are incorrect. Since K_xf = K_x ⊗ K_f, the correct derivatives are ∂(K_xf+S)/∂θ_i = (∂K_x/∂θ_i) ⊗ K_f and ∂(K_xf+S)/∂L_{i,j} = K_x ⊗ (E_{i,j} L^T + L E_{j,i}). The paper instead writes ∂K_x/∂θ_i ⊗ I_M and I_N ⊗ (E_{i,j} L^T + L E_{j,i}), which are not equivalent in general and are only valid under a special structure. Since Eq. (22) for the partial-observation case relies on these formulas, the error propagates to that extension.
  3. [Section 4.1.3, Step 2] The Monte Carlo EM objective is written as N log( (1/K) Σ_k |(F~(k))^T K_x^{-1} F~(k)| ), but the correct Monte Carlo estimate of the expected complete-data log-likelihood term is N (1/K) Σ_k log|(F~(k))^T K_x^{-1} F~(k)|. The logarithm must be inside the average; the stated form is the log of an average determinant, which is not an unbiased estimate of the expected log-determinant and will bias the EM update.
minor comments (5)
  1. [Section 4.1.2, Eq. (7)] The intermediate expression '-N/M Σ_m ∂log σ_m/∂σ_i' appears to be a typo; the final result agrees with the correct derivative -N/σ_i, but the displayed formula should be corrected to avoid confusion.
  2. [Eq. (11)] The line 'log|F^T K_x^{-1} F| |N^{-M}|' is confusing; since |F^T K_x^{-1} F / N| = |F^T K_x^{-1} F| N^{-M}, the constant should be written additively as -M log N. The final expression is unaffected because constants are discarded, but the presentation needs fixing.
  3. [Section 4.1.2] The phrase 'By taking the derivative of zero' should read 'By setting the derivative to zero'.
  4. [Notation] The symbol Q is used for K_f^{-1} in Section 4.1.2 and for (K_xf+S)^{-1} in Section 4.2.1; this overloading is confusing and should be replaced with distinct symbols.
  5. [Section 5] Several self-citations to the author's TPE papers appear tangential to the derivation content, and Watanabe et al. (2022) and (2023) share the same arXiv identifier; the references should be cleaned up.

Circularity Check

0 steps flagged · score 0.0 of 10

No circularity: the MTGP derivations are self-contained algebra from the stated model; the author's self-citations are contextual and not load-bearing.

full rationale

The paper is a self-contained derivation reference rather than an empirical prediction paper. Its starting point is the explicitly stated MTGP model from Bonilla et al. (2007): the prior N(f_{1:N}|0, K_xf) and the likelihood N(y|f, S), given in Section 4 and Eq. (3). Every subsequent formula is obtained by ordinary linear algebra and matrix calculus: the complete-data log-likelihood in Eqs. (3)-(6), the EM stationary-point updates in Eqs. (7)-(12), the marginal log-likelihood in Eq. (13), its derivative in Eqs. (14)-(16), and the Kronecker diagonalization in Eqs. (25)-(26). None of these steps assumes the formula being derived, and no parameter is fitted to a subset of data and then renamed as a prediction. The author's self-citations (Watanabe 2023; Watanabe and Hutter 2022, 2023; Watanabe et al. 2022, 2023) appear only in footnotes and in the practical-considerations section as pointers to alternative TPE-based Bayesian optimization methods; they are not used as evidence for any MTGP identity and therefore are not load-bearing. Any possible algebraic errors, such as a sign issue in Eq. (16) or the Monte Carlo EM averaging in Section 4.1.3 Step 2, are correctness concerns rather than circularity: the derivation chain would be wrong, not self-referential. Accordingly, no circular step can be exhibited from the paper's own equations, and the fair circularity score is 0.

Assumptions & free parameters 0 free parameters · 5 assumptions · 0 invented entities

The central derivations rely only on standard GP assumptions and elementary linear algebra facts. No free parameters are introduced and no new entities are postulated. The paper's mistakes are algebraic errors, not hidden assumptions.

assumptions (5)
  • domain assumption Kernel matrices K_x and K_f are positive definite and symmetric, hence invertible.
    Stated at the end of Section 3 (Notations) and used for inverses and log-determinants throughout Section 4.1 and 4.2.
  • domain assumption The joint prior over latent functions is N(f_{1:N} | 0, K_x⊗K_f) and the likelihood is N(y | f, I_N⊗Σ).
    This is the MTGP model from Bonilla et al. (2007), introduced in Section 4. The derivations of EM and gradients assume this model.
  • standard math Kronecker product identities: |A⊗B|=|A|^M|B|^N, (A⊗B)^{-1}=A^{-1}⊗B^{-1}, and (A⊗B)(C⊗D)=(AC)⊗(BD).
    Used in Eqs. (4), (5), and (25). These are standard results cited in Appendix A.1.
  • standard math Matrix calculus identities: d log|A|/dz = tr(A^{-1} dA/dz), d(A^{-1})/dz = -A^{-1} dA/dz A^{-1}.
    Used in Eqs. (14) and (15). Standard and stated in Appendix A.1.
  • standard math The diagonalizations K_x = U_x Λ_x U_x^T and Σ^{-1/2} K_f Σ^{-1/2} = V_f Λ_f V_f^T exist with unitary U_x, V_f.
    Used in Eq. (25). Guaranteed for symmetric matrices; requires Σ to be full rank so Σ^{-1/2} is defined.

how reviews work

0 comments
Cite this review

Pith. "Pith review of Derivation of Output Correlation Inferences for Multi-Output (aka Multi-Task) Gaussian Process." pith.science (2026). https://pith.science/paper/HISQE6LE

@misc{pith2026250107964,
  author       = {Pith},
  title        = {Pith review of: Derivation of Output Correlation Inferences for Multi-Output (aka Multi-Task) Gaussian Process},
  year         = {2026},
  howpublished = {\url{https://pith.science/paper/HISQE6LE}},
  note         = {Machine review of arXiv:2501.07964}
}
read the original abstract

Gaussian process (GP) is arguably one of the most widely used machine learning algorithms in practice. One of its prominent applications is Bayesian optimization (BO). Although the vanilla GP itself is already a powerful tool for BO, it is often beneficial to be able to consider the dependencies of multiple outputs. To do so, Multi-task GP (MTGP) is formulated, but it is not trivial to fully understand the derivations of its formulations and their gradients from the previous literature. This paper serves friendly derivations of the MTGP formulations and their gradients.

Discussion (0). Continue with ORCID to comment.

Reference graph

Works this paper leans on

39 extracted references · 31 canonical work pages

  1. [1]

    and Lawrence, N

    Alvarez, M. and Lawrence, N. (2008). Sparse convolved G aussian processes for multi-output regression. Advances in Neural Information Processing Systems

    work page 2008

  2. [2]

    Balandat, M., Karrer, B., Jiang, D., Daulton, S., Letham, B., Wilson, A., and Bakshy, E. (2020). BoTorch: A Framework for Efficient Monte-Carlo Bayesian Optimization . In Advances in Neural Information Processing Systems

    work page 2020

  3. [3]

    Bergstra, J., Bardenet, R., Bengio, Y., and K \'e gl, B. (2011). Algorithms for hyper-parameter optimization. Advances in Neural Information Processing Systems

    work page 2011

  4. [4]

    Bonilla, E., Chai, K., and Williams, C. (2007). Multi-task G aussian process prediction. Advances in neural information processing systems

    work page 2007

  5. [5]

    Brochu, E., Cora, V., and de Freitas, N. (2010). A tutorial on Bayesian optimization of expensive cost functions, with application to active user modeling and hierarchical reinforcement learning. arXiv:1012.2599

    arXiv 2010

  6. [6]

    Daulton, S., Balandat, M., and Bakshy, E. (2020). Differentiable expected hypervolume improvement for parallel multi-objective B ayesian optimization. Advances in Neural Information Processing Systems

    work page 2020

  7. [7]

    Daulton, S., Eriksson, D., Balandat, M., and Bakshy, E. (2022). Multi-objective B ayesian optimization over high-dimensional search spaces. In Uncertainty in Artificial Intelligence

    work page 2022

  8. [8]

    Dong, K., Eriksson, D., Nickisch, H., Bindel, D., and Wilson, A. (2017). Scalable log determinants for G aussian process kernel learning. Advances in Neural Information Processing Systems

    work page 2017

Show all 39 references
  1. [9]

    and Poloczek, M

    Eriksson, D. and Poloczek, M. (2021). Scalable constrained Bayesian optimization. In International Conference on Artificial Intelligence and Statistics

    2021

  2. [10]

    Falkner, S., Klein, A., and Hutter, F. (2018). BOHB : Robust and efficient hyperparameter optimization at scale. In International Conference on Machine Learning

    2018

  3. [11]

    Feurer, M., Letham, B., Hutter, F., and Bakshy, E. (2018). Practical transfer learning for Bayesian optimization. arXiv:1802.02219

    2018 arXiv

  4. [12]

    Gardner, J., Kusner, M., Xu, Z., Weinberger, K., and Cunningham, J. (2014). Bayesian optimization with inequality constraints. In International Conference on Machine Learning

    2014

  5. [13]

    Garnett, R. (2022). Bayesian Optimization . Cambridge University Press

    2022

  6. [14]

    Gelbart, M., Snoek, J., and Adams, R. (2014). Bayesian optimization with unknown constraints. arXiv:1403.5607

    2014 arXiv

  7. [15]

    Kandasamy, K., Dasarathy, G., Oliva, J., Schneider, J., and Poczos, B. (2019). Multi-fidelity G aussian process bandit optimisation. Journal of Artificial Intelligence Research , 66

    2019

  8. [16]

    Kandasamy, K., Dasarathy, G., Schneider, J., and P \'o czos, B. (2017). Multi-fidelity B ayesian optimisation with continuous approximations. In International Conference on Machine Learning

    2017

  9. [17]

    Lin, J., Ament, S., Balandat, M., and Bakshy, E. (2024). Scaling G aussian processes for learning curve prediction via latent K ronecker structure. arXiv:2410.09239

    2024 arXiv

  10. [18]

    Lindauer, M., Eggensperger, K., Feurer, M., Biedenkapp, A., Deng, D., Benjamins, C., Ruhkopf, T., Sass, R., and Hutter, F. (2022). SMAC3 : A versatile B ayesian optimization package for hyperparameter optimization. Journal of Machine Learning Research , 23

    2022

  11. [19]

    Maddox, W., Balandat, M., Wilson, A., and Bakshy, E. (2021). B ayesian optimization with high-dimensional outputs. Advances in neural information processing systems

    2021

  12. [20]

    and Krishnan, T

    McLachlan, J. and Krishnan, T. (2008). The EM algorithm and extensions . John Wiley & Sons

    2008

  13. [21]

    Nguyen, T., Bonilla, E., et al. (2014). Collaborative multi-output G aussian processes. In Uncertainty in Artificial Intelligence

    2014

  14. [22]

    Ozaki, Y., Tanigaki, Y., Watanabe, S., Nomura, M., and Onishi, M. (2022). Multiobjective tree-structured P arzen estimator. Journal of Artificial Intelligence Research , 73

    2022

  15. [23]

    Ozaki, Y., Tanigaki, Y., Watanabe, S., and Onishi, M. (2020). Multiobjective tree-structured Parzen estimator for computationally expensive optimization problems. In Genetic and Evolutionary Computation Conference

    2020

  16. [24]

    Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J., Chanan, G., Killeen, T., Lin, Z., Gimelshein, N., Antiga, L., et al. (2019). PyTorch : An imperative style, high-performance deep learning library. Advances in neural information processing systems

    2019

  17. [25]

    Shahriari, B., Swersky, K., Wang, Z., Adams, R., and de Freitas, N. (2016). Taking the human out of the loop: A review of B ayesian optimization. Proceedings of the IEEE , 104

    2016

  18. [26]

    Song, J., Chen, Y., and Yue, Y. (2019). A general framework for multi-fidelity B ayesian optimization with G aussian processes. In International Conference on Artificial Intelligence and Statistics

    2019

  19. [27]

    Stegle, O., Lippert, C., Mooij, J., Lawrence, N., and Borgwardt, K. (2011). Efficient inference in matrix-variate G aussian models with iid observation noise. Advances in Neural Information Processing Systems

    2011

  20. [28]

    Swersky, K., Snoek, J., and Adams, R. (2013). Multi-task Bayesian optimization. In Advances in Neural Information Processing Systems

    2013

  21. [29]

    and Bishop, C

    Tipping, M. and Bishop, C. (1999). Probabilistic principal component analysis. Journal of the Royal Statistical Society Series B: Statistical Methodology , 61

    1999

  22. [30]

    Titsias, M. (2009). Variational learning of inducing variables in sparse G aussian processes. In Artificial Intelligence and Statistics

    2009

  23. [31]

    Wackernagel, H. (2003). Multivariate geostatistics: an introduction with applications . Springer Science & Business Media

    2003

  24. [32]

    Watanabe, S. (2023). Tree-structured P arzen estimator: Understanding its algorithm components and their roles for better empirical performance. arXiv:2304.11127

    2023 arXiv

  25. [34]

    Watanabe, S., Awad, N., Onishi, M., and Hutter, F. (2023). Speeding up multi-objective hyperparameter optimization by task similarity-based meta-learning for the tree-structured P arzen estimator. arXiv:2212.06751

    2023 arXiv

  26. [36]

    and Hutter, F

    Watanabe, S. and Hutter, F. (2023). c- TPE : Tree-structured P arzen estimator with inequality constraints for expensive hyperparameter optimization. arXiv:2211.14411

    2023 arXiv

  27. [37]

    and Seeger, M

    Williams, C. and Seeger, M. (2000). Using the N ystr \"o m method to speed up kernel machines. Advances in Neural Information Processing Systems

    2000

  28. [38]

    Wilson, A., Knowles, D., and Ghahramani, Z. (2011). G aussian process regression networks. arXiv:1110.4411

    2011 arXiv

  29. [39]

    Wilson, J., Borovitskiy, V., Terenin, A., Mostowsky, P., and Deisenroth, M. (2020). Efficiently sampling functions from G aussian process posteriors. In International Conference on Machine Learning

    2020

  30. [40]

    Wistuba, M., Kadra, A., and Grabocka, J. (2022). Supervising the multi-fidelity race of hyperparameter configurations. Advances in Neural Information Processing Systems

    2022

  31. [41]

    Yang, K., Emmerich, M., Deutz, A., and B \"a ck, T. (2019). Multi-objective B ayesian global optimization using expected hypervolume improvement gradient. Swarm and evolutionary computation , 44

    2019

Pith tools

Reviewed August 10, 2026 · model on record in the stance chip above.