Graphical Multioutput Gaussian Process with Attention
Authors: Yijue Dai
Wenzhong Yan
Feng Yin
Venue: ICLR 2024 spotlight (2024)
Keywords: Gaussian process regression
Multioutput Gaussian process
Attention mechanism
Primary area: probabilistic methods (Bayesian methods, variational inference, sampling, UQ, etc.)

Abstract: Integrating information while recognizing dependence from multiple data sources and enhancing the predictive performance of the multi-output regression are challenging tasks. Multioutput Gaussian Process (MOGP) methods offer outstanding solutions with tractable predictions and uncertainty quantification. However, their practical applications are hindered by high computational complexity and storage demand. Additionally, there exist model mismatches in existing MOGP models when dealing with non-Gaussian data. To improve the model representation ability in terms of flexibility, optimality, and scalability, this paper introduces a novel multi-output regression framework, termed Graphical MOGP (GMOGP), which is empowered by: (i) Generating flexible Gaussian process priors consolidated from dentified parents, (ii) providing dependent processes with attention-based graphical representations, and (iii) achieving Pareto optimal solutions of kernel hyperparameters via a distributed learning framework. Numerical results confirm that the proposed GMOGP significantly outperforms state-of-the-art MOGP alternatives in predictive performance, as well as in time and memory efficiency, across various synthetic and real datasets.

Published as a conference paper at ICLR 2024
GRAPHICAL MULTIOUTPUT GAUSSIAN PROCESS WITH ATTENTION
Yijue Dai, Wenzhong Yan & Feng Yin" School of Science and Engineering The Chinese University of Hong Kong, Shenzhen, China y yinfengecuhk. edw..eduuccoo cn, yijuedai @link. aeu.hk.edu.cm
ABSTRACT
Integrating information while recognizing dependence from multiple data sources and enhancing the predictive performance of the multi-output regression are chal- lenging tasks. Multioutput Gaussian Process (MOGP) methods offer outstanding solutions with tractable predictions and uncertainty quantification. However, their practical applications are hindered by high computational complexity and storage demand. Additionally, there exist model mismatches existing MOGP models when dealing with non- Gaussian data. To improve the model representation abil- ity in terms of flexibility, optimality, and scalability, this paper introduces a novel multi-output regression framework, termed Graphical MOGP (GMOGP), which is empowered by: (i) Generating flexible Gaussian process priors consolidated from identified parents, (ii) providing dependent processes with attention-baut based graphical representations, and (iii) achieving Pareto optimal solutions of kernel hyperparameters via a distributed learning framework. Numerical results confirm that the proposed GMOGP significantly outperforms state-of-the-ar MOGP al- ternatives in predictive performance, as well as in time and memory efficiency. across various synthetic and real datasets.
INTRODUCTION
For inference on functions, the Gaussian Process (GP) serves as the most prominent Bayesian non- parametric model, offering not only function values but also predictive distributions. Capturing the input-output relationship with a multivariate Gaussian at finite locations grants tractable pre- dictions and predictive uncertainty (Rasmussen & Williams 2006). The scenarios where resource constraints impede the collection of ample data samples motivate the development of Multioutput GP (MOGP) models, as the MOGP can correlate data from multiple diverse sources and obtain better predictive accuracy than isolated models. Practical applications of the existing MOGP meth- ods are widespread, such as time-series analysis (Dürichen al.) 2014). water resource forecasting Pastrana-Cortés et al. 2023), and group-structured data clustering (Leroy et al. 2023).
The popular MOGP models are mainly rooted in the Linear Model Coregionalization (LMC) pro- posed by Goulard & Voltz ( (1992). Modeling all outputs as the linear combinations of shared latent GPs, the LMC correlates the diverse data by covariance matrix. There also exist other MOGP models, such (I) transformation-bur ing the MOGP into a sequence of single output GPs (Requeima et al. 2019); (2) discrepancy- models, which focus transferring information from low-fidelity outputs to progressively enhance the performance of high-fidelity target outputs (Perdikaris et al. [217 Requeima et al., 2019). How- ever, the challenges of providing understandable output correlations and a probable output hierarchy have impeded their advancement. To date, improved learning performance of the LMC is achieved by deriving representative covariance functions (Bonilla et al. 2007 Wilson et al.) 2012) Dai al.) 2020 Chen et al., (022). by adapting common mean process [eroy et (012)) and by introducing neural embeddings or efficient frameworks (Liu al. 2022 Chung & Al Kontar 2023). However, the predictive performance of the LMC- models may be inferior to that of isolated GPs when sufficient training data exists for each output (Bruinsma et 2020). See the results in Section 5
models that treat outputs augmented inputs, transform-
"The corresponding author Feng Yin. Our code is available http:///tttuh..cm//lllaaaaaaaa

---

Published as a conference paper at ICLR 2024
In this paper, we propose a novel multi-output regression framework, termed graphical MOGP (GMOGP). Unlike the LMC, which interprets the output correlations by covariance based on a set of shared GPs, the GMOGP, built upon probability graphical model, can directly learn the condi- tional dependence and imply graphical representations of the multiple outputs. Moreover, learning the kernel hyperparameters of the GMOGP can be formulated as multi-objective optimization (MOO) problem with Pareto optimal solutions. Further extensions by non-linear transformations make the transformed GMOGP (TGMOGP) available fit non- Gan-Gan data. Enhanced predictive accuracy and efficiency validate the efficacy of the GMOGP across various synthetic/real datasets.
2
BACKGROUND
This section concludes basic concepts and notations in Gaussian process regression (GPR), serving as the bedrock of the upcoming GMOGP framework detailed in Section 3
2.1
GAUSSIAN PROCESS REGRESSION
GP characterizes a distribution over functions fully by mean function m(x) and a kernel function k(x, x'; 0), i.e.,
f(x) ~ GP (m(x), k(x, x);
(1)
The commonly used squared exponential (SE) kernel, lis (x, x'; oj exp( |x
x'I
is parametrized by length-scale and a signal variance of (Rasmussen & Williams 2006). Given a dataset D : {X.y} = {x.: yn | N n= the GPR model can be described as:
f(xx). + cn. n =
N,
(2)
where X, E Rd denotes an input location, and the noise en is assumed i.idd with e, ~ NN d²). In particular, the scalar-valued observations {ym, E K] direct single-output GP (SOGP) models.
2.1.1 MULTIOUTPUT GAUSSIAN PROCESS
For multi-output regression, an MOGP can be derived as:
f(x) GP(mm (x) 0ar)),
where the vector-valuee function f(x) (x), (x),.. (x)((s)))) is determined by the mean function mas (x) E SR and the matrix-valuus kernel function KM (x, x') E SRS. Note that the hyperparameters 0 are searched in a vector-valued reproducing kernel Hilbert space (RKHS).
At finite locations X E Nd the prior and the likelihood of the MOGP has the following forms:
(Prior): p(f(X)) (m (X), KM XX
X; 0mx))
(Likelihood):
PY X, = -N (((())) @ In).
(5)
where (1) ..... y()) (S) E R. and (+) (1 (6) ER are the labels of the 8th output.
n The noise terms are assumed to have zero mean and a covariance matrix n diag , oz).
The kernel function for the MOGP prior evaluated at X is derived follows:
(x;
XX,
KM (x,
:
:
R SNNSN
(6)
(X,
Xss 0.0)
Generally, we train the MOGP models by minimizing their negative log marginal likelihood (NLL):
[( x T (X, (x - Σ) + log K (X, + E|,
(7)
with Y mm (X), and Σ = Σ @ IN. Integrating over the parameter space differentiates the marginal likelihood methods from other non -By-esysi counterparts. The automatic incorporation of a trade-off between model fit and model complexity renders the marginal likelihood a valuable metric for model selection (Schölkopf et al. (2002) Micchelli & Pontil| (005). Conditioning the joint Gaussian prior the observations, the predictive distribution for test imput x, turns out to be:
pff. |x., X, Y. Ass)
(8)
2

---

Published as a conference paper at ICLR 2024
with (omitting kernel hyperparameters)
{ (
f. = (x.,,xx.))) (X,
Σ)
(9)
(10)
In order to infer the underlying functions and correlate the multiple outputs, the LMC models tailor distinct coefficients to each output via Q shared independent GPs. The elements of their kernel functions (Eq. 6)) are formulated as (x,))) = > aigarighig x'), where i, 1 E I=
(1.2...., S). and different modes of the coefficients Aiggisg correspond to varied LMC variants (Alvarez et al., 2012: Liu et al.) 2018). Although the LMC can generate dependent processes and measure output correlation, their applications computational and storage burdens associated with the Ain-dimensional correlation matrix, and (ii) model mismatch/inflexiblliyy when the underlying likelihood (Eq. (5) deviates from Gaussian distributions. Given these challenges, an alternative framework that can leverage the advantages of MOGPs, generate flexible prior/predictive distributions, and freely accommodate various efficient learning methods considered Section 3
narrowed by two significant impediments: (i) The
GRAPHICAL MULTIOUTPUT GAUSSIAN PROCESS
Ubiquitously, there exists interplay among the multiple outputs. giving rise varied research, e.g.. graph construction, transfer learning, and multi-task learning (Qiao et al.) 2018 Vandenhende et al.) 2021 For instance, the dependence (denoted by arrows) between a larget output (marked yellow) and the other six outputs can be described by graphical representation shown in Figure ((a) In the multi-output regression, the strong bond between the trend of the target CO2 and Temperature has been well discovered (solid link), while the other bonds are vague (dashed link) and worth exploring.
Art
Temperature
Education
CO2
174
cover
GDP
Vehicle
(a)
(b)
Figure 1:
Illustrations of: (a) graphical representation among the target CO2 and the other outputs, (b) an example of a directed acyelic graph (the cross means that there are no directed cycles).
As directed graphical models (a.k.a. Bayesian network) can express the conditional dependence structure among random variables and make probabilistic statements for broad class of distribu- tions based specific graph (Koller & Friedman, 2009), propose graphical MOGP model to learn the output dependence jointly with the regression task.
Mathematically, graph is built with nodes and links. In the directed graphical model, each node can represent random variables, while the links convey probabilistic dependence and have directionality indicated by arrows. For example, given nodes (U1, U2,..., ... 06) and graph G, shown in Figure [(b) the joint distribution can be decomposed by repeatedly applying the product rule of probability:
p(11, 02, U6) p(U6 U1, U2, Us, 14, Us (aspp(ss/o (1, U2, u3, v4) ... (...-pp(2ll)nn
(11)
G. p(vo W1, U2 v3, Us (t))((ssvvr /12, ())p(r)/u, vg) (2)))(2())()(()))))
(12)
The Eq. 12 is derived according to the conditional independence, e.g., p(ue Vi, 02 Us, U4 Us) p(ua P1, U2, Us Us), as there no link from node 24 to UG the graph G.. In more general terms, the joint distribution of multiple nodes v = fuieer defined by a graph can be expressed as:
=
S
p(v)= II p(u) [-1
(13)
where pa, denotes the parents of G; e.g., pag {ui, 02, va } in Eq. 12). Note that the joint distribution only admits directed acyclic graphs; see more details in (Bishop & Nasrabadi 2006).
If there exists link going from node i node j, the node regarded the parent of the node j.
3

---

Published as a conference paper at ICLR 2024
Accordingly, for S outputs, model each output SOGP, and generate multivariate Gaussian random variables evaluated X (represented by node 1& j E Z). Following the result in Eq. 13 the joint distribution defined by the specific graph structure, with the target node f (representing the (O2) connected to the heads of arrows (see Figure [(a), can be derived as follows:
p(fx fX
fx ( p pa.) II p(52')
(14)
SEPa,
Since the links are unknown in practice, the parents of the target node are assumed to contain all
other nodes (except the target one) at first, i.c., pa; = { f? 1.... xx......... 3/L J. and can be adjusted via attention mechanism detailed in Section 3.2 Alternately, we can treat each output the other
as the target node (represented by ( i Σ), and learn the conditional dependence nodes. The diagrammatic illustration of the output bonds and the graph with respect to different target at the first stage (all vague) is given in Figure 2 Moreover, conditioning on the target node, the parents are dependent (see Remark 3.1. which paves the way to generate dependent processes.
Illustration of the bonds among six outputs and the dependence structure between the different target (marked in yellow) and its parents (marked in purple) at the first stage.
Figure 2:
Remark 3.1 (Conditional Dependence). Defined by the specific graph with joint distribution form: $2 p( (())
((k)) p(fx
(()) (
)((x)) >p(xx) (R)) if conditioning on the target fx) the condi-
tional distribution p(S2) 1xx 15x = p((X () fx Dppppp()1)) (ou2s11k).. (. 1 s/(()))), yielding the par-
ents p and (()) are dependent, ie., 1 M. (00 15x
3.1 GRAPHICAL MULTIOUTPUT GAUSSIAN PROCESS PRIOR
Bayesian methods are invaluable within the machine learning community due to the ability to in- corporate prior knowledge. Appropriate prior beliefs represent detailed inductive bias, improve interpretability, and prevent overfitting (Theodoriis) [015] Lotfi et al. (022). the GMOGP, we model the conditional distribution in (4.) as Gaussian with aggregated information, namely,
p x pa. (pa)paa)-00 -N ( (1
Σ
a + m,, ko, (X, i
(15)
jEPa
where at, E R and m, are bias parameters added into the mean, and ke, (X, XX) NxN is covariance function. Conditioning on the states of its parents, each target node is of the form:
JijiX + m, + thi
(16)
jEPa
with m. and tb; (N, ko. (X) characterizing the 5th output. Given the parents with p( p()2))
N(mj, ke (X, X)). can derive the GMOGP prior via (E((() Σjem, (n)B(r()) m, and cov(fl)
(uX) ELKDD {E,epa an - EISKD) + w/ (see Proposition 3.1
Proposition 3.1 (GMOGP Prior). In the multi-outpul regression, defined by the specific graph with a target node and given its parents, the GMOGP prior for each target node (i E followss p(f N EjEm, m, m;, Ejem, ko, (X ke, (x,x)) The proof is given in Appendix A.1
Following the noise setting in Section
2.1.1 the posterior of the 7th output is proportional to the product of the prior and the likelihood, with hyperparameters θ (0,bie and a (aep)
(0) .5 O O, ov) -. ax)pp (0) oz)
(17)
4

---

Published as a conference paper at ICLR 2024
Enhanced representation ability provided by the graph structure in the input space has been demon- strated in recent graph GP methods, such as Ng et al. (2018) and Fang et al. (202I). in which the graph structures among the imput features structure is unclear and intriguing. To learn the graph while improving the performance of the multi- output regression,
well discovered. For the multiple outputs, the graph
propose to incorporate the attention mechanism into the GMOGP to modify the parents, as the attention-base based graph construction methods have shown impressive performance, such /Vaswani al. 2017 Veličković et al., 2018 Kim & Oh 2022
3.2 LEARNING PARENTS WITH ATTENTION
Centering on the CO2 level (Figure I[a). the dependence on the number of vehicles, GDP level, and the degree of education is obscured. Intuitively, highly conditional dependent parents deserve large coefficients {at,) jepu, Without knowing the graph adjacency, such coefficients can be learned by using an attention mechanism J = exp(es,)/(L + Ejepe, exp(ei,)) and a scoring function:
e = Lea (IF ( Way + tas).
(18)
where the weights Wij and bias Cıj are learning parameters. The inner product-based scoring func- tion widely used, as seen in works such as Vaswani et al.] (2017) and Devlin et al. (2018), and
can reveal geometric constructions involving angles, lengths, and distances (Schölkopf et al. [2002).
An alternative scoring function employs the concatenation of inputs (Veličković et al., 2018 Brody et al., (022). However, existing scoring functions parametrized by common parameter matrices showed static attention - a fixed node consistently achieves the highest attention (Brody ct al. 2022).
In Eq. 18), we introduce modified scoring function, which can achieve dynamic attention with the aid of disparate learning parameters and non-linear function. The learned coefficients of various real datasets are listed in Section 5 and Appendix B In practice, values (y(0),')) at the initial learning phase. When oi.j R 0, the parents set pa, is adjusted by unlinking the node j (j E I.j * i). For classic MOGPs, the dependence between outputs is measured by symmetric covariance matrices. In comparison, our GMOGP can not only learn more flexible asymmetric dependence measure but also capture the covariance from cov(I') sK3
substitute the observation
Ejer,, oi. skk)) + Liy (X, where Iúj is the i, element of the identity matrix (see Appendix A.I). In addition, the GMOGP can handle heterotopic data (X + (2) A...., (S) * X) by introducing an extra weight matrix to Eq. 18) for aligning extracted input features.
3.3 MODEL LEARNING AND INFERENCE
In Section 2.I.I all kernel hyperparameters are learned through the marginal likelihood Lio,e Both compulational/storagg requirements and learning challenges, brought about by the high- dimensional covariance matrix and multiple vector- RKHSs, impede the practical implemen- tation of classic MOGPs Yang et al. 2020). For the GMOGP model, each target node i, E I has its own prior and likelihood (Eq. 17p). Accordingly, the type maximum likelihood can be conducted for every output, and the parameters corresponding to the 7th output, i.c., 7(0)) == {O, a,, J.). can be updated via minimizing:
(6) 900 (15(0) ( (Ke (x) 3(0) + log KG (X,x) }.
(19)
where g(0)= = y 1(1) (EjeR, am.m and (0) (X, x) ni (X, X) + kX, X). Note that the size of the correlation matrix k()) (X, X) is much smaller than the (X, in Eq.(6), and the hyperparameters θ are shared among the multiple outputs for information exchange. Instead of solving with a high-dimensional problem, the separate objective with the shared kernel hyperparameters can be modeled by a multi-objective optimization (MOO) problem, i.e., F(O) ((8) (O))"" Moreover, applying the weighted sum method with the objective function: Σ W: (()) (O), wi > 0, to solve the MOO problem provides sufficient condition for Pareto optimality] of the kernel hyperparameters (Marler & Arora 2010). In the GMOGP, the
(O), ((3)
2 A solution point is Pareto optimal if it is not possible to move from that point and improve at least one objective function without causing detriment to any other objective function (see details in Appendix A.5)
5

---

Published as a conference paper at ICLR 2024
weights be set equal, since there is priority among the outputs, and each marginal likelihood c() 27700, E I follows Gaussian that can be normalized so that no one dominates the objective.
To summarize, the proposed GMOGP empowers the vanilla MOGPs with: (i) The attention- asymmetric dependence measure (asj # aji), (ii) less computational and storage requirements, and (iii) Pareto optimal solutions of the kernel hyperparameters and extra graphical representations. Figure 3 presents the detailed workflow of the GMOGP hyperparameter optimization scheme. As separable model, each target node possesses distinct training samples and parents, undergoing an in- dependent inference procedure (O(N*)). It contributes individual information via 0; and aggregates dependent information through shared 8j jEpa, prompting the application of diverse distributed learning schemes. Key distinctions between the distributed GMOGP and the classic distributed GP (DGP) (Deisenroth & Ng. 2015) are summarized in Remark 3.2
Global update step for @
(ept "
(E6,, oD) (02"1,0"" " (er, on (9th) vict
(erll of (V (ar))maae
vellertt of
(ert) aag"))
af al (000)((a))len
pay
pay
pag
Inference
Prediction
(X.)
(.
Figure 3: Illustrating the GMOGP workflow, each target output transmits individual knowledge while aggregating optimized information from its parents, enabling separable training and inference.
Remark 3.2 (Distributed scheme). The distributed framework can learn all parameters jointly and does not require adjustments, such as (i) simplifying approximations, e.g., the block diagonal ap- proximation commonly employed in DGP approaches; (ūi exira prediction fusion methods, as the prediction x. can be separately inferred from k(@) (x., x)(kk)) (X, + ( (iii) data segmentation and (iv) additional client selection methods, as we can apply the identified parents.
3.4
NON-GAUSSIAN GMOGP PRIOR
The Gaussian priors on functions are widely used in Bayesian machine learning literature. benefiting from posterior consistency, tractable inference, and closed-form solutions (Rasmussen & Williams
2006 Whereas, the case of non-Gaussian prior/likelihood is often overlooked and ubiquitous in reality, e.g., heavy-tailed, bounded, and multimodal data (Snelson et al. 2003 Sendera et al. 2021
Typicalemmedie remedies include deep GP (DGP) (Damianou & Lawrence 2013) and transformed GP (TGP) RRios& & Tobar, 2019). The performance of the TGP, transformed by marginal flows, was tested in (Maroñas et al., 2021) and demonstrated competitiveness with the DGP while incurring lower computational costs. However, the naive application of marginal flows to MOGP models involves untraceable learning procedures and inefficient high-dimensional approximation approaches, as ID quadrature cannot be utilized. Lately, TGP model, tailored for multi-class classification, has been proposed. It can generate dependent processes with different marginal flows from a single sample of latent process (Maroñas & Hernández-Lobato 2023). In the context of the multi-output regression problem addressed in this paper, the single latent function neither sufficiently flexible nor accom- modating distinct prior beliefs. Also, the strong dependence imposed by the single GP sample can degrade regression performance, especially when dealing with low-quality or biased outputs.
In contrast, the GMOGP model is unconstrained in applying non-linear transformations (termed flow). Concretely, extra model flexibility can be extended by compositing K element-wise invert- ible transformations {Go, F + K-K-1 Correspondingly, each target node defined in Eq.(6
aggregates non-linear relations, i.e., G (/ (1() = GG EEepa, aigs (() + m: + +bbb) Applying

---

Published as a conference paper at ICLR 2024
the inverse function theorem and the change of variable formula iteratively (Rezende & Mohamed 2013), we can derive the transformed GMOGP prior as:
K-
OG det
P,0 oi (fR GG x) = -P-co (10))
II K=0
vk E 10. K --1 13. (20)
where f0x = f& and fitlx ( (fke) with parameters +; (80) 01 KK-1 such
((1)
(2)
as // (Sk, Pk. Au in Sinh-Archsinh (SAL) flow: G
Gg (1x)) Ca sinh(ps aresinh Aa).
We refer readers to Rios (2220) and Maroñas et al. (2021) for available flows and validity analysis.
Consolidating with the non-Gaussian prior, the transformed GMOGP (TGMOGP) can be easily implemented by adapting the off-the-shelf sparse variational inference algorithms in Hensman et al. 2015) or Maroñas et al. (2021) through optimizing the negative evidence lower bound (NELBO):
p( (3)
min (7 2..0 K
Σ
E (1622)
102
+ E 4 (u!"))
log q(u (0)
(21)
with inducing points un E M M << N and q(u))) = ((m KC The detailed derivation
of the NELBO and variational gap are elaborated in Appendix A.2 For inference, by the law of the
unconscious statistician (L.OTUS) the predictions of output i, E can be estimated from the first
moment of predictive distribution: ply=") > p(y GG. O ()((6 )ath (Snelson et al. 2003
In total, the dependencies of the variables in the GMOGP are concluded in the tollowing Figure
Invertible transformation
2/(1)
W (0
VJ
of
Exs
y!
Wy
Attention
Figure 4: Visualizing variable dependencies in the GMOGP. Each target output has their own par- ents, transformations, and samples with knowledge exchanged by kernel hyperparameters 0j, jem,,-
4 RELATED WORK
The graph-based GP models have been widely developed in local structured feature spaces and rela- tional learning discipline, resulting in the development of graph GP (Ng et al. 2018) and relational GP (Chu al., 2006 Silva et al. 2007). Without knowing neighbor graph structures or a com- plete Laplacian matrix, attention-bass shown be successful in supervised learning problems (Veličković al.. 2018). The idea of finding neighborhoods using non-negative kernel (NNK) regression (Shekkizhar & Ortegal (2020) explores geometric understanding and connections to the classic matching pursuit approach (Tropp & Gilbert) 2007). In the GMOGP framework, we can draw parallels with the kernel matching pursuit objec-
graph construction methods have sprung up and have been
tive Vinneet Bengio 2002). where y!!) 3k (X,, X) can be regarded as matching through constructed kernel dictionary and coefficient 3 (1) (X,X) + -ffx)--yyd Moreover, attention mechanisms have been utilized in the multi-task learning community, such as the multi-task atten- tion networks (MTAN) in (Liu al. 2019). In addition, the weighted sum objective function has been well-established for balancing task contributions (Vandenhende et al. 2021).
Regarding to the graph among the multiple outputs, the GMOGP introduces the probability graph- ical model to the MOGP regression, leading to an efficient MOGP alternative with attention-based Soxible/non-Gaussian priors and various graphical representations.
7

---

Published as a conference paper at ICLR 2024
Table I: The average test RMSE of the synthetic experiments described in Section 5.1 All metrics
are compared against the baselines: [1] Isolated SOGPs, [2] LMC (de Wolff et al. 2021). [3] free-
form task similarity model (FICM) (Bonilla al. 2007). [4] Gaussian process regression network
(GPRN) (Wilson et al. 2012), and [5] convolution process (CMOGP) (Alvarez & Lawrence 2011
(Np The flow parameters, Vn The variational parameters.)
Average RMSE
Test NLL
Kaim
Number of Parameters
1
SOGP
0.5653++0.003
0.4891 +0.0043
N
4$ (2+ 25
2
LMC.
0.5917+-0.006
0.5543++00.006
S X N
3
FICM
0.5544110.0444
0.4798110.0176
S X N
5)
GPRN
0.581910.0207
0.5787110.0445
X
2(S + +)+3
5
CMOGP
0.5539110.0089 0.4689110.014
X N
(2 (2+ S)Q 25
66
GMOGP
0.0044
0.1636110.0133
S(S + 5) +1
7
TGMOGP
0.5343:10.0023
0.635411.0023
S(S + + 3/xF + Vn
5
EXPERIMENTS
In this section, we benchmark the predictive performance of the proposed (transformed) GMOGP against baseline MOGP methods and isolated SOGP using various synthetic and real-world datasets.
5.1 SYNTHETIC DATA EXPERIMENTS
multi-output regression task with non-Gaussian noise and different function compositions is eval- uated. Five outputs (S 5) are generated by the following functions specialized at X R x800
y (1) = f1(X) + 61,
(22)
y- f1(X) f2(X)
(23)
y (3) sinh(2aresinb(f f2(X)) 63);
(24)
( = (tan (X) +
(25)
(5) -5/s(X)/4 + E5,
(26)
where fi(x) = 2 c1 + x2). f2(x) = (21 (2) = == f((x) = log(21 + +3),
and 61, ........ 6s are i.i.d. Gaussian noise with common standard deviation 0.2. The detailed experiment settings and available marginal flows are elaborated in Appendix B.
Table shows the compared predictive performance 600 test points. The GMOGP is competitive with the FICM but exhibits a smaller test NLL., and surpasses other baselines both test RMSE and NLL. Given the test error for each output shown in Figure 5(c) the TGMOGP achieves improved performance, especially in the 3rd output, confirming the enhanced representation ability w.r.t. the non-Gaussian data. The inferior performance of the classic MOGP methods (even worse than the isolated SOGP) indicates inefficient model learning and correlation measures in high- nigh-dimensional space. Moreover, increasing the number of independent latent processes Q LMC variants yields little performance gain (shown in Figure 5(b).. Other results are obtained with Q S. By contrast, our TGMOGP achieves superior results constantly among all different training sizes (Figure
5057
50GP
NDN
UMC
N=100
LMC
TICH
GPAN
RANSE
T24D07
GWOSP
TGHOGP
Thaning Sampiees
Number of
0
05
(a)
(b)
(c)
Figure 5: The sub figures show (a) the average RMSE changes with the number of training samples, (b) the RMSE versus the number of latent independent GPs, and (c) the test error of each output.
8

---

Published as a conference paper at ICLR 2024
Table Comparison of test RMSE on real datasets, where SGPRN (Li et al. 2021) and variational
LMC (V-LMC) are tested. The shadowed results are learned with two distributed computing units.
Datasets
SOGP
V-LMC100
FICM
SGPRN
GMOGP
TGMOGP 100
JURA
0.605=0.01
0.443+0.01
0.394+0.05
0.438=0.02
0.376=0.01
0.382:+0.01
ECG
0.24
0.22911.01
0.222:10.01
0.232:10.02
0.219110000
0.217:10.00
EEG
0.343-0.05
0.207110003
0.147:10.03
0.261-10003
0.082-10011
0.117:10.00
SARCOS;
1.139=0.01
1.06340.001
0.792+40.04
0.844-1.004
0.64
0.558:10.00
KUKA
0.05110.01
0.14100.01
0.03100.01
0.12.10.02
0.02/ 0.02 10.00
0.04/ 0.04 1.0011
Test NLL
-0.25+0.01
-0.51+0.01
-0.65 +0.02
-0.55=0.01
-1.81/ -1.76 +0.02
-3.49/
-3.48 +.001
SARCOS2
0.26=0.05
0.29+0.04
0.33+0.03
0.21/ 0.22 #0.02
0.16/ 0.16 #0.01
Time/Iter
20.75(s)
370.2(s)
419.3(s)
22400(s)
32.41/ 21.07 (s)
4.65/ 3.43 (s)
Table 3: Instances of the learned attention coefficients values from the GMOGP (O)j 1).
Coefficient| 01,2 01,3
01,4
ar2,
022,3
0r2,4
Cr3,1
Or3,2
0,3,4
0x4,
04,,2 0.4,3
SARCOS, 1.9c-4
0,998
2.9c-4
5.4c-4
6.9c-4
0.898
7.3c-5
5.1e-5
0.989
1.3e-3
0.966 6.7e-4
5.2 REAL-WORLD DATA EXPERIMENTS
In this section, we evaluate the predictive performance of the GMOGP alongside scalable MOGP models in real-world applications. The data description and tasks are detailed in Appendix BB3 where we select the first 2k/20k training samples (denoted as SARCOS from the SARCOS dataset to investigate how performance varies with different training sizes. The results Table 2 show that our GMOGP outperforms bascline models on all real datasets in accuracy and efficiency. For variational inference with 100 inducing points, the TGMOGP100 100 achieves lower predictive error than the V-LMC100 with largely reduced time cost, and even surpasses the models trained with full data. One probable reason is that the variational distribution generated from a free-form Gaussian in the TGMOGP 1o0 is also transformed by the SAL flow, which improves the possibility of learning sufficient statistic. Another reason is that the TGMOGP transformed by the SAL flow can describe asymmetry and oon-aussian tailed distributions (Jones & Pewsey| 2019), and the SARCOS data is negatively-skewer with mean median < mode. The trend illustrated in the Figure 5(a] is also shown by the different rankings across the two SARCOS data. As for the two larger datasets (KUKA (12k) and SARCOS2), we learn the GMOGP-based models through the distributed learning framework using two distributed computing units (the shadowed results), where comparable results achieved with reduced time consumption. In comparison of the test negative ligllikelihood, the outstanding result of the TGMOGP model implies it has a more appropriate model complexity.
The graphical representations learned along with the improved predictive performance the multi- output regression tasks (KUKA and SARCOS)) can be indicated by the attention coefficients (Table identify the parent of the node representing the 2nd output
3). In the SARCOS (4-outputs), is the node corresponding to the 4th output (02.4 = 0.8978), since we unlink the other two nodes with coefficients approaching O. Similarly, the dependence between the 3rd and 4th outputs can be implied by 0334 0.9998. To understand the graphical representations, we calculate the Pearson correlation coefficients between the 2ddd and 4th outputs, which give high values of 0.7744 and 0.9657 coincidently. Other interpretations of a traffic prediction task are explored in Appendix B.3
6 CONCLUSION
We propose graphical multioutput Gaussian process (GMOGP) model, an innovative framework tailored for efficient multi-output regression and graphical representation learning. Noteworthy model flexibility and optimality are verified by superior predictive performance achieved across various synthetic and real datasets. Distributed learning frameworks and sparse variational infer- ence methods can be directly applied to the proposed GMOGP framework, giving chance to deal with large datasets and a great number of outputs. Moreover, varied graphical representations and conditional dependence empower the vanilla MOGP models with more representation ability.

---

Published as a conference paper at ICLR 2024
7 ACKNOWLEDGMENTS
This work was supported by the NSFC under Grant No. 62271433, and part by Shenzhen Science and Technology Program under Grant No. JCYJ2022053014388000 and No. RCJC202106091044 481 14.
REFERENCES
Mauricio Alvarez and Neil D Lawrence. Computationally efficient convolved multiple output Gaussian processes. The Journal of Machine Learning Research, 12:1459-1500 2011.
Mauricio A Alvarez, Lorenzo Rosasco, and Neil D Lawrence. Kernels for vector-valued functions: A review. Foundations and Trends() in Machine Learning, 4(3):195-266, 2012.
Christopher M Bishop and Nasser M Nasrabadi. Pattern recognition and machine learning, vol- ume 4. Springer, 2006.
Edwin V Bonilla, Kian Chai, and Christopher Williams. Multi-task Gaussian process prediction. In Advances in neural information processing systems, volume 20, 2007.
Shaked Brody, Uri Alon, and Eran Yahav. How attentive are graph attention networks? In Interna- tional Conference Learning Representations, 2022. URL https::/ net/
Forum? d========F2xxxxxxxx
Wessel Bruinsma, Eric Perim, William Tebbutt, Scott Hosking, Amo Solin, and Richard Turner. Scalable exact inference in multi-outpin Gaussian processes. In International Conference Machine Learning, pp. 1190-1201. PMLR, 2020.
Yair Censor. Pareto optimality in multiobjective problems. Applied Mathematics and Optimization, 4(1):41-59, 1977.
Kai Chen, Twan Van Laarhoven, Elena Marchiori, Feng Yin, and Shuguang Cui. Multitask Gaussian process with hierarchical latent interactions. In ICASSP 2022-2022 IEEE International Confer- ence on Acoustics, Speech and Signal Processing (ICASSP), pp. 4148 4152. IEEE, 2022.
Wei Chu, Vikas Sindhwani, Zoubin Ghahramani, and S Keerthi. Relational learning with Gaussian processes. Advances in Neural Information Processing Systems, 19, 2006.
Seokhyun Chung and Raed Al Kontar. Federated multi-output Gaussian processes. Technometrics, pp. 1-14, 2023.
Yijue Dai, Tianjian Zhang, Zhidi Lin, Feng Yin, Sergios Theodoridis, and Shuguang Cui. An inter- pretable and sample efficient deep kernel for Gaussian process. In Conference on Uncertainty in Artificial Intelligence, pp. 759-768. PMLR, 2020.
Andreas Damianou and Neil D Lawrence. Deep Gaussian processes. In Artificial intelligence and statistics, pp. 207-215. PMLR, 2013.
Taco de Wolff, Alejandro Cuevas, and Felipe Tobar. Mogptk: The multi-output Gaussian process toolkit. Neurocomputing, 424:49-53, 2021.
Mare Peter Deisenroth and Jun Wei Ng. Distributed Gaussian processes. In International Conference on Machine Learning, pp. 1481-1490 July 2015.
Jacob Devlin, Ming- Wei Chang. Kenton Lee, and Kristina Toutanova. Bert: Pre-training of deep bidirectional transformers for language understanding. arXiv preprint 110iv:1810.04805, 2018.
Robert Dürichen, Marco AF Pimentel, Lei Clifton, Achim Schweikard, and David A Clifton. Mul- titask Gaussian processes for multivariate physiological time-series analysis. IEEE Transactions Biomedical Engineering, 62(1):314-333 2014.
Jinyuan Fang, Shangsong Liang, Zaiqiao Meng, and Qiang Zhang. Gaussian process with graph convolutional kernel for relational learning. In Proceedings of the 27th ACM SIGKDD Conference Knowledge Discovery & Data Mining, pp. 353-363, 2021.
10

---

Published as a conference paper at ICLR 2024
Michel Goulard and Marc Voltz. Linear coregionalization model: tools for estimation and choice of crois-variogram matrix. Mathematicaal Geology, 24(3):269-28) 1992.
James Hensman, Alexander Matthews, and Zoubin Ghahramani. Scalable variational Gaussian pro- cess classification. In Artificial Intelligence and Statistics, pp. 351-360. PMLR, 2015.
Chris Jones and Arthur Pewsey. The sinh-arcsinh normal distribution. Significannee 16(2):6-7, 2019.
Dongkwan Kim and Alice Oh. How to find your friendly neighborhood: Graph attention design with self-supervision. arXiv preprint arXiv:2204.0453 04479, 2022.
Daphne Koller and Nir Friedman. Probabilistic graphical models: principles and techniques. MIT press, 2009.
Arthur Leroy, Pierre Latouche, Benjamin Guedj, and Servane Gey, Magma: inference and prediction using multi-task Gaussian processes with common mean. Machine Learning. 111(5):1821-184,, 2022.
Arthur Leroy. Pierre Latouche, Benjamin Guedj. and Servane Gey. Cluster-specifie predictions with multi-task Gaussian processes. Journal of Machine Learning Research, 24(5):1-49 2023.
Shibo Li. Wei Xing, Robert M. Kirby, and Shandian Zhe. Scalable gaussian process regression networks. In Proceedings of the Twenty-Ninth International Joint Conference on Artificial Intel- ligence, 2021.
Haitao Liu, Jianfei Cai, and Yew-Soon Ong. Remarks on multi-output Gaussian process regression. Knowledge-Beud Basedd Systems, 144:102-121) 2018.
Haitao Liu, Jiaqi Ding, Xinyu Xie, Xiaomo Jiang, Yusong Zhao, and Xiaofang Wang. Scalable multi-task Gaussian processes with neural embedding of coregionalization. Knowledge-Bau Based Systems, 247:108775, 2022.
Shikun Liu, Edward Johns, and Andrew J Davison. End-to-end multi-task learning with attention. In Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, pp- 1871-1880, 2019.
Sanae Lotfi, Pavel Izmailov, Gregory Benton, Micah Goldblum, and Andrew Gordon Wilson. Bayesian model selection, the marginal likelihood, and generalization. In International Con- ference on Machine Learning, pp. 14223-14244 PMLR, 2022.
R Timothy Marler and Jasbir S Arora. The weighted sum method for multi-objective optimization: new insights. Structural and multidisciplinary optimization, 41:853-862, 2010.
Juan Maroñas and Daniel Hernández-Lobato, Efficient transformed Gaussian processes for non- stationary dependent multi-class classification. In International Conference on Machine Learn- ing. pp. 24045-2408 PMLR, 2023.
Juan Maroñas, Oliver Hamelijnck, Jeremias Knoblauch, and Theodoros Damoulas. Transforming Gaussian processes with normalizing flows. In International Conference on Artificial Intelligence and Statistics, pp. 1081-1089. PMLR, 2021.
Charles A Micchelli and Massimiliano Pontil. On learning vector-valued functions. Neural compu- tation, 17(1):1772044 2005.
Yin Cheng Ng, Nicolò Colombo, and Ricardo Silva. Bayesian semi-supervised learning with graph Gaussian processes. Advances in Neural Information Processing Systems, 31, 2018.
Julián David Pastrana-Cortés, David Augusto Cardenas-Peña. Mauricio Holguín-Londoñoe Germán Castellanos-Domine Domiinuez, and Álvaro Angel Grozco-Gutiérree Multi-output variational Gaussian process for daily forecasting of hydrological resources. Engineering Proceedings, 39( 88:3, 2023.
Paris Perdikaris, Maziar Raissi, Andreas Damianou, Neil D Lawrence, and George Em Karniadakis. Nonlinear information fusion algorithms for tita-efffcient multi-fidelite modelling. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 17301198):00100511, 2017.
11

---

Published as a conference paper at ICLR 2024
Lishan Qiao, Limei Zhang, Songean Chen, and Dinggang Shen. Data-driven graph construction and graph learning: A review. Neurocomputiing. 312:336-351, 2018.
Carl Edward Rasmussen and Christopher KI Williams. Gaussian processes for machine learning-
MIT press Cambridge, MA, 2006.
James Requeima, William Tebbutt, Wessel Bruinsma, and Richard E Turner. The Gaussian pro- cess autoregressive regression model. In Proceedings of the 22nd International Conference on Artificial Intelligence and Statistics, pp. 1860-1869. PMLR, 2019.
Danilo Rezende and Shakir Mohamed. Variational inference with normalizing flows. Interna- fional conference on machine learning, pp. 1530-1538. PMLR, 2015.
Gonzalo Rios. Transport Gaussian processes for regression. arXiv preprint arXiv:2001. 2020.
Gonzalo Rios and Felipe Tobar. Compositionally-wann Gaussian processes. Neural Nenworks, 118:235-246, 2019.
Bernhard Schölkopf, Alexander J Smola, Francis Bach, al. Learning with kernels: support vector machines, regularization, optimization, and beyond. MIT press, 2002.
Marcin Sendera, Jacek Tabor, Aleksandra Nowak, Andrzej Bedychaj, Massimiliano Patacchiola, Tomasz Trzcinski, Przemystaw Spurek, and Maciej Zieba. Non-Gaussian Gaussian processes for few-shot regression. In Advances in Neural Information Processing Systems, volume 34, pp-
10285-10298, 2021.
Sarath Shekkizhar and Antonio Ortega. Graph construction from data by non-negative kernel re- gression. In ICASSP 2020-2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pp. 3892-3896. IEEE, 2020.
Ricardo Silva, Wei Chu, and Zoubin Ghahramani. Hidden common cause relations in relational learning. Advances neural information processing systems, 20, 2007.
Edward Snelson, Zoubin Ghahramani, and Carl Rasmussen. Warped Gaussian processes. In Ad- vances in neural information processing systems, volume 16, 2003.
S. Theodoridis. Machine Learning: Bayesian and Optimization Perspective. 05 2015. ISBN 978-0- 108-0-12-80111 12-801522-3.
Michalis Titsias. Variational learning of inducing variables in sparse Gaussian processes. In Artificial intelligence and statistics, pp. 567-574. PMLR, 2009.
Joel A Tropp and Anna Gilbert. Signal recovery from random measurements via orthogonal matching pursuit. IEEE Transactions on information theory, 53(12):4655 2007.
Simon Vandenhende, Stamatios Georgoulis, Wouter Van Gansbeke, Marc Proesmans. Dengxin Dai, and Lue Van Gool. Multi-task learning for dense prediction tasks: A survey. IEEE transactions on paltern analysis and machine intelligence, 2021.
Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Łukasz Kaiser, and Illia Polosukhin. Attention is all you need. In Advances in neural information processing systems, volume 30, 2017
Petar Veličković, Guillem Cucurull, Arantxa Casanova, Adriana Romero, Pietro Liò, and Yoshua Bengio. Graph attention networks. In International Conference Learning Representations, 2018. httlltttt:///peeei
net/ fwweee//oommmmddeeeeeee
Pascal Vincent and Yoshua Bengio. Kernel matching pursuit. Machine learning, 48:165 18:165-187, 2002.
Andrew Gordon Wilson, David A. Knowles, and Zoubin Ghahramani. Gaussian process regression networks. In Proceedings of the 29th International Coference on International Conference on Machine Learning, pp. 1139-1146, 2012.
Qiang Yang, Wei-Neng Chen, Tianlong Gu, Hu Jin, Wentao Mao, and Jun Zhang. An adaptive stochastic dominant learning swarm optimizer for high-dimensional optimization. IEEE Trans- actions Cybernetics, 12600 2020.
12

---

Published as a conference paper at ICLR 2024
A MATHEMATICAL APPENDIX
In the part A, we will give some insights and detailed derivation of the graphical multioutput Gaus- sian process (GMOGP) prior, the evidence lower bound (ELBO) of the transformed GMOGP, etc.
A.1 GMOGP PRIOR DERIVATION
In the main paper, we can derive the target node in the GMOGP as the following aggregated form, i.e.,
Σ
aijsx ( + m, + th.
(27)
jEPa;
Recalling that each node represents the random variables generated by a single output Gaussian pro- cess, i.e., we have p(Fk)) Nmm,, Ko, (X, x)), E pa, and tb ~0 ki, (X, Therefore,
the GMOGP prior of target node p(5x²) follows Gaussian with mean:
1 = E
(E( + m;
(28)
EE
Σ Epp,
mi m:
(29)
and covariance:
( ) (1 - 11) '
cov(f2
E
(30)
Σ ai, ( - + tb }
11 E(
(31)
/EPa;
ai.coov (1) : 114
+ ki, (X,X),
(32)
Epp,
where,
cov(T (k)) = cov( 11
(33)
(1P - B((X)) (1 - =
E
(34)
E (1L - B(S(1) { Σ
aik (5L - -B(((>)) + a/s H
(35)
AEPa;
CLL (K )(0)
(1+
(36)
kEPa;
Notably, under the specific graphical representation (such the decomposition in Figure 2 there are no parents of other nodes j. i.e., pajuzi = {0). Then, we have Σxepy, 0s, ( (()) cov 5k²) = oi.y (X, X). As a result, substituting this result to the above Eq. 32 we can derive the covariance form for our target node i, i.c.
1 (X, cov c
Σ jEPa
(X, + (x,
(X, XX.
(37)
Remark A.1. In the context of the probabilistic graphic model, only directed acyclic graphs (DAGs) can be generated and represented. Given a well-defined ordering w.n.t the joint distribution (see Remark ) the covariance can be evaluated recursively starting from the lowest numbered node.
Remark A.2 (Ordering). The decomposed product term of the joint distribution is asymmetrical, i.e. an implicit order must be chosen, such as (Vi, U2...., v8g) in E2. Following a different ordering, we can obtain different decomposition and hence different graphical representation.
13

---

Published as a conference paper at ICLR 2024
A.2 ELBO OF THE TRANSFORMED GMOGP
Following a similar derivation as Titsias 2009) and (Maroñas et al. 2021), the log marginal like- lihood of output i ..2,...,5 S admits the following inequality: (We denote specific locations
with additional subseripts, such that f( Go. (10g. = G 2 o G (6) G (Jo Km.
)
with fox
P py(()) 1 P
(12 ux
log p(y(()()) 2 EE
a (112 uxk
log
(38)
(6)
k 2uK
II P (yl)
(1) (()
(())
Kxx
p
[xx uk
E
(1 u4k)
log
(39)
(1)
Kx uK
(
(8)
[log (x( IJ xxxx
)
E
( xxx
+ +
(1 [xx 2)
log
(6)
1xx .uk
ELL
KKT.0
(40)
where we define the approximate posterior as:
q( () = piSk' 1 ()
(41)
OG ok
(uk (8)
= =p(fky 2 )q(2b)) (0)
det
(42)
Oug (5)
Here,
pisky (2 = pllky. u
1/(
(43)
GG U
8G G
%
x
K
pl
(1)
II k=
det
ouk
/ plu
(44)
0G
aG (n)
05
Bu
K-1
OG
(u1 (5)
pl (/0)
(1)
p(ub
II
det
(
auk
k=0
K-
aG 61 (122)
JG 100 (JER)
OG
(1)
(us
(
-1
JG
(uk
(
GG
of
(u (
del
del
05(0
Duk
Ou
asky
Ou (1)
k=0
evaluated to zero for marginal flows
(45)
K-1
GG a
= =P(f
(1)
II k=
det
(46)
Moreover, applying marginal flows (coordinate-wist we can derive a more simplified form:
aG 19 (F(?)
N II
dG
det
(47)
01
dfk xm
14

---

Published as a conference paper at ICLR 2024
and suppose
(u)))) Nm (m))
(48)
with m E AM and KP MMM being the variational parameters. The detailed derivation
(4)
and simplification of the expected log-likelihodd (ELL) and KL divergence in Eq. (40) are specified separately in the next subsections.
A.2.1 KL DIVERGENCE FOR A TARGET OUTPUT
Following the approximate posterior in Eq. 42 the conditional term in the KL can be canceled, namely:
P ( (
(1) K
KT, (4)
log
Kx
(49)
(1
(8) (1)
[ ukk
80 (u
pl (8)
ta K )p(
K
det-
x
k=
log
(50)
80 (1)
pl xx
a K )a(u6)"
111 K-
del
8u
log plag (8)
-F
(51)
("
(((e)))
log p(ug (3)
2( (1)))
g( q(u))) (0))
(
- KL (q(u))) p(
(52)
where we first marginalized fky and then applied the law of the unconscious statistician (LOTUS). Concretely, given an invertible transformation GG.) (any transformation implies valid stochastic processes, we refer the readers (Rios) 2020) for more details), the expectation of any function g(-)
(() under the distribution p(ux) K admits:
E plu
(8) (
E p(c(f))
[g(Goo (m)
(53)
x
A.2.2 EXPECTED LOG-LIKELIHOON UNDER MARGINAL FLOWS
After integrating out the uk and applying the LOTUS rule over the expectation with regard to ( ), we can simplify the ELL ( term in Eq.44 as:
N
[ogp (y!0) (ypp(yh0)) (12).)
ELL
(1 xx
(54)
Σ
N
/ [logp (x10) 11 ){
(fh (() [2
(2
(ux)
du drkk
(0)
(55)
N
Σ
E
(1k
logg (y(²) 15kx
(56)
[logp (y10) Go (f0f) ))]
(ysp(y0)
E
(1022
(57)
15

---

Published as a conference paper at ICLR 2024
The distribution q ( can be calculated by:
al (l66))) / p u0 12 (a(ub
(1) ydu (C)
(58)
= -N (x, (X, 2) (2. !
(x, (x. 2) (2, (2) | (2, KTa (2. -k6 (z))x))
(59)
with
plu!)) N (o, k0 (2,2))
(60)
1 14b = N N( ( (x, () (2, up -
(x, x) 18 (x, () (2, -)) (2, x)):
(61)
where R Mxd denotes the inducing points locations, and X Nxd are data inputs.
In conclusion, the resulting ELBO of the log marginal in Eq. 40) for output i is given by:
ELBO ELL KL
(62)
log plac (1)
E (162)) (lgppy(") Go, [102)]] (( + (usi))
(63)
g(ue))) (6)
(0) log
Σ
E (100 |logp (y(h))
Co. Co. 10 ))]
a(
(1))0)
uu
(64)
Note that the observation value yn only depends on the function evaluation at position Therefore, minimizing the objective function in the main paper (Eq, 19 with zero mean equals to:
P1:
max (
Σ ('g
(65)
i=.,.....S
with ( = { au, o.). According to the above inequality in Eq. Eq,40), every output can derive the corresponding lower bound ELBO (i) = 1,2,.... S. Applying the following logic: a a. d: + < + d, (a, b, c, IR),
(66)
we can formulate the variational lower bound of the transformed GMOGP as follows:
S
log p(u))) (3)
log logp(y (0) 2
E
logp(y(") Go (6 + (w)))
(67)
q( ((ug))
(2)
OLBO
Then, instead to solve the optimization Pi in Eq. 65
Eq.65)
we can maximizing its lower bound, namely:
logr(y Go. (1)2) + o (1
p (
P2
max
log
( 4.. (l
7(
1-1.2....
(68)
8
log pl p(u!))
(0)
min (,,0) an ml! k()); 2........
E
(1
logr(y(0) Go (ogp(y")GG,,(0))
(G, +B (a(0))
q(()))) (1)
(69)
which recovers the negative ELBO (NELBO) in the main paper, and the objective function can be further decomposed with locations as below:
8
N
log p(ug)) (1))
logp (y Go. GG. })] +B
(70)
q((u))) (8)
16

---

Published as a conference paper at ICLR 2024
A.2.3 VARIATIONAL GAP
The gap between the evidence and the ELBO for each output can be measured follows:
log llg (LBO
p (x((0)
(10
_K
E 0 ukk
[ogp(s(1)) E
-B, (1 x)
log
(71)
q (5xxx u4x
(6)
(y(t) 1
(1)
(0)
logg (fxx ux )
Kxx
uK
E
a (12 u12)
-E (1R xx
log
p(y(1)
(72)
a (fR
E ( wl?)
log
(73)
P (J 18(1)
=KL (L u
ll
Kx
(74)
The gap is yielded as the Kullback-Leibler divergence between gf Kg This fact forms the basis of the variational inference algorithm for approximate Bayesian inference.
2 and p( p(fky xx ux y()))
A.3 GMOGP PREDICTION COMPARISON
For the GMOGP model proposed in Section joint Gaussian:
we can derive the conditional distribution from the
ply (2 p(f O a.)((
(75)
Correspondingly, the estimates of the output i at a test data can be calculated from:
GMOGP:
f(0 = k(0 (x., X)
[ (X,x)+
+
(0)
(76)
The counterpart in the classic LMC model admits:
Q (x X) - 1I +1:NJ] j=1 9=1
LMC:
[1]
(77)
N-dimensional vector
where [i] denotes the ith element of the LMC prediction vector, and β = (X,(XX, represents the S-dimensional vector calculated via the high-dimensional gram matrix. Regarding to the output the prediction only counts the elements from index (j 1)N + 1 to Nj.
-
In contrast, it is noticeable that the first part on the right-hand side of Equations (76) and 777 are both derived composite kernel functions of the test point x, and support training points. As shown in the experiments on synthetic data, we find out that little predictive performance gain can be brought by increasing the number of latent independent processes (see Figure 5(b)). Moreover, the LMC with constant coefficients Ojj may lack of flexibility and representation ability.
In the context of kernel methods, the predictions are determined by coefficients and non-linear mapping function in a kernel reproducing Hilbert space. For the above two models, the coefficient / and Ø both convey the correlated information from other observed output values. Distinct from the LMC, the GMOGP admits less computational complexity and graph learning.
Another alternative view of GMOGP is that it can remedy the cancellation of inter-output transfer in the FICM proposed in Bonilla et al. (2007) with noiseless observations. Since the noiseless case with block design, the predictions for output depend only on the observations ((), In other
17

---

Published as a conference paper at ICLR 2024
words, there is cancellation of information transfer among other outputs. Specifically, given the kernels, the prediction at a test data x, for output can be concluded as:
( (xx)][ = (Kar(x 6 l(x., X)) (RM(x,,) @ (),
[1]
(78)
{ {( (Kas(x, (KM (x) & X) k(x,, X) X)) ())
(79)
(80)
Here, the [f (x.) d represents the Tt element of the vector-valued prediction. Compared with the GMOGP in ( the noiseless version, the information other outputs can be implied from the inner product based attention coefficients enrolled in the aggregated kernel function.
A.4 JOINT DISTRIBUTION OF MULTIPLE OUTPUTS
From the main paper, the joint distribution of all S outputs defined by the specific graph structure can be formulated as:
(1) (2) p(Xx Ix {(0 p(fx (pa))p(paa
(81)
=
( pa.) II
(82)
jEPa,
In the GMOGP model, we model the target output with the following conditional distribution:
p(Xx pa.) = N (N
Σ
anss& + m, ko (X.x)) i
(83)
jEpa
and ) N(m, (X, x)).j E pa.. After the model training, given the learned parameters, obtain the corresponding joint distribution from Eq 82 While in the TGMOGP, since transform the GMOGP prior with the marginal flow, the TGMOGP prior becomes:
=
aG II del 22
= P0
(84)
Under the equality derived in Eq. 27). the invertible transformations on the target node i can be represented as:
Co. (1
G#.
Σ
Eassll + m, + t/i
(85)
jEpa,
Given the parent nodes pa, we have
K
O of (112)
PG.. (1 pa.) p(xx
pa.)
det
(86)
k=0
Therefore, the resulting joint distribution can also deduced with learned coefficients and hyper- parameters.
A.S PARETO OPTIMALITY
Pareto optimality is a foundational concept in the optimization community. In single objective prob- lems, the Pareto optimal solution unique. In the context of the multi-objective optimization, the multiple objectives need to be optimized simultaneously, looking for set of Pareto optimal so- lutions. This process is called multi-objective optimization (Censor 1977). In short, the Pareto optimal solution is set of 'non-inferior' solutions in the objective space, and defines boundary beyond that none of the objectives can be improved without sacrificing at least one of the other objectives.
18

---

Published as a conference paper at ICLR 2024
The benefits of solving the multi-ouiput regression in multi-objective problem (MOO) instead of the classic high-dimensional optimization considered in the MOGP methods are discussed here. Since the kernel hyperparameters è are shared among all different outputs to achieve information exchange, we can build an MOO problem,
F()) (4 (0), (())),..,((((()))
(87)
with S objectives the following form:
Siio (1)10)) ( (X.x) + ofIN) g(0)) log KG [X, X) of }
(88)
There exists competition among the multiple objectives, as the outputs correspond to multiple data sources and distinct hyperparameters. The illustration of the competition between two objectives is given in Figure 6 where the Pareto optimal solutions represent the points at the Pareto front.
Without loss of generality, the Pareto optimal solutions of the kernel hyperparameters are better than other feasible solutions, and can obtain more stable performance compared to those learned by high-dimensional optimization problem built in state-of-the-an (SOTA) MOGP methods. Further- more, the objective functions for target outputs are parameterized by unique attention coefficients, all of which are simultancously minimized zero.
Objective
Feasible solutions
Parcto front
Infeasible
Objective 2
Figure 6: The description of the Pareto optimal solutions with respect to the multi-objective problem.
B ADDITIONAL DETAILS ON EXPERIMENTS
In this section, we explicit the detailed experiment settings, learning parameters for all competitive models, and extra experiment results of various real-world datasets.
B. LEARNING PARAMETERS
In the competing methods, the primitive kernel function is selected to be the Squared Exponential (SE). The hyperparameters contain a length-scale an oupptt-scalessignal variance c, namely
kse(x, x') c. exe (- ||x (-1Ix-x'1 /1²)
x'1
(89)
The detailed learning parameters specified in Table 4 where NF denotes the number of stacked transformations in the flow (specific parameters of the flow are listed section B.2). and the Vm, represents the quantity of variational parameters (including the inducing points. the mean and covari- ance in the variational distribution). Also, the Q denotes the number of latent independent Gaussian processes used in the LMC variants, and the S represents the number of outputs.
19

---

Published as a conference paper at ICLR 2024
Table 4: The detailed learning parameters for competing models.
Acronym
Learning parameters (style and size)
[1]
SOGP
noise: S, length-scale: S, output-scale: S, mean: S
2 LMC
noise: S+1, mean: S, coefficients: SQ, length-scale: Q. output-scale: Q
3 FICM
noise: S+1, mean: S, raw-variance: S(S 1)/2, length-scale: 1
4 GPRN
noise: S+1, mean: S length-scale: Q + output-scale: Q ++1
5 GMOGP
noise: S+1, mean: S, attention: S + S. length-scale: S, output-seale: 5
6 TGMOGP
noise: S+1, mean: S, attention: S + S, length-scale: S, output-scale: S, NF Vm
B.2 MARGINAL FLOWS
In the experiments, stack K = 3 or 4 layers of composite Sinh-Archsinh flow with Affine flow (SAL), which can be formulated as below:
R2 c! (sinh (b()). arcsinh ( + d(i).
(90)
For the output i, E I, each layer in the flow has four free parameters, i.e., a (()) (0) d E R. Some commonly used and available formulas of the flow are listed the following Table More discussions and validity analysis can be found in (Rios & Tobar 2019).
Table 5: The list of common marginal flow types with input p( The parameters are constrained so that each individual flow is strictly increasing/ddcreasingg functions.
Flow
Forward
Inverse
Parameters
Log
log
exp(F
Softplus
log
exp b
J0x
+1)
log exp
Kxx
-1)
Affine
(5()
Kx
-a)/b
a, R
Sinh-Archsinh sinh arcsinh
sinh
(t
aarssinh
[
(1)
+aa)
a E R
Boxcox
sgn (
fox
sgn
a A + 1)
A (()
+
A>
TanH
tanh ( (50g + +c))
+
0, c, b,
B.3 REAL WORLD EXPERIMENTS
This section shows more detailed experiment results of the multi-output regression tasks across various real datasets. The data descriptions and specified tasks are demonstrated Table
Table 6: Description of the applied real datasets and multi-output regression tasks.
Dataset
Nuai
Nwa
Dx
S
Tasks
JURA
249
100
3
Cadmium, nickel, and zinc concentration prediction in Jura Extrapolation sampled signals from electrodes on scalps Traffic time- series prediction for three adjacent cells Extrapolation of signals from electrodes navel and head Torques prediction of seven joints for an lightweight am
EEG
256
100
7
TRAFFIC
500
172
6
ECG
6000
2000
5
KUKA
12235
5325
7
SARCOS;
2000
4449
21
Torques prediction of the 2,3,4,7 joints for anthropomorphic arm
SARCOS2
20000
4449
21
4
B.3.1 EEG DATA FITTING FOR THE FIRST OUTPUT
The main results of the proposed GMOGP models on the EEG dataset are shown in the main paper. Here, we want to draw the predictions of the output regarding to the electron Fg and show more data
20

---

Published as a conference paper at ICLR 2024
(a)
(b)
-10
(c)
(d)
Figure 7: The predictive performance of the EEG data with respect to the signal of electron Fo- Sub-figures show the fitting results from different methods: (a) GMOGP, (b) FICM, (c) LMC, (d) SOGP, where the cross items represent the training data points, and the inverted triangles are the locations of the 100 test points.
fitting details. In Figure 7 lower predictive uncertainty.
the GMOGP outperforms the other competing methods with
B.4 KUKA PREDICTIVE PERFORMANCE FOR THE MULTIPLE OUTPUTS
In Figure 8 show the test RMSE of each output, and compare the average RMSE the seven outputs. The standardized KUKA data has close mean, median, and mode, which may be better described by the Gaussian predictive distributions.
F1CM
T1EM
TIH
VLNC2
V.LMG
GHOSP
GNOGP
GMOGP
15MD0P.0
1GMOUP
Outpuct1
OUtDUE
Outbu:
500P
FICM
FICM
FXM
V-MMCa>
V-NGa
VLMG
V-LMCu
VMCuu
GNOGP
TENOGP
TGMOGP
Ouipuc4
Outeu.
12M
HCH
LMC
V-NE
VLMG
Y4l Waae
V-IMu
V-
VA
TGKOGP
TGNOSPax
TGHOGP g0
Durpit?
Average
Figure 8: The sub-figures show the test RMSE for output to output The last two show the average RMSE and time cost per iteration.
21

---

Published as a conference paper at ICLR 2024
B.4.1 TRAFFIC DATA EXTRAPOLATION
In this section, we test our proposed GMOGP models on the real time-series dataset:
Traffic data: The data containing 28 days of traffic measurements of 3 cells are collected by com- munication base stations in China. For each cell, gathered two traffic time-series data (the amount of traffic that has been accumulated per hour in 28 days). Each traffic data contains 672 samples. We use the first 500 samples for training and the last 172 for testing. This multi-output regression task can be modeled as GMOGP with six outputs, and the data are drawn in Figure
M A hm
0.5
1.0
0.5
1.0
0,
1.0
0.5
0.
0.75
0.25
0.
100
200
Figure 9: The description of the traffic data from neighbored cells in 28 days. The lines paired with black and red from top to bottom represent the traffic data in the same cell of different orientations.
1151
0.1700
1.0
0.506
0.099
0.098
0.094
0.1
0.1675
0.1850
2
0.129
0.117
0.116
RMSE
9.1625
6
0.122
0.357
1.0
0.123
0.126
0.145
0.1600
0.1575
0.119
0.399
0.122
1.0
0.121
0.12
0:1550
151
0.147
0.225
0.165
0,146
1.0
0.167
0.1525
0.171
0.157
SOGP V-LMC
FICM SGPRN GMOGPTGMOGP
(a)
(b)
Figure 10: Traffic data predictive results, namely, (a) average RMSE of different methods, and (b)
learned attention coefficient values
Figure 1O(a) shows the average test RMSE on the traffic dataset over five runs. The results demon- strate that the GMOGP-barn based models outperform other MOGP methods and the SOGP model. Since real data non-smooth and contains zero values, the isolated SOGP model can hardly extrapo- late well with insufficient model flexibility brought by primitive kernel. According to the attention coefficients learned by the GMOGP model Figure [0bb we can interpret graphical representa- tions and construct the joint distribution of all nodes. Intuitively, we can find out the output has strong conditional dependence with the output 2, and also can noticed by the coefficient 01,2. In fact. the st and the 2nd outputs are collected in the same cell with similar trends in practice. How- the other output with different envelope and regularity can provide little useful information when predicting the first output. Comparing to the classic MOGP models that correlate all outputs without selection, GMOGP can provide an efficient alternative for multi-output regression.
22

---

Published as a conference paper at ICLR 2024
C LIMITATIONS, EXTENSIONS, AND FUTURE WORK
The main limitations the GMOGP can be mitigated through various avenues:
Interpretation of the attention coefficients: The attention coefficients signify the depen- dence among nodes and graph structures in this paper. However, interpretations of the coefficients learned by the TGMOGP and SOTA graph attention networks are vagued and anticipated in Explainable Artificial Intelligence (XAI) field.
2. Dynamic dependence measure: The interaction or dependence dynamically changes with inputs.
3. Extensions of online learning frameworks and graph construction schemes: Opportunities exist for extending the GMOGP to online learning frameworks and enhancing graph con- struction schemes in an online setting.
4. Adaptation of attention mechanism/scoring function/ddeenddecce measure for different tasks: Consideration of diverse attention mechanisms, scoring functions, and dependence measures tailored to specific tasks in practice.
5. Exploration of different aggregation models (Eq. (16)..
6. Exploration of different objective functions and weighting schemes tailored for varied ap- plications.
In the future, exploration the of GMOGP with various types of adaptive/inuut-deeeeddnt covariance structures and flows would be instructive. Additionally, extending the GMOGP to equip data-dependent mean function, model multiple fidelity outputs, and facilitate sequential forecasting would be of great interest. We hope the GMOGP will inspire further research into explainable networks, bridging benefits from learning models and statistics.
23


Reviews:
[{"review_id": "jc1LxIVvb0", "reviewer": "Reviewer_9P5Q", "summary": "This work proposes a multi-output regression framework, termed graphical MOGP (GMOGP). The GMOGP is built upon a probability graphical model that can learn the conditional dependence and imply graphical  representations of the multiple outputs. The framework allows learning of the model hyper-parameters via a multi-objective optimization problem. The work also provides an extension by non-linear transformations in order to fit non-Gaussian data.", "questions": "---Specific comments---\n\n-In the Introduction, where it reads \"the input-output relationship with multivariate Gaussian...\", better \"...with a multivariate Gaussian...\"\n\n-I did not see any literature related to MOGPs with convolution processes which might indeed be a more general approach than the LMC. See for instance works like:\n\n\"P. Boyle and M. Frean, Dependent Gaussian processes (2005)\", \"J-J. Giraldo et al, Correlated Chained Gaussian Processes for Modelling Citizens Mobility Using a Zero-Inflated Poisson Likelihood (2022)\", \"T. M. McDonald et al, Nonparametric Gaussian Process Covariances via Multidimensional Convolutions (2023)\", \"M. A. Alvarez and N. D. Lawrence, Computationally Efficient Convolved Multiple Output Gaussian Processes (2011)\"\n\n-I do not follow what the authors mean by \"on a set of shard GPs\". Is shard an adjective? Was it shared GPs maybe?\n\n-In section 2: \"...is parametrized by: length-scale $l$ and signal variance...\", better to use \"...by: a length-scale $l$ and a signal variance...\"? \n\n-After Eq. (3), change \"vector-valued function... are determined by...\", by \"vector-valued function... is determined by...\"\n\n-After Eq. (4), with the current notation $Y$ might be understood as a matrix instead of a vector, it might be better to clarify the dimensionality of $Y = [{\\mathbf{y}^{(1)}}^\\top,...,{\\mathbf{y}^{(i)}}^\\top,...{\\mathbf{y}^{(S)}}^\\top]^\\top \\in \\mathbb{R}^{NS}$.\n\n-After Eq. (7), shouldn't it be $K_{i,i^{\\prime}}(\\mathbf{x},\\mathbf{x}^{\\prime})=\\sum_{q=1}^{Q}a_{iq}a_{i^{\\prime}q}k_q(\\mathbf{x},\\mathbf{x}^{\\prime})$?\n\n-In section 3.3, lower case \"applying\" in \"In addition, Applying\". Also, include \"a Gaussian\" for \"...follows Gaussian that can be normalized...\"\n\n-After Eq. (18), the functions $f_1(x), f_2(x),f_3(x),f_4(x)$ should receive a vector $\\mathbf{x}=[x_1,x_2]^\\top$ as argument, i.e., $f_1(\\mathbf{x}), f_2(\\mathbf{x}),f_3(\\mathbf{x}),f_4(\\mathbf{x})$.\n\n-In the figure 5(c), correct the legend \"GRPN\" to \"GPRN\".\n\n---Other Questions---\n\n-In the GMOGP I only noticed a derivation for isotopic data, but I did not see any mention that the method would allow the use of heterotopic data, would the method allow its application to heterotopic data? \n\n-In line with the last question, how is the behaviour of the GMOGP in the context of notable unbalanced data per output?\n\n-It is mentioned that the weights $w_i$ are set equal since there is not priority among the outputs, is this assumption still valid in the context of markable unbalanced data along the outputs?\n\n-From the practitioner perspective, what would it be an appropriate way to initialise the weights $w_i$? Is there much effect on initialising them quite small (say $(0,1]$) or quite big ($>10$)?\n\n-It might be useful to also show the performance comparison with a MOGP with Convolution processes a more general and powerful approach than the LMC. \n\nIt would be important to also discuss about the possible limitations of the work, there was nothing related to it in the paper.", "limitations": "", "rating": 8, "confidence": 4, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "Originality: To the best of my knowledge the work seems original, built upon different ideas from the state-of-the-art that the authors appropriately merged to compliment each other.  \n\nQuality: The experiments, both with synthetic and real world data, are useful to provide an idea of the approach. \n\nClarity: The work is very well written in general with neat notations and explanations of the work. Quite consistence in the mathematical notation. There were just a few typos to easily fix.\n\nSignificance (importance): The application of Multioutput Gaussian processes tend to be prohibited in the practice when the number of data observations and  outputs grows. Therefore this work is significant due to allowing the MOGP computational complexities to be considerably reduced and additionally permitting the non-conjugate likelihood/prior scenario to be applied in a practical manner.", "weaknesses": "The experiments, both with synthetic and real world data, are useful to provide an idea of the approach. Nonetheless, the work lacked of discussion, citation and experimentation of the Gaussian process models based on convolution processes, a more general and powerful approach than the Linear Model of Coregionalisation (LMC).\n\nThere is not discussion in the paper about the limitations of the work.", "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": "", "extra_scores": {}}, {"review_id": "EqzsvW4EaR", "reviewer": "Reviewer_sgU9", "summary": "The paper introduces a new construction of multioutput Gaussian processes (MOGP) by leveraging probability graphical models to capture dependencies among outputs through an attention mechanism. The Graphical MOGP (GMOGP) can be seen as a generalization of the classic MOGP-LMC, wherein dependencies (i.e., coefficients of coregionalization) are determined based on the output's parent nodes. The evidence lower bound of the GMOGP is introduced for the joint estimation (in a variational fashion) of the graphical structure and kernel parameters. The proposed framework is tested on several examples, and results allow assessing its competitivity w.r.t. the state-of-the-art.", "questions": "**Questions**\n\n- Although the diversity of the examples, it is not completely clear that the model properly estimated the dependencies between the outputs. Can this be verified in a toy example? I suggest considering a toy example with 3 outputs where one of them is highly correlated to another, and the remaining one has no influence on the other ones. For instance, one may simulate $f_1, f_2$ from a GMOGP and $f_3$ from a SOGP. For the GMOGP, I suggest considering $\\alpha_{1,2}^\\ast \\gg \\alpha_{2,1}^\\ast$. The other parameters of the graphical representation and kernels can be chosen arbitrarily. Once an illustrative sample is chosen, a new GMOGP can be trained considering data from the 3 outputs and initial values of cross-attention coefficients (non-diagonal elements) all equal (to 0.5). The framework is expected to recover the values of $\\alpha_{1,2}, \\alpha_{2,1}$ close to the ones used to generate the dataset and to estimate small values for $\\alpha_{1,3}, \\alpha_{2,3}, \\alpha_{3,1}, \\alpha_{3,2}$. It is also possible to compare the kernel parameters estimated by the GMOGP with the ones used to generate the data.\n- When considering systems with a significant number of \"potentially correlated\" outputs (e.g. $S \\geq 20$), one may think of promoting sparsity in the estimation of the attention coefficients. This may allow focusing the analysis on a subset of \"highly correlated\" nodes (outputs) while isolating the ones exhibiting a lesser influence on the others (i.e. attention coefficient close to zero). For instance, if $\\alpha_{i,j} \\approx 0$ for all $j = 1, \\ldots, S$ and $j \\ne i$ (there is no contribution to the other outputs $i$), and $\\alpha_{j,i} \\approx 0$ for all $i = 1, \\ldots, S$ and $i \\ne j$ (there is no relation with the other outputs $i$), then the output $j$ can be considered independent and it can be excluded from the graphical model. In the aforementioned example with 3 outputs, one may consider excluding the third output. Have the authors considered this experimental setup?\n- Promoting sparsity in the attention coefficient values may simplify the complexity of the resulting model but implies overheads in the inference procedure since all the coefficients are jointly estimated. Is there a possibility to adapt the GMOGP inference scheme to sequentially estimate those coefficients? For instance, is it possible to start the GMOGP model with a minimal number of nodes (1-2) and add new ones using a criterion based on the level of the attention values?\n- Can the authors provide further insights about the limitations of the proposed framework and potential solutions?\n\n**Other minor remarks**\n- Punctuation marks in the equations need to be double-checked throughout the paper (e.g. first line Eq.(3), expression after Eq.(6), Eq.(17)).\n- Page 2, Section 2.1: \"A Gaussian process...\" $\\to$ \"A GP ...\"\n- Page 3, Footnote: To begin with capital letter\n- Page 4, after Eq (12): $k_{\\theta_i}$ is a **Co**variance function\n- Page 5, after Eq (16): \"In addition, Applying... \"\n- Page 6, Section 3.4: \"1-d quadrature\" $\\to$ \"1D quadrature\"\n- Page 9, Table 3: to display all the attention coefficients in the same table or to add pictures instead (see Figure 10b, Appendix B.4.1)\n- In the references: \n    - (Brody et al., 2021): already published in the ICLR 2022 (https://openreview.net/forum?id=F72ximsx7C1). The reference needs to be updated.  \n    - To promote uniformity in the references (styles, capital letters, names of the journals and conferences, ...). \n    - (Hensmann et al., 2015): gaussian $\\to$ Gaussian\n    - Capital letters after \":\"\n    - (Williams and Rasmussen, 2006): the order of the authors is not correct.\n    - (Velickovic et al., 2017): already published in International Conference on Learning Representations (https://openreview.net/forum?id=rJXMpikCZ). The reference needs to be updated.", "limitations": "", "rating": 8, "confidence": 4, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "A new construction of MOGP based on probability graphical models is introduced. The graphical model allows capturing dependencies through the definition of the parents of a given output. This makes it possible to define correlations between outputs that are not \"bidirectional\" (i.e. $\\alpha_{i,j} \\ne \\alpha_{j,i}$ for $i \\ne j$), which is not the case for the MOGP-LMC where a symmetric coregionalization matrix is required to promote invertibility. Both the graphical structure (attention coefficients) and kernel parameters are estimated in a variational context where the evidence lower bound of the GMOGP is introduced. The GMOGP inference involves less computational and storage ressources compared to the MOGP-LMC, allowing a better scaling of the model. In addition, the GMOGP can be performed in a distributed learning framework. The transformed GMOGP (TGMOGP) is introduced as well to deal with non-Gaussian priors/likelihoods.\n\nThe proposed GMOGP and TGMOGP are tested on several numerical examples, and results allow assessing their competitivity (and in many cases their superiority) w.r.t. other models from the state-of-the-art. Open-source codes based on pytorch/gpytorch are provided.\n\nThe paper is well-written and well-organized. The main theoretical developments are readable and notions well-motivated. The proposed discussions allows to easily understand how the GMOGP is placed w.r.t. the literature.", "weaknesses": "- The authors have claimed that their framework ensures pareto optimality but there is no numerical evidence. I suggest considering an example with only 2 outputs and adding plots allowing them to validate their claim. They can consider a well-known MOO toy example where the Pareto front has been already studied.\n- Limitations of the proposed framework are not discussed throughout the paper. I believe that scaling the model to large datasets and/or to systems with tens/hundreds of outputs remain a challenge. I refer to **Questions** for further concerns.", "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": "", "extra_scores": {}}, {"review_id": "sD6kJIsBGH", "reviewer": "Reviewer_pYyL", "summary": "This paper delves into the multiple output regression problem, a significant challenge in both machine learning and statistics. It focuses on predicting multiple target variables concurrently. Multiple output Gaussian process regression (MOGP) is an extension of Gaussian process regression (GP) tailored for handling multiple target variables. While MOGP methods provide excellent solutions with uncertainty quantification, they are hindered by substantial computational complexity and storage requirements. The paper introduces an innovative approach to address these complexities and enhance prediction quality by incorporating graphical models and an attention mechanism.", "questions": "Please see the discussions in the Weaknesses section.", "limitations": "", "rating": 8, "confidence": 4, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "In general, the paper has been well-written, particularly in the first two sections where the authors aptly define the problem, establish motivation, and present existing baselines. The foundational concept of addressing the MOGP problem through the lens of graphical models is intriguing, as it offers a means to elucidate the interplay between various target variables.\n\nThe primary contribution of this paper, which involves leveraging distributed learning and attention mechanisms, is commendable given the capabilities of these models. The experiments and numerical results underscore a significant enhancement in prediction quality and execution time when compared to other baseline methods.", "weaknesses": "The current manuscript still exhibits several weaknesses. Notably, in certain sections, the paper's coherence is lacking, making it challenging for readers to follow. Mathematical concepts, in particular, suffer from inadequate definitions. For instance, the paper introduces $f^{(i)}_X$ as a latent variable representing the GP functions for a specific target variable i. However, it later employs $f^{(i)}_X$ to denote the distribution of a GP conditioned on its parents (as seen in Eq. 11 and Eq. 13). Given that prior and conditional distributions hold paramount importance in this work, it is imperative that related variables are precisely defined, and vague expressions are avoided.\n\nAdditionally, a considerable portion of the explanations provided in Section 2 concerning directed graphs appears to be redundant. The authors delve into the intricacies of directed graphs and the definition of a Bayesian network, yet the second paragraph on page 4 reveals that the graph relating the GPs in this paper is a fully connected undirected graph (as indicated by the definition of $pa_i$ for all output functions). It begs the question of why the authors did not emphasize this type of graph from the outset. There are several undirected graph models that can be pertinent when dealing with Gaussian variables, including Gaussian graphical models, pararnormal/nonparametric graphical models, and functional graphs. \n\nIndeed, the description of the distributed Gaussian process (DGP) in the paper is lacking. The conventional DGP model, as proposed by Deisenroth and Ng in 2015, operates under two fundamental assumptions. Firstly, it assumes that all GP experts share the same hyperparameters during the training phase, which serves as a mechanism to address overfitting. Secondly, it assumes perfect diversity between the target variable GPs. However, this second assumption does not align with the scope of MOGP, as interactions between GPs are a crucial aspect of the model. The paper also falls short of providing a clear explanation of how the distributed framework depicted in Figure 4 functions.\n\nTo mitigate computational costs, the model introduces DGP, which can also be applied to conventional MOGP. However, the paper falls short in discussing why this solution is faster than the other baselines. Both the attention mechanism and non-Gaussian transformation inherently raise computational costs, making it essential to have a comprehensive and detailed discussion on this issue within the paper. This is particularly important because it is one of the claims made in the paper. Regrettably, the authors only touch upon this matter in the experiments section, which is not sufficient to provide a thorough understanding of the computational efficiency of their approach.\n\nA similar issue pertains to the prediction quality. The proposed solution integrates a combination of existing methods within a unified framework. However, the question arises: why does this method enhance prediction quality and surpass other baselines? While the experiments demonstrate a significant improvement, the underlying reasons have not been sufficiently elucidated within the paper.", "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": "", "extra_scores": {}}, {"review_id": "tEfKAYgYZV", "reviewer": "Reviewer_X7C1", "summary": "To alleviate the computational complexity and expressiveness of MOGP models, this article proposes a multioutput regression model where each output is modelled as a (marginal) single-output GP, then, these variables are \"tied\" together as nodes of a graph (eq 11). The adjacency of the nodes is then learnt via an attention mechanism (Sec 3.2). The authors also proposed transformed (non-linear) and sparse variants of their method. The article features synthetic and real-world experiments against some MOGP models.", "questions": "Please refer to the comments in the previous section", "limitations": "", "rating": 6, "confidence": 4, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "The idea of connecting GPs with graphs and attention networks is of interest for the community, a number of researchers will be interested in this interface.", "weaknesses": "Although the general idea of building a graph where nodes are GPs is attractive and to some extent promising, the concept is not properly exploited in the article. In this regard, the most relevant weaknesses of the paper are: \n\n- The proposed dependence structure between an output node and its parents is Gaussian (eqs 12 and 13), which in turn results in an overall (MO)GP. Therefore, the proposed model could be understood as a (standard) MOGP with a specific kernel. In addition to the conceptual impact of this fact, there is also a practical consequence: It is unclear if the exhibited experimental superiority of the proposal comes from the introduced graph idea or the fact the there is an implicit, more expressive, kernel being used. From the paper is not possible to clarify this because the proposal is only compared empirically to classic kernels for MOGP (e.g., LMC), leaving behind the MO spectral mixture, deep kernels, cross spectral mixture, etc. \n\n- The proposed model introduced the (hyper) parameters \\alpha as weights learned via attention. These hypers have a direct effect on the linear combination of the mean of the marginal GP which models each node. Though the intuition of attention is that it represents cross relationships among variable of interest, there is no reason to assume that they _linearly_ affect the output. In general, attention networks have an MLP stage after the attention stage, where the MLP \"discovers\" the role of the attention weights. \n\n- The authors claim that standard MOGP have _enlarged searching space with countless local optima_. This is general not true, in GPs (and MOGPs), hyperparameters are usually of very low dimensionality (in particular in the LMC considered in this work as benchmark).\n\n- Diagrams in Figures 3 and 4 are not explained, neither in the body of the paper nor in the Figures' captions. Therefore, they do not contribute to a better understanding of the paper.\n\n- In one of the experiment, MOGP performed worse than single-output GP. This is clearly a training artifact: If SOGP are a restriction of MOGP (same kernel), the latter has to perform at least equal or better than the former, unless MOGP is incorrectly trained. \n\n- The purpose of Experiment 5.1 is not clear: synthetic data are not Gaussian, meaning that the measured performance of the methods considered is not indicative of their applicability for the general case or even for the Gaussian case (a hypothesis shared by all models).\n\n- Given the proposed model builds upon a direct effect among  nodes, how can one deal with the missing data case?\n\n- There are relevant works in the literature that are not mentioned in this paper. A model that is particularly similar is the Semiparametric latent factor model [SLFM](https://proceedings.mlr.press/r5/teh05a.html), which builds an MOGP by linearly combining a set of single-output GPs (rather than relating them through a graph as done in this paper). Also, there have been a number of expressive kernels proposed in the last decade, but this paper only considers LMC (one of the first approaches) and GPRN.", "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": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (spotlight)