Gaussian Processes for Shuffled Regression
Authors: Masahiro Kohjima
Venue: NeurIPS 2025 poster ()
Keywords: shuffled regression
gaussian processes
Primary area: probabilistic_methods
TL;DR: We propose a probabilistic model using Gaussian processes for learning regression functions from data without correspondence.

Abstract: Shuffled regression is the problem of learning regression functions from shuffled data where the correspondence between the input features and target response is unknown. This paper proposes a probabilistic model for shuffled regression called Gaussian Process Shuffled Regression (GPSR). By introducing Gaussian processes as a prior of regression functions in function space via the kernel function, GPSR can express a wide variety of functions in a nonparametric manner while quantifying the uncertainty of the prediction. By adopting the Bayesian evidence maximization framework and a theoretical analysis of the connection between the marginal likelihood/predictive distribution of GPSR and that of standard Gaussian process regression (GPR), we derive an easy-to-implement inference algorithm for GPSR that iteratively applies GPR and updates the input-output correspondence. To reduce computation costs and obtain closed-form solutions for correspondence updates, we also develop a sparse approximate variant of GPSR using its weight space formulation, which can be seen as Bayesian shuffled linear regression with random Fourier features. Experiments on benchmark datasets confirm the effectiveness of our GPSR proposal.

Gaussian Processes for Shuffled Regression
Masahiro Kohjima
NTT, Inc. masahiro. masahiro.kohjimaene komiiommtt.mmmmmccccmo
Abstract
Shufiled regression is the problem of learning regression functions from shuffled data where the correspondence between the input features and target response is unknown. This paper proposes probabilistic model for shuffled regression called Gaussian Process Shuffled Regression (GPSR). By introducing Gaussian processes as prior of regression functions in function space via the kernel function, GPSR can express a wide variety of functions a nonparametrie manner while quantifying the uncertainty of the prediction. By adopting the Bayesian evidence maximization framework and a theoretical analysis of the connection hetween the marginal iikelihood/rredictive distribution of GPSR and that of standard Gaussian process regression (GPR), we derive an GPSR that iteratively applies GPR and updates the input-output correspondence. To reduce computation costs and obtain closed-form solutions for correspondence up- dates, we also develop sparse approximate variant of GPSR using its weight space formulation, which can be seen as Bayesian shuffled linear regression with random Fourier features. Experiments on benchmark datasets confirm the effectiveness of our GPSR proposal.
inference algorithm for
Introduction
The purpose of shuffled regression (SR) is to learn regression functions from shuffled data where the correspondence between the input features and the target response is unknown [1 3 4 This
situation often arises with data collected independently from multiple devices/viewpoints 4 5 6]
or when personal information must be anonymized for privacy reasons, such as in cases involving clinical history 7 Since traditional supervised regression methods, collectively referred to as coupled regression (CR), rely data with input-output correspondences (the pairs of features and responses) and so cannot deal with shuffled data, SR has become the focus of recent studies
9 10 12].
In SR literature, main focus is the use of linear models [ 4 9
(an exception is 12 which uses
neural networks as explained in $2). Accordingly. the use of Gaussian processes (GPs)
[3] has not been well investigated. Considering the known benefits of GPs such as nonparametric flexibility and uncertainty quantification, efforts are needed to realize their benefits in SR.
This study proposes probabilistic model based on GPs for SR called Gaussian Process Shuffled Regression (GPSR). By introducing GPs as a prior of regression functions in function space via the kernel function, GPSR can express wide variety of functions in a nonparametric manner while quantifying prediction uncertainty shown in Table T
By use of the Bayesian evidence maximization framework and a theoretical analysis of the connection between the marginal likelihood/predictive distribution of GPSR and that of standard Gaussian process regression (GPR), we derive an mas-too-mpllemet inference algorithm for GPSR using GPR subroutine; it iteratively applies GPR for the optimization of the kernel's hyperparameters and solves the quadratic assignment problem (QAP) for the updates of the input-output correspondence. Although this algorithm is valid for scall-scale data, it seems prohibitive for handling medium/large-
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

Table I: Comparison of the proposed method (GPSR) with Gaussian process regression (GPR) and the existing methods for shuffled regression (SR). Item (a) indicates which meth- ods can handle shuffled data, while item (b) indicates which methods can express non-linear regression functions. Finally, item (c) shows the methods that can quantify the uncertainty of their prediction.
Table 2: Classification of GPSR, GPR, Bayesian Linear Regres- sion (BLR) and Bayesian Shuf- fled Linear Regression (BSLR), We establish the connection between these four methods.
Problem Type SR
GPR 13
SLR [4]
SDR 12
GPSR
CR
(ours)
function space weight
GPR
GPSR
(a) shuffled reg.
Model
(b) nonlinear
BLR
BSLR
uncertainty
space
scale data due to the difficulty of solving QAP (QAP NP-hard in general) and the computation cost of GPR, which scales cubically with the number of samples. Therefore, we also develop sparse approximate variant of GPSR called Sparse Spectrum GPSR (SS-GPSR) using the weight space formulation of GPSR. This can be seen as Bayesian Shuffled Linear Regression (BSLR) with random Fourier features, and allows us to reduce both memory and computation cost and to use sorting operations to obtain closed-form solutions of correspondence updating. Note that both developing an algorithm for BSLR and clarifying the relations among four methods in Table are also contributions. Experiments on benchmark datasets confirm the effectiveness of our proposals.
2 Related Works
This work is positioned within the literature of shuffled regression and that of Gaussian processes. We describe below prior works in these research lines.
Background of shuffled regression (SR). SR 3 which is also called "regression without "permuted and unlinked regression" [ arises in various fields, including flow cytometry for measuring chemical characteristics cells [4], image/point cloud registration [5] 0 and linkage of health records [7]. Thus various theoretical aspects and problem settings of SR have been studied. For example, Pananjady et al. show that recovering the correspondence between input and output is NP-hard in general (3) the conditions
correspondence" [4, "uncoupled regression" [15 16
under which unique solutions can be achieved have been explored [17]. Some existing works consider
the setting with additional data such pairwise comparison data [10 15 16 18]
We focus on the
setting where no additional data is, only the shuffled data.
Models and algorithms for SR. Earlier works on SR mainly use linear models categorized as shuffled linear regression (SLR). Only in recent years have approaches that utilize more flexible models, such as the neural-network-baseee method called shuffled deep regression (SDR) [12). been considered, and our study using GPs follows this research line. For comparison, see Table [ For parameter estimation, variants of the pppectationmmaximizatiii (EM) [19] are frequently used, 20 In contrast, our inference algorithm is derived from the Bayesian evidence
4 8 and
e.g., 4 maximization framework [2] 22 which has, up now, not been used in the context of SR. Since the evidence penalizes overly complex models by integrating over parameters, this algorithm offers the way of hyperparameter tuning without validation data while still avoiding overfitting.
Gaussian processes (GPs). For supervised regression, Gaussian process regression (GPR) provides nonparametric flexibility and the quantification of uncertainty. However, naive GPR has limited practicality because its computation cost scales to the cube of the number of data points. For handling datasets with more than a few thousand points, various scalable approaches such as using inducing points 23
25 26), grid points [27 28), and spectral points or called random Fourier 29 30 have been proposed. These so-called sparse GPs are still an active part of
24
features (RFF)
GP research e.g., 32
33 Our study adopts the approach using spectral FF (2 30 to
develop SS-GPSR, not only to reduce computation cost, but also because equations for updating the input-output correspondence can be obtained closed-form solutions via a sorting operation.

---

3 Preliminaries
In this section, we provide definitions of shuffled data and its connection coupled data. We also review Gaussian process regression (GPR) for "coupled" data used in supervised regression. All symbols are listed in the Appendix A
Definition of Shuffled Data. Let A' C R and y S R be the imput space and response space, respectively. Shuffled data DsD are defined by pairs of input-set X, {x11 X,2, ... X,L } and response-set set V/ (Ve,, V/22. ). i.e., DSD = (Xx, ( where X E X, Vi y, T. is set size | and M is the number of pairs. Unlike "coupled" data, the correspondence between input and response is unknown in shufßld data; the response corresponding input X is one of the elements in response-set y; - (Wil, 3.r.). but it is not known which one. For simplicity, we define the total number of inputs/responses in DSD as N & ML and denote all inputs and fi=1= ) respectively. Also we equate symbol y; with all responses as X {X,),M, 1= and y {y:
L-dimensional column vector y; (yil- ... M/rt) T and symbol y with N-dimensional column vector y (yT .Y)) T whose {(i-1)I + f)-th element is V/e-
Connection between Shuffled Data and Coupled Data. Shuffled data DSD are regarded as a general representation of coupled data DeD defined as pairs of input x; X and response y; ), DCD = [x,, v:) This is because shuffled data with set size = 1 are equivalent to coupled data. To distinguish it from shuffled data, bold symbols/vectors, etc. related to couple data will henceforth be underlined; we denote all inputs and all responses in DCD as X
1x: ] and
= (y1..... yn) respectively. Note that early works on shuffled regression such as 22 focus on shuffled data where the number of pairs is M = 1. In practical scenarios of data collection, data are often collected multiple times by repeated measurements or by changing the target group of users; even if the input-response correspondence is obscured we can know from which measurement/user- groups the data was collected, so we consider shuffled data where the number of pairs is M > 1.
Gaussian Process Regression (GPR) for Coupled Data [13]. GPs are collection of random variables [f(x)]x x E } for which any finite subset follows a multivariate Gaussian specified by mean function m(x) X IR and kernel/covariance function k(x, x') X' X - TR. GPR for coupled data uses (zero-mean) GPs with kernel ky depending on hyperparameter as priors over function f. which denoted as GP(O, (x, x')). Each response yi assumed to be obtained by adding Gaussian noise to f; = f(x)), i.e., y/ = f + e, where ~ 0,, and B is the precision. Introducing function values vector f = (f,, ... fN)T, this data generative process is summarized as where K is the gram matrix whose (i,j)-th element is Ko ( x,) and IN is the N x N identity matrix. Figure Th shows the graphical model representation of GPR. We get the aaaalytically-taetabbl latent function values f, X(y
X) [f; 0, K) and f(y]f)
f-
marginal Tikelihood by integrating over f p(y(f) ( (f (fff Cg) where C K + 5 II
Thus hyperparameter is estimated by the evidence framework 21 22 of marginal likelihood (evidence) defined as
that maximizes the logarithm
N log(2r). (1)
(GPR. Marginal L.) LGPR (6; DCD) - logply = log C
2 ' y 2
2
Similarly, we get the predictive distribution response y/- for test input x, as follows.
(GPR Prediction) PGPR(y.lx,, DCD) = Nyy=: morn(x.)).. o (xPR(x.) k ' (x.) = k(xx,(xx x.) + - - C k.
where k an N-dimensional vector defined by k- = (I),,x),, (x(, Thus GPR can make predictions with information on confidence levels.
Spectral Representation of Shift Invariant Kernel. One bottleneck of GPs is that, in direct implementation, the memory and computation requirements for inverting covariance C scale as O(N2) and ( respectivel.. Accordingly, Lazaro et al. [29] use the spectral representation or the
'This could be extended to a setup where I, different for each sample but that is omitted for simplicity.
2 Also called empirical Bayes, Type-II maximum likelihood or generalized maximum likelihood 34 35
3

---

(8) Gaussian Process Regression (GPR)
(b) Gaussian Process Shuffled Regression (GPSR)
Responses
1/4
Responses
V/11
V/12
V12
g+
Me1
BKL
yx
Shuffle by Permulation
11
T
Gaussian Field
f.
fz
Gaussian Field
f12
fiz.
Inpuls
X1
X2
X;
Inputs
X11
X12
XIL
X=1
X.L
X
Figure 1: Graphical models of (a) Gaussian process regression (GPR) and (b) proposed Gaussian process shuffled regression (GPSR). Shaded nodes represent observed variables. The thick horizontal bar represents set of fully connected nodes. Dependency on hyperparameters is omitted. GPSR includes GPR as a special case since these two are identical when set size L
so called random Fourier features (RFF) vb :' + F& S IR that approximate the (shift-invariant) kernel function,
k(x, x x')x=(((x))))> (xx), (cxs sin (2 x). () sin xi) (3)
where H' H/2 and 8y. is D-dimensional vector sampled from the kernel's spectral density, e.g., (s; 0p: (4.T (4T² ") for the (ARD) Gaussian kernel defined by k(x, x') = (f
x') A(x x')) This, together with the weight-space-view formulation of GPR, enables us to create sparse-spectrum approximation of GPR through Bayesian linear regression (BLR) on feature space Fe C Rt This reduces the memory requirement to O(N since no longer store the N covariance matrix Cg. As a result, the predictive distribution and the logarithm of marginal distribution for BLR, denoted by PBLR (3/ x.., DCD) and LBLR can be computed with lower memory usage and computation cost. See Appendix B for thorough definitions and detailed explanations. Although any positive definite kernel can be adopted in our method, assume that kernel k: is shift-invariam when deriving the approximate variant of GPSR using RFF.
Gaussian Process Shuffled Regression (GPSR)
This section presents our probabilistic model based on GPs for shuffled regression (GPSR).
4.1
Generative Process and Marginal Likelihood of GPSR
GPSR is constructed by the following generative process of shuffled data DsD {(x,, y:) y:) -
Similar to GPR [33 in the previous section, we put GPs prior on f. ~ GP (x( x )))) and
denote function values at i-th input-set X; = {x:1. ... } as f; = (f(x.)),. fxxxt) So
N (=ML)- dimennional vector f = (f] ... f T is subject to the multivariate Gaussian given by
X) = ON,
(4)
element is ko x.e). For each sample = .... M, using I. x L latent permutation matrix II, response-set y: (ya,,...
where K is the N x N gram matrix whose (('
" " ie., y; follows ply: [f.) = (y; II,f,, 3 -I,). This yields the joint probability distribution of all
is determined by shuffling latent function values and adding Gaussian noise,
responses y = (yT ... . given f as
p(ylf) Ny:: II, I-IN,,, II diag([I, II, M.).
(5)
From Eq. (4) and 5 we get the marginal likelihood for GPSR by integrating over function value
p(y[X)= ply (lylf)p((f(X)) NV Ow, S IIKIT + = II [K I IT
(6)
C
A permutation matrix is square binary matrix that contains exactly one in each row and column, with all
other elements being 0.

---

Figure DP shows the graphical model representation of GPSR. Note that the order of generating function values and shuffling by II can be exchanged: the identical marginal likelihood is derived even if we first shuffle inputs X and then generate the function values of shuffled inputs.
For further analysis of the marginal likelihood, denote the logarithm of 66 as symbol Lersn.
N
(GPSR Marginal L.) CGPSR PsD) log p(y|)) log S
S
(2T), (7)
complexity
data Att.
and define pseudo coupled data (PCD) as follows:
Definition 4.1. (Pseudo Coupled Data) Given shuffled data DsD = (((xx)/=11 (1e) ) and permutation matrices {II, I i== we define the inversely permuted response y ( :.. . M as
(Inversely Permuted Response) ý, " ys. y y).
(8)
We also denote the l-th element of y, as Yit. Pseudo coupled data DeCD is defined as a set of M.L.
N=ML pairs of input X4t with inversely permuted response lit, DPCD = {(xie, I 1,t=1
Marginal likelihood of GPSR is equivalent that of GPR ( with pseudo coupled data. Proposition 4.2. (Equivalence of Marginal Likelihood) Larsn (0; DsD) = Carn (0; DrcD) holds.
Proof. By the multiplicativity of determinant and the properties of permutation matrix (II = II-1 and II II 1), we get
N 1 TIIC, 'II yy log(2
LGPSR(#; (sD) log( II
Cgl
II
(9)
2
1
N log(2 (o:R(0;DPc))) 7 2
log Cel
II
y) (II y
2
y ST
This result allows us to use GPR as subroutine of the inference detailed in the next subsection.
4.2 Inference
we adopt the Bayesian evidence maximization framework [1 for estimating kernel hyperparameter 0 and permutation matrix II. To effectively utilize the Proposition 4.2 adopt the alternating optimization scheme that iteratively updates @ and II.
Similar to GPR 13
Optimization of 0. Proposition 4.2 allows us to employ exactly same approach GPR when estimating hyperparameter 8. We update # by applying an optimization routine such as scaled conjugate gradient (SCG) and L-BFGS using the following partial derivatives of LGPSR w.r.t. the hyperparameters,
1 (CO @LoPSR (0; DSD) @LCPR (0; DrCD) 00
8CB 00
8Ca
-- C
(C, (10) 36
00
Optimization of II. Considering that permutation matrix is a discrete variable, we cast this as known combinatorial optimization problem. By extracting the terms related to II from LGrSR (Eq. (9h) with multiplication by -2 for simplicity, we get
(IT y) T = tr(C, II y(II" y)T) tr(C, 'I' yy yy III UQAP III,
(11)
where we used the trace trick (z.TAz Az ir( (zz T)). We can regard UQAP as the objective function of quadratic assignment problem (QAP), and so we can employ one of the existing QAP solvers (36)
fl In the later experiments, we use simulated annealing, which is easy to implement and closely related to Markov Chain Monte Carlo (MCMC) used existing shuffled regression methods 20). For
more details on the inference procedure, see Appendix C Note that QAP is a NP-hard problem and this procedure seems prohibitive when handling medium/large scale data; we derive an approximation algorithm for GPSR that can avoid solving QAP in $5
Here II has block diagonal structure and minor modification of constraint conditions may be needed.
5

---

(a) M
(b) M
(c) M
Figure 2: Predictive mean and variance of GPSR for various sample sizes, M= to 3, with set size L=3. Marker A indicates the elements in the input-set and response-set in shuffled data. Same color means same index i. GPSR captures the true function (y = sin .)) as the sample size is increased.
4.3 Making Predictions
We derive the predictive distribution of y/ for input, in GPSR by considering the joint distribution of the [N++)-dimensiond vector y' = (y y-) Note that (x., B-) is "coupled" pair not included in DSD as shown in Fig. Db. We define the latent function value at x, as f. = =f(x)) and = (fT,f.)
(f''; 0 N-1),, K') where K' is (N+1) (N+1) gram matrix computed using X' X x.. From Eq. (5). and assuming p(y. [f.) My, lf., 8-1) similar to GPR, y' is determined by following p(y' = (yy;; II'', (N+1)) where II' diag 1).
From GPs prior on p(f')'))
Then, the marginal distribution y' is given by
p(y' ) /
or WN'; ( where S'-II'K'II' +B I (N++). (12)
Note that S' can be cast in the following block form:
II
ON
Ca
k
II
S ON ) - ko(x, x.) x.)+ (13)
IIk
S'=
ON
k
xx)
ON
where k the N-dimensional vector whose ((i-1)L + f)-th element is ke (xae, x.). The predictive distribution of y. given shuffled data is the conditional Gaussian of (12;;
; it is given by
(GPSR Prediction)
PGrsn x., DsD) = -N mGrsH (x, @SPRR (x.)),
(14)
mGrrn((.) Ss IITH (ICC 'II' Jy = k Cp (II' y) = C
@CPRR (x.) ki (z., x.. 8- (Ilk) S IIk (x. x.) 3- k C k.
Note that IIT II- is used in the above derivation. It is obvious, similar to Proposition 4.2 that the predietive distribution of GPSR is equivalent to that of GPR (2) given pseudo-coupled data.
Proposition 4.3. (Equivalence Prediction) PGPSR (y- [x. DSD) ) = PGPR (y- x., DPCD) holds.
Proof. The proof is completed by replacing y in Eq. 2 by response vector of PCD y.
Figure 2 demonstrates the effectiveness of GPSR with the derived inference algorithm and predictive distribution for the simple toy problem. GPSR can output prediction with its uncertainty and capture the true function as the sample size M is increased. However, the prediction and inference procedure presented in this section is prohibitive for handling medium/large scale data since it has two difficulties: (i) need to compute the inversion of the covariance matrix Cş at the cost of ( operations (same as GPR), and (ii) need to solve QAP for estimating permutation matrix II. The next section provides a sparse variant of GPSR that can resolve these two difficulties at once.
Sparse Spectrum GPSR (SS-GPSR): Approximation of GPSR using RFF.
The key idea behind our sparse variant of GPSR, sparse spectrum GPSR (SS-GPSR), is to use weight-space-viers of GPSR and random Fourier features (RFF) 29 30 Here we start with a model
with a finite-dimensional feature map Φ X' + F, whose inner product corresponds to kernel function We then consider using RFF defined in (3) to approximate GPSR with shift-invariant kernel that may not have corresponding diiee-dimensionll feature map.

---

5.1 Weight Space View of GPSR: Bayesian Shuffled Linear Regression (BSLR)
We derive the predictive distribution and inference algorithm for SS-GPSR from the following Bayesian shuffled linear regression (BSLR) models defined as
p(w)= (W; 0, a IH), p(y, w) = Nyy B.
(15)
where +; is an I H design matrix whose b-th is ((xx)) We see that BSLR is weight- space-view of GPSR as follows. Using the equivalent expression of (15) y; = I, w + & where E, ~N((E; 3 Ir), we have the following mean and covariance of y; and yz:
Ew,e [y]] E Byy,lyyy] 1 = II, 2 Ew ww
II, $ Ew w| + E [c,] OL,
(16)
Φ II + E. [E II - P + + I
So the joint distribution of response-set set y1. Yas is given by zero-mean Gaussian with N x N covariance R II((= $ + 8 INI & is N x H design matrix + (Φ .... $ j/
Thus we can confirm that BSLR is a weight-space-view of GPSR since R is equivalent to the covariance S in GPSR (6) if the kernel is defined by inner product of the feature map, k(x, a ((x) 6(x).
5.2 Prediction and Bayesian Evidence Maximization Framework for BSLR
We derive the predictive distribution and inference algorithm in BSLR. From Bayes rule, p(w (DsD) xx(w w, / -)- (w) M ply, (X,, w) holds and thus a posterior distribution of w is given by
( = (W; W. A )) w 3A-1 (II y) 3A- y, A +34
T d. (17)
So we obtain the following predictive distribution PBSLR (v. x., DSD) = f ply, . X. w () dw
and marginal likelihood LESTR DsD)= p(y, w XX Jow (see Appendix x
D for derivation):
(BSLR Prediction)
(B.LR DSD) = N Ny.: MBLLR ), J.SLR (x.)),
(18)
(x,ta( w Tx(x), B 8-1 + o(x,)
T A
(xx,).
(BSLR Marginal L.) LBST. DsD) 11
N
B
log +
log
log |A| E(w),
2
2m
(19)
E(w) = ow T w + Bly 11 = By Ty y + w Aw 23
Σ tr IIIII, wy
(LAP(IL))
Similar to GPSR, the predictive distribution and marginal likelihood of BSLR are equivalent to those of BLR given pseudo-coupled data (proofs are given in Appendix F Proposition 5.1. (Equivalence of Prediction) PESLE (v- Xx, Dsn) PRL. (v.lx, DPcD) holds.
Proposition 5.2. (Equivalence of Marginal Likelihood) LBSLR (w; DsD)= (w; DPCD) holds.
We can also show the predictive distribution BSLR (8) is equivalent to that of GPSR 14 using the connection between GPSR and GPR, GPR and BLR (Appendix B)) BLR and BSLR.
Proposition 5.3. (Equivalence of Prediction) PGPSR( |x., DsD) = PBSLR (y= (y. DaD) holds if the kernel function k is given by k(x, x') = a (--1((x) (xx).
and 5.1 sequentially, we get PGpSR(y.|x. Xx,, DsD) PGPR (y, X,, DPCD) PELR y, x. DPCD) PESla (y. xx., Dsp).
Proof. Applying Proposition
4.3
Proposition 5.3 states the predictive distribution of GPSR can be obtained via BSLR for kernels
that admit dimeessoiiniii feature maps (regardless of whether the kernel is inffttinvaarant).. For shift-invariant kernels that do not admit dimee-nimeniional feature maps, an approximate predictive distribution can constructed using RFF / 3 as feature map Φ. Therefore, we can more efficiently construct the (approximate) predictive distribution using BSLR when the number of samples N exceeds the dimension size H since the expression (T8) involves the inverse of H H matrix A instead of N x N gram matrix Cg.
7

---

5.3 Inference
By leveraging the result of Proposition 5.2 a simple inference algorithm that iteratively updates the posterior of w and II can be derived, similar to GPSR. Thanks to the equivalence of marginal likelihood, we can estimate posterior of w by exactly same way of BLR (update w and A -1 by (7)) Also we can update II by minimizing the trace term UL. in (19) that can be seen as the objective of linear assignment problem (LAP). This is a simple instance of LAP, so the optimal permutation matrix II, can be found by sort operation in lo L) without using tailored algorithms for LAP 37 38, 39 This is because the objective function ULAP (TT) = vII, yy is maximized by setting the permutation matrix such that the E-th largest element of response vector y; and that of the predictive mean vector v; - P,w w match each other for all L. For more details on the inference procedure, see Appendix F
Thus our SS-GPSR using the predictive distribution (18) and the above inference procedure reduces memory and computation costs associated with GPs and avoids solving QAP. Although existing studies also use sort operations 4 [2], the best of our knowledge, this is the first show that the permutation matrix can be estimated by sorting, even in the evidence maximization framework.
6 Experiments
6.1
Setting
Data. We evaluated GPSR and SS-GPSR using four publicly available data sets found in the UCI 5l airfoil data (Airfoil), concrete compressive strength data (Concrete), Boston housing data (Housing), auto-MPG data (MPG). We also made Housing- 1 with input feature dimension of D = 1 by extracting RM (average number of rooms) features for visualization purposes. We prepared 5 data sets by randomly dividing the data and using 60% for training, 20% for validation, and 20% for testing. We made the shuffled data used for training/validationn by randomly dividing the training/validationn data into M sets each with L elements (i.e., M x 0.6xoriginal datasize/ L) and shuffling the indices in each set. Note that validation data were used only for neural network-basee baselines as explained in next paragraph. The test data in the evaluations were coupled data. We used test mean squared error (test MSE) as the performance metric. We ran 5 trials using 5 sets of training, validation and test data.
machine learning repository
Baselines. We compared (SS-)GPSR with kinar-model-basee method (SLR) 4 and the state-of- the-art method that uses neural networks (SDR) [2]. Both methods use (Gaussian-basedd MSE loss and were trained by the (stochastic) sparse EM algorithm. SDR adopted early-stopping with validation data (which were also shuffled data) 42). As oracle baselines, we also examined linear- regression (LR), GPR, and deep-regression (DR) using coupled data with data size of N LM. Both SDR and DR used a one-hidden--ee feedforward network with ReLU activation function.
Hyperparameters. GPSR (and Oracle GPR) used Gaussian kernels as in $3, and SS-GPSR used its RFF approximation with dimension size H 100 in common for all datasets. The other hyperparameters such as precision a and | in (SS-)GPSR are estimated when applying GPR or BLR in inference procedure (see Alg. [ and 2 in Appendix). So GPSR and SS-GPSR do not use validation data unlike SDR. Moreover, due to their slow convergence, both GPSR and SDR used warm initialization similar to 12). i..., the initial value of permutation matrix is set using the estimated result of SLR. For SS-GPSR, ten different runs were performed with various initial permutation matrices, and the solution that maximized the marginal likelihood selected. For more details of implementations and hyperparameters such as optimization setting for SDR, see Appendix G
6.2 Results
Table 3 shows the results of the experiments. We observe that SDR, GPSR and SS-GPSR which can express nonlinear functions, outperform SLR in almost all settings. When comparing SDR and
'This problem also can be seen as an optimal transport problem between two 1-D point clouds consisting of the same number of elements having unit mass 44) https: /ttpp//////rc.e..eee icc.. uci wdw/em/mmmmae hit//indexxai php
8

---

Table 3: Results benchmark datasets with various sizes L. Average and standard deviation of test MSE are shown. Bold/underline means the Ist/2nd best (lowest) MSE among the methods for shuffled regression. Note that the performance of oracle methods including linear regression (LR), deep regression (DR) and GPR trained using coupled data are also displayed in the Dataset column.
Dataset
Existing Methods
Proposed Methods
Source
SLR
SDR
GPSR
SS-GPSR
Airfoil
23.66 + 1.39
19.57 + 1.73
9.05 + 0.99
12.48
1.36
LR: 22.33 1.72
22.82 + 1.62
20.15
1.91
14.26
1.46
14.11
2.07
DR: 17.79 2.49
8
21.05 + 1.66
21.73 + 2.02
19.14 + 2.14
17.21
2.64
GPR: 6.42 + 0.47
16
25.18 + 2.01
23.83 + 1.99
22,44 = 2.28
24.22 + 1.78
Concrete LR: 119.06 + 7.00
2
139.79 + 4.13
101.62 + 26.15
45.05 = 7.64
52.20 + 2.64
4
119.66 + 4.96
101.04 + 27.47
67.22 + 8.79
58.57 4.31
DR: 82.66 + 37.06
128.62 + 4.62
121.67 20.89
97.79 5.54
68.19 + 3.84
GPR: 39.21 + 11.08
16
161.95 + 30.28
165.32 + 33.61
146.74-334.5
137.26+20.75
Housing-1D LR: 41.09 7.54
51.40 6.78
37.08114.47
33.60 8.29
34.45 8.37
40.95 7.69
35.19 9.59
33.64 8.26
34.85 8.63
DR: 41.25 + 23.52
41.53 8.42
35.72 + 8.74
35.22 8.45
36.24 + 8.68
GPR: 33.70 + 8.01
16
43.07 8.88
37.78 + 8.61
37.18 9.19
37.93 8.51
Housing LR: 23.86 5.69
2
33.68 6.10
12.31 1 4.71
11.19 4.51
11.70 1 4.95
24.57 5.11
11.73 1 3.53
11.44 4.91
12.70 5.17
DR: 11.03 3.38
2.11
18.25 7.42
18.84 8.14
17.11 8.04
GPR: 9.78 = 4.22
16
45.03 + 16.43
39.11
41.25=17.86
46.88=23.72
MPG LR: 11.92 1.10
12.31 + 1.30
8.96 + 1.55
8.21 + 1.20
7.60 1.29
13.24 + 1.59
10.54 + 1.75
8.71 + 1.83
8.47 = 1.42
DR: 7.69 + 1.00
8
15.15
+ 3.25
11.55
+ 3.75
10.94 + 3.21
9.05 = 1.80
GPR: 7.75 0.93
16
21.50
+ 6.49
15.38 + 4.23
6.72 18.43
10.65 3.03
Predictive Hean
400
Oradle GPR
GPSR (L-2)
(MEDV)
GPSRicurs
Test Data
(MEDV)
40
Test NL.1
350
GPSR
GPSR (L-4)
20
20
300
250
Housing:
Hluusin
(RM)
(RM)
Dataset
(a) Learned regression functions
(b) Predictive uncertainty
(c) Negative log- le-liiklihhod
Figure 3: (a) Learned regression function (predictive mean) and (b) predictive uncertainty of GSPR for Housing- (L = 8). The shaded regions represent 11 +2, and +3 standard deviations, respectively. (c) Comparison of negative log-like ligeliikl on test data (Test NLL).
(SS-)GPSR, either of our proposals (GPSR or SS-GPSR) demonstrates the best performance across all datasets, except for Housing (L = 16), confirming our methods' effecttiveness. Note that GPSR and SS-GPSR perform similarly well across comparable number of datasets, so we cannot judge which is superior in terms of prediction performance. From Fig. 3a we can see that GPSR well captures the non-linear structure thanks superior performance comparable to that of oracle methods for Housing-1 D. Figure 3b illustrates that GPSR effectively captures the inherent uncertainty the data, while Fig. 3e] the negative log-likelihode on the test data, demonstrating the method's effectiveness in uncertainty quantification.
nonparametric flexibility, and this contributes to its
We also investigated the detailed behavior of SS-GPSR as illustrated in Fig. In Fig. 4al we observe that SS-GPSR with the random initialization scheme can identify better solution than warm initialization (in terms of marginal likelihood) and converges within about iterations. Given that SDR requires longer time to converge without warm initialization 12 these findings highlight
We confirmed that GPSR has similar convergence behavior to SDR.

---

Log- Log-Marginal Likellhood
17.5
-540
Test MSE
15.0
-560
12.5
-580
Randam
-600
Init by SUR
5
10
15
20
10
10
10
the number of iterations
H: Dimension Size of RFF
(a) Convergence behavior
(b) Effect of RFF
Figure 4: (a) Convergence behavior of SS-GPSR with different initialization schemes (10 different runs for random init.) for MPG (L = 8), and (b) effect of the dimension size of RFF on SS-GPSR for MPG (L 2~8).
the effectiveness of SS-GPSR. Figure furtther further explores the impact of the dimension size of RFF, H. It is shown that although MSE performance is degraded when H is smaller than 100, it stabilizes when Hl exceeds 100. Thus. SS-GPSR utilizing RFF works well without requiring sensitive tuning of the dimension size. These results support the usefulness of SS-GPSR.
7 Conclusion
In this study, we proposed GPSR and SS-GPSR for learning regression functions from shuffled data. By theoretically analyzing the connection between GPSR, GPR and their weight-space formula- tions, we derived the easy-to-implement inference algorithms for (SS-)GPSR, and confirmed their effectiveness by experiments on public benchmark datasets.
Limitation and Future Work. The proposed methods are the first to apply GPs for shuffled regression, and promising direction for future work is extend their applicability to higher- dimensional and more complex settings. This may involve incorporating advanced techniques such as deep kernel learning [33, which combines neural networks to learn task-specific feature representations, as well sparse GPs described in $2 which offer scalable solutions for large-scale problems. Additionally, exploring fully Bayesian approaches that introduce prior distributions over the kernel parameters and permutation matrix offers another promising avenue for future research, as similar techniques have been successfully applied related domains 20
Societal Impact. Shuffled regression methods, including ours, useful for analyzing e.g., in- dependently collected data as mentioned in $ However, it's important to note that coupled (test) data with input-output correspondence is necessary evaluate the regression performance (That's why the experiments in this paper are based on benchmark data). When applying the method to real-world applications without coupled (test) data, it's crucial to carefully assess the performance of the estimated model when making predictions or decisions. Furthermore, the interpretation of input-output correspondences via latent variable representations requires careful consideration, as such usage falls outside our primary scope and may pose risks in scenarios involving privacy concern..
Acknowledgments and Disclosure of Funding
We would like to thank the anonymous reviewers for their insightful comments and constructive feedback, which have helped improve the quality and clarity of this paper. This research was supported by the Human Informatics Laboratories at NTT, Inc. as part of the authors' regular employment responsibilities. The authors declare no competing interests related this work.
References
[1] Alexandra Carpentier and Teresa Schlüter. Learning relationships between data obtained independently. In Artificial Intelligence and Statistics, pages 658-666, 2016.
10

---

[2] Daniel J Hsu, Kevin Shi, and Xiaorui Sun. Linear regression without correspondence. Advances in Neural Information Processing Systems, 30, 2017.
[3] Ashwin Pananjady, Martin Wainwright, and Thomas A Courtade. Linear regression with shuffled data: Statistical and computational limits of permutation recovery. IEEE Transactions on Information Theory, 64(5):328--3300, 2017.
[4] Abubakar Abid and James Zou. A stochastic expectation-maxinine approach to shuffled linear regression. In Annual Allerton Conference on Communication, Control, and Computing, pages 470-477, 2018.
[5] Philip David, Daniel Dementhon, Ramani Duraiswami, and Hanan Samet. Softposit: Simul- taneous pose and correspondence determination. International Journal of Computer Vision, 59:259-284, 2004.
[6] Feiran Li, Kent Fujiwara, Fumio Okura, and Yasuyuki Matsushita. Generalized shuffled linear regression. In IEEE/CVF International Conference on Computer Vision, pages 6474 6483, 2021.
[7] Xu Shi, Xiaoou Li, and Tianxi Cai. Spherical regression under mismatch corruption with application to automated knowledge translation. Journal of the American Statistical Association, I (16(536))1153-1999 2021.
[8] Martin Slawski, Mostafa Rahmani, and Ping Li. A sparse representation--asee approach to linear regression with partially shuffled labels. In Uncertainty in Artificial Intelligence, pages 38 48. 2020.
[9] Guanhua Fang and Ping Li. Regression with label permutation in generalized linear model. In International Conference Machine Learning, pages 9716-9760, 2023.
[10] Ikko Yamane, Yann Chevaleyre, Takashi Ishida, and Florian Yger. Mediated uncoupled learning and validation with bregman divergences: Loss family with maximal generality. In Artificial Intelligence and Statistics, pages 4768-4801, 2023.
[11] Martin Slawski and Bodhisattva Sen. Permuted and unlinked monotone regression in r' d: an approach based on mixture modeling and optimal transport. Journal of Machine Learning Research, 25(183)11557, 2024.
[12] Masahiro Kohjima. Shuffled deep regression. In AAAI Conference Artificial Intelligence. pages 19238-13245 2024.
[13] Christopher KI Williams and Carl Edward Rasmussen. Gaussian processes for machine learning. MIT press Cambridge, MA, 2006.
[14] Manolis C Tsakiris, Liangzu Peng. Aldo Conca, Laurent Kneip, Yuanming Shi, and Hayoung Choi. An algebraic-geomeeric approach for linear regression without correspondences. IEEE Transactions on Information Theory, (6(8))5130-511 2020.
[15] Liyun Xu, Junya Honda, Gang Niu, and Masashi Sugiyama. Uncoupled regression from pairwise comparison data, Advances in Neural Information Processing Systems, pages 3992- 4002, 2019.
[16] Ikko Yamane, Junya Honda, Florian Yger, and Masashi Sugiyama. Mediated uncoupled learning: Learning functions without direct outut- correspondences. In International Conference on Machine Learning, pages 11637-11647, 2021.
[17] Jayakrishnan Unnikrishnan, Saeid Haghighatshoar, and Martin Vetterli. Unlabeled sensing with random linear measurements. IEFE Transactions on Information Theory, 64(5)3337-3253,, 2018.
[18] Masahiro Kohjima, Yuta Nambu, Yuki Kurauchi, and Ryuji Yamamoto. General algorithm for learning from grouped uncoupled data and pairwise comparison data. International Conference Neural Information Processing, pages 153-164, 2022.

---

[19] Arthur P Dempster, Nan M Laird, and Donald B Rubin. Maximum likelihood from incomplete data via the em algorithm. Journal of the Royal Statistical Society: Series B (Methodological), 39(1):1-22, 1977.
[20] Zhenbang Wang, Emanuel Ben- Den-David, and Martin Slawski. Regularization for shuffled data problems via exponential family priors on the permutation group. Artificial Intelligence and Statistics, pages 2939-2959, 2023.
[21] David JC MacKay. Hyperparameters: optimize, or integrate out? In Maximum Entropy and Bayesian Methods: Santa Barbara, California, USA, 1993, pages 43-59. Springer, 1996.
[22] Christopher M Bishop. Pattern recognition and machine learning. Springer, 2006.
[23] Edward Snelson and Zoubin Ghahramani. Sparse gaussian processes using pseudo-inputs. Advances in Neural Information Processing Systems, 18, 2005.
[24] Joaquin Sattonerr-Candelle and Carl Edward Rasmussen. A uniffying view of sparse approx- imate gaussian process regression. Journal of Machine Learning Research, 6:1939-1959, 2005.
[25] Michalis Titsias. Variational learning of inducing variables in sparse gaussian processes. In Artificial Intelligence and Statistics, pages 567-574, 2009.
[26] James Hensman, Nicolo Fusi, and Neil D Lawrence. Gaussian processes for big data. In Uncertainty in Artificial Intelligence, pages 282-290,10 2013.
[27] Andrew Gordon Wilson. Covariance kernels for fast automatic paltern discovery and ex- trapolation with Gaussian processes. PhD thesis, University of Cambridge Cambridge, UK. 2014.
[28] Andrew Wilson and Hannes Nickisch. Kernll interpolation for scalable structured gaussian processes (kiss-gp). In International Conference on Machine Learning, pages 1775-1784, 2015.
[29] Miguel Siizar--Gredilll Joaquin Quinonero-Candele Carl Edward Rasmussen, and Aníbal R Figueiras- Vigueiras-Vidal. Vidal. Sparse spectrum gaussian process regression. Journal of Machine Learning Research, 11: 1018865-1888 2010.
[30] Ali Rahimi and Benjamin Recht. Random features for large-scale kernel machines. Advances in Neural Information Processing Systems, 20, 2007.
[31] Rishit Sheth, Yuyang Wang, and Roni Khardon. Sparse variational inference for generalized gp models, In International Conference on Machine Learning, pages 1302-1311, 2015.
[32] James Hensman, Nicolas Durrande, and Arn Solin. Variational fourier features for gaussian processes. Journal of Machine Learning Research, 18(151):1-52, 2018
[33] Vidhi Lalchand, Wessel Bruinsma, David Burt, and Carl Edward Rasmussen. Sparse gaussian process hyperparameters: optimize or integrate? Advances in Neural Information Processing Systems, pages 16612-16623, 2022.
[34] Grace Wahba. A comparison of gev and gml for choosing the smoothing parameter in the generalized spline smoothing problem. The annals of statistics, pages 1378-1402. 1985.
[35] José M Bernardo and Adrian FM Smith. Bayesian theory. John Wiley & Sons, 2009.
[36] Eliane Maria Loiola, Nair Maria Maia De Abreu, Paulo Oswaldo Boaventura-Netto, Peter Hahn, and Tania Querido. A survey for the quadratic assignment problem. European Journal of Operational Research, 17602):657-690, 2007.
[37] Harold Kuhn. The hungarian method for the assignment problem. Naval Research Logistics Quarterly, 2(1-2):83-97, 1955.
[38] Dimitri P Bertsekas. Auction algorithms for network flow problems: tutorial introduction. Computational Optimization and Applications, 1:7-66, 1992.
12

---

[39] Marco Cuturi, Olivier Teboul, and Jean- -han-Philippe Vert. Differentiable ranking and sorting using optimal transport. Advances in Neural Information Processing Systems. 32, 2019.
[40] Julien Rabin, Gabriel Peyré, Julie Delon, and Marc Bernot. Wasserstein barycenter and its application to texture mixing. In International Conference on Scale Space and Variational Methods Computer Vision, pages 435 446, 2011.
[41] Gabriel Peyré, Marco Cuturi, et al. Computational optimal transport: With applications to data science. Foundations and Trends@ in Machine Learning. 11(5-6)3555-077, 2019.
[42] Ian Goodfellow, Yoshua Bengio, and Aaron Courville. Deep Learning. MIT Press, 2016.
http:://ww..ceepleeannnnnooooo org
[43] Andrew Gordon Wilson, Zhiting Hu, Ruslan Salakhutdinov, and Eric P Xing. Deep kernel learning. In Artificial Intelligence and Statistics, pages 370-378, 2016.
[44] David JC MacKay. Bayesian interpolation. Neural Computation, 4(3):115-447 1992.
[45]
Diederik Kingma and Jimmy Ba. Adam: A method for stochastic optimization, 2014.
[46]
Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Kopf, Edward Yang, Zachary Vetito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala. Pytorch: An imperative style, high-performance deep learning library. Advances in Neural Information Processing Systems, pages 8024 8035, 2019.
13

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification:
Our contributions are clearly explained in the abstract and introduction.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made in the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or
NA answer to this question will not perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results
be expected to generalize other settings.
is fine to include aspirational goals as motivation as long is clear that these goals are not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification:
See the Conclusion section.
Guidelines:
The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those not discussed in the paper.
The authors are encouraged to create separate "Limitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification. asymptotic approximations only holding locally). The authors should reflect how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution low images are taken in low lighting. Or speech-to-text system might be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach address problems of privacy and faimess.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't 't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory assumptions and proofs
Question:
For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
Answer: [Yes]
14

---

Justification: We provide the assumptions and proofs in main text. Some propositions used in the proofs are explained in Appendix.
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide short proof sketch provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental result reproducibility
Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that it affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Question:
Answer: [Yes]
Justification: Detailed explanation of the experimental setting are provided $6 and Appendix G.
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is novel architecture, describing the architecture fully might suffice, or if the contribution a specific model and empirical evaluation, it may be necessary either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., in the case of large language model), releasing of a model checkpoint, or other means that appropriate to the research performed.
While NeurlPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) the contribution is primarily a new algorithm, the paper should make clear how reproduce that algorithm. (b) If the contribution is primarily new model architecture, the paper should describe the architecture clearly and fully.
(c) the contribution is a new model (e.g., a large language model), then there should either be way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome deseribe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should possible for other researchers to have some path to reproducing or verifying the results.
5. Open access to data and code
15

---

Question: Does the paper provide open access to the data and code, with sufficient instruc- tions faithfully reproduce the main experimental results, as described in supplemental material?
Answer: [No]
Justification: We do not provide the code but provide implementation details of proposed methods and baselines in $ and Appendix G. Also, the experiment is conducted open- access data (provided in UCI machine learning repository).
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurIPS code and data submission guidelines (https://nips.co cc/
for more details.
While we encourage the release of code and data, we understand that this might not be possible, "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for a new open-source benchmark).
The instructions should contain exact command and environment needed to run to reproduce the results. See the NeurlPS code and data submission guidelines (https
y /nips. /nips.cc//ppp http:ppbbicgggieeeeeeeeeeee
for more details.
The authors should provide instructions on data access and preparation, including how access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended to the paper) recommended, but including URLs to data and code permitted.
6. Experimental setting/details
Question:
Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary understand the results?
Answer:
[Yes]
Justification: Detailed explanation of the data and hyperaparameters are provided in $6 and Appendix G.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper to level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, as supplemental material.
7. Experiment statistical significance
Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Question:
Answer:
[Yes]
Justification:
We provide both mean and standard deviations in Table 3.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by emoor bars, confi- dence intervals, or statistical significance tests, least for the experiments that support the main claims of the paper.
16

---

The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g.. Normally distributed errors).
It should be clear whether the error bar is the standard deviation or the standard error of the mean.
is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that out of range (e.g. negative error rates).
If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text.
8. Experiments compute resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time execution) needed to reproduce the experiments?
Question:
Answer: [Yes]
Justification: See Appendix G.
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental runs as well estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary failed experiments that didn't make into the paper).
9. Code of ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurlPS Code of Ethics
Answer:
[Yes]
Justification:
We have reviewed NeurlPS Code of Ethics
Guidelines:
The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.
If the authors answer No. they should explain the special circumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g., there a special consid- eration due to laws regulations in their jurisdiction).
10. Broader impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [Yes]
Justification:
See the Conclusion section.
Guidelines:
The answer NA means that there no societal impact of the work performed.
If the authors answer NA or No, they should explain why their work has no societal impact why the paper does not address societal impact.
17

---

Examples of negative societal impacts include potential malicious unintended uses (e.g., disinformation, generating fake profiles, surveillance), faireess considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional unintentional) misuse of the technology.
If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms monitor how system learns from feedback over time, improving the efficiency and accessibility of ML).
11. Safeguards
Does the paper deseribe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification: Our paper poses no such risks.
Guidelines:
The answer NA means that the paper poses no such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Question:
Answer: [Yes]
Justification: We provide URLs of the used code and dataset.
Guidelines:
The answer NA means that the paper does not use existing assets.
The authors should cite the original paper that produced the code package dataset.
The authors should state which version of the asset used and, possible, include URL.
The name of the license (e.g., CC-BY 4.0) should be included for each asset.
For scraped data from a particular source (e.g., website). the copyright and terms of service of that source should be provided.
If assets are released, the license, copyright information, and terms of use in package should be provided. For popular datasets, paperswithcone com/ datasets has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
18

---

For existing datasets that are re-packaged, both the original license and the license of the derived asset (if has changed) should be provided.
If this information is not available online, the authors are encouraged to reach out to the asset's s creators.
13. New assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [NA]
Justification: Our paper does not release new assets.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details of the dataset/code/moddl as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset used.
At submission time, remember to anonymize your assets (if applicable). You can either create anonymized URL or include anonymized zip file.
14. Crowdsourcing and research with human subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well as details about compensation (if any)?
Answer: [NA]
Justification: Our paper does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Including this information in the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included in the main paper.
According the NeurIPS Code of Bthies, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional review board (IRB) approvals or equivalent for research with human subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or equivalent approval/review based the requirements of your country or institution) were obtained?
Answer: [NA]
Justification: Our paper does not involve crowdsoureing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Depending on the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurlPS Code of Ethics and the guidelines for their institution.
For initial submissions, not include any information that would break anonymity (if applicable), such the institution conducting the review.
19

---

16. Declaration of LLM usage
Does the paper describe the usage of LLMs if it is an important, original, or non-standard component of the core methods in this research? Note that if the LLM is used only for writing, editing, or formatting purposes and does not impact the core methodology. scientific rigorousness, originality of the research, declaration is not required.
Question:
Answer: [NA]
Justification: The core method development in our research does not involve LLMs as any important, original, or non-standard components.
Guidelines:
The answer NA means that the core method development in this research does not involve LLMs as any important, original, or non-standard components.
Please refer to our LLM policy battps://neure wwwparippssc.coooeeeeeeeee cc/Conferencess///// for what should or should be described.
20

---

A
Symbol Definitions
Table shows the symbols used in our paper. Bold symbols//ectors related couple data underlined.
Table 4: Notation summary for the paper
Symbol
Description
X
input space C _R where D is the dimension size of input space
Y
response space y 'C R
Fv
feature space constructed by given feature map V on X, F, S R
kp
kernel function with hyperparameter & used for Gaussian process, Kg A' x X + R random Fourier features approximating (shift-invariant) kernel, 1 : X Fv feature map used for linear regression model, ó X F+
w
H-dimensional vector parameter (regression coefficients) of linear regression model hyperparameter (precision) in prior distribution on w
precision parameter of observation noise (commonly used regardless of the model)
the size of input/response- set in shuffled data (which reduces to coupled data when L=1) the number of pairs in shuffled data the total number of inputs/responses N ML
M
N
DSD
shuffled data (pairs of input-set X, and response-set ys), DSD y:))
X,
i-th input-set X; = [xa," ... X.L] where X5 EX i-th response-set y: (U-1,". Mii) where yie E L-dimensional vector representation of y: y: (111,... (ViL) all inputs in shuffled data X = { M
y
yi
X
N N-dimensional vector representation of all responses y = (y YYM)
L-dimensional vector rep. of function values of i-th input set, (f.1, fiL)
f
N N-dimensional vector representation of all function values, f (f,
f
II;
L X L permutation matrix used for shuffling function values f, (or w))
II
N X N permutation matrix defined by II = . I....... II
K
N N gram matrix whose ((i-1)
(1'1) L+++")-th element is ko (o Xye)
Co
N N covariance matrix C K + 8 Iy (IN is the N x N identity matrix)
N X N covariance matrix defined by IIKII + IN = IIC
4;
L X H design matrix whose ε-th row (6 ((x:))
+
N H design matrix whose ((-1))++((-th row is d(xi))
DCD
coupled data (pairs of input x; E A' and response yi E J). DCD = {(x, Wil all inputs coupled data X = [x:
X
Y
N-dimensional vector representation of all responses, y (y,,.... ... UN) N-dimensional vector representation all latent function values, f (f,,... [NN N x N gram matrix whose (i, j)tth element is (o(x,, xx)
K
C.
N N covariance matrix C K + I- (IN is the x N identity matrix)
N H design matrix whose i-th row is ((x;)
B Details of Bayesian Linear Regression for Coupled Data
Bayesian linear regression (BLR) with parameter w R and feature map Φ: Fs CRRR is derived by considering the following prior distribution and model,
p(w)= (ww 0, -I p(u: w) = Ny w 8-1) where a and 3 are hyperparameters. From Bayes theorem, p(w (DeD) x ply, w X = p(w) II p(M[x, w) holds and thus a posterior distribution of w given coupled data Dcp is
(20)
written
(BLR Posterior) p(w[PeD) = N(w; W,W W. A W 3A $
A = aI 4 (21)
21

---

where + is the N x H design matrix whose i-th row is o(x;)" T. We get the predictive distribution of y test point x. by integrating over parameter w:
(BLR Prediction)
PBLR(U+|E Xx, DcD) = / p(yo|xe, (eppw/DDDDdd = NN MBLR BLR (22)
(x. = w ((x)) (BLR(x))= = B- + o(x))
A (xx)).
Thus the memory required for Bayesian linear regression is H( to store design matrix 4. Also, the inverse of H x H matrix A can be computed in O(H)) by direct computation.
Remark. Note that w in 21) is both the mean and the mode of the posterior distribution p( ((w/ (DCD), so W is equivalent to the penalized maximum likelihood/maximen too-Bayyssian/requentttt schemes, however, the predictive distribution is constructed by plugging the estimator into the model (without integration over a posterior distribution). i.e., p(y- xx, w) (yy w ((xx),3-1) )) and so the predictive variance is different from (2)
a posteriori (MAP) estimator. In
It is known that the predictive distribution of BLR (22) and that of GPR are equivalent when the kernel function & defined by the dot product of feature map Ø.
Proposition B.1. (Equivalence of Prediction, c.g.,
22 PGPR (ya [x,, DCD) PBLR (v. x,, DCD)
holds if the kernel function k is given by ((x, x')
[xx) 6(x).
We define the logarithm of marginal likelihood for BLR by symbol LBLK as follows:
(BLR Marginal L.) GBLR(W; DCD) A log ply X = ply, w W/ww.
(23)
Note that LBLR is an analytically tractable quantity computed
H log a + log 2x
N 3
1
(W DCD)
log A
W 2
w
y
Dw 1
(24)
2
C Inference Procedure of GPSR
We detail the inference procedure of GPSR explained in $4.2 As stated, simulated annealing (SA) is used for solving QAP since it is easy to implement and closely related to Markov Chain Monte Carlo (MCMC) used in existing shuffled regression methods 4 20). Although not explicitly men- tioned in the main text, can estimate both the precision parameter and the kernel parameter 0 in the GPR subroutine. The equivalence of marginal likelihoods between GPSR and GPR (shown in Prop. 4.2) enables the use of existing GPR subroutines, making the algorithm easier implement.
Algorithm shows the pseudo code that implements the inference procedure of GPSR. At initializa- tion (line 1-3), permutation matrix II is initialized. Note that, reported in [12), convergence can be accelerated by warmm-iitiaalizatii
i.e., the initial setting of II is set by using the result of e.g., SLR [4). Then, the response vector of pseudo-coupled data (Definition 4.1) DPCD. ý II y, is set. At each (algorithmic) time step f, GPR-step (line 5-6) and Correspondence Update (CU)-Step (line 7-14) are iteratively conducted. As the GtR-step, the parameters of 0 optimized by applying GPR using pseudo-coupled data DrCD represented by X and y. Next, as the CU-step, permutation matrix II is updated by randomly swapping its rows following SA. Let denote the permutation matrices before and after random swapping as (old) and II (new), respectively. SA accepts this swapping following probability q:
exp(-//T)) A ((I(ww() - ((I(ol)))
(25)
where T' is the temperature parameter. From the definition of q, the random swapping is accepted with probability when the solution is improved. Otherwise, SA accepts the swapping with probability exp(-A/T). i... the solution is updated with certain probability even when the solution becomes worse, which helps to escape from local optima solutions.
Note that the above algorithm can be extended or modified in various ways, such as using the sparse for the GPR-step and using optimization methods
GP methods deseribed $2 23 24
26 28 33
other than SA shown in 336 for the CU-Step.
22

---

Algorithm Inference Procedure of GPSR
Input: kernel function (with hyperparameter 0, 3), shuffled data Dsn, optimization parameter (e.g..
# of iterations/steps T mazz, Smex, initial temperature Timit, and cooling rate ?)
Output: hyperparamettr 0, 33 response vector of pseudo coupled data y, permutation matrix IT
/* Initialization */
2: Randomly initialize II (Optional: II could be initialized by solving SLR).
3: Set the response vector of pseudo coupled data (Definition 4.1) as y II y.
4: for t = 1 to T maz do
5: 1* GPR-Step (Optimization of hyperparameter using partial derivative Eq. (TO)
6: Update 9 and / by SCG/L-BFGS (This could be done by e.g., CPRRfit(X. y) scikit-leam)
7: 1* Correspondence Update-Step (Optimization of II) by simulated annealing*/
8: Set temperature T <- Timit
9:
for s S Smar do
Randomly select {1.... ,N1 and l, {1,.... L} (e ('). Swap {(i-1)L + f)-th row and {(i-1)L+ + f') row of II with probability (Eq. 25 12: Update temperature T r''.
10:
11:
13:
end for
14: Update the response vector of pseudo coupled data, ý Iy
15:
end for
D Details of Bayesian Shuffled Linear Regression (BSLR)
Here we derive BSLR's ' predictive distribution [8 and marginal likelihood
19 We start with the
parameter posterior distribution as follows:
Derivation of posterior distribution (sD) distribution p(w.y] X) is written as
From the BSLR model
15
the joint
M ply: Xi, w)
p( w. x) p(w)
4
3
3
2 (y IIlww) T(y:))y----r)) (y - -l(bw)
exp
W
W
exp
2T
2
2T
a
4
4
exp
a
W
B
w
T
Φ T
II कw + 3w Φ II
B
2 yr'y)
2m
2m
2
2
cancel out
Q
2
expo
w
(II y))
w
al + 36
w
2 y
B
Tyy
2m
2m
A
#
exp
1
Aw 2w 3
Ts (II")) y) + By })}
(26)
2
2
By the trick of adding 0 = w A
w w ---wwand I A A Eq. 266 is expanded into
H
N
p
pW,N
exp
w Aw 2w AA- 8$ III y) +3y T f)}
2T
2T
w
2
8
N
exp
w Aw - 2w Aw + By y w
Aw w Aw)
2T
2T
H
w) A(w { 08 Ty T Www)]]
exp
(27)
2m
2T
Thus we get the posterior distribution w) A(w - w)).
from p(w DsD) & p(w.y|X) -
23

---

Derivation of predictive distribution (US [X., DSD) and conditional Gaussian distributions, we get
18
From the property of marginal
PBSLR(V.(x., DgD) / ww
Tw, "J (w; W, A d)w
(N.x:;(xx.) "w, + x(x.) ) A A- (xx)).
(28)
Derivation of logarithm of marginal distribution CBSLR (w: DSD) (= log p(y, d(w)dw)
(19)
By the trick of adding 3A - 'I' and A aIn + _ T. we can expand the final term of the argument of the exponential function in (7 follows.
A 1-1ww w 'w and IN II II with the definitions
(By [By] y Aw))
By 2w Aw W Aw)
(ByT
[8yT: y 2w AA
B II y w (aIN 36 x)w)
cancel out
[0y] y 2w 8Φ II
y
Sw
II II tw aw ww]
IN
1 B(y Hlw) T
3 2
=
(y
1lww) aw
lly
11
w.
(29)
Substituting (2) into 27
we get the following expression of joint distribution:
p(w y(w,y|X)
B
exp
! (w w) A( w) 2 fly 1l 1 - w
3
2T
2T
From the property used for normalizing constant of Gaussian distributions, w) A(w w))dw = (2m) # ||| holds and thus we get
-expl-
P(y|X)= = / ply, w/W)dw
-
a
#
exp
W W) A(w - ww) 2 -1yy 1l
12- a T [wT]]ct
zw
2m
2T
a
3
(2t))
A
1/2 exp
9
2 Ily
11T 12
a
2 w w
(30)
2x/
2x/
We get LBSLR (w; DSD) (19) by taking the logarithm of (30). The expression of E(w) using trace function is obtained by the trace trick (z z = tr(zz) )).
E(w) aw w + B
I
aw w + 3 i=1
(y: II,
(y; W.)
aw w + 3 Σ
(y, y: 2 (I, + (II,
=a w 28
tr( 3 y; w
Il p.w) TN
23 Σ M Ir(II By
By y + (aI + 30 T ④)) w.
E
Connection Between BLR and BSLR
Here we show the connection between BLR and BSLR. Similar to the relationship between GPSR and GPR discussed in $ we can prove the equivalence their predictive distributions and marginal likelihoods using the definition of PCD. We begin by presenting proof about predictive distributions.
24

---

Proof of Proposition 5.7 We need to show the equivalence of the posterior mean of BSLR w and that BLR w and posterior variances BA
17
when pseudo coupled data PCD are given. The design matrices ( +
A, A) are equivalent from the definition. So we get w = 3A-1 II
Marginal likelihood of BSLR 19 is also equivalent to that of BLR (23 with pseudo coupled data.
we show the equivalence of the design matrices, posterior means and variances. Equivalence of the remaining terms can be shown as follows:
Proof of Proposition 5.2 Similar to the proof of Prop. 5.1
lly H1 -yyalm
(I IIII II | ww/
118
Pw 1|
Iy
F Inference Procedure for SS-GPSR (BSLR with RFF)
We detail the inference procedure of SS-GPSR explained in $5.2 Although not explicitly mentioned in the main text, we can estimate both the parameter posterior and hyperparameters a. 3 in the BLR subroutine (22 [4 The equivalence of marginal likelihoods between BSLR and BLR (shown in Prop. 5.2) enables the use of existing BLR as subroutines, making the algorithm easier to implement.
Algorithm shows the pseudo code for inference procedure. Overall structure is analogous to Alg. At initialization (line 1-5), permutation matrix II and the response vector of pseudo coupled data (Definition 4.1) ý is set. Also, vector Sp, used in RFF is sampled from the kernel's spectral density, At each (algorithmic) time step I, BLR-step (line 7-8) and Correspondence Update (CU)-Step (line 9-13) are iteratively conducted. As the BLR-step, posterior of the parameters w and hyperparameter a, 3 is update by applying BLR using pscudo coupled data represented by X and y. Next, as CU-step, permutation matrix y is updated by sort operation.
G Implementation Details of Experiments
This appendix explains the details of our experimental setting and implementation.
Dataset: As stated in the Experiments section, we used four publicly available data sets, Airfoil, Concrete, Housing, and MPG. The description of each dataset and the pre- rrepprocessing applied are as follows. Airfoil records 1503 samples obtained from a series of aerodynamic and acoustic tests of two and three-dimensional airfoil blade sections with 5 features including angle, velocity and so on. Concrete records 1030 samples of the compressive strength of concrete with 8 features including the amount of cement, water and so on. Housing records 506 samples of home prices (MEDV: Median value of owner-occupied homes) in Boston with 13 socioeconomic and environmental features such as average number of rooms per dwelling (RM). MPG records 398 samples of the city-cycle fuel consumption in miles per gallon of automobiles with 7 features including weight and horsepower. We excluded samples with missing values and converted the categorical feature into one-hot vector for MPG. Input features were standardized for all data sets.
Oracle method (LR, GPR and DR): For oracle methods other than DR, used scikit- learr Specifically, LR uses sklearn.linear halpiiiiearmmll_lel_ LinoreeeiieaaReeesssion and GPR uses sklearn.gaussine process.
with default settings (we adopted Gaus- sian kernel for GPR). Details of DR are presented the next paragraph.
Neural satwork-based methods (DR and SDR): As stated in the Experiments section, we used one-hiiden-llayr feedforward neural network with the ReLU activation function for neural network- based methods (DR and SDR). Hyperparameters of these methods were set following [12]. The number of units was set to 20 for all problems. The parameters of DR were optimized using Adam 45 with learning rate of 0.001. The parameters of SDR were optimized by the stochastic sparse EM algorithm using Adam with a learning rate of 0.001. The mini-batch size of DR and that of SDR were 32 and 32/ L, respectively. The maximum number of epochs was 2000 in common. We used
https://scikitt-eese
25

---

Algorithm 2 Inference of SS-GPSR (BSLR using Random Fourier Features)
Input: dimension size of RFF H. shuffled data Dsn, maximum number of iterations T Tma Output: parameter posterior (w, A 11), hyperparameter a; B, and response vector of PCD y
1: /* Initialization */
2: Randomly initialize II (Optional: could be initialized by solving SLR).
3: Set the response vector of pseudo coupled data (Definition 4.1) as y + II'
4: Sort the elements of y; = (ya, ... Vyi)) Ascending sort simplify CU-Step.
5: Sample from the kernel's spectral density.
6: for 1 to T maz do
7: /* BLR-Step (Optimization of parameter's posterior distribution and hyperparameters) */
8: Update w, A -1 at, (This could done by e.g., BaysianRidge.fite y) in scikit-leam).
9: 1* Correspondence Update-Step (Optimization of II by sorting operation) */
for to .M do
11: ids f- argsort(v:) // Get the indices sorting elements of vi (= w) in ascending order.
12: y; + yi (idss //Re-orderring of y; by (argsort of) ids provides ý; <- II Iyy,, y;.
end for
14: end for
carly-stopping with validation data (coupled data for DR and shuffled data for SDR) 42) The above was implemented using PyTorch [46]. Experiments were run on a computer with Apple MI.
Proposed methods (GPSR and SS-GPSR): Here, provide details of the proposed methods (GPSR and SS- AS-GPSR). Their inference algorithms are summarized Algorithm [ (Appendix [c] and Algorithm (Appendix F respectively. As mentioned in the Experiments section, GPSR utilized warm initialization, i.e., the initial value of the permutation matrix is set using the estimated result of SLR. For SS-GPSR ten different runs were performed with various initial permutation matrices, and the solution that maximized the marginal likelihood was selected. For GPSR, the maximum number of iterations Tmaz
2
set to 100. The parameter for simulated annealing was configured with Smaz = 10 10-1 'N (Lg,(L), an initial temperature of T = 1.0 and a cooling rate y = 0.99. In the case of SS-GPSR, the maximum number of iterations Mea was set to 20. The dimension size of RFF H was fixed at 100 for all datasets. Note that for GPR and BLR steps in Algorithm T and [2] we used the seikit-learn analogous to oracle methods (BLR uses sklearn.lines hyperparameters such as precision a
model. Badoel.Byyssi So and 3 in (SS-) GPSR are estimated when applying GPR or BLR. So GPSR and SS-GPSR do not validation data unlike SDR.
26


Reviews:
[{"review_id": "wRVa9BFN6k", "reviewer": "Reviewer_fcj5", "summary": "The paper studies the shuffled regression problem from a Gaussian Process (GP) perspective. By introducing the GP, they additionally provide the uncertainty measure while improving the SOTA performance. Moreover, to make the method scalable, they introduce the spectral representation of the kernel to reduce the time complexity and memory consumption.", "questions": "The questions are listed in the weaknesses.\n\nBesides, I wonder what the computational complexity is with the spectral decomposition method. The authors only provided the memory complexity, but the computation complexity is missing.\n\nLastly, it seems that a deep version of the proposed method could be developed by replacing the kernel with a deep kernel, following the approach of Wilson et al. [1]. This extension may further boost performance and help adapt the method to high-dimensional or more complex data distributions. Exploring this direction, or at least discussing its feasibility, would enrich the contribution and position the work more clearly within the broader GP literature.\n\n[1] Wilson, Andrew Gordon, et al. \"Deep kernel learning.\" Artificial intelligence and statistics. PMLR, 2016.", "limitations": "The scalability of the proposed method remains uncertain when applied to complex, high-dimensional datasets. In particular, it is unclear how the number of spectral features ($H$) must grow and how computational cost scales in such settings.", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 2, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strengths:\n\n1. The paper clearly presents the idea of Gaussian Process Shuffled Regression (GPSR) and is overall well-written and easy to follow.\n2. It is the first to introduce Gaussian Processes (GPs) to the shuffled regression setting, addressing a novel and relevant problem.\n3. A scalable variant, BSLR, is proposed to address the computational limitations of GPSR.\n4. The paper demonstrates that GPSR and BSLR produce equivalent results under exact spectral approximation of the kernel, providing theoretical justification for the scalable method.\n5. The proposed method is evaluated on four real-world datasets and consistently outperforms existing state-of-the-art baselines.\n\nWeaknesses:\n1. Although the optimization details of GPSR are clearly described, the derivation or computation of the feature map $\\phi$ in BSLR is omitted, which negatively affects the clarity and completeness of the scalable method. It would also be helpful to include error bounds for the spectral approximation to better understand the conditions under which BSLR remains reliable.\n2. The SOTA baseline SDR is claimed to be a deep method, yet the implementation in the experiments uses only a single hidden layer. This may raise fairness concerns in the experimental comparison.\n3. The paper does not report training or inference time, despite the known scalability challenges of GP-based methods.\n4. There is no empirical scalability analysis. It would be valuable to see how the proposed method scales in terms of runtime, memory, and performance with increasing numbers of samples ($N$) or spectral features ($M$).\n5. All experimental datasets are low-dimensional, where a spectral dimension $H=100$ suffices. However, it remains unclear how $H$ needs to scale with higher-dimensional input, which is important for understanding the general applicability of the method.\n\nI would be inclined to raise my score if the authors address the following key concerns: (1) ensuring a fairer comparison by using deeper architectures for the SDR baseline as claimed; (2) reporting training and inference time as well as memory consumption; and (3) providing a clearer analysis or empirical evidence regarding the scalability of the proposed method, particularly in high-dimensional or large-scale settings.", "was_revised": true, "final_justification": "I have considered the rebuttal responses and other reviews. All my questions are resolved. Therefore, I'd like to raise my score.", "rebuttal": "We thank the reviewer for the positive comments and suggestions. Your comments will be incorporated into the camera-ready manuscript. Detailed responses to each comment are provided below. We believe that we have addressed all of the reviewer's concerns and would be happy to clarify anything that remains unclear.\n\n**1. Although the optimization details of GPSR are clearly described, the derivation or computation of the feature map \\phi in BSLR is omitted, which negatively affects the clarity and completeness of the scalable method. It would also be helpful to include error bounds for the spectral approximation to better understand the conditions under which BSLR remains reliable.**\n\nThe spectral approximation (RFF) for supervised learning, described in Eq. (3) in Sec. 2, can be directly applied in our method. Therefore, the theoretical support remains valid, including the uniform approximation bound by Rahimi & Recht (2007) and Sutherland & Schneider (2015). \n\n**2.The SOTA baseline SDR is claimed to be a deep method, yet the implementation in the experiments uses only a single hidden layer. This may raise fairness concerns in the experimental comparison.**\n\nWe designed the network architecture of SDR based on the original SDR paper [12].  Although the reviewer seems to have concerns about fairness, the current experimental setup is already advantageous for SDR using neural networks. This is because, as explained in the Baselines and Hyperparameters sections of Section 6, SDR utilizes early stopping based on validation data, while our proposed methods do not use any validation data.\n\n**3. The paper does not report training or inference time, despite the known scalability challenges of GP-based methods.**\n\nAs stated in our paper, the proposed method is developed by clarifying its connection to existing methods, Gaussian Process Regression (GPR) and Bayesian Linear Regression (BLR). Consequently, the training and inference times can be directly derived from these established methods.\n\nThe specific runtime for each method can be written as follows:\n- **Training time of GPSR**: (number of iterations) × {(GPR's hyperparameter optimization) + (simulated annealing)} \n- **Training time of SS-GPSR**: (number of iterations) × {(BLR's parameter and hyperparameter optimization) + (sorting operations)} \n- **Inference time of GPSR**: identical to that of Gaussian Process Regression (GPR) \n- **Inference time of SS-GPSR**: identical to that of Bayesian Linear Regression (BLR)\n\nTherefore, the investigation of training/inference time essentially boils down to investigating existing GPR and BLR algorithms, which is not the essential contribution of this study.\n\nNote that we report the convergence behavior of the objective function with respect to the number of iterations in Figure 3b.\n\n**4. There is no empirical scalability analysis. It would be valuable to see how the proposed method scales in terms of runtime, memory, and performance with increasing numbers of samples (N) or spectral features (M).**\n\nAs explained previously, the computation time and memory capacity are linked to the performance of the existing algorithms used as optimization subroutines (such as GPR and BLR), rather than being directly related to the contributions of our study.\n\n**5.All experimental datasets are low-dimensional, where a spectral dimension H=100 suffices. However, it remains unclear how H needs to scale with higher-dimensional input, which is important for understanding the general applicability of the method.**\n\nWe believe that this criticism is a result of a misunderstanding of current status of shuffled regression studies. Unlike supervised learning, shuffle regression studies primarily use low-dimensional data (and linear models). \n\nFor example, in a state-of-the-art work on linear models [11] (published in JMLR2024) and the SDR paper [12] (published in AAAI2024) focusing on nonlinearities, the input dimension size of their experiment is relatively small. Specifically, the experiments in [11] used data with input dimension size d=2, 4 (see Fig. 5 in [11]), and [12] used the UCI dataset similar to ours. \n\nWe respectfully ask that this gap between supervised learning and shuffled regression be taken into consideration during the discussion phase. \n\n**Besides, I wonder what the computational complexity is with the spectral decomposition method. The authors only provided the memory complexity, but the computation complexity is missing.**\n\nRandom Fourier features (RFF) approximate the spectral properties of a kernel using random sampling, rather than performing spectral decomposition directly. The computational cost of BLR using RFF described in Appendix A. \n\n**Lastly, it seems that a deep version of the proposed method could be developed by replacing the kernel with a deep kernel, following the approach of Wilson et al. [1]. This extension may further boost performance ...  Exploring this direction, or at least discussing its feasibility, would enrich the contribution and position the work more clearly within the broader GP literature.**\n\nWe appreciate your suggestions. We also believe the deep kernel is a promising approach for future shuffled regression studies, particularly in addressing high-dimensional problems. We will include this point in the camera-ready version."}, {"review_id": "7G3DtrefB9", "reviewer": "Reviewer_TTAf", "summary": "This paper introduces Gaussian Process Shuffled Regression (GPSR), a probabilistic model for learning regression functions when input-output correspondences are unknown, a setting referred to as *shuffled regression*. GPSR places Gaussian process priors over functions in the input space and iteratively recovers the correspondence using a Bayesian evidence maximisation framework. To address computational scalability, the authors propose a sparse variant (SS-GPSR) using random Fourier features, enabling closed-form correspondence updates via sorting. Theoretical connections to standard GPR and Bayesian linear models are established, leading to a principled and efficient inference scheme. Empirical results on benchmark datasets provide comparisons between (SS-)GPSR, classical shuffled regression methods (linear and deep shuffled regression), and oracle baselines using unshuffled data.", "questions": "The paper carefully describes and discusses the key aspects we could expect in this framework. I only have two remaining questions after reading: \n\n- Could you provide an empirical comparison of training/prediction times for all methods involved in the experiments? As SLR is used as initialisation in GPSR, I suspect it's orders of magnitude quicker (which is fine, but it would be valuable to know the exact time/accuracy tradeoff). For problems with large values of L, we can observe that the difference between (SS-)GPSR and SLR decreases; therefore, discussing this aspect may be more fair and would provide more nuanced takeaways to readers. \n\n- Although it is a key desirable property of GP-based frameworks, no evaluation of the uncertainty quantification (UQ) is provided. I understand that competing methods do not provide off-the-shelf uncertainty measures for comparison. However, it would be valuable at least to evaluate whether UQ is well-calibrated for (SS-)GPSR, for instance, by computing the number of testing points effectively lying within a credible interval. If a correct UQ is preserved in the (SS-)GPSR framework, it would be a significant benefit compared with existing approaches, and this feature would deserve to be highlighted accordingly. On the other hand, if UQ is somewhat damaged during the inference procedure, it could be presented as a natural and valuable path for future work, which could be mentioned in the conclusion.", "limitations": "Yes", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 4, "significance": 3, "originality": 3, "strengths_and_weaknesses": "The article is well-written and didactic, making it easy to follow throughout successive derivations. All involved methods are carefully detailed and re-explained in the Appendix for the most basic ones. The problem being addressed, extending GP regression to the case of shuffled measurements, is clearly positioned within existing literature. As a GP expert, I was unaware of this unusual *shuffled* context, which made it interesting to discover. Although the proposed method mainly builds on existing ideas, they are combined in a sensible way. In particular, Random Fourier Features (RFF) tackles both GP scalability and the quadratic assignment problem (QAP) at the same time, which is an elegant proposal. The experimental setting is convincing and effectively highlights the key properties of the methods (except for two important aspects, which are detailed in the *Questions* section). \n\nOverall, I believe that this article meets the threshold for publication and constitutes an incremental yet valuable contribution to both the GP and the shuffled regression literatures.", "was_revised": true, "final_justification": "Thank you for answering all my queries carefully and updating your future manuscript accordingly. I've read the other reviews and your replies to the different comments, and now I better understand the overall contribution. I believe this article provides a valuable contribution and expands the range of regression problems that can be addressed. I've upgraded my rating to 5.", "rebuttal": "We thank the reviewer for the positive comments and suggestions. Your comments and suggestions will be incorporated into the camera-ready manuscript. Detailed responses to each comment are provided below.\n\n**Could you provide an empirical comparison of training/prediction times for all methods involved in the experiments? As SLR is used as initialisation in GPSR, I suspect it's orders of magnitude quicker (which is fine, but it would be valuable to know the exact time/accuracy tradeoff).**\n\nAs stated in our paper, the proposed method is developed by clarifying its connection to the existing methods, Gaussian Process Regression (GPR) and Bayesian Linear Regression (BLR). Consequently, the training and inference times can be directly derived from these established methods.\n\nThe specific runtime for each method can be written as follows:\n- **Training time of GPSR:** (number of iterations) x {(hyperparameter optimization of GPR) + (simulated annealing)}\n- **Training time of SS-GPSR:** (number of iterations) x {(parameter and hyperparameter optimization of BLR) + (sorting operations)}\n- **Inference time of GPSR:** identical to GPR\n- **Inference time of SS-GPSR:** identical to BLR\n\nThus, comparing the training and prediction times essentially boils down to evaluating existing algorithms and optimization methods, including GPR, BLR, and Adam for SDR, which is not the primary contribution of this study.\n\nAdditionally, the reviewer's comment seems to show the reviewer is more interested in the proposed method's convergence behavior than in computational time. We would like to highlight that, as shown in Figure 3b, the proposed method achieves convergence in approximately 10 iterations, meaning that a large number of iterations is not required.\n\n**For problems with large values of L, we can observe that the difference between (SS-)GPSR and SLR decreases; therefore, discussing this aspect may be more fair and would provide more nuanced takeaways to readers.**\n\nThank you for your suggestions. Here are some clarifications: L represents the size of the set, which refers to the number of inputs in the input set as well as the number of responses in the response set. A larger L makes the shuffled regression problems more challenging. While we can confidently state that SS-GPSR significantly outperforms SLR for MPG data, regardless of the value of L, we will also mention in the camera-ready version that a larger L may reduce the performance difference when compared to SLR.\n\n**Although it is a key desirable property of GP-based frameworks, no evaluation of the uncertainty quantification (UQ) is provided. I understand that competing methods do not provide off-the-shelf uncertainty measures for comparison. However, it would be valuable at least to evaluate whether UQ is well-calibrated for (SS-)GPSR, for instance, by computing the number of testing points effectively lying within a credible interval … if UQ is somewhat damaged during the inference procedure, it could be presented as a natural and valuable path for future work, which could be mentioned in the conclusion.**\n\nThank you for your suggestions. Based on the reviewer's recommendation, we conducted a preliminary analysis of uncertainty quantification (UQ) for the data visualized in Figure 3a. The analysis indicates that approximately 33%, 55%, and 80% of the test samples fall within ±1, ±2, and ±3 predicted standard deviations from the predictive mean, respectively. While we cannot draw definitive conclusions due to the presence of outliers in the test data illustrated in Figure 3a, it seems that the uncertainty quantification (UQ) is somewhat damaged. (ideally, we would expect 68.3%, 95.4%, and 99.7% of the samples to fall within ±1, ±2, and ±3 standard deviations, respectively, for a normal distribution).\n\nOne possible reason for the damage to UQ is that our current approach relies on a \"point estimate\" of the permutation matrix without incorporating a prior distribution. To address this, we may consider using regularization techniques for the permutation matrix [20], and exploring Lalchand's fully Bayesian approach to Gaussian Process Regression (GPR) for supervised regression [32] as a future direction. We will include additional analysis of UQ and the above discussions in the camera-ready version of our paper."}, {"review_id": "LWvmSG9VXO", "reviewer": "Reviewer_HE9T", "summary": "This paper proposes Gaussian Process Shuffled Regression (GPSR) and its scalable variant SS-GPSR for the problem of learning regression functions from shuffled data, where the input–output correspondence is unknown. The authors formulate a Bayesian evidence maximization framework in which both the regression function and permutation matrix are inferred. GPSR uses a Gaussian process prior, and inference alternates between Gaussian process regression and solving a quadratic assignment problem. To address the scalability challenges of GPSR, the authors develop SS-GPSR, a sparse approximation based on random Fourier features and a weight-space formulation. A key technical insight is the introduction of a latent variable $w$, which yields a marginal likelihood that matches that of GPSR at optimality but depends only linearly on the permutation matrix, enabling efficient updates via sorting. Experimental comparisons against baselines show that SS-GPSR achieves competitive accuracy with improved scalability.", "questions": "1. **Interaction of kernel and permutation optimization**\n   In Figure 2, GPSR is applied to a toy sine wave task. Can you provide insight into how optimization of the kernel hyperparameters and the permutation matrix interact? For example, if the data is highly shuffled, will the GP favor a short length scale and overfit? Could this lead to poor local optima?\n\n2. **Lack of kernel optimization in SS-GPSR**\n   SS-GPSR does not optimize the kernel hyperparameters $\\theta$. Is this the main price paid for scalability? In settings where tuning $\\theta$ is important for performance, how does this limitation affect results?\n\n3. **Convergence behavior of alternating updates**\n   Since $w$ and $\\Pi$ are each optimized globally in SS-GPSR, is there any convergence guarantee for the alternating procedure? Could the method oscillate between solutions with the same likelihood? Empirical observations or theoretical discussion here would be helpful.\n\n4. **Runtime and scalability experiments**\n   A central claim is that SS-GPSR is more efficient than GPSR and SDR. To support this, it would be helpful to include runtime-to-convergence plots for varying $L$, and wall-clock comparisons across methods.\n\n5. **Larger-scale experiments**\n   Current experiments are limited to relatively small datasets. If feasible, then I think evaluating performance at larger $N$ and especially larger $L$ would better illustrate the scalability of SS-GPSR and its competitiveness with neural methods like SDR.\n\n\nIf the authors can provide meaningful insight into the optimization dynamics (especially around kernel/permutation interactions), clarify the limitations of not tuning $\\theta$ in SS-GPSR, and supply runtime or larger-scale evaluations, I would be inclined to raise my score.", "limitations": "Yes", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 2, "strengths_and_weaknesses": "#### **Strengths**\n\n* The paper presents a clear, principled formulation of shuffled regression using Gaussian processes and evidence maximization.\n* The derivation of SS-GPSR is elegant, with the introduction of the auxiliary variable $w$ allowing efficient optimization through a linear dependence on the permutation matrix.\n* The paper is clearly written and accessible, even through dense derivations.\n* The closed-form inference and uncertainty quantification offered by the GP-based models set them apart from purely neural alternatives.\n\n#### **Weaknesses**\n\n* While the mathematical development is thorough, the core approach ultimately reduces to kernel ridge regression with a permutation matrix estimated via alternating optimization. From this perspective, the methodological contribution may feel incremental, though well-executed.\n* Experiments are relatively limited in scale and do not fully evaluate how SS-GPSR behaves as dataset size or set size $L$ increases.\n* SS-GPSR does not optimize kernel hyperparameters ($\\theta$), and the impact of this design choice on performance is not explored.\n* No runtime comparisons are provided, even though efficiency is one of the main selling points of SS-GPSR.", "was_revised": true, "final_justification": "Thank you for your responses to all my questions. My main concern was regarding kernel hyperparameter optimization, as I'm not sure how the basis functions can be optimized once sampled. However, afterreading the rest of the reviews and rebuttal, I don't think this should bar the paper, and so I've raised my score.\n\nI believe the authors adequately addressed my comment", "rebuttal": "We thank the reviewer for the insightful comments and suggestions. Your comments and suggestions will be incorporated into the camera-ready manuscript. Detailed responses to each comment are provided below. We believe that we have addressed all of the reviewer's concerns and would be happy to clarify anything that remains unclear.\n\n**1. Interaction of kernel and permutation optimization: In Figure 2, GPSR is applied to a toy sine wave task. Can you provide insight into how optimization of the kernel hyperparameters and the permutation matrix interact? For example, if the data is highly shuffled, will the GP favor a short length scale and overfit? Could this lead to poor local optima?**\n\nIn cases with a small amount of data and a poor initial value of the permutation matrix, it is possible to experience overfitting, as the reviewer pointed out. However, with a sufficient amount of data, such overfitting is unlikely to occur because the \"bad\" permutation matrix necessitates a substantial variation in the predicted mean or a large value of the predicted variance to accurately represent the Pseudo Coupled Data (PCD) generated by the “bad” permutation. The use of the evidence maximization framework also contributes to avoiding overfitting.\n\n**2. Lack of kernel optimization in SS-GPSR: SS-GPSR does not optimize the kernel hyperparameters \\theta. Is this the main price paid for scalability? In settings where tuning θ is important for performance, how does this limitation affect results?**\n\nIn SS-GPSR, optimization of alpha, used in a prior distribution of parameter w, is performed instead of kernel hyperparameters. Indeed, the current optimization of scalar-alpha may seem to be less expressive, but our method can also be applied to more flexible settings, such as automatic relevance determination using a vector-alpha \\{\\alpha_h \\}^H_{h=1} having different values for each dimension h. Therefore, we believe that this point does not necessarily represent a limitation of the proposed method. \n\n**3. Convergence behavior of alternating updates: Since w and Π are each optimized globally in SS-GPSR, is there any convergence guarantee for the alternating procedure? Could the method oscillate between solutions with the same likelihood? Empirical observations or theoretical discussion here would be helpful.**\n\nEmpirically, as is shown in Figure 3b, the objective function value converges to a certain value for each initial value setting. We have checked the numerical values up to 15 decimal places, but no fluctuations in the objective function values are observed, and no results suggesting that the solution oscillates have been obtained so far.\n\n**4. Runtime and scalability experiments: A central claim is that SS-GPSR is more efficient than GPSR and SDR. To support this, it would be helpful to include runtime-to-convergence plots for varying L, and wall-clock comparisons across methods.**\n\nAs stated in our paper, our methods are developed by clarifying its connection to the existing methods, Gaussian Process Regression (GPR) and Bayesian Linear Regression (BLR). Therefore, the training and inference times are immediately derived from these existing methods. \n\nThe specific runtime for each method can be written as follows::\n- **Training time of GPSR:** (number of iterations) x {(hyperparameter optimization of GPR) + (simulated annealing)}\n- **Training time of SS-GPSR:** (number of iterations) x {(parameter and hyperparameter optimization of BLR) + (sorting operations)}\n- **Inference time of GPSR:** identical to GPR\n- **Inference time of SS-GPSR:** identical to BLR\n\nTherefore, the comparison of computation time essentially boils down to evaluating existing GPR and BLR algorithms and optimization methods for deep learning such as Adam, which is not the essential contribution of this study. Note that we report the convergence behavior of the objective function with respect to the number of iterations in Figure 3b.\n\n**5. Larger-scale experiments: Current experiments are limited to relatively small datasets. If feasible, then I think evaluating performance at larger N and especially larger L would better illustrate the scalability of SS-GPSR and its competitiveness with neural methods like SDR.**\n\nThe dataset sizes of our experiment are almost similar to that used in the experiments of e.g., Xu’s work [15] and the SDR paper [12] (these two papers also used the UCI-ML dataset similar to ours). \n\nWe believe that this criticism is a result of a misunderstanding of the current status of shuffled regression studies. Unlike supervised learning, shuffled regression studies primarily use small-scale, low-dimensional data and linear models (see, for example, state-of-the-art work on linear models [11] published in JMLR2024, Xu’s work [15] published in NeurIPS2019, and the SDR paper [12] in AAAI2024). \n\nWe respectfully ask that this gap between supervised learning and shuffled regression be taken into consideration in the discussion phase."}, {"review_id": "Jeqx8JMqT0", "reviewer": "Reviewer_aYMq", "summary": "This paper introduces a novel Gaussian Process-based framework for shuffled regression, termed Gaussian Process Shuffled Regression (GPSR), along with its sparse variant, SS-GPSR. The goal of shuffled regression is to learn a regression function from \"shuffled data\" where the correspondence between input features and target responses is unknown. GPSR leverages Gaussian Processes (GP) as a functional prior, employing kernel functions for nonparametric modeling and providing principled quantification of predictive uncertainty. The authors establish the theoretical equivalence between the marginal likelihood and predictive distribution of GPSR and those of standard Gaussian Process Regression (GPR). Building on this, they develop an alternating optimization algorithm that iteratively updates kernel hyperparameters (via standard GPR subroutines) and input-output correspondences (by solving a Quadratic Assignment Problem, QAP). To improve computational efficiency, the authors further propose SS-GPSR, which approximates GPSR using Random Fourier Features (RFF) to map the model into weight space, resulting in a Bayesian Shuffled Linear Regression (BSLR) formulation. This enables closed-form updates of the correspondence via sorting operations. Empirical results on multiple UCI datasets demonstrate that GPSR and SS-GPSR outperform existing shuffled regression approaches, particularly in modeling nonlinear relationships and quantifying uncertainty.", "questions": "The authors are encouraged to address the following critical issues through additional analysis or experiments, as they directly impact the paper’s overall evaluation:\n\n(1) SS-GPSR updates the correspondence matrix $\\Pi_i$ via sorting of the predicted means $v_i$ and observed responses $y_i$ (\\S5.2). However, when the underlying regression function is non-monotonic or the noise level is high, this procedure may yield incorrect pairings. A theoretical analysis is needed to specify the optimality conditions for sorting---e.g., under what smoothness or monotonicity assumptions Proposition 5.1 holds. In addition, synthetic experiments should be conducted to quantify how the sorting error rate and regression MSE vary with increasing noise levels, similar to Figure 2.\n\n(2) Current experiments are limited to low-dimensional tabular datasets. The effectiveness of GPSR and SS-GPSR in higher-dimensional feature spaces (e.g., $D > 100$) or under non-stationary conditions (e.g., covariate shift) remains unexplored. Please include experiments on at least one high-dimensional dataset (e.g., CIFAR-10 image embeddings) or temporally shuffled sensor data collected from heterogeneous devices. Comparative analysis with SDR in these settings would strengthen the empirical claims.\n\n(3) While SS-GPSR is claimed to be more computationally efficient than GPSR, no quantitative comparisons of runtime or memory usage are provided. Please include such metrics and analyze the computational trade-offs, particularly contrasting the theoretical complexities $\\mathcal{O}(N^3)$ (GPSR) and $\\mathcal{O}(H^3)$ (SS-GPSR).\n\n(4) In GPSR, the QAP subproblem is solved via simulated annealing (\\S4.2), a stochastic optimization approach that lacks direct theoretical grounding in Bayesian inference. This introduces a potential disconnect with the marginal likelihood maximization framework. Please discuss how simulated annealing affects convergence and posterior estimation—possibly in comparison with MCMC methods~\\cite{ref20}—or provide a theoretical or empirical guarantee of the monotonic increase of the marginal likelihood $\\mathcal{L}_{\\text{GPSR}}$ during the iterative process.", "limitations": "While the authors briefly acknowledge certain limitations and social implications in the conclusion section, several critical aspects require further elaboration:\n\n(1) GPSR assumes the underlying regression function follows a Gaussian Process (GP) prior, yet the paper does not discuss potential failure cases when this assumption is violated—for example, in the presence of heavy-tailed or non-Gaussian data distributions. It is recommended to include a synthetic experiment in the appendix that compares GPSR with quantile-based shuffled regression methods when the ground-truth function f deviates significantly from the GP prior.\n\n(2) SS-GPSR relies on Random Fourier Features (RFF), which inherently assume a stationary (translation-invariant) kernel (§3). However, the applicability of the method to non-stationary kernels is not examined, such as deep kernels [42]. The paper should clarify this kernel restriction and discuss possible extensions or limitations regarding kernel choice.\n\n(3) Although the conclusion mentions that prediction performance should be evaluated with caution, it does not address the concrete risks of incorrect input-output alignment, such as patient-label mismatches in medical records. It would be beneficial to discuss whether a confidence-based rejection mechanism is necessary for privacy-sensitive or high-risk applications, e.g., imposing an uncertainty threshold to discard low-confidence predictions.", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strengths:\n\n(1) The proposed method is rigorously designed with clear theoretical underpinnings. In particular, Propositions 4.2 and 4.3 establish a theoretical equivalence between GPSR and standard GPR, and the experimental evaluation is comprehensive, involving multiple datasets and varying set sizes.\n\n(2) The work is the first to introduce Gaussian Processes into the shuffled regression setting and further develops an efficient sparse variant, SS-GPSR.\n\n(3) Shuffled regression has practical relevance in scenarios involving privacy-preserving data (e.g., medical records) and multi-source data integration. The adoption of nonparametric modeling enhances the expressive capacity of existing approaches in these settings.\n\nWeaknesses:\n\n(1) SS-GPSR updates the input-output correspondence via a sorting operation (§5.2), which is optimal only under the assumption of a strict rank alignment between the predicted means v_i and responses y_i. The paper does not analyze the robustness of this assumption under noise or model misspecification, nor does it quantify the impact of sorting errors on regression performance. Notably, the derivation in Appendix F omits consideration of distributional shifts in the noise.\n\n(2) The experiments are limited to low-dimensional UCI datasets (with a maximum feature dimension D=13), and the method's robustness in high-dimensional settings (e.g., image embeddings) or under complex noise structures (e.g., heteroscedastic noise) remains untested. Moreover, SS-GPSR exhibits noticeable performance variance at L=16 in Table 3 (e.g., MPG dataset shows a standard deviation of 3.03 in MSE), but the source of this instability is not discussed.\n\n(3) Implementation details of the simulated annealing algorithm in Appendix C are incomplete—for example, the temperature schedule T and the acceptance probability q are unspecified. Furthermore, the paper does not explain how the QAP solver avoids convergence to suboptimal local minima (§4.2).\n\n(4) Although SS-GPSR is claimed to reduce computational complexity, the paper does not provide runtime comparisons with existing baselines (e.g., SDR), nor does it analyze the joint impact of the RFF dimension H and the set size L on computational cost (§6).", "was_revised": true, "final_justification": "I have read the other relevant comments and counterarguments. All my questions have been answered. I am very grateful to the authors for their detailed rebuttal responses, and therefore I am willing to increase my rating.", "rebuttal": "We thank the reviewer for the insightful comments and suggestions. Your suggestions will be incorporated into the camera-ready manuscript. Detailed responses to each comment are provided below.\n\n**(1) SS-GPSR updates the correspondence matrix via sorting of the predicted means vi and observed responses yi (\\S5.2) ... A theoretical analysis is needed to specify the optimality conditions for sorting---e.g., under what smoothness or monotonicity assumptions Proposition 5.1 holds. In addition, synthetic experiments should be conducted to quantify how the sorting error rate and regression MSE vary with increasing noise levels, similar to Figure 2.**\n\nThere may be some confusion regarding our method and the focus of our study. Let us clarify.\n\n*Regarding the optimality of the sorting operation in our method:*\nThe term \"optimization\" (or \"optimality\") used for the permutation matrix in our paper refers specifically to the optimization of the log-marginal likelihood, which acts as the loss function given other parameters (w, alpha, and beta). This interpretation differs from the definition of \"optimality\" mentioned by the reviewer. Furthermore, Proposition 5.1 does not require any additional assumptions, contrary to the reviewer's critique.\n\n*About the position of our Study:*\nWe recognize that the analysis of permutation recovery is performed in some theoretical papers for shuffle regression. However, the main contribution of our research is the proposal of a new GP-based method for shuffled regression. There are various methodological studies, such as the SLR paper [4], Xu's work [15] (presented at NeurIPS 2019), and the SDR paper [12] (presented at AAAI 2024), that do not provide the analysis of permutation recovery.\n\nFurthermore, our primary objective is to learn a regression function rather than to enhance permutation recovery performance. As the reviewer noted, when true function values (noise-free values) are very similar, it is likely that the assignments may not reflect true pairings due to alterations in relationships caused by noise. However, the impact on function estimation performance tends to be minimal if the observed noisy values are similar. In fact, we have confirmed through our experiments on several real noisy datasets that our proposed method achieves performance that is either superior to or competitive with state-of-the-art techniques.\n\n**(2) Current experiments are limited to low-dimensional tabular datasets. The effectiveness of GPSR and SS-GPSR in higher-dimensional feature spaces (e.g., D>100) or under non-stationary conditions (e.g., covariate shift) remains unexplored. Please include experiments on at least one high-dimensional dataset (e.g., CIFAR-10 image embeddings) or temporally shuffled sensor data collected from heterogeneous devices. Comparative analysis with SDR in these settings would strengthen the empirical claims.**\n\nWe believe that this criticism is a result of a misunderstanding of current status of shuffled regression studies. Unlike supervised learning, shuffle regression studies primarily use low-dimensional data (and linear models). \n\nFor example, in a state-of-the-art work on linear models [11] (published in JMLR2024) and the above mentioned papers [12][15], the input dimension sizes used in their experiments are relatively small. Specifically, the experiments in [11] utilized datasets with input dimensions of d=2 and d=4 (see Fig. 5 in [11]), while [12][15] employed the UCI ML dataset, similar to ours.\n\nWe kindly request that this *gap between supervised learning and shuffled regression* be considered during the discussion phase.\n\n**(3) While SS-GPSR is claimed to be more computationally efficient than GPSR, no quantitative comparisons of runtime or memory usage are provided. Please include such metrics and analyze the computational trade-offs, particularly contrasting the theoretical complexities O(N3) (GPSR) and O(H3) (SS-GPSR).**\n\nAs stated in section 4.2, the optimization of the permutation matrix in GPSR is the quadratic assignment problem (QAP), which is NP-hard. In contrast, the permutation optimization in SS-GPSR is a simple instance of the linear assignment problem (LAP) that can be solved by sorting operations~(see Sec. 5.2). This clearly shows the efficiency of SS-GPSR. \n\nIf the reviewer's concerns focus on the parameter optimization steps, our theoretical analysis enables the use of existing optimization methods for Gaussian Process Regression (GPR) and Bayesian Linear Regression (BLR). Consequently, the runtime (both training and inference time) and memory usage can be directly derived from these established methods. \n\nThe specific runtime for each method can be written as follows:\n- **Training time of GPSR:** (number of iterations) x {(hyperparameter optimization of GPR) + (simulated annealing for QAP)}\n- **Training time of SS-GPSR:** (number of iterations) x {(parameter and hyperparameter optimization of BLR) + (sorting operations for LAP)}\n- **Inference time of GPSR:** identical to GPR\n- **Inference time of SS-GPSR:** identical to BLR\n\nTherefore, investigating training/inference time essentially boils down to evaluating existing GPR and BLR algorithms, which is not the essential contribution of our study. \n\n**(4) In GPSR, the QAP subproblem is solved via simulated annealing (\\S4.2), a stochastic optimization approach that lacks direct theoretical grounding in Bayesian inference. This introduces a potential disconnect with the marginal likelihood maximization framework.**\n\nThank you for your insightful comment. It is important to note that in our (SS-)GPSR formulation, the permutation matrix is treated as one of the hyperparameters, similar to the kernel's parameter. We acknowledge that our method is not fully Bayesian; however, it is still a Bayesian approach, similar to the classical GPR which does not put a prior on the kernel parameter. This formulation enables us to develop an easy-to-implement algorithm and avoids the need for Markov Chain Monte Carlo (MCMC) methods. We agree that extending GPSR to a full Bayesian framework could be a promising direction for future work, and the research by Lalchand on supervised regression [32] could be beneficial in this regard.\n\n**(1) GPSR assumes the underlying regression function follows a Gaussian Process (GP) prior, yet the paper does not discuss potential failure cases when this assumption is violated—for example, in the presence of heavy-tailed or non-Gaussian data distributions.**\n\nIn the limitations paragraph of the Conclusion section, we clearly state that our methods rely on the assumption of Gaussian Processes (GP). We will add the examples provided by the reviewer (heavy-tailed or non-Gaussian data distributions) in the camera-ready. \n\n**It is recommended to include a synthetic experiment in the appendix that compares GPSR with quantile-based shuffled regression methods when the ground-truth function f deviates significantly from the GP prior.**\n\nThank you for your suggestion, but we have not been able to find the quantile-based shuffle regression methods mentioned by the reviewer. Since shuffle regression is still in its infancy compared to supervised regression, we have some doubts about whether such a method actually exists. If the reviewer provides us with more information on some quantile-based shuffled regression methods, such as the paper title, we will consider the possibility of adding it as a supplement to the camera-ready version. \n\n**(2) SS-GPSR relies on Random Fourier Features (RFF), which inherently assume a stationary (translation-invariant) kernel (§3)...The paper should clarify this kernel restriction and discuss possible extensions or limitations regarding kernel choice.**\n\nWe will clearly describe that the current SS-GPSR focuses on transition-invariant kernel but the further extensions would be possible for handling non-stationary kernel such as deep kernel in the camera-ready. \n\n**(3) Although the conclusion mentions that prediction performance should be evaluated with caution, it does not address the concrete risks of incorrect input-output alignment, ... It would be beneficial to discuss whether a confidence-based rejection mechanism is necessary for privacy-sensitive or high-risk applications.**\n\nThank you for your suggestion. Although the main focus of our study is on learning regression functions, not permutation recovery, we agree that caution should be exercised if the reader were to use the estimation results of the permutation matrix. We will include this point in the camera-ready version."}]

Author rebuttal:

Meta-review: The authors propose a GP-based model for shuffled regression, i.e. regression in which sets of data points map to sets of responses, but the correspondence is not known in advance.  Existing methods for the problem include a linear regression approach and a neural network based approach; the authors propose a GP together with a permutation chosen by evidence maximization.  In particular, they propose two methods: one that uses a function space perspective and computes the permutations via a (hard) quadratic assignment problem; and one that uses Bayesian linear regression with random Fourier features and computes the permutations via a (much simpler) linear assignment problem.  The accuracy of the new methods is compared to the existing methods on a number of benchmarks.

I would summarize the paper as having three main pieces:

- The idea of using a GP to solve a shuffled regression problem at all (and using an evidence maximization approach, though this is done in a pretty standard way)
- The idea of converting to a feature-vector formulation and computing the permutation by a linear assignment problem.
- The idea of integrating all this with a low-rank approximation to the kernel and solving as a weight space form / regularized least squares problem.

The first two pieces seem to be new; the third piece is standard.  The exposition lumps the second two pieces together (though one could use a Cholesky factorization of the kernel matrix to get to a linear assignment problem, so the use of a low-rank approximation is not really fundamental here) and only considers a single flavor of low-rank approximation via random Fourier features (though an SVD, Nystrom approximation, or other approach would work without any real modification).

The paper is well-written, and the authors generally liked the idea.  However, there were two broadly-shared concerns.  First, several reviewers commented that the input space was rather low-dimensional in all the examples, and asked what would happen for higher-dimensional problems.  The authors replied that most of the shuffled regression literature deals with fairly low-dimensional inputs (unlike in more general supervised learning settings), and asked that this be taken into account.  This response seemed to basically satisfy the reviewers (or at least it did not prompt further discussion).  Second, the reviewers all pointed out that despite the claim that one selling point of the methods is scalability -- particularly for the random-feature-based method -- there are no timing experiments.  The authors essentially replied that the new algorithms use known building blocks with known performance, that one can reason about the performance in terms of those building blocks, and that "investigating existing GPR and BLR algorithms [...] is not the essential contribution of this study."  The reviewers again did not comment further on this, though I still find the lack of concrete timing numbers somewhat unsatisfactory.

Two other significant points came up in the reviews and were addressed by the authors.  First, two reviewers (fcj5 and aYMq) pointed out that the SS-GPSR method as described seems to restrict one to stationary kernels; the authors replied that this is not a fundamental restriction, and it will be addressed briefly in the camera-ready.  Second, reviewer TTAf mentioned that no uncertainty quantification was actually shown in the paper, though this is one of the benefits of using a GP-based method.  The authors replied that they had done some preliminary experiments in response to this question (which they shared); the results were not entirely satisfactory, but they said they would say more in the camera ready.

The key reasons to accept would be that this is a new approach to an interesting class or problems (shuffled regression), and it seems to achieve pretty good accuracy compared to competitors.  The key shortcoming is a lack of experimental study related to performance, though I also think separately treating the conversion to a linear assignment problem and the approximation via a low-rank method (of which RFF is just one example) might clarify some of the problems and help with reader questions about applicability to non-stationary kernels similar to what the reviewers have.
Decision: Accept (poster)