Revisiting Nonstationary Kernel Design for Multi-Output Gaussian Processes
Authors: Qiaochu Xu
Zi Yang
Ying Li
Michael Minyi Zhang
Pablo M. Olmos
Venue: ICLR 2026 Poster (2026)
Keywords: Nonstationary kernel
Multi-ouput Gaussian Process
Bayesian non-parametric
Primary area: probabilistic methods (Bayesian methods, variational inference, sampling, UQ, etc.)

Abstract: Multi-output Gaussian processes (MOGPs) provide a Bayesian framework for modeling non-linear functions with multiple outputs, in which nonstationary kernels are essential for capturing input-dependent variations in observations. However, from a spectral (dual) perspective, existing nonstationary kernels inherit the inflexibility and over-parameterization of their spectral densities due to the restrictive spectral–kernel duality. To overcome this, we establish a generalized spectral–kernel duality that enables fully flexible matrix-valued spectral densities — albeit at the cost of quadratic parameter growth in the number of outputs. To achieve linear scaling while retaining sufficient expressiveness, we propose the multi-output low-rank nonstationary (MO-LRN) kernel: by modeling the spectral density through a low-rank matrix whose rows are independently parameterized by bivariate Gaussian mixtures. Experiments on synthetic and real-world datasets demonstrate that MO-LRN consistently outperforms existing MOGP kernels in regression, missing-data interpolation, and imputation tasks.

Published as a conference paper at ICLR 2026
REVISITING NONSTATIONARY KERNEL DESIGN FOR MULTI-OUTPUT GAUSSIAN PROCESSES
Qiaochu Xu
Zi Yang
Ying Lit:1.2
Michael Minyi Zhang
Pablo M. Olmos
School of Computing and Data Science, The University of Hong Kong 2 Department of Statistics and Data Science, National University of Singapore
3 Department of Signal Theory and Communications, University Carlos III de Madrid
{qi Iqiaochu. xu, ziyang2023, lymnli 98 cconnect.hku.hie mzhangl 58hku.hk ying- 8nus. edu. sg remartineings 3c3m.es
ABSTRACT
Multi-oulput Gaussian processes (MOGPs) provide a Bayesian framework for mod- eling non-linear functions with multiple outputs, in which nonstationary kernels are essential for capturing input-dependent variations in observations. However, from a spectral (dual) perspective, existing nonstationary kernels inherit the inflexibility and over-paramettrizzaionn of their spectral densities due to the restrictive spec- tral-kernel duality. To overcome this, we establish generalized spectral-kernel duality that enables fully flexible matrix-valued spectral densities - albeit at the cost of quadratic parameter growth the number of outputs. To achieve linear scaling while retaining sufficient expressiveness, rank nonstationary (MO-LRN) kernel: by modeling the spectral density through a low-rank matrix whose rows are independently parameterized by bivariate Gaus- sian mixtures. Experiments on synthetic and real-world datasets demonstrate that MO-LRN consistently outperforms existing MOGP kernels in regression, missing-date interpolation, and imputation tasks. Code is publicly available at https:/ / / thub. /ate///nnte////////0/
propose the multi-output low-
INTRODUCTION
Gaussian processes (GPs) (Williams and Rasmussen, 2006) offer principled Bayesian non- parametric framework that both flexible and interpretable for modeling complex nonlinear functions. Extending GPs to the multi-output setting leads multi-output Gaussian process (MOGPs) (Alvarez et al., 2012), which, rather than modeling each output independently as in standard GPs, explicitly capture statistical dependencies across multiple outputs. This capability is particularly valuable in diverse application domains, such as early sepsis detection (Futoma et al., 2017), traffic speed predic- tion (Rodrigues et al., 2018), financial risk modeling (Nguyen et al., 2014), and beterogeneouu-ooututt modeling (Moreno-Muñoz et al., 2018).
Formally, MOGP defines priori
vector-valued function f DRD RV: f(x)
with zero mean and matrix-valued covariance (kernel) function K(x1.x2) =
[x)((),,2))]1 x2)]k, E VVV (Alvarez et al., 2012), i.e.,
f(x) ~ GP (0, ((x,, X2))
Given a set of inputs X corresponding outputs
X NXD (distinct from the generic x1, x2 above), the F= ==(((x)): f(x])) NNv
follow joint Gaussian distribution:
p(F) = N(ve)
Equal contribution. Corresponding author.

---

Published as a conference paper at ICLR 2026
where
x)XX = [ X, X)
(,=1 E RNVXXN
is block-structured covariance matrix where each element of the block X (xp: Xz)- Intuitively, the diagonal blocks k. (x,XXX
(X,X) X) E R capture output-specific patterns, while off- diagonal blocks k,y (X. X) E RNX, * j. capture cross-output dependencies. Note that when V 1, the above definition reduces to the standard GP (named single-output GP (SOGP) for convenience) with scalar kernel (x1, X2). Using a standard GP for V-dimensional outputs assumes independence across dimensions, i.e., p(F) 11 I (FF 1,3 0, k, (x, X)) which is equivalent to
XCXX, X) being block .digonna..
key aspect GPs the choice of the covariance function, as it encodes prior assumptions about the latent functions f, which in turn determine the patterns the model can capture in the data. Stationary kernels are the most prominent choice (Ulrich al., 2015; Parra and Tobar, 2017; Alvarez and Lawrence, 2008), where covariance depend only on relative distances X2. Consequently, an MOGP model with such kernel cannot capture input-dependent variations the observations which are prevalent in real-world data.
Thus, designing nonstationary kernels for MOGPs is essential. A straightforward approach is to extend the linear model of coregionalization (LMC) by replacing its SOGP base kernel with nonstationary variant (Bayle and Frean, 2004; Álvarez et al., 2010). However, from a spectral (dual) perspective, this extension still restricts the form cross-output covariances, thus limiting flexibility (Parra and Tobar, 2017). Recently, Altamirano and Tobar (2022) introduced the multi-output harmonizable spectral mixture (MOHSM) kernel, which constructs a two-level mixture spectral density and maps to the kernel domain via existing duality to overcome this problem. Nevertheless, because the duality enforces structural conditions on the spectral density, its design ends up over-parameterized but still inflexible. These limitations carry over the kernel, reducing its ability model nonstationarity and creating redundant parameter space that is difficult to optimize.
In this paper, we propose an expressive and paramettr--fficieet nonstationary kernel for MOGPs, Our contributions are as follows:
1. We establish a new duality between the spectral density and the kernel, removing con- ventional restrictions and thereby enabling, in theory. fully flexible matrix-valued spectral densities, though with quadratic parameter growth V.
2. We propose both a parameter-efffiient and sufficiently expressive multi-output low-rank nonstationary (MO-LRN) kernel by first specifying a low-rank spectral density with in- dependent bivariate Gaussian -ixtture factors reduce parameter growth linear in the output dimension, and then mapping it back to the kernel space via the new duality.
3. We provide experimental results covering tasks such as regression, missing-dee interpola- tion, and imputation, both synthetic datasets and diverse real-world benchmarks, to show that MO-LRN achieves state-of-the-aa (SOTA) performance compared with all existing MOGP kernels.
2 BACKGROUND AND PROBLEM STATEMENT
In this section we examine the design of the MOHSM kernel by taking the next-gen spectral mixture (NG-SM) kernel, which is currently the leading nonstationary kernel for SOGPs, reference point. We first present NG-SM and MOHSM ($ 2.1, $ 2.2), and then revisit the MOHSM design through the lens of the NG-SM's principles, highlighting its key shorteomings via theoretical analysis and empirical validation.
2.1 NEXT-GEN SPECTRAL MIXTURE (NG-SM) KERNEI.
The NG-SM kernel is single output nonstationary kernel that is theoretically capable of approx- imating any continuous kernel arbitrarily well (Yang al., 2025). Its core design principle is to specify dense spectral density and subsequently derive the corresponding kernel function via the well-established duality between kernels and spectral densities. Formally, this duality is characterized by the following theorem:
Theorem 1 (Universal Bochner's Theorem). A complex-valued bounded continuous kernel k k(x), X2) on TR. D is the covariance function of mean square continuous complex-valued random process on
2

---

Published as a conference paper at ICLR 2026
(a) Periodic function
(b) Piecewise function
MOHSM
NG-SM
MOHSM
NG-SM BM-BMM C.I.
Training cata
Figure 1: MOHSM vs. NG-SM fits on scalar observations. Panel (a) shows periodic function, and panel (b) piecewise function. In both panels, the blue line denotes the MOHSM posterior and the red line the NG-SM posterior mean. shaded light blue and pink regions indicate the corresponding 95% confidence intervals. Black dots mark the training data.
RDD if and only if
1(x1, x2) 4 [xxp(iuT>x -ww x2) + exp(iw))x x1 [w! xx))
(1)
exxp(iw] (x1 x2)) + (ip( (x1 xx))]] u(dwi, dwz2,
where u is the Lebesgue-Stiene Stieltjes measure associated with some function p (W1, w2). When u(w1,w)) = 8(wi W2) u(w) (with w1 = w2), this theorem reduces to Bochner's 's theorem.
By specifying the spectral density as mixture of bivariate Gaussians, i.e., p(w),w2) Eg (11 w2), where each (.(w1,,))) w2) a bivariate Gaussian component, the NG-SM kernel can be derived via Eq. (1). Leveraging the fact that Gaussian mixtures are dense in the space of continuous functions (Plataniotis and Hatzinakos, 2000), the NG-SM kernel can approximate any continuous kernel arbitrarily well.
2.2 MULTI-OUTPUT HARMONIZABLE SPECTRAL MIXTURE (MOHSM) KERNEL
The current exemplar for nonstationary nulti-output kernels the MOHSM kernel (Altamirano and Tobar, 2022). It follows the general paradigm of specifying a matrix-valued spectral density and mapping back to the kernel domain via multivariate duality result known as Kakihara's theorem (see Appendix A.I and Theorem 2 in Altamirano and Tobar (2022)). Specifically, by defining = W1 W2 and w (w1 w2)/2. 2,, the (i.j)-th element matrix-valuee spectral density [S is:
Sig = Σ w SE()) (w) SE (w),
where ((9)) are scalar weights and SE SE are squaredl-eponentiall functions (Rasmussen and Nickisch, 2010). Applying Kakihara's theorem yields the (i.j)-th element of corresponding kernel kiý (x1 X2)
Σ g=1
p 1 exp 1x - 0 T (1) (x + cos (* (x T (1)
+ ((5)) exp Li || (2)
X2. x x2) and (0 M Σ .4,,4435 2 4=1 are the kernel hyperpa- rameters. Equation (2) represents a locally stationary obtaine by multiplying a stationary kernel with non-negative modulation function-whin which generally limits its ability to capture strong nonstationarities (Altamirano and Tobar, 2022). To overcome this limitation, Altamirano and Tobar
where &
(0)

---

Published as a conference paper at ICLR 2026
Table I: Simulation settings and runtime comparison for NG-SM and MOHSM
No. Mixture
No. Positional Shifts (P)
Parameter Size
Model
No. Iterations
Runtimes [s]
Reference
Components (Q)
MOHSM
2000
2 N/A
29
63
(Altamirano and Tobar, 2022)
NGSM
2000
11
21
(Yang et al., 2025)
treat (2) as a base kernel ke and localize it by replacing the global Gaussian envelope in x with window centered set of shift points (xp)>= i.e., (x)( X2)
(1) xx xxxxxpll (3) exp (- (xx T EnS (x++))))))) cos ( (a) + exp (-0p-222
The final MOHSM kemnel is then the superposition of the base kernels, (x,, Ep (x), x2), which models nonstationarity by representing the data as a union of locally stationary regimes. Notably, the resulting spectral density exhibits two-level mixture structure: mixture over shift points Xp, P l,..., P, with each shift point further associated with mixture of Q components.
2.3 REVISITING THE MOHSM DESIGN
Despite targeting MOGPs, MOHSM still faces intrinsic limitations in modeling nonstationarity in contrast to NG-SM. To illustrate this limitation clearly and fairly, consider scalar observations, under which the MOGP reduces to a standard GP, and compare the two kernels through theoretical analysis and empirical validation.
Theoretically, when V = 1, the two-level mixture spectral density of MOHSM is neither guaranteed be dense (see analysis in 3.1) in the dual space nor efficient in parameter size (Altamirano and Tobar, 2022). Consequently, MOHSM cannot approximate arbitrary nonstationary SOGP kernels, despite requiring at least times more parameters than NG-SM. In contrast, NG-SM employs dense Gaussian mixture, enabling it to approximate any nonstationary kernel with a relatively small number of parameters.
Empirically, we compare GP regression models with MOHSM and NG-SM kernels on periodic and data with input-dependent frequency transitions (see Figure simulation settings and runtimes are in Table 1). The MOHSM-based model fails to represent global input-dependent transitions, whereas the NG-SM-based model accurately models the data with significantly fewer parameters and substantially lower runtime.
3 METHODOLOGGY
This section begins by establishing a relaxed duality theorem between the spectral and kernel spaces, which loosens the conventional constraints on the spectral density. Building this result, introduce novel nonstationary kernel for MOGPs in & 3.2. Finally, $ 3.3 outlines the procedure for performing MOGP regression inference, later used in the experiments to evaluate the proposed kernel.
3.1 ADVANCED KAKIHARA THEOREM
As discussed in 2.3, the core limitation of the MOHSM design is due the restrictions arising from the original Kakihara theorem (Appendix A.1). ). In addition, the spectral density of NG-SM kernel is constructed using a dense bivariate Gaussian mixture, which cannot be directly adopted in MOGP setting as it would also be subject to the same limitation. We now examine these restrictions in detail. For convenience, we define the measure associated with the spectral density Sij as Fij.
1) F5 must be positive remi-definite (PSD) measure. This requires the diagonal elements of spec- tral density Siz to be PSD function, excluding fully flexible choices such as dense bivariate Gaussian mixtures.
2) Exchangeable within Hermitian symmetry: Fr, (W), W2) = F (2 W1). In MOHSM, ensure the exchangeability, the input of the spectral density set to & = wi - We and w (w) + w2)/2, which leads to both reduced flexibility and over-parumeterzzation. To observe this, recall that
4

---

Published as a conference paper at ICLR 2026
S - ((w) p(w, Xp), where omit its specific functional form for brevity. First, + and w modeled independent, which removes W, interactions and prevents representing general bivariate spectral densities. Sccond, as discussed in $ 2.2, aver-parameterized two-level spectral mixture is required to model nonstationarity, since single mixture layer yields only locally stationary kernel. This local stationarity arises because, after the Fourier transform, @ = (w1+w2)/2 maps (xi + x2) /2 in the spatial domain, producing stationary kernel, while w w1 W2 maps to X X2. which only appears as non-negative weighting function.
To address these constraints, propose advanced version of Kakihara's theorem that relaxes the above conditions, enabling more flexible choices of spectral densities.
funtiiitiion (kig (x1, X2) ) Theorem 2 (Advanced Kakihara Theorem). A family of complex-valued functions
on R serves as the covariance functions of a harmonizable multivariate stochastic process on R if and only if they can be represented as:
(15(X1,X2)
1
// RxR XRD
eexp(i(w]xx w2xx)) + exp(i(w3] w x2))
+ +exp [ww] (x1 - x2)) + exp( (x1 x2)) dFig(W1, w2).
(4)
where F(w1, = Fer( w2 )b is the matrix-valued Sebbeguu--ttelttee bimeasure associated with some matrix-valued function P(W1,w>) = [Pis (Wi, w2)lk which satisfies the Hermitian symmetry condition:
Psy(w1, = Pji(w (2)
Proof. See the proof in Appendix A.I.
We use the symbol P instead ofS S to denote the spectral density in order to distinguish MOHSM and MO-LRN. Moreover, even though using the same notation, the bimeasure F in the two theorems differs in its functional form and associated integral representation.
Remark 1.
With both constraints removed, Theorem 2 implies that Hermitian symmetric matrix-valuee bivariate spectral density determines nonstationary MOGP kemel.. For flexibility, entry P. using a bivariate Gaussian mixture with complex-valued coefficients, while enforcing Hermitian symmetry (Pij Pxx). Diagonal terms are restricted to real nonnegative coefficients, Since the linear span of Gaussian atoms dense (Plataniotis and Hatzinakos, 2000), this parameterizzaion is, in principle, universal over admissible Hermitian symmetric spectral densities and thus can approximate any nonstationary MOGP kernel. However, realizing this flexibility requires O (DV?Q) parameters, rendering the approach impractical.
approximate each
3.2 MUL.TI-UUTPUT LOW-RANK NONSTATIONARY (MO-LRN) KERNEL
Instead of parameterizing each spectral entry with the bivariate Gaussian mixture under Hermitian symmetry, we introduce a both garameter-efficient and sufficiently expressive design of spectral density. Specifically, we first assign a latent vector r; E Fº each output and define
Puy(w), w2) = -r,",,
(5)
for reducing the number of parameters from O(DV²Q) to O(DVQ) while automatically enforcing
symmetry (Pi, = P..). This construction closely resembles latent factor models (e.g. Koren et al., 2009: Mnih and Salakhutdinov, 2007), where pairwise interactions are captured via inner products of low-dimenssioal embeddings. From this perspective, the latent embedding encodes output-specific spectral characteristics, while the relative positions of r, and r, in the embedding space determine the strength and sign of their cross-output spectral interactions.
To retain flexibility, we parameterize the q-th component in the embedding, r; = (1( as:
1Q)TT
1 w (1) .V
W1
Σ
ic
(0)
(6)
w22
(g)
:2
For brevity, the input (w) is sometimes omitted.
5

---

Published as a conference paper at ICLR 2026
where ur (0) > 0 is the component weight, M41 14/ E RD are the mean vectors, and the covariance
(a) (a)
blocks are diagonal matrices Σ/1 diag((o
(9) 11 }²),
(9)
D diag((o (1)))
:2 )²), with
(g) dil o2 E RD.
(g) (q)
The cross-covariance is defined as Σ/0) = Pi diag( (g) diag(o ( where Pj E [-1, 1] denotes
the correlation coefficient.
Substituting (6) into (5) yields:
Piylwi, w2) = r.
Q
(a) N
((1r) w2
(w))
(a) mij
which remains bivariate Gaussian mixture. Here, 2ij mi RRD and S E R2Dx2D are the scaling factor, mean vector, and covariance matrix of the g-the bivariate Gaussian density. These quantities are entirely determined by the paired component parameters 0 (q) in Appendix A.2.
(g) (g)
{wl) (a) M2g k1 k2 Pk (ke(".)) Further details on the parameter derivations are given
To ensure the resulting kernel is real-valued, we set each (wy(w1,w)) = 1 Piy(w1 w2) + (6(-w1--w2)) By applying the duality relation in Eq. (4) of Theorem 2, we obtain real-valued MOGP kernel, whose explicit form is given below.
Definition (Multi-output Low-Rank Nonstationary Kernel). The (1,j)-th element of the MO-LRN kernel (x1, X2) is defined as
Q
(1) 20
cos
151) (m)
(e) T
X1
(9) T (m
192)
X2)
exp (-2
( 1/1 x1 - 2x (5 15c)
x2 + X2 S
9=1
+ cos
(ml X1 (m!)))
x2)
exp -- s bjix x2 2x s x S!
+ o08
(mmlg))T (xx xa)
exp
[ 1 x2)² | (x1 - x2)
+ cos ((mlg) [x1 xe) exp (- [x1 x2) 0 xxx
where each component has parameters (zi) m. (4) ((() with partitions mij (4)
(g m. '32
(a) (9)T
IT
s(s) (s
and SQ9)
all analytically determined from hyperparameters 0, with the explicit
derivation steps given in Definition in Appendix A.2.
Remark 2. Although each element of the spectral density P.j appears to be modeled as a bivariate Gaussian mixture, denseness can only be guaranteed for the diagonal terms (i (i= j}. since only these are independently parameterized. Nevertheless, this construction remains more expressive than MOHSM, whose diagonal spectral terms are not even guaranteed to be dense, as discussed in 2.3.
3.3 MULTI-OUTPUT GAUSSIAN PROCESS REGRESSION
The training procedure for MOGP regression parallels that of the SOGP-by maximizing the marginal log-likelihood of the observed data with respect to the hyperparameters. Formally, we represent each data point as a pair (x,, ) where x,, E D and yn RV. Let ly . NN VV denote the vertically stacked outputs.
We assume that the observations y are generated from the inputs X via noisy MOGP mapping. Specifically.
y = voc(F) tc, C~ C~N(0, Em). vec(F) ~ NO (X,,
(8)
where F) denote the vectorized latent function evaluations over all outputs and inputs, X(X, X) the MOGP prior covariance matrix, and Σn, is the observation noise covariance, taking the form E, = IN & (aag(a?..., with of > 0 is the noise variance specific to the i-th output.
Due to the conjugacy of the Gaussian prior and Gaussian likelihood, we can marginalize out vec(F) to obtain the marginal log- bigliiilhhoo::
1
N log 2n,
logp(y () = - (K+E,,)
log K. + Σn
(9)
2
Maximizing log p(yp(yl] () yields the maximum likelihood estimates of hyperparameters 0.

---

Published as a conference paper at ICLR 2026
Mean Aloolute Error Performance
W
Tralning Comparisents)
Parameter Size Campurisn
MG-EN MOBSM
MI-LRRN CLL.
Test 6a0
MORSM 2P
MORS (ina))
MC-LRN
MOHSM
(a) MOGP regression plot
(b) Comparison on model performance
Figure 2: MOGP regression results comparing MO-LRN and MOHSM on synthetic data. Panel (a) shows regression plots, while Panel (b) reports MAE, training time, and parameter size.
Table 2:
Comparison of various MOGP kernlls. An expressive kernel refers whose spectral density is flexible in the dual space. Here, denotes the index in the LMC summation over latent processes, with V.
Able to model non-stationary pattern
model stationary patter
Expressive kemne
Model
#Parameters
Reference
MOSM CONV
O(DQV)
(Purra and Tobar, 2017)
O(DQV)
(Alvurez and Lawrence, 2008)
LMC-SM
O(IDQ)
(Wilson and Adams, 2013a)
LMC-NGSM
O(IDQ)
(Yang et al., 2025)
MOHSM
O(PDQV)
(Altumirano and Tobur, 2022)
MO-LRN
O(DQV)
This Work
EXPERIMENTAL RESULTS
This section first presents detailed comparison between our MO-LRN kernel and the MOHSM kemel on MOGP regression using synthetic datasets ($ 4.1). We then show that MO-LRN achieves superior performance on real-world MOGP regression tasks 4.2), well on missing-data interpolation and imputation tasks ($ 4.3) across diverse time-series datasets. Comprehensive experimental setups are provided in Appendix B, and benchmark implementations Appendix B.2 and B.3. Appendix C.2 reports additional regression experiments on other real-world datasets.
4.1 SYNTHETIC EXAMPLE
In this subsection, we conduct detailed comparison between the proposed MO-LRN kernel and the MOHSM kernel on MOGP regression task. To this end, we construct two-dimensional output Yn by combining the scalar periodic signal and the piecewise- -erewiii---riodicee signal introduced in $2.3 (see more detail in Appendix B.I.I), Collecting J forms the complete dataset, on which we perform MOGP regression with both kernels.
For MOHSM, performance is highly sensitive to the choice of (P, Q) and requires carefully designed initialization strategies (Altamirano and Tobar, 2022). We therefore compare our kernel with MOHSM across multiple settings, evaluating three metrics: test set accuracy measured by the mean absolute error (MAE), training time, and parameter count. As shown in Figure 2b, our kernel consistently outperforms MOHSM while requiring substantially fewer parameters and significantly shorter training time. By contrast, MOHSM not only delivers inferior accuracy but also remains dependent costly initialization and large parameter budget.
To visualize the performance gap between the two kernels, Figure 2a shows the regression fits obtained using our kernel and the MOHSM kernel under its lowest-MAE configuration (P = 4, Q 4 with the initialization strategy). The results clearly indicate that the MOGP model with our kernel can precisely fit the data and effectively capture nonstationary patterns. In contrast, although the MOHSM-based model outperforms a baseline that applies SOGP regression independently to each dimension (as shown in Figure 1) by leveraging cross-output correlations, still fails to fully capture nonstationary patterns due to its limited spectral density design, thus resulting in inferior overall performance.

---

Published as a conference paper at ICLR 2026
80000
CONV
LMC-NGSM
70000
LMC-SM
MO-LRN (Ouns)
60000
MOHSM
MOSM
Runtime (8)
50000
40000
30000
20000
4
100
150
150
400
450
500
Number of Data Points
Figure 3: Runtime comparison on the ETT dataset against varying number of data points.
Table 3: Comparison of MAE, NMAE, RMSE and NLPD
kernels the ETT dataset. Each column (OT, HUFL. HULL. MUFL, MULL, LUFL, LULL) reports the reconstruction accuracy one output variable, while the last column (Overall) summarizes the average performance across all outputs. Mean and standard deviation are computed over five runs. The best result per column is bolded, and the second-best boott is underlined.
TARGET
OT
HUFL
HULL
MUFL
MULL
LUFL
LULL
OVERALL
METRIC
MAE
CONV
0.276 1. 0.003
0.439 a: 0.140
0.214
0.047
0.455 1: 0.136
0.781
0.119
0.348
0.091
0.220 0.039
0.390 1. 0.073
0.004
MOHSM
0.289
0.005
0.328
0.003
0.251
0.006
0.365 :: 0.004
0.389
0.007
0.359
0.002
0.290 =: 0.001
0.325 4: 0.001
MOSM
0.295
0.024
0.309
+ 0.0%6
0.216
0.027
0.333
0.099
0.459
0.117
0.379
0.022
0.210 0.016
0.314 + 0014
0.372
MO-LRN
0.148
0.003
0.256
0.017
0.125
0.005
0.287
0.023
0.008 0.208
0.003 0.248
0.133
0.003
0.201 1: 0.006
METRIC
NMAE
CONV
0,294
0.038
0.424 =: 0.131
0.229 1 0.033
0.457 1: 0.118
0.659
0.335
0.468 1: 0.090
0,303 :: 0.119
0.409 1: 0.058
LMC-SM
0.319
0.008
0.369
: 0.003
0.266
0.002
0.409
0.005
0,771
0.057
0.488
0.005
0.331
0.040
0.422 1 0.016
MOBSM
0.428
MOSM
0.344
0.028
0.346
0.041
0.258
0.033
0.382
0.044
0.546
0.020
0.452
0.026
0.227 + 0.017
0.365 + 0.016
LMC-NGSM
0,238
0.021
30303 = 0.008
0.199 0.012
0.336 0.002
0.361
0.039
0.016 0.444
0.194 :: 0.002
0.296 :: 0.013
METRIC
RMSE
CONV
0.317 0.018
0.453 2: 0.107
0.252
0.014
0.190 1 0.091
0.754 0.282
0.558
0.169
0.346 0.151
0.153 1: 0.020
LMC-SM
0.340
0.030
0.389 :: 0.035
0,280 1: 0.033
0.439
1: 0.045
0.796
0.079
0.549
1: 0.075
0.376 :: 0.048
0.453 4: 0.045
MOHSM
0.376
1: 0.005
0.407
2: 0.004
0.333
0.005
0.467
0.002
0.005 0.574
0.540
0.001
0.379 0.001
0.440 1: 0.002
0.049
0.431
LMC-NGSM
0.261
+ 0.019
0.344
= 0.000
0.236
0011
0.392
+ 0.003
0.469
+: 0.034
0,322
0003
0.234 :: 0.001
0.350 + 0.013
MO-LRN
0,210
0.003
0.334
0,023
0.004 0.188
0,393
0J030
0,412
0.007
0.353
0.003
0.178
0,002
0.295 1. 0.008
METRIC
NI.P
0.045
0.661
LMC-SM
0,286 0.021
0.429
0.009
0.066
0.012
0.554
0,008
1,283
0.033
0.889
0.002
0.506
0.198
0.573 1: 0.031
MOHSM
0.733 2: 0.022
0.757
= 0.019
0.702
1 0.023
0,809 1 0.016
0.920
0.011
0.881
1: 0.011
0,740
:: 0.021
0.792 1: 0.018
MOSM
0.473
0.073
0.467
0.118
0.208
0.118
0.586
0.116
1,088
0.038
0.846
0.047
0.103
0.075
0.533 4: 0.043
LMC-NGISM
0.146 0.166
0.365
0.046
.8009
0.133
0.020
0,866
0.362
0.271
0.049
0.172
0.382 0.1.53
MO-LRN
00159 - 0.012
0.334 = 0.087
-0.271 1: 0.038
0.505 : 0.097
0.691
0.092
0.385 1 0.007
-0.33
0.014
0.166 2: 0.031
4.2 ELECTRICAL TRANSFORMER TEMPERATURE DATA
We further evaluate our kernel on MOGP regression task using the real-world electricity transformer temperature (ETT) dataset (Zhou al., 2021), which contains one-week records of oil temperature and six load-related features We treat time as the input and the oil temperature together with the six load-related features as the observation 2 The dataset is randomly split into 70% for training and 30% for testing.
For fair comparison, we benchmark against both stationary and nonstationary kemels. The stationary baselines include: (i) the multi-output spectral mixture (MOSM) kernel (Parra and Tobar, 2017), (ii) the convolution (CONV) kernel (Alvarez and Lawrence, 2008), and (iii) the LMC with spectral
See Appendix B. 1.2 for detailed description of the dataset and preprocessing.
8

---

Published as a conference paper at ICLR 2026
PH2.5
2000
1000
TEMP
PRES
0041
0290
5327
1204
00-45
2200
0207
1.902
15-
90% Confidence
Thaining Data
105
Ramoved
Removed interval
Duttu (March 20133
Dabe (March 2013)
Figure 4: Interpolation and Imputation plots on air quality dataset
mixture base kernel (LMC-SM) (Álvarez et al., 2010; Wilson and Adams, 2013b). The nonstationary baselines include: (i) the LMC with a NG-SM base kernel (LMC-NGSM) (Yang et al., 2025; Álvarez al., 2010), and (ii) the MOHSM kernel (Altamirano and Tobar, 2022)². More detailed comparison can be found in Table 2.
Each kernel is evaluated on the test set over five independent runs, and Table 3 reports the mean and standard deviation of MAE, normalized MAE (NMAE), RMSE, and negative log predictive density (NLPD). It can be seen that our proposed MO-LRN kernel outperforms all others across the four metrics, while LMC-NGSM almost achieves the second-best performance. This implies that the ETT dataset naturally contains nonstationary patterns. Consequently, stationary kernels such as MOSM, LMC-SM, and CONV exhibit inferior performance compared to LMC-NGSM, despite it being relatively simple nonstationary MOGP kernel.
For nonstationary kernels, both LMC-NGSM and MOHSM perform worse than our proposed kernel, as they are theoretically limited in expressiveness from spectral perspective. From the viewpoint of spectral density design, MOHSM should, in principle, outperform LMC-NGSM kernel. However, this not observed in practice, likely because its large number parameters makes training substantially more challenging, potentially leading to convergence to suboptimal solution and preventing it from realizing its theoretical advantages.
Figure reports the runtime of the benchmark MOGP kernels on the ETT dataset as the number of data points increases. Our proposed MO-LRN kernel exhibits computational efficiency comparable to existing kernels, with the exception of the MOHSM kernel, whose runtime grows rapidly and quickly approaches the upper limit.
4.3 AIR QUALITY DATA
In this subsection, evaluate kernel robustness to missing data in two tasks: imputation (recon- structing large continuous missing segments) and interpolation (predicting at isolated held-out points). For this purpose, we use data from the Aotizhongxin station of the Beijing multi-site air-quality dataset (Zhang et al., 2017) for March 1-8, 2013, comprising 184 hourly records of eight outputs: six pollutants (PM2.5, PM10, SO2, NO2, CO. 03) and two meteorological variables (TEMP, PRES), with no original missing values (see more detail in Appendix B.1.3). For imputation, we remove different continuous intervals from each pollutant output; for interpolation, we randomly drop 20% of isolated points from all outputs. We then train an MOGP regression model on the resulting datasets containing both types of missing data.
See Appendix B.2 for benchmark implementation details.

---

Published as a conference paper at ICLR 2026
Table 4: MAE for interpolation and imputation tasks across different kemels the air quality dataset. Each column (PM2.5, PM10, SO2, NO2, CO. O3. TEMP. PRES) corresponds to one target variable, and the last column (Overall) summarizes the average performance across all targets. Mean and standard deviation are computed over 5 runs. The best (lowest) values are highlighted in bold, and the second-best are underlined.
TARGET
PM2.5
PM10
S02
NO2
Interp.
Imput.
Interp.
Imput.
Interp.
Imput.
Interp.
Imput.
MOHSM
0.237 + 0.024
0,411
2.0. 108
0.244
0.035
0.331 + 0.103
0:449 + 0.020
0.375
0.116
0.315 = 0.056
0.345 :: 0.153
0.011
0.1002
0.173
0.077
0.160
0.448
0.093
0.332
0.082
LMC-SM
0.116 + 0.006
0.062
# 0.030
0.152
+: 0.004
0
:: 0.003
0.381
+ 0.109
0,684
0.299
0.317
= 0.077
0.796
:: 0.354
MOSM
0,241 + D.41
0.486
0.357
0.213
0451
0.161 + 0.097
0.345
+ 0. 103
0.619
0.180
0.285
0.116
0.434
+ 0.123
LMC-NGSM
0.093
0.006
0.136
0,051 0.000
0.462
0.116
0.397
MO-LRN
0.123 + 0.019
0.109
4 0.079
0.137
0.012
0.103
4 0.016
0.273
0.035
0.484
0.115
0.187
0.028
0.325 0.041
TARGET
CO
TEMP
PRES
Overall
Interp.
Imput
Interp.
Impot.
Interp.
Interp.
Interp.
Imput.
MOHSM
0.310 + 0.022
0.465
0.048
0.385
0.088
0,982 + 0.065
0.347 + 0.015
0.118
0.024
0.301 # 0.011
0.485 + 0.036
CONV
0.328 + 0.044
0.526 4: 0.111
0.421
1: 0.009
0.674 :: 0.015
0.676 +: 0,297
0.041
a: 0.049
0.316 4: 0.042
0.387 :: 0.043
LMC-SM
0.196
0.085
0.483
0.170
1,188
0.110
0.507
0.446
0.296
0.496 +
MOSM
0.314 1 0.018
0.497
2 0.155
0,318
0.040
0.714 1: 0.656
0,450 1: 0.211
0.098
0.022
0.289
:: 0.024
0.485 4: 0.065
NGSM
0.278 4 0.009
0.175
1 0.006
0.449
1 0.186
0.965
1 0.138
0.137
1: 0.007
0.020
0,003
0.250
4: 0.007
0.372
:: 0.037
MO-ERN
0.220 0.034
0.360
0.118
0.328
0.016
0.723
0.054
0.132 0.016
0.022 0.005
0.178 0.008
0.351 0.032
Figure 4 presents the results of imputation and interpolation with associated confidence intervals using an MOGP model equipped with our MO-LRN kernel. For interpolation, the model accurately estimates isolated missing values (orange triangles) with narrow confidence bounds. For imputation, the reconstruction of large continuous missing segments (red squares) yields broad but accurate confidence intervals that fully cover the ground-truth values. Overall, these results indicate that the MO-LRN kernel effectively captures both intra- and cross-output correlations.
Table 4 compares the MAE for interpolation and imputation across different kernels, averaged over five runs. The model with MO-LRN attains the lowest errors in both tasks, excelling in reconstructing large continuous missing segments and isolated missing values. Stationary SM, and CONV cannot adapt to input-dependent patterns, resulting in weaker performance. Among non-stationary kernels, LMC-NGSM lacks the capacity capture cross-output correlations, and MOHSM is constrained by its spectral design, leading to inferior performance in both tasks compared with MO-LRN.
LMC-
5 CONCLUSIONS
In this paper, we introduce the MO-LRN kernel, a novel nonstationary MOGP kernel grounded in new spectral-kernel duality that eliminates restrictive structural constraints and, in principle, permits fully flexible matrix- spectral densities. To avoid the quadratic parameter growth of this unconstrained form, MO-LRN adopts low-rank spectral density with independent bivariate miiusseanmiittue
Gaussian mixture factors, reducing complexity to linear in the number of outputs while preserving sufficient modeling expressiveness. Extensive experiments on synthetic and real-world datasets for regression, interpolation, and missing-data imputation tasks demonstrate that MO-LRN consistently outperforms existing MOGP kernels.
6 ACKNOWLEDGMENTS
Michael Zhang was partially supported by the University of Hong Kong Seed Fund for PI Research #2402101367.
Pablo M. Olmos has received funding from the EU Horizon Europe research and innovationprogram under the Marie Sklodowska- Curie Grant agreement No. 101226456: MLCARE, Machine Learning Computational Advancements for peRsonalized mēdicine. Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the EU or the European Commission (EC). Neither the EU nor the EC can be held responsible for them.
Pablo M. Olmos was also supported by the Comunidad de Madrid IND2024/TIC-34728, IDEA- CM project (TEC-2024/COM-89) the ELLIS Unit Madrid (European Laboratory for Learning and Intelligent Systems), the 2024 Leonardo Grant for Scientific Research and Cultural Creation from the BBVA Foundation, and from the Ministerio de Ciencia, Innovacion y Universidades. Spain, DEITESSAN,,ID222-1121111
under Grant Agreements No. PID2024- 157856NB- EPiCENTER (funded by MICIU/AEV 10. 13039/50 10001 1033 and by ERDF/UE).
10

---

Published as a conference paper at ICLR 2026
Ethics Statement: This work is primarily theoretical, introducing the multi multi-output low-rank non- stationary (MO-LRN) kerell for multi-output Gaussian processes. In the spirit of scientific integrity and transparency, the details of all experiments conducted and the open access datasets used are thoroughly described the paper. Despite its theoretical focus, we recognize that any powerful predictive tool, once deployed, carries potential risks. A key concern is the potential for model using our kernel to generate harmful or misleading information, such as perpetuating societal biases trained on flawed data. The responsibility thus lies with the practitioner audit their data and deployment context, ensuring the application is vetted for fairness and does not lead to negative societal consequences.
Reproducibility Statement: We provide the source code for toy example with corresponding output in the supplementary materials to demonstrate the core functionality of our proposed method. To replicate main experimental results, the Appendix contains all necessary details, including the specific hyperparameters, descriptions of the datasets, and our data processing pipeline. We commit releasing the full, documented codebase on public repository upon acceptance of the paper to ensure complete verification and to support future research.

---

Published as a conference paper at ICLR 2026
REFERENCES
Christopher KI Williams and Carl Edward Rasmussen. Gaussian Processes for Machine Learning. volume 2. MIT Press Cambridge, MA, 2006.
Mauricio A Alvarez, Lorenzo Rosasco, Neil D Lawrence, et al. Kemelss for vector-valued functions: A review. Foundations and Trends in Machine Learming, 4(3):195-266, 2012.
Joseph Futoma, Sanjay Hariharan, Katherine Heller, Mark Sendak, Nathan Brajer, Meredith Clement, Armando Bedoya, and Cara O'brien. An improved multi-output Gaussian process RNN with real-time validation for early sepsis detection. In Machine learning for healthcare conference, pages 243-254. PMLR, 2017.
Filipe Rodrigues, Kristian Henrickson, and Francisco C Pereira. Multi-output Gaussian processes for crowdsourced traffic data imputation. IEEE Transactions on Intelligent Transportation Systems, 20(2):594 603, 2018.
Trung V Nguyen, Edwin V Bonilla, et al. Collaborative multi-output Gaussian processes. In UAI, pages 643-652, 2014.
Pablo Moreno-Muñoz, Antonio Artés, and Mauricio Alvarez. Heterogeneous multi-output Gaussian process prediction. Advances in Neural Information Processing Systems, 31, 2018.
Kyle R Ulrich, David E Carlson, Kafui Dzirasa, and Lawrence Carin. GP kernels for cross-spectrum analysis. Advances in Neural Information Processing Systems, 28, 2015.
Gabriel Parra and Felipe Tobar. Spectral mixture kernels for multi-output Gaussian processes. Advances in Neural Information Processing Systems, 30, 2017.
Mauricio Alvarez and Neil Lawrence. Sparse convolved Gaussian processes for multi-output regres- sion. Advances in Neural Information Processing Systems, 21, 2008.
Phillip Boyle and Marcus Frean. Dependent Gaussian processes. Advances in Neural Information Processing Systems, 17, 2004.
Mauricio Álvarez, David Luengo, Michalis Titsias, and Neil D Lawrence. Efficient multioutput Gaus- sian processes through variational inducing kernels. In Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics, pages 25-32. JMLR Workshop and Conference Proceedings, 2010.
Matías Altamirano and Felipe Tobar. Nonstationary multi-output Gaussian processes via harmonizable spectral mixtures. International Conference on Artificial Intelligence and Statistics, pages 3204- 3218. PMLR. 2022.
Zi Yang, Ying Li, Zhidi Lin, Michael Minyi Zhang, and Pablo Olmos. Multi-view oriented GPLVM: Expressiveness and efficiency. In Advances in Neural Information Processing Systems, 2025.
Kostantinos N Plataniotis and Dimitris Hatzinakos. Gaussian mixtures and their applications to signal processing. Advanced Signal Processing Handbook, pages 89-124, 2000.
Carl Edward Rasmussen and Hannes Nickisch. Gaussian processes for machine learning (GPML) toolbox. The Journal of Machine Learning Research, 1 11:001 1-3011, 2010.
Yehuda Koren, Robert Bell, and Chris Volinsky. Matrix factorization techniques for recommender systems. Computer, 42(8):30-37, 2009.
Andriy Mnih and Russ R Salakhutdinov. Probabilistic matrix factorization. Advances in Neural Information Processing Systems, 20, 2007.
Andrew Wilson and Ryan Adams. Gaussian process kernels for pattern discovery and extrapolation.
In International Conference on Machine Learning, pages 1067-1075. PMLR, 2013a.
12

---

Published as a conference paper at ICLR 2026
Haoyi Zhou, Shanghang Zhang, Jieqi Peng. Shuai Zhang, Jianxin Li, Hui Xiong, and Wancai Zhang. Informer: Beyond efficient transformer for long sequence time-series forecasting. Proceedings of the AAAI conference artificial intelligence, volume 35, pages
2021.
Andrew Wilson and Ryan Adams. Gaussian process kernels for pattern discovery and extrapolation. In International Conference on Machine Learning, pages 1067-1075. PMLR, 2013b.
Shuyi Zhang, Bin Guo, Anlan Dong, Jing He, Ziping Xu, and Song Xi Chen. Cautionary tales on air-quality improvement in Beijing. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 473(2205)20110457, 2017.
Yûichirô Kakihara. Multidimensional second order stochastic processes, volume 2. World Scientific, 1997.
Taco de Wolff, Alejandro Cuevas, and Felipe Tobar. MОGPTК: The multi-output Gaussian process toolkit. Neurocomputing. 424:49-53, 2021.
Pierre Goovaerts. Geostatistics for natural resources evaluation. Oxford university press, 1997.
Taco de Wolff, Alejandro Cuevas, and Felipe Tobar. Gaussian process imputation of multiple financial series. In ICASSP 2020-2020 IEEE International Conference Acoustics, Speech and Signal Processing (ICASSP), pages 8444 8448. IEEE, 2020.
13

---

Published as a conference paper at ICLR 2026
APPENDIX
Advanced Kakihara Theorem and Multi-output low-rank nontationary kernel
15
A.I Proof of Theorem 2
15
A.2 Derivation the Multi-output Low-Rank Nonstationary Kernel
17
Experimental Details
20
B.1 Dataset Description
20
B.1.1
Synthetic Dataset
20
B.1.2
Electricity Transformer Temperature (ETTT)
20
B.1.3
Air Quality
B.2
Benchmark Methods
B.2.1
Linear Model of Coregionalization
................................
22
B.2.2
Convolution Process
22
B.2.3
Multi-output Spectral Mixutre (MOSM) Kerel
22
B.3
Hyperparameter Settings
23
Additional Experiments
23
C.I Dataset Description
23
C.2 Experiment Results
24
D The Use of Large Language Models (LLMs)
24
Future Work
24
14

---

Published as a conference paper at ICLR 2026
A ADVANCED KAKIHARA THEOREM AND MULTI-OUTPUT LOW-RANK NONTATIONARY KERNEL
A.1 PROOF OF THEOREM 2
Before proving Theorem 2, we first recall Kakihara's Theorem (Kakihara, 1997) and introduce lemma that reformulates its conditions the measure Fij in terms of the spectral density.
Theorem 3 (Kakihara's Theorem (Kakihara, 1997)). A famity of complex-vaiie evaludd covariance functions (kig (x1, {2) on R corresponds to a harmonizable multivariate stochastic process if and only each function admits the represenlation
/1 XRD c(w] x: xa]]]x--jaa xx) dFig(wi, W2);
Xx1,
(10)
where F (F) is a matrix-valued spectral bimeasure such that:
(1) For all i, Fú is positive semi-definite (PSD) measure;
(2) For all j. (A, B) = (j,(B, A) for all measurable sets A, B S TRP.
Lemma 1. Let F.g (A, SA JB Sty ky(l1, w2) dwz dwo,, where Siz is the spectral density. The two conditions of Kakihara's theorem, namely
for all Fú is PSD measure;
(2) for all j and measurable A.B.R' F) hold if and only if the spectral matrix satisfies:
A) Sä is PSD for each diagonal entry;
B) Sky(wi, = w1), for all i,
Now, we prove Theorem 2 by showing necessity (=) and sufficiency (+=) as below.
Necessity (=): From Valid Kernel to the Theorem's 's Form.
By Theorem 3 and Lemma 1, (kig (x1, x2)) form the covariance functions of a multivariate harmonizable process if and only if the associated spectral density matrix S satisfies two conditions: (A) its diagonal entries Su w2) are PSD functions, and (B) its off-diagonal entries satisfy the Hermitian symmetry Siy(w1, (w) = Sji(wy w1). Accordingly, we construct as follows:
1 Sij(w,, (z)) = [Pis(w), + P.j (w w) + Piş(w )0(w) w1 + (u)(()))( - w2)l.
(11)
where 8 is the Dirac delta, (1)(w) Pfilw, joint density function, and Py((w))PPy(4n) are marginal distributions. Morcover, we denote the corresponding Lebesgue- Stieltjes measure as F.P. By replacing the Eq, (11) into the duality given by Eq. (10) and Lemma 1, we have
kis (X1, / cllwl xi -w] Sty (wس w2m dw dww?
1 = // Pú Poy(w), We) x wi x2) duy جwي + //
Pip (102, w/ dell xi -nj x2) dw dwo
1 // Ply(w1)((w2 w/ w(w x w x2) dw duez
/f Po)(w2)0(wr www x - x2) dw dwez
+
4
15

---

Published as a conference paper at ICLR 2026
(Swap wi t+ we in 2nd term, integrate in 3rd and 4th terms)
1
// Piy(wr, Psy(wr,wy)elle wy) x1 w x²) dwodwwy +
1 4
/f Pay(wr wyjetlus x1 x2) dw 1dw2
4
1 Pry(wi)e (y(wi)elus أn (x1-x)) dwl + 1 / Peslwa)elus -xa جحلل
(Replacing the marginal density by the joint density for factorization)
= 1 # [Wl]] X1 [] + cllog x1 xsx + oww (x1 x2) + juo (x1 x2) x Pis(wi, w2) dw dwaz.
Finally, by setting the Fy as Sabesgue-SStiettjer measures associated with the spectral density functions (ij(wi, W2), we can express the kernel as:
ki x2) - w2x2)) + exp(i(w221 - w w1x)))
1 1 //
(12)
+ exp(iwe(X1 - xxx) + exp(iwi(Xx x2)]]dFij (w)
(13)
where Pij(w), W2) = (ri(wi, W2). This result is aligned with Theorem 2.
Sufficiency (=): From the Theorem's Form to a Valid Kernel.
Suppose we construct the (i.j)-th element of kernel (X1, x²) using the formula:
ky xx) 1 l/ [/wo x1 w/ x2) + ellW]] x1 w| x2) + Tor (x1 xx) + e w (x1 xa)] dF! (wi, w2)
where the measure F. is associated with a density Ps such that d (W1,W2) = Pij(wi, du d ,dw9, where Pi,(wi, w2)
To establish that kij defines valid kernel for harmonizable multivariate stochastic process, it
suffices to verify that the constructed kernel satisfies conditions (A) and (B) in Lemma 1.
Step 1: Proof of Conditions (A).
By reversing the factorization from the necessity proof, we can write ki as the Fourier transform of spectral density Sij.
1 // P.y ( W2) (ww xi w/ x2) dw dwz
1
(x1, x2)
=
1
// P.y Puy w/ )ell] xi - x2) dw dwws
4
1 f/ Pry(w))((22 wI wwl x wi x2) d dwyr
1 // Piz(wz)ō(w) w2)e "w x1 -m x2) du dwz-
We can then express kig (x,( X2) = |J (ww x1 - xx) Siy(wi (wi) w2) dwzz, where the spectral density is:
1
7 [P((w)),z) w2) + Pij (wس, wI) + (w)))(W> wi) + Pi)(w2)((w) w2)], which is a PSD function for Sa, and thus condition (A) is satisfied.
Siy (wi,
Step 2: Proof of Conditions (B).
16

---

Published as a conference paper at ICLR 2026
It remains to verify the Hermitian symmetry Sty(w1, w2) Sys(w2, w1). We start from the right- hand side:
Sfi(w2, wi) + Pss(wi, w2) + ()i(w2))(w) w2) + Pr((w1)((wz wi)] (
Using assumption that Pej(w, w2) (ji W2). this implies for the marginals that Piy(w) = Pji(w). Substituting these into the expression for (j:(w2, w1):
1 Sty w1) 4
(www, Piy(w2,ws)) Piy(w2)0(w1 w2) + Pi)(w1)8(wz w1)l
= = w1).
The Hermitian condition (B) holds.
Since the constructed kernel Ki possesses spectral density Siy that satisfies the conditions given by Lemma 1, ki is valid covariance kernel for a a harmonizable multivariate stochastic process. This completes the proof.
A.2 DERIVATION OF THE MULTI-OUTPUT LOW-RANK NONSTATIONARY KERNEL
Before detailing the kernel construction, we first introduce two mathematical identities that are essential for the derivation.
Identity 1 (Product of Gaussian). Let Pi( (w) and (w) be two 2D-dimensional multivariate normal distributions:
1
Pi(w)
exp
- (w- m1)) m1)) Σ '.:(w... (ww-))))
(1
T (w 14,)) 2
P,(W)=
( 1/ cxp
(w
Then, the product of these two densities proportional to another multivariate normal distribution:
Ps (w (w) 2ιj ((w miy. Sty)
where the resulting mean, covariance, and normalizing constant are
Sij (2, + " mij Stg E :4 + Σj Ht))
exp {
T,)TT 1 -m; m1) (E + Σ,) - -1,y) }
2/j 2mT) Σ Σg
Identity 2 (Characteristic Function). For a random vector z ~ (u, E), its characteristic
<(t) Ele E 9= 6 ((2l,, exp i41 1 Σt
1
This formula enables us to bypass direct integration and obtain closed-fome analytical result.
Step Spectral Density Construction.
As discussed in $ 3.2, we propose spectral density design that is both varameter-efficient and expressive. For each output we assign a latent vector r, E TRQ and define (w) w2) r r,- And each component ((4) of r, is modeled as bivariate Gaussian distribution:
(0) (1س) w N
((w)
Ww
where (0))) E R is positive weight, and the mean vector and covariance matrix are partitioned as:
(c)
(g)
(g)
11
Σ n
ie
17

---

Published as a conference paper at ICLR 2026
The covariance matrices are parameterized by 2/1 diag( )2) 11
x/ =diag((
(0) (0) and the
12
cross-covariance by (g)) = p diag(o (1)) diag(o 12
(4)), where o il :2
(4) (q)
E RD, and ( E 1,
denotes the correlation coefficient.
By Identity 1, the spectral density has the following form:
Q Σ
Poy(w1, W2) -r
9=1
Σ
p(e) ij
(w1
=
Σ
ij
(2)
N
(wi)
W2
m.y
g=1
q=1
where the parameters of each mixture component derived as:
5/ (1x (2
(4) ( (E (10)) + ) )/))))
ij mij
2j w w (2t) (2T) M(10 Σ - exp { ( (44 (0) - (0)
(4)
T
(1)
+ ) - - 20))
Since we only focus on real-valued kernels, we eliminate the imaginary part by enforcing a symmetric spectral density:
P.j (11,w))
p pyymm.()) (w) W2)
9=1
Q
p(1) (W), w2) + p w/2
2
Q
ij
W
W1
m (
+ +N
((w)
m 1j
((9)
)
2
W2
W2
-=1
Step 2: Transform Spectral Density to the Kernel Domain.
We now transform spectral density into the kernel domain by Theorem 2, yielding the (t.j)-th element of kernels ki (x1,X))
= 1 //
[ wxx + cllee] x1 -wi x2) + e! (xi -x2) + cww (x1- xxa] pyym dwprios
4
Q
Σ
//
[[w]] xi تس x2) + قساا x1 mn x2) إسسي (xi - x2) + وسسي (xi - xa)] poy (symm. dw dwo
4
Q
Σ
c [(w] xi - fm (ex p psymm.(a) dw dww + /f ((-ww
x/tw] xetwk x1) paymm.(a) dw dww
Torm 1 (1(00)
Term (1
+ 1
[n (xi x2) (4) du dww2 /f jn (xi xx) p paymm.(a)) dw ddwz
Term (1!11)
(1f1))
18

---

Published as a conference paper at ICLR 2026
In fact, terms 1-4 correspond to the characteristic function of Gaussian distribution. Specifically, we can derive their general analytic form by setting (w)) and t.
(4)
/ Nw'' (2ymm.(9) (w) dw = e2 ut
/
t
Rig
[(wlm (g) ( + +Wwl m. ) du
2
(4)
1j
f w
* (wwlm (2) /dw ydw +
f -mij (0) )dw )dw
2
(9)
1 ---t st + +exp (-i(mL)) 1 s ))
exp ( Tt
(By Identity 2)
2
(n)
cos (ml))
(
1 (-2e- - s
(0)
2ij
exp
Thus, we can use this result to derive the analytic form of terms 1-4 as following.
Term 1: By setting
= =( The result is:
xx]
(0)
((90
3 cos (m/
(m) x1 (m.))) x2) exp
1 [x 2x ( x2 + x s s/12x2)
Term 2: By setting t = (~x;)) The result is:
24 cos (mm
(0)) T x1 - (m
(1) ( T x2) exp (-2(xx (x] 1/1 - 2x] |/5x + x S(( (
Term 3: By setting t= = (x1=xx)) The result is:
100
#) cos (m/)) (2) T (xi xe) exp (-2(x) x2) S -xa))
(
Term 4: By setting t = =xxxxxx00 (x1-x)) The result is:
(n)
10
cos
(m)
(1)
(xi x2) exp
- 2
-x1 (x1 - x2) xx s -xa))
4
where m (g)
(9) and
((()
S(9) ffl
(S
ije (g)
m.
sje
s
ij2
Thus, the explicit form the q-th component of kiy (X1, X2) can be derived by combining four terms. Finally, summing over all Q components yields the MO-LRN kernel as below.
Definition (Multi-output Low-Rank Nonstationary Kernel), The MO-LRN kernel for Multi-output
(q) s.()) ) q {l..... Q}. The parameters Gaussian Process with component parameters (2()) .m are partitioned as:
(
(f (
m (g)
and 11 =
81
Gje)
S
S/g
tje
The resulting cross-covariance kernel kig (x1, X2) is given by the explicit formula:
Q (9)
kig (x1, X2)
2gj
cos
(4) (ai
exp
3/1())
1
+ +cos
a2
(g)
exp
2
9=1
(g)
1
a(4)
(4) a
g(1)
+ +os
ag
exp
2 Dg
+ cos
exp
2
19

---

Published as a conference paper at ICLR 2026
(9) where the auxiliary functions Oik and { are defined for each component g as:
a1 (m (g) T X1
(m (g)
t)1l
192) X2:
8 x 191 X1 2 (S(1) (1))
s(9)
X2 + X 5(1)
ijc
22
(9) 10) (0) T (mif1) X2:
22 (m
3 x2 1i]
s(9) - 2x S (19)
iic X2 + x 5
43 m (0) T [xx -x2),
(
94 (x1 x2) [ (x1 -x2),
( (4) T (x1 X2),
01 (m)
[x X2). 8/00 (x1 x2)) T 122
The parameters of the kernel are ultimately determined by the set of hyperparameters 0 (q)
{wg Pil 22 (g)
(g) (4) Q.V
1g=1.=11
B EXPERIMENTAL DETAILS
B.1 DATASET DESCRIPTION
This subsection details the datasets employed in our experiments and the comprehensive preprocessing method.
B.1.1 SYNTHETIC DATASET
The synthetic data is generated from two functions that exhibit
variation:
1) Piecewise Function:
(2nx) + (x) sin (a + (sin(sin66TT)) + (in())))) + if 2
if -1
if ff=10
(14)
2) Periodic Function:
(x) s(n))inn(2)) f(z) +
(15)
where A(x) 0.5 + 1.5 E the amplitude function that increases with x. f(x = - the frequency function that decreases with x, and ~ (0, 0.1) represents Gaussian noise.
The final dataset comprises 400 ouput pairs over the domain æ E [0, 5]. where each input xn associated with the two-dimensional output y, (xx),,,,, f (xx)) We randomly split the dataset into 70% for training and 30% for testing, and standardize both output dimensions before model fitting.
B.1.2 ELECTRICITY TRANSFORMER TEMPERATURE (ETT)
The electricity transformer temperature (ETT) dataset' (Zhou et al., 2021) records transformer data over two years (July 2016-June 2018) 15-minute intervals, yielding 69,680 samples with seven variables: oil temperature and six power load measurements. For our experiments, one-week subset (July 1-7, 2016), giving 670 samples across all seven variables. We randomly split the data into training (70%) and test (30%) sets with shuffling, All variables were standardized to mean and unit variance. The timestamps were converted into minutes and rescaled to [0, 600] using min max normalization, preserving relative time intervals while reducing scale. Each standardized variable was treated as an output channel, enabling multi-output modeling and joint prediction of all seven transformer metrics.
selected
hhttps:// thub. com/ wtt..com//houuaaaayoyoooa// httph/ohoohooooooooo Brtoouuaoy///lbattte... /Pitatasee
20

---

Published as a conference paper at ICLR 2026
Table 5: Hyperparameter settings for the MO-LRN kernel.
Kernel Sctup
Parameters
Synthetic Data
ETT
Air Quality
Mixture densities (Q)
2
Optimizer Setup (Adam)
Parameters
Synthetic Data
ETT
Air Quality
I Learning rate
0.10
0.02
0.01
Weight decay
10-5 5
N/A
# Iterations
500
4000
4000
Table Consolidated hyperparameter settings for the ETT and Airquality datasets across various models.
ETT Dataset
Model
Q
Ir
# Iterations
LMC-NGSM
4
N/A
4
0.02
4000
LMC-SM
4
N/A
4
0.02
4000
CONV
N/A
N/A
4
0.02
4000
MOSM
N/A
N/A
4
0.02
4000
MOHSM
N/A
4
0.02
4000
Airquality Dataset
Model
Ir
# Iterations
LMC-NGSM
4
N/A
4
0.01
4000
LMC-SM
4
N/A
4
0.01
4000
CONV
N/A
N/A
4
0.01
4000
MOSM
N/A
N/A
4
0.01
4000
MOHSM
N/A
4
0.01
4000
B.1.3 AIR QUALITY
The air quality dataset (Zhang et al., 2017) contains hourly measurements from Beijing monitoring stations between 2013 and 2017. used one week data (March 1-8, 2013) from Aotizhongxin station, yielding 168 hourly samples with eight features: PM2.5, PM10. SO2, NO2, CO, 03, temperature, and pressure. Timestamps were converted to hours since March 1, 2013, and rescaled to [0, 160) using min-max normalization. All features were standardized to zero mean and unit variance. For interpolation, we randomly removed 20% the data points across all variables. For imputation, removed contiguous intervals in different channels: PM2.5 (0-20h), PM10 (25-45h), SO2 (50-70h), NO2 (75-95h), CO (100-120h), and 03 (125-145h).
B.2 BENCHMARK METHODS
We present detailed discussion of the theoretical foundations of the benchmark methods summarized Table 2. For implementation, we rely on the MOGP Toolkit (MOGPTK)* (de Wolff et al., 2021).
'https://arch" ics ies.uci......................... www...ewwww.ttteeeeeeeeeeeeaee /bei /ngmjingemmul t +aitetai quality+data httpps://ii thub. wtt.:/c.o//ammcoa.e.... 1e/mogptk
21

---

Published as a conference paper at ICLR 2026
B.2.1 LINEAR MODEL OF COREGIONALIZATION
The linear model of coregionalization (LMC) (Goovaerts, 1997) represents each output fx(x) as linear combination of latent Gaussian processes:
(x)
d ug (x).
(16)
Here:
ug (x) independent latent Gaussian processes within the same group share the same
covariance function k; (.,).
d/ are the scalar coefficients that mix the latent functions create the observed output [x(x).
I is the number of latent process groups (each associated with a different covariance function), and J; the number of latent processes within each group
From the perspective of MOGP, the kernel function of LMC is given by:
K(x1, x2) = fov(((x), = Σ (,,
(17)
where B, a x V positive semi-definite matrix determined by the coefficients o captures the correlations between the outputs that are explained by the i-th group of latent functions. The rank of B, is Jr, which controls the complexity of the correlation structure for that component.
B.2.2 CONVOLUTION PROCESS
The process convolution method (Boyle and Frean, 2004) offers alternative way generate correlated outputs. Each output (x( is obtained by convolving a shared latent GP u(x) with a smoothing kernel (x). Intuitively, this corresponds taking single underlying random process and "smoothing" it differently to produce the observed outputs.
Formally, the model for the e-th output is expressed convolution integral:
fx(x)= = / G - z)u(z)dz,
(18)
where
u(x) is a shared latent with kernel function k(-,)). Gy the smoothing kernel specific the w-th output.
Based on this convolutional structure, the cross-covariance between fo, (xp) and Joy (Xg), (i.e. kij (Xp, x,)) is given by
J [ / (x, -2)Gey 2')k(2, d'add''.
lig(X xx)
(19)
B.2.3 MULTI-OUTPUT SPECTRAL MIXUTRE (MOSM) KERNEL
The multi-output spectral mixture (MOSM) (Parra and Tobar, 2017) kernel constructed by designing the spectral density and then mapping it back to the kernel space via the well-established duality for stationary MOGP kernel. form is given by:
(Ty =
exp (- +
T + -10)) cos ((T ((9) +
(4) 2
(20)
9=1
where 0g w i
(g)
(2x)
m/2 ( 1/2 and the superindex (-) denotes the parameter of the 9th
component of the spectral mixture.
22

---

Published as a conference paper at ICLR 2026
Model Performance and Runtime
1000
0.019
Runtime (s)
0.018
800
Runtime (s)
NMAE
600
0.016
0.015
200
trade-off
Number of Spectral Components Q
Figure 5: Model performance and computational cost with varying numbers of spectral components Q in MO-LRN. The left y-axis shows the NMAE, and the right y-axis shows the corresponding runtime.
B.3 HYPERPARAMETER SETTINGS
To select a reasonable value for Q. we evaluate both runtime and NMAE on the GOUN dataset as Q varies (see Appendix C.1 for details), and present the results in Figure 5. As Q increases, the model achieves slightly better predictive accuracy (lower NMAE), but at the cost of substantially longer computation time. Particularly, the accuracy gain beyond
4 marginal, whereas the runtime grows rapidly, suggesting that Q 4 provides favorable balance between accuracy and efficiency.
In addition, to ensure a fair comparison, we the number mixture components the same value across all competitors, Moreover, in the toy example, MOHSM is also optimized for 500 iterations for consistency. The configurations for our model and the baselines are summarized in Table 5 and Table 6, respectively.
C ADDITIONAL EXPERIMENTS
In this section, evaluate the proposed MO-LRN kernel MOGP regression tasks using additional datasets, all obtained from the MOGPTK package (de Wolff et al., 2021). We further include runtime analysis on the ETT dataset with varying sample sizes demonstrate the efficiency of our model.
C.I DATASET DESCRIPTION
We first provide brief description of the datasets and then present the corresponding results.
GONU dataset: real-world financial dataset consisting weekly prices for gold, Brent crude oil, the NASDAQ composite index, and a broad USD index (de Wolff et al., 2020). Details are follows:
Oil (USD): Europe Brent spot price from the U.S. energy information administration (EIA).
USD Index: Federal Reserve trade-weighted nominal broad U.S. dollar index (FRED).
Gold (USD): London bullion market association (LBMA) gold price at 10:30 a.m. London time, distributed via FRED.
NASDAQQ Adjusted daily closing prices of the NASDAQ Composite index (symbol iXIC, Yahoo Finance).
EEG dataset: Electroencephalograaly (EEG) recordings from human neonates. Multi-channel EEG was recorded from 79 term neonates admitted the neonatal intensive unit (NICU) at the
23

---

Published as a conference paper at ICLR 2026
Helsinki University Hospital. The median recording duration was 74 minutes (IQR: 64 to 96 minutes). For our experiments, we use eight of the available twenty-two sensors as outputs: Fpl, Fp2. Fz, Cz, T3, T4, 01, and O2. These correspond to standard electrode positions in the international 10- EEG system, covering frontal, central, temporal, and occipital regions. This results in an eight-dimensional multi-output regression task, enabling the study of cross-correlations between EEG channels.
Currency dataset: Daily currency exchange rates with respect to the U.S. dollar, collected over the two-year period 2017-2018. We consider ten major currencies in addition the U.S. dollar reference: Canadian dollar (CAD), Euro (EUR), Japanese yen (JPY), British pound sterling (GBP), Swiss frane (CHF), Australian dollar (AUD), Hong Kong dollar (HKD), New Zealand dollar (NZD). South Korean won (KRW), and Mexican peso (MXN).
Bramble dataset: A real-world environmental dataset of tidal height measurements from four coastal weather stations in South England: Bramblemet, Cambermet, Chimet, and Sotonmet. The stations have been continuously recording since April 2012, with measurements every five minutes. For our experiments, we use the tidal height data from one-week period in June 2020. To improve trainability, the data are detrended before model fitting, and 90% of the data points are randomly removed to accelerate training and reduce memory requirements.
C.2 EXPERIMENT RESULTS
All datasets are randomly partitioned into 70% training and 30% testing splits. Evaluation is conducted on the test set over five independent runs. Table 7 reports the NMAE results on the GONU dataset, with corresponding regression plots shown in Figures 6-11. Results on the EEG dataset are summarized in Table 8 with regression plots in Figures 12-17. For the currency dataset, NMAE scores are presented in Table 9, 10 and regression plots in Figures 18-23. Finally, Table 11 provides the NMAE results on the Bramble dataset, with regression plots in Figures 24 24-29.
The proposed MO-LRN kernel consistently outperforms all baselines on the cross-domain multi- output datasets, attaining the lowest overall NMAE and delivering superior or comparable accuracy across individual outputs. The LMC-NGSM kernel typically emerges the second-best performer. These findings indicate that non-stationary patterns are widespread in time-series data, and capturing them is essential for effective regression. In contrast, stationary kernels such MOSM, LMC-SM, and CONV show inferior performance, as they lack the ability to capture non-stationary patterns.
Among the non-stationary kernels, both LMC-NGSM and MOHSM fall short compared to proposed approach, as their spectral formulations inherently restrict expressiveness. In particular, LMC-NGSM constructs kernels linear combinations of NG-SM components, forcing auto- and cross-covariances to share the same structure. This restriction limits its ability to capture complex cross-output interactions. Although MOHSM is theoretically more flexible than LMC-NGSM, its excessive parameterization makes training difficult, often leading to suboptimal convergence and weaker empirical performance.
D THE USE OF LARGE LANGUAGE MODELS (LLMS)
We used large language models (LLMs) only for language polishing and readability improvements.
E FUTURE WORK
In flexible MOGP kernels (Parra and Tobar, 2017: Altamirano and Tobar, 2022), scalability remains an open challenge. One potential direction is the use of inducing point methods, which offer substantial computational savings. However, key limitation is that they often overlook redundancy from repeated input locations across outputs, leading to inefficient representations and possible numerical instability. Another promising approach random Fourier features (RFFs). While RFFs provide a straightforward route to linear-time approximations, applying them to complex MOGP kernels poses difficulties, particularly in ensuring the positive definiteness of the joint kernel matrix. In future work, plan to address these issues to enhance the scalability of flexible MOGP kernels including MO-LRN and extend their applicability to a broader range of tasks.
24

---

Published as a conference paper at ICLR 2026
Table 7: Comparison of NMAE across kernels on four outputs (Oil, Gold, NASDAQ. USD) of GONU dataset.
Mean and standard deviation are computed over five runs. Best (lowest) bold, second-best is underlined.
MODEL
Oil
Gold
NASDAQ
USD
OVERALL
METRIC
NMAE
CONV
0,02333 :: 0.00035
0.06135 = 0.00092
0.03221 = 0.00050
0.008239 + 0.00012
0.03128
0.00047
LMC-SM
0,02057 = 0.00031
0,05301 = 0.00080
0.01966 = 0.00029
0.009844
0.000155
0.02577 = 0.00039
MOHSM
0,04836
0.00073
0.12455 0.00187
0.04644
0.00070
0.020996
=0000332
0.06059
0.00091
MOSM
0,01363
+ 0.000200
0,05293 + 0.00079
0.02515
0.00038
0.005669 + 0.00009
0.02485
0.00037
LMC-NGSM
0,01747
0.00026
0.02519 + 0.00038
0.01559
0.00023
0.005906
0.000099
0.01654
0.00024
MO-LRN
0.01218
0.0M023
0.02513 = 0,00195
0.01439
0.0M2278
0.004996
0.000632
0.01467
0.00723
Table 8: Comparison of NMAE across kernels on the eight outputs (Fpl. Fp2, Fz. Cz, T3, T4, 01. 02) of EEG dataset. Mean and standard deviation five runs. Best (lowest) in bold; second-best underlined.
MODEL
Fp!
Fy2
Fr
Cx
T4
01
02
OVERALL
METRIC
NMAE
QUNV
0005
10.85
00013
0,502
0.742
(0,00)
0.736
0.520
1MC-SM
0.846
0.006
0.876
0.007
0.777 00011
0.986
0.005
0.593
0.003
0.704
0.001
0.928
0.009
0.613
0.006
0.791
000006
MOHSM
0.006
0.820
0812-
0.002
0.680
0.001
0.670
0.001
0735
0.008
0.670
0.005
MOSM
0.576 0.006
0.83
1.000
0.770
1 0.013
0.10
1 0.005
0.702
0.00
0.654
0.003
0750
0,008
0.652
0.005
0.731
0.000
LMC-NGISM
IL548
+ 00004
0.551
+ 5..0
0.14
= 0.010
1313
= 0.309
0.390
0.083
0.302
0
03003
0.008
0.267 = 11003
0.61
= -.0005
MO-LRN
I.502
0.084
0.503
2 0.084
0.579
0.009
0,468
= 0.004
0.527
0.004
0.354
= 0.003
0.517
0.005
0.259 = 0.003
0.464
0.005
Table 9: Comparison of NMAE across kemelss five curreney outputs (EUR, CAD, JPY, GBP, CHF). Mean and standard deviation over five runs. Best (lowest) in bold; second-best underlined.
MODEL
||
EUR/USD
CAD/USD
JPY/USD
GBP/USD
CHF/USD
METRIC
NMAE
CONV
0.00563
0.00017
0.00683 + 0.00020
0,00721 + 0,00022
0.00770
0.00023
0.00593 = 0,00018
LMC-SM
0.00653
0.00020
0.00857 1 0.00026
0,00742 1 0.00022
0.00732
0.00022
0.00702 1 0.00021
MOHSM
0.29028
0.00871
3.51393 + 0. 0.10542
0.01158
0.00035
0.07513
0.00225
0.08108 0.00243
MOSM
0.0041
0.00012
0.01326
0.00040
0.01415
0.00042
0.00917
0.00028
0.00629 + 0.00019
LMC-NGSM
0.00548
+ 0.00016
0.00445
= 0.00013
0.00713 = 0.00021
0.00921
: 0.00028
0.00386 = 0.00012
MO-LRN
0.00239
0.000077
0.00989
0.00030
0.00923
0.00028
0,00731
0.00022
0.00221
0.00007
Table 10: Comparison of NMAE across kernels five currency outputs (AUD, HKD, NZD, KRW, MXN) plus OVERALL. Mean and standard deviation are computed over five runs. Best (lowest) in bold; second-best underlined.
MODEL
AUD/UDD
HKD/USD
NZDHUSD
KRWAISD
MXNASD
OVERALJ
METRIC
I|
NMAE
CONV
0.00672 # 0.00220
0.000336
= 0.000011
0.00669 +: 0,00020
0.00635
0.00019
0.01344 =: 0.00040)
0.00669
0.00007
LMC-SM
DO0733
DO0022
0.001366
-0.00844
0.00836 + 0,00025
0.00834
D.0KI25
9.00525
+ -0.00016
9.00675
0.D00077
MOHSM
0.02322
0.00237
0.00007
0.01785
0.48130
=..00472
MOSM
0.00530
0.00166
0.00026
0.00001
0.00595
1 0.00018
0.003304
0.00009
0.00686
0.00021
0.00684
0.000077
LMC-NGSM
0.00543
1 0.00016
0.00034
0.00001
0.00812
1 0.00024
0.00774
1: 0.00023
0.006R5
0..0021
0.00586
0.00005
MO-LRN
0,00513
0.00024
0.00012
0.00003
0.00326
0,00010
0.01084
0.00033
0.00515
0.00015
0.00335
0.09006
Table 11: Comparison of NMAE across kernels on the four outputs (Bramble, Camber, Chi, Soton) of bramble dataset. Mean and standard deviation computed over five runs, the best (lowest) value is bold, and the second-best underlined.
MODEL
Bramble
Camber
Chi
Soton
OVERALL
METRIC
NMAE
CONV
0.0825
0.0025
0.0621
0.0019
0.0599
0.0018
0.1071 0.0033
0.0779 0.0024
LMC-SM
0.0791 + 0,0024
0.0651 + 0.0020
0.0590
0.0018
0.1055 = 0.0032
0.0772 = 0.0024
MOHSM
211.530 + 2.115
4.027 + 0.040
0.609 = 0.006
1.269 + 0.013
54.858 + 0.556
MOSM
0.0831
0.0025
0.0617 0.0019
0.0594 = 0.0018
0.1082 1 0.0032
0.0781 0.0024
LMC-NGSM
0.0504
0.0015
0.0612 0.0018
0.0690
0.0021
0.0642 0.0019
0,0612 0.0018
MO-LRN
0.0497 0.0015
0.00977 + 0.00029
0.0150 + 0.00045
0.0183 = 0.00055
0.0237 + 0,00071
25

---

Published as a conference paper at ICLR 2026
NASDAQ
1400
1300
1200
1300
130
115
2015/07 2016/11 2017/67 2017/01
2017/07 2018/01 2018/07
2015/07 2016/01 2017/07 2017/01 2017/07 2018/01 2018/07
Train data
Test data
Prediction
95% Confidence
Figure 6: MOGP regression plots on the GONU dataset with the MO-LRN kernel.
Gold
NASDAQ
1300
7000
6200
1100
5200
Oll
USD
140
100
2015/07 2016/01 2017/07
2015/07 2017/07 2017/0L
Train data
Test data
Prediction
95% Confidence
Figure 7: MOGP regression plots on the GONU dataset with the MOHSM kernel.
Gold
NASDAQ
1400
1300
6200
1300
5200
4200
130
115
2015/07 2016/01 2017/07 2017/01 2017/07 2018/01 2018/07
3015/07 3015/01 2017/07 2017/011 2017/07 2018/01 2018/07
Train data
Test data
Prediction
95% Confidence
Figure 8: MOGP regression plots the GONU dataset with the MOSM kernel.
26

---

Published as a conference paper at ICLR 2026
Gold
NASDAQ
1400
1300
1200
6200
5000
4000
130
125
120
Train data
Test data
Prediction
95% Confidence
Figure 9: MOGP regression plots on the GONU dataset with the LMC-NGSM kernel.
1400
1000
1300
1000
5000
1300
O1
USD
00
120
2015/97 2016/01 2012/97 2012/01 2017/07 2018/01 2018/77
2015/07 2016/01 2017/07 2012/01 201 2017/07 2018/01 2018/07
Train data
Test
Prediction
95% Confidence
Figure MOGP regression plots on GONU dataset LMC-SM kernel.
Gold
NASDAQ
1300
7000
1200
1300
5000
1207
Oil
130
315
310
2015/07 2015/01 2017/07 2017/01 2017/07 2015/01 2018/07
2015/07 2016/01 2017/07 2017/01 2017/07 2015/01 2015/07
Train data
Test data
Prediction
95% Confidence
Figure 11: MOGP regression plots the GONU dataset with the CONV kernel.
27

---

Published as a conference paper at ICLR 2026
Bed
04
0.
-0.
-008
EGG O1-Ber
ZEG 02-Net
850 T3 Bef
EEG T4-Ret
Train data
Predict on
Tist cotta
95% Co fffeeccc
Figure 12: MOGP regression plots on the EEG dataset with the MO-LRN kernel.
BEG
data
Figure 13: MOGP regression plots on the EEG dataset with the MOHSM kernel.
EEG Cadl
SEG Fp2-RAt
EIG Fp2-der
eto rantel
03-Btl
#80 02 Bef
12-RM
Prediction
Confidence
Figure 14: MOGP regression plots on the EEG dataset MOSM kernel.
28

---

Published as a conference paper at ICLR 2026
Bedi
Spt.tee
100
140
AM
EEG 01-Bet
ZEG 02-Net
Train data
Prodiction
95% Co confiiencc
Figure 15: MOGP regression plots on the EEG dataset with the LMC-NGSM kernel.
Fol-kef
Fp2-Ref
036
0.
EEG T3-3at
R60 T6.
Prediction
TAcc datta
Figure 16: MOGP regression plots on the EEG dataset with the LMC-SM kernel.
reG Ceffitt
EE0 Pal.Bet
Pa2.fer
#60 Pzchel
BBG O2-Bet
B1G Th-bel
EEC T4-Ret
Train data
Prediction
55% Confideee
Figure MOGP regression plots on EEG dataset with the CONV kernel.
29

---

Published as a conference paper at ICLR 2026
135
1.0
120
AUN/USO
G8PUUSO
4.0
0.09
JPYNSD
KRWASS
1.00
NZDAUSO
140
data
55% Confidance
Figure 18: MOGP regression plots on the currency dataset with the MO-LRN kernel.
USD
P7/05D
m133
NZD.SD
95% Confidence
Figure 19: MOGP regression plots on the currency dataset with the MOHSM kernel.
AUD/USO
117
115
CBPNUSD
HKDNS5
3.0
JPYUUSD
MOONUGO
1186
MZD/USD
10
Prodiction
95% Confidence
Figure 20: MOGP regression plots on the currency dataset with the MOSM kernel.
30

---

Published as a conference paper at ICLR 2026
AUDURD
CADUSO
CHFAUSO
179
EURJUSO
OBANUD
0.0
230
6.0
JPEAUSD
1130
M20
45% Confidence
1.0
Figure 21: MOGP regression plots on the currency dataset with the LMC-NGSM kerml.
CADUUD
LB
Ph Rack
135
L20
EURUSO
GSP/USD
370
JPY/US
1055
1100
Text
90% Confidence
Figure 22: MOGP regression plots on the currency dataset with the LMC-SM kernel.
AUDAUSO
CAGAICD
CHEAUSO
115
0E2./50
H5DU50
JPYSSD
KRWIUSD
MHNUND
120
130
140
Tranndtt
Test data
Prediction
Confidance
Figure 23: MOGP regression plots on currency dataset with the CONV kemel.
31

---

Published as a conference paper at ICLR 2026
Bramble
Camber
Chi
Soton
05401
06/07
96/01
06/05
06/00
D6/08
00/01
06/11
05/04
06/96
06/0a
Train data
Test data
Prediction
95% Confidence
Figure 24: MOGP regression plots on the bramble dataset with the MO-LRN kernel.
Camber
40
-200
-000
Soton
10
-5
-10
-15
06/01
05/022
06/28
06/04
06/05
0t/C6
06/b7
06/02
06/01
06/022
06/t0
06/26
06/07
Train data
Test data
Prediction
Confidence
Figure 25: MOGP regression plots on bramble dataset with MOHSM kernel.
Bramble
Camber
Chi
Soton
06/07
Train data
Test data
Prediction
95% Confidence
Figure 26: MOGP regression plots the bramble dataset with the MOSM kerne..
32

---

Published as a conference paper at ICLR 2026
Bramble
Chi
Soton
06/01
06/02
06/03
06/04
06/05
06/37
06yna
06/01
06/32
06/03
06/04
05/05
06/06
06/17
06/03
Train data
Test data
Prediction
95% Confidence
Figure 27: MOGP regression plots on the bramble dataset with the LMC-NGSM kernel.
Bramble
Camber
00/01
06/02
D0/01
08y04
oxyos
06/05
06/07
06/06
06/04
06/02
06/03
06/04
06/05
06/06
DK/117
06/08
Train data
Test data
Prediction
95% Confidence
Figure 28: MOGP regression plots on the bramble dataset with the LMC-SM kernel.
Bramble
Camber
Chi
Soton
Train
Test data
Prediction
95% Confidence
Figure 29: MOGP regression plots the bramble dataset with the CONV kemel.
33


Reviews:
[{"review_id": "5Xs6qHs7x9", "reviewer": "Reviewer_5TkQ", "summary": "This paper proposes the Multi-Output Low-Rank Nonstationary  kernel for multi-output Gaussian processes. By introducing a generalized spectral–kernel duality and modeling the matrix-valued spectral density via a low-rank Gaussian mixture, MO-LRN achieves linear scalability and superior expressiveness, outperforming existing MOGP kernels on regression and imputation benchmarks.", "questions": "I do not have quetions.", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "The paper’s strengths lie in its clear theoretical and practical contributions.\n\n1 It establishes a new spectral–kernel duality that removes conventional restrictions and enables fully flexible matrix-valued spectral densities for multi-output Gaussian processes.\n\n2 This paper introduces the MO-LRN kernel, a parameter-efficient yet expressive nonstationary kernel that reduces parameter growth from quadratic to linear through a low-rank spectral density with independent Gaussian-mixture factors. I really like this idea.\n\n3 Extensive experiments validate  its effectiveness and scalability in modeling complex multi-output, nonstationary processes.", "weaknesses": "I really enjoy reading this paper but I am not an expert in GP. From a general perspective, I cannot find obvious flaw in this paper. I will read other reviewers' comments along with authors' feedback and ajust my score.", "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": "tuMAtyppJl", "reviewer": "Reviewer_kt5w", "summary": "The manuscript presents a novel formulation of a multi-output non-stationary spectral kernel. A low rank approximation is adopted in order to maintain a manageable number of trainable parameters. The model is benchmarked against synthetic and real world data, and shown to be competitive for point prediction tasks. The approach also demonstrates favourable computational efficiency.", "questions": "Under what conditions should we expect this formalism to fail? It might be helpful to construct a somewhat adversarial synthetic example where we know the model cannot reproduce the cross-covariances, and see how it performs. \n\nAs pointed out in 4.1, spectral kernels tend to be highly sensitive to the initialisation strategy, but it was not clear to me what initialisation is used for MO-LRN? And how should practitioners decide on a suitable configuration such as Q for a given dataset?", "limitations": "", "rating": 4, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The manuscript is well written and clearly presented\n\nThe approach is a novel solution to tackling a challenging problem in the literature.\n\nThe experimental results are clearly presented along with a suitably broad selection of benchmarks and a good set of baselines.", "weaknesses": "My primary concern with this work is the absence of quantitative probabilistic performance metrics. The metrics shown, such as MAE and RMSE, test only the quality of the point predictions of the model. But here we are working with probabilistic models, so without metrics such as NLPD to verify the uncertainty calibration, it is difficult to recommend acceptance. This concern would hold for any probabilistic model, but is particularly acute here where we have many degrees of freedom (due to non-stationarity and multi-output), as it becomes more challenging to regulate the full posterior.\n\nIt is unclear what we are sacrificing when moving to the low rank setting. \n\nScalability remains a significant concern that is only briefly discussed in the Appendix. For multi-output datasets that are long enough to reveal detectable non-stationary features, we quickly hit the upper limits of Figure 4, where the computational cost becomes prohibitive.", "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": "sAmvjw02UZ", "reviewer": "Reviewer_Ls9o", "summary": "The paper introduces a multi-output version of the next-gen SM kernel from Yang (2025).", "questions": "* Table 2: \"$I$ denotes the index in the LMC summation over latent processes, with $I < V$.\"\n   Should this be $I \\le V$?\n * Table 2: what does \"Expressive kernel\" mean?\n * Eq (1) is missing parentheses around the 4 terms in the integral.\n * \"When ω_1 = ω_2, this theorem reduces to Bochner’s theorem.\"\n    ω are not parameters of the theorem, so this statement makes no sense. Do you mean to say that the theorem reduces to Bochner's theorem when u(ω_1,ω_2)=δ(ω_1-ω_2)u(ω), so when u is non-zero only when ω_1 = ω_2?", "limitations": "", "rating": 8, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "* The proposed kernel gives an elegant solution to the problem.\n* The method has good results.", "weaknesses": "* The appendix contains the hyperparameter values, but only for MO-LRN, not for the other kernels.\n  * The appendix also specifies how the kernel was optimized for MO-LRN (500 iterations of Adam), but not for the other kernels. Is the same procedure used for all kernels?\n  * It is not clear how the hyperparameters like learning rate were chosen.\n* The source code is not released (yet)", "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": "5KBtXZHPuE", "reviewer": "Reviewer_ZUMe", "summary": "I applaud the authors for revisiting multi-task GPs through non-stationary kernels. I fully agree with the paper's motivation. The manuscript is well written and a delightful, helpful read. My comments and questions mainly seek clarification.", "questions": "- Are there restrictions or bounds on the m_ij, z_ij, S_ij^(q), and other hyperparameters? \n- What is the reasoning behind the choice of competitor methods?  \n- Why is there no comparison to a baseline coregionalization?", "limitations": "", "rating": 10, "confidence": 4, "soundness": 4, "presentation": 3, "contribution": 4, "strengths": "- The paper proposes a new, very flexible kernel for multi-task GPs, which is an important area of research\n- The manuscript is well-written and comprehensible without too much effort.\n- Sound theory.\n- A reasonable number of test cases.\n- A good choice of competitor methods.", "weaknesses": "- Some extra implementation details would be good. I found it hard to reproduce the kernel from Definition 1. For instance, one could describe typical bounds for the hyperparameter.", "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 (Poster)