Shuffled data problems refer broadly to situations in which the goal is to perform inference for a functional of the joint distribution of a pair of random variables (x, y) (such as, e.g., their covariance) based on separate samples {x i } n and {y i } m i=1 that involve matching pairs {(x ⇡ ⇤ (i) , y i )} m i=1 pertaining to the same statistical unit, where the map ⇡ ⇤ : {1, . . . , m} ! {1, . . . , n} may only be observed incompletely. This is a rather common scenario in data integration problems in which different pieces of information about a shared set of entities reside in multiple data sources that need to be combined in order to perform a given data analysis task. The process of identifying matching parts across two or more files is often far from trivial in the absence of unique identifiers, and has thus grown into a vast and active field of research known as record linkage, e.g., Binette and Steorts (2020). The above shuffled data model represents a direct approach to account for mismatches in record linkage and the impact on downstream data analysis. Shuffled data problems were first systematically discussed in DeGroot et al. (1971), with little progress until a few years ago given advances in computation (Gutman et al., 2013). Recently, shuffled data problems have generated widespread interest, fueled by applications in signal processing (Unnikrishnan et al., 2018;Pananjady et al., 2018), correspondence problems in computer vision (Pananjady et al., 2017;Li et al., 2021) and NLP (Grave et al., 2019;Shi et al., 2021), biomedical data analysis (Ma et al., 2021a;Abid and Zou, 2018), and data privacy (Domingo-Ferrer and Muralidhar, 2016;Gordon et al., 2021).
On the theoretical side, several papers have investigated the statistical limits of signal estimation and permutation recovery in unlabeled sensing in which the goal is to recover a signal ✓ ⇤ from n noisy linear measurements
, where ⇡ ⇤ is an unknown permutation, i.e., m = n and ⇡ ⇤ is one-to-one (Unnikrishnan et al., 2018;Pananjady et al., 2018;Hsu et al., 2017;Abid et al., 2017;Tsakiris and Peng, 2019). Another line of research has studied the setting in which x and y are scalar and related by a monotone map (Carpentier and Schlüter, 2016;Rigollet and Weed, 2019;Flammarion et al., 2019;Ma et al., 2020;Balabdaoui et al., 2021).
A common conclusion from these works is that shuffled data problems are generally plagued by both statistical and computational challenges. First, the combinatorial nature of ⇡ ⇤ makes it hard to devise computationally tractable approaches with provable guarantees. Existing algorith-mic "solutions" involve integer programming (Tsakiris and Peng, 2019;Peng and Tsakiris, 2020;Mazumder and Wang, 2021), the EM algorithm (Gutman et al., 2013;Abid and Zou, 2018;Tsakiris et al., 2020), sampling and approximate inference (McVeigh et al., 2019;Steorts et al., 2016;Klami, 2012). Regardless of the computational challenges, shuffled data problems tend to be highly susceptible to noise and prone to overfitting. In fact, statistical guarantees typically involve unrealistically stringent signal-tonoise requirements (Pananjady et al., 2018;Hsu et al., 2017). Loosely speaking, this issue results from the fact that the set of permutations grows rapidly in size with n. This observation suggests that suitable forms of regularization hinging on prior information on ⇡ ⇤ are needed to constrain the size of the parameter space under consideration. Several papers consider partial shufflings in which varying fractions of (x i , y i )-pairs are already observed with the correct correspondence (Slawski and Ben-David, 2019;Slawski et al., 2019Slawski et al., , 2020;;Zhang and Li, 2020;Peng et al., 2021), and only the remaining portion of the data is subject to shuffling. Another constraint commonly encountered in record linkage is that ⇡ ⇤ is block-structured with known composition of the blocks based on auxiliary variables that are required to agree for matching records. In domains such as signal processing and computer vision, ⇡ ⇤ is often constrained to act locally in the sense that indices are shuffled only within small time windows or image regions (Ma et al., 2021b;Abbasi et al., 2021).
The goal of the present paper is the development of a regularization framework for shuffled data problems that integrates those and other constraints in a unified way. To that end, we introduce an exponential family prior on the permutation group that is flexible enough to accommodate any kind of prior information that can be expressed solely in terms of index pairs (i, j). This prior turns out to be conjugate for canonical shuffled data problems in which the likelihood conditional on a fixed permutation can be expressed as product over the corresponding (x, y)-pairs. Inference is based on the MC-EM algorithm considered in Wu (1998) and Abid and Zou (2018). We show that for exponential family likelihood, the resulting M-step is particularly scalable since it involves n instead of n 2 terms. Moreover, computation of the MAP estimator of ⇡ ⇤ with the remaining parameters fixed reduces to a linear assignment problem, and hence remains computationally tractable. Theoretical results and a collection of experiments for various shuffled data setups demonstrate the usefulness of regularization based on the proposed prior in comparison to the unregularized counterpart and other baselines.
Notations. We denote by D = {(x i , y i )} n i=1 the observed merged data, subject to shuffling. We use (x, y) for a generic pair of matching records. We use X and Y for the row-wise concatenation of {x i } n i=1 and {y i } n i=1 , respectively. We let p(•) denote the density (PDF) of a list of variables in (•), and accordingly p(• | •) is used for conditional PDFs. We write u ⇠ p to express that random variable u has density p. The symbol E (...) [•] is the expectation w.r.t. (. . .). The Hamming distance on the permutation group P(n) of [n] = {1, . . . , n} is denoted by d H . The symbol tr is used for the matrix trace, and I n denotes the identity matrix of dimension n. The cardinality of a set is denoted by | • |, and I denotes indicator function.
Conventions. We often refer to a permutation via the underlying map ⇡ and the associated matrix ⇧ = (⇡ ij ) in an interchangeable fashion, and accordingly P(n) and subsets thereof may refer to both maps and matrices. Asterisked symbols such as ⇡ ⇤ , ✓ ⇤ , ⇤ etc. refer to ground truth parameters; non-asterisked symbols such as ⇡, ✓, etc. refer to generic elements of the associated parameter spaces.
Our approach will be presented as follows: we start with a brief motivation, followed by a more formal systematic introduction, and conclude with technical details pertaining to computation and model fitting.
) n i=1 and corresponding amplified slope b ML . R: Re-paired data (x b ⇡(i) , y i ) n i=1 based on the proposed Hamming prior.
Consider the simple linear regression setup y i =
x ⇡ ⇤ (i) ⇤ + ⇤ ✏ i , where x i and ✏ i are independent standard normal random variables, 1 i n, and ⇡ ⇤ permutes 10% of the indices uniformly at random. Suppose that the sign of ⇤ is known to be positive. Then the ML estimator of ⇡ ⇤ (or equivalently, the MAP estimator under a uniform prior over P(n)) is given by the permutation b ⇡ ML that matches the corresponding order statistics in {x i } n i=1 and {y i } n i=1 :
As shown in Figure 1, b ⇡ ML performs rather poorly. The scatterplot of the matching of corresponding order statistics is far from that of the underlying correct pairing. In fact, b ⇡ ML is associated with massive overfitting. Specifically, let
denote the resulting ML estimators of ⇤ and 2 ⇤ , respectively. It is straightforward to show that
) based on the proposed prior.
in probability as n ! 1. This is alarming since it implies that the least squares fit absorbs all the noise.
Figure 1 shows that the ML estimator is too aggressive in forming "corrected" pairs (x b ⇡ML(i) , y i ) given that only 10% of the observations are actually mismatched, and among those 10%, only a fraction contributes substantial mismatch that exceeds the noise inherent in the problem. Sparsity of ⇡ ⇤ is often a reasonable assumption in post-linkage data analysis (Chambers and Diniz da Silva, 2020;Slawski and Ben-David, 2019), where sparsity here means that set of mismatches {i 2 [n] : ⇡ ⇤ (i) 6 = i} has significantly smaller cardinality than n. Given an upper bound on the number of mismatches, say k, it is appropriate to consider the following constrained ML estimator of ⇡ ⇤ :
where id is the identity map on [n] and d H (⇡, ⇡ 0 ) = P n i=1 I(⇡(i) 6 = ⇡ 0 (i)) denotes the Hamming distance on P(n). To the best of our knowledge, there is no efficient algorithm for computing the maximizer directly. However, there exists a Lagrangian multiplier > 0 such that (3) is equivalent to the optimization problem
which is a linear assignment problem with cost matrix C = I(i 6 = j) x j y i , which is computationally tractable according to the discussion following (7) below. The right panel of Figure 1 highlights the improvement that can be achieved by the resulting estimator which here only makes a small number of re-pairings of (x, y) capturing pairs that correspond to massive mismatch error in the left panel.
Figure 2 illustrates scenarios in which ⇡ ⇤ is not sparse (with a mismatch rate exceeding 80%), but constrained to be a "local shuffling" in the sense that max i2[n] |⇡ ⇤ (i) i| r, i.e., the corresponding permutation matrix is a band matrix with bandwidth at most r. This scenario is particularly relevant when the data is recorded sequentially (e.g., over different time points) or across a spatial domain endowed with a notion of distance, and it is known that ⇡ ⇤ can only mix up the order of data inside a specific time window or within a local neighborhood. There are numerous applications in which ⇡ ⇤ is locally constrained such as genome se-quencing (Abid et al., 2017), signal processing (Balakhrisnan, 1962;Abbasi et al., 2021), or computer vision (Ma et al., 2021b). The illustrative example in Figure 2 can be thought of a regression problem in which the signal is a sine with known frequency but unknown (positive) amplitude ⇤ , i.e., y ti = ⇤ sin(10⇡t i )+0.1✏ i , 1 i n (left panel). However, the observed data is of the form (y t ⇡ ⇤ (i) ) n i=1 for some unknown (local) permutation ⇡ ⇤ (middle panel). If
⇤ is known to be positive, then the (unconstrained) ML estimator b ⇡ ML of ⇡ ⇤ matches the order statistics {µ (i) } n i=1 and {y (i) } n i=1 , where
In order to improve over the ML estimator using the prior knowledge of local shuffling, we impose the constraint that the alternative estimator b ⇡ does not pair any indices more than r = 3 apart. This estimator can be obtained as solution of the optimization problem max
where c ij = 0 if |i j| r and c ij = +1 otherwise. As in (4), the problem in (5) side is a linear assignment problem and hence computationally tractable, and corresponds to MAP estimation under the family of priors considered below. The corrected, i.e., repaired data
based on this approach are depicted in Figure 2 (R).
The priors discussed in the two examples of the previous subsection can be understood as specific instances of a more general family of prior distributions over P(n). Specifically, we consider the family of priors
where > 0 is the concentration parameter, and the matrix M (which is not required to have any specific properties) defines the mode(s) of the distribution argmax ⇧2P(n) h⇧, Mi, where h•, •i here represents the trace inner product on the space of matrices. In the same vein, the mode(s) of the distribution correspond to the set of matrices closest to M with respect to the same norm. Moreover, the distribution specified by ( 6) is of exponential family form with respect to the trace inner product (Wainwright and Jordan, 2008).
Linear Assignment Problems (LAPs). Linear assignment problems are a class of optimization problems for computing optimal one-to-one matchings of two sets of items (Burkard et al., 2009). LAPs are of the form
where C is a given cost matrix. By the Birkhoff-von Neumann theorem (Ziegler, 1995), the minimum over P(n) can be replaced by the minimum over DS(n), the set of n-by-n doubly stochastic matrices. Therefore, (7) reduces to a linear program in n 2 variables and n 2 + 2n linear constraints. This implies that computing a mode of ( 6) is tractable.
Specific examples. Below, we consider a few examples of interest that are special cases of ( 6).
(I) Hamming prior. Consider the choice M = I n . For any ⇧ 2 P(n), we then have h⇧,
denotes the Hamming distance on P(n). Since n does not depend on ⇡, Eq. ( 6) reduces to
which appeared in the first example of the preceding section, cf. ( 4), in which the goal was to take into account the underlying low rate of mismatches. The prior ( 8) is a specific Mallows' prior p(⇡) / exp( d(⇡, ⇡ 0 )) for a base permutation ⇡ 0 and a metric d on P(n) (Mallows, 1957).
(II) Local shuffling prior. As in the second example in §2.1, suppose we want to have the prior p place most of its mass on permutations that move indices within small windows, i.e., |⇡(i) i| tends to be small. This can be achieved by choosing the entries of M in (6) as M ij = (|i j|) for a non-decreasing function . The choice (u) = 0 if |u| r for a positive integer r and (u) = +1 otherwise yields the approach (5) that underlies the example in Fig. 2 above.
(III) Block prior. In record linkage, it is common that
) is block diagonal with known block composition. For example, suppose that gender, ethnicity, and age group are used as matching variables, and that these three categorical variables are free of errors. In this case, mismatches can only involve pairs (i, j) falling into the same block corresponding to a specific combination of the above variables. Such known block structure can be encoded via prior ( 6) by choosing M ij = 1 if (i, j) are not contained in the same block and M ij = 0 otherwise. This corresponds to a uniform prior for each block, i.e., p(
The prior for each block does not have to be uniform; e.g., a Hamming prior as in Example (I) above can be used instead. Moreover, the hard block constraint can be relaxed.
(IV) Lahiri-Larsen prior. In their seminal work on linear regression in the presence of mismatch errors, Lahiri and Larsen (2005) and Chambers (2009) 6). An example for Q is the so-called exchangeable linkage model (Chambers, 2009;Zhang and Tuoto, 2021)
In this case, the resulting prior is equivalent to the Hamming prior considered in Example (I). More complex priors are obtained depending on the structure of Q.
We now outline how the above prior can be integrated into generic shuffled data problems. The proposed Monte-Carlo EM (Wei and Tanner, 1990) framework builds upon the paper by Wu (1998) that has been rediscovered in the more recent work Abid and Zou (2018). The MC-EM scheme in Wu (1998) was further developed in Gutman et al. (2013) based on the concept of data augmentation (Tanner and Wong, 1987). None of Wu (1998); Abid and Zou (2018); Gutman et al. (2013) consider informative priors for ⇡.
Conditional & Integrated Likelihood. Suppose we are given data D = {(x i , y i )} n i=1 potentially contaminated by mismatch error. Let p(x j , y i ; ✓) be the likelihood (depending on a parameter ✓) for the pair (x j , y i ), (i, j) 2 [n] 2 . The likelihood for ✓ resulting from D conditional on a specific ⇡ 2 P(n) is given by
Conjugacy. It is worth noting that under (9), the posterior p(⇡|D, ✓) is a member of the family of distributions specified by p(⇡), i.e., the latter is a conjugate prior. This follows from the observation that
with M D,✓, = log(p(x j ,
The (conditional) likelihood ( 9) can be maximized with respect to both ✓ and ⇡ as, e.g., in Pananjady et al. (2018); Abid et al. (2017); Slawski and Ben-David (2019). Alternatively, ✓ is considered as the quantity of primary interest, which suggests the integrated likelihood
As seen in §2.1, maximizing the conditional likelihood is prone to overfitting, prompting a need for regularization. The use of the integrated likelihood mitigates that problem at best slightly, but not substantially (cf. supplement for details), hence regularization remains relevant.
MC-EM scheme. The Expectation-Maximization (EM) algorithm (Dempster et al., 1977) is an established heuristic for minimizing the negative log-likelihood `(✓) = log L(✓) corresponding to (11) via a sequence of surrogates { e `(t) (✓; ✓ (t) )} t 0 that are minimized successively:
where e `(t) (✓;
the so-called expected complete data negative loglikelihood. The surrogates { e `(t) (• ; ✓ (t) )} tend to be easier to minimize since they are linear combinations of standard likelihood terms as encountered for fixed and known ⇡. Surrogates are updated according to
Here, the main challenge of this scheme is the E-step, i.e, the calculation of the expectation on the right term in (12). For any pair (i, j), we have
Since the summation over P(n) is not computationally tractable, the expectation needs to be approximated, e.g., via Monte Carlo simulation. Since for the same reason, the posterior p(⇡|D, ✓ (t) ) is only accessible up to an unknown constant (cf. ( 10)), it is appropriate to resort to Markov Chain Monte Carlo (MCMC) (Gelman et al., 2013). The Metropolis-Hastings (MH) algorithm can be used to generate a Markov Chain {⇡ (k) } k 1 whose stationary distribution equals p(⇡|D, ✓ (t) ). This yields the approximation
where b denotes the length of the "burn-in" period, and m denotes the total length of the Markov chain. Substituting (13) into ( 12) then yields what is known as MC-EM scheme, cf. Algorithm 1. Conveniently, there is a proposal distribution for the MH algorithm that is easy to work with, known as Fisher-Yates sampling: it generates a new permutation from the current one by swapping the assignments of a pair of indices (cf. Algorithm 2).
Initialization. The choice of the initial iterate ✓ (0) can critically impact the quality of the solution that is returned by EM schemes given that the latter is a local strategy that finds a stationary point of a (in general) non-convex objective near the initial iterate. Several consistent initial estimators are known for regression setups depending on the structure of ⇡ (Lahiri and Larsen, 2005; Chambers and Diniz da Silva, 2020; Slawski et al., 2021;Peng et al., 2021), and those naturally lend themselves as initial iterate.
Careful initialization of the MH subroutine is important in order to ensure that p(⇡|D, ✓ (t) ) is explored well given that |P(n)| = n! while the number of MCMC iterations m is limited. Fortunately, under the prior (6), computing the mode argmax ⇡ p(⇡|D, ✓ (t) ) reduces to an LAP of the form (7) in virtue of (10). Initialization via the mode has the advantage that the Markov chain is started in a high density region. The hope is that the resulting iterates (generated according to a localized proposal distribution) will pick up most of the mass of p(⇡|D, ✓ (t) ) so that (13) will well approximate the underlying expectation.
Algorithm 1 Monte Carlo EM (MC-EM) algorithm
r(e ⇡, ⇡ (k) ) min n p(e ⇡|D,✓; )
Reduction under exponential family likelihood. For a variety of exponential family models, the expected complete data negative log-likelihood (12) involves n instead of n 2 terms. Specifically, (12) will be
Examples of interest are presented in the sequel.
(i) Least squares regression. In this case, we take log p(x, y; , 2 ) = (y x > ) 2 /(2 2 ), which corresponds to the negative likelihood of a linear regression model with i.i.d. zero-mean Gaussian errors with variance 2 . This yields the following expression for the expected complete data negative log-likelihood:
which is identical to a least squares objective with design matrix X and response vector E[⇧ > |D, ✓ (t) ]Y.
(ii) Generalized linear models. In this case, we have log p(x, y; ,
, where a, and c denote scale, cumulant, and partition function, respectively. Similar to above, one shows that 1 a( )
While the canonical link is assumed above, this is not necessary to achieve the reduction from n 2 to n terms.
where z = [x > y > ] > ; zz > has diagonal blocks xx > , yy > and off-diagonal blocks xy > , yx > . Thus tr(⌦zz > ) = tr(⌦ xx xx > ) + tr(⌦ yy yy > ) + 2tr(⌦ yx xy > ), where ⌦ xx , ⌦ yy etc. are the corresponding sub-matrices of ⌦. Hence,
where
Computational complexity of Algorithm 1. For exponential family models benefitting from the above reduction, updating ✓ in the M-step involves n terms, and is computationally equivalent to a standard estimation problem. Apart from the initialization of the Markov chain, the approximate E-step has complexity O(m), where m denotes the length of the Markov chain. Computing the acceptance probability, updating ⇡ (k) , and keeping track of b E[⇧|D, ✓ (t) ] > Y within Algorithm 2 can be done in time O(1) since the proposal distribution only changes ⇡ (k) at two positions. However, m is recommended to be of the order ⌦(n), heuristically justified by the fact that in the worst case a permutation is the product of n 1 transpositions.
Remarks. (i) Following Tanner and Wong (1987) and Gutman et al. ( 2013), the MC-EM approach can be converted into a fully Bayesian approach: given that MC-EM involves sampling from p(⇡|D, ✓), one can as well sample from p(✓|D, ⇡) in an alternating fashion, which yields a Gibbs sampler for the joint posterior p(✓, ⇡|D). (ii) The framework herein is not limited to permutations: we may be given D x = {x i } N i=1 and D y = {y i } n i=1 , N > n (w.l.o.g.), and then consider maps ⇡ : [n] ! [N ] represented by a matrix ⇧ 2 {0, 1} n⇥N with unit row sums. Priors of the form (6) given a mode M 2 R n⇥N as well as conditional and integrated likelihoods can be defined analogously to (9) and ( 11). (iii) We think of the sampling scheme as a template rather than an efficient approach; improving efficiency, e.g., along the lines of Zanella (2020); Grathwohl et al. (2021), is left for future work.
In this section, we present some analysis of the proposed prior from a regularization perspective and provide guidance on the choice of the hyperparameter . Data-driven selection of based on Empirical and Hierarchical Bayes approaches are detailed in the supplement.
Hamming prior. Our first results concerns the MAP estimator of ⇧ ⇤ under the Hamming prior (8). Specifically, we consider the linear regression setup
as considered in prior work on shuffled linear regression (Pananjady et al., 2018;Hsu et al., 2017). The theorem be-low considers the sparse setting in which the underlying ⇡ ⇤ satisfies the constraint d H (⇡ ⇤ , id) k for k "small enough" as made precise below. For simplicity, it is assumed that ⇤ and 2 ⇤ are known; various estimators for the regression parameter in this scenario have been proposed (Zhang and Li, 2020;Slawski et al., 2021;Peng et al., 2021).
with probability at least 1 2/n and 1 3/n, respectively, where µ = (µ i (x)) n i=1 and SNR = k ⇤ k 2 2 / 2 ⇤ . Theorem 3.1 implies that if is chosen larger than the threshold 0 , the MAP estimator b ⇧ will be 2k-sparse, which matches the sparsity of ⇧ ⇤ up to the factor 2. By the triangle inequality, d
⇧ and ⇧ ⇤ will be close in Hamming distance. Moreover, for values such that 3 0 C 0 for C > 3, we obtain
The dependence on the signal-to-noise ratio SNR is improved compared to the naive estimator b ⇧ 0 = I n , whose error scales as ⇤ p k log(en/k)SNR 1/2 ; for small SNR, one cannot hope for improvements over b ⇧ 0 in general. In light of the discussion in §2.1, the improvement over the maximum likelihood estimator b ⇧ ML whose error scales as ⇤ p n, is substantial as long as k is small relative to n. The next result yields a lower bound on ensuring a prescribed level of sparsity k based on the prior only.
Proposition 3.2. Suppose that ⇡ follows the Hamming prior (8) with parameter . Then for all 2 k < n
Proposition 3.2 asserts that the hyperparameter of the prior (8) should be chosen proportional to log(n k) ⇠ log n as n gets large in order to ensure that the prior places essentially no mass outside the Hamming ball {⇡ : d H (⇡, id) k}. The threshold ⇠ log n is sharp in the sense that if log(n k), the prior will place at least mass ⌦(1) = 1 4 !k/k! ⇠ 1 4e for not too small k outside that Hamming ball. The likelihood p(D|⇡) favors permutations with best fit to the given data, so that the posterior mass P ⇡|D ({⇡ : d H (⇡, id) k}) will be less than the prior mass. It is thus natural to consider ⇠ log n as initial point.
Local shuffling prior. The next statement addresses the scenario in Fig. 2 for Lipschitz functions. In particular, the level of penalty needed for the MAP solution to satisfy the condition max i |b ⇡(i) i| r is provided.
⇠ N (0, 1), where
, with probability at least 1 exp( ( p 2 1) 2 /2) 2/n. Under the same event,
In theory, the assertion of the above proposition can be achieved by setting = 1. Established solvers of LAPs require the entries of the cost matrix to be finite. In addition, solver accuracy can degrade with the magnitude of the entries (Bernhard, 2021).
In this section, we present the results of experiments conducted with synthetic and real data. In the supplement, we also show an example demonstrating the use of the data augmentation approach discussed at the end of §2.3 as an alternative to the Monte-Carlo EM scheme.
Synthetic data. We consider data generation according to the following three models ( 1 i n):
Linear Regression (LR):
i=1 and {✏ i } n i=1 are i.i.d random samples from the N (0, I d ) and N (0, 1) distributions, respectively. The regression parameter ⇤ is sampled uniformly from the sphere { 2 R d : k k 2 = 3}, and ⇤ 0 ⇠ N (0, 1). For MVN, we let µ ⇤ = 0 and ⌦
Finally, ⇡ ⇤ is a permutation selected uniformly at random from one of the following constraint sets:
(ii) r-Banded:
where P(B) denotes the set of block-structured permutations corresponding to B blocks of uniform size n/B, i.e., {1, . . . , n/B}, . . . , {(B 1)(n/B) + 1, . . . , n}. Note that in (i), k refers to the number total of mismatches, whereas in (iii) k refers to the number of mismatches per block. We fix n = 1, 000, d = 20, ⇤ = 1, ⇢ ⇤ = 0.8, p = q = 5, B = 50. The mismatch rates k/n and k • B/n in (15)(i) and (15)(iii), respectively, are varied between 0.2 and 0.5 in steps of 0.05, and the bandwidth r in (15)(ii) is varied between 3 and 10. For each setup and each value of k and r, 100 independent replications are performed. The following approaches are compared:
(I) naive. Standard maximum likelihood estimation as used for parameter estimation in the absence of mismatches, which corresponds to fixing ⇡ = id as the identity.
(II) oracle. The unknown permutation ⇡ ⇤ is considered as known, and standard maximum likelihood estimation is used for parameter estimation after fixing ⇡ = ⇡ ⇤ .
(III) robust [for setting k-Sparse only]. For setup LR, the regression parameter is estimated on the robustfit function in (MATLAB, 2019). For setup GLM, the regression parameter is estimated based on the robust GLM estimation method (Wang et al., 2020) that uses observationspecific dummy variables and penalization. For setup MVN, ⌦ 1 ⇤ is estimated according to the robustcov function in MATLAB which implements the minimum covariance determinant estimator (Rousseeuw and Driessen, 1999).
(IV) EM, EMH, EML, EMB. Algorithm 1 using uniform, Hamming, local shuffling, and block-Hamming prior, respectively, which reflect the constraint sets (i) to (iii) in ( 15). The EM iterations are initialized by setting ⇡ = id, and the number of EM iterations is limited to 400. The number of MCMC iterations per EM iteration is set to 8k, half of which are counted towards the "burn-in period". We note that a modified MH algorithm is used for EML (cf. supplement); for EMB, the MH scheme in Algorithm 2 is applied blockwise. For the Sparse and SparseBlock settings, the hyperparameter is chosen based on Proposition 3.2, which suggests / log(n). For the Banded setting, the choice = 1 was found to yield good performance. 2009) and Lahiri and Larsen (2005) with the choice Real data. We consider seven benchmark data sets for shuffled data problems. The data sets are preprocessed versions of their original counterparts (details on data processing can be found in the supplement). Even though the data sets themselves are real, the permutations that scramble the given matching pairs {(x i , y i )} n i=1 are synthetic; for each data set, we consider 100 independent random permutations depending on the underlying setting. We consider the same list of competitors and associated configurations as for the synthetic data experiments. Asterisked ground truth parameters here refer to oracle estimates based on knowledge of ⇡ ⇤ , and relative estimation error (REE) is defined accordingly in terms of those ground truth parameters.
Hamming & Block prior. As shown in Fig. 4, the proposed approach consistently improves over naive least squares once the fraction of mismatches exceeds 0.2, and yields substantial improvements as that fraction increases. The regularized EM approach based on the priors in EMB, EMH, EML noticeably reduces error induced by shuffling.
Local shuffling prior. As shown in Table 1, the EM approach with local shuffling prior achieves significant error reductions compared to the naive approach and the EM approach without regularization.
In this paper, we have proposed a framework for regularized estimation in shuffled data problems by means of an exponential family prior on the permutation group. The exponential family form is convenient for computational purposes yet sufficiently rich to incorporate common forms of prior knowledge. The proposed prior is not tailored to specific data analysis problems, but can be applied generically. The results in this paper confirm the importance of regularization in shuffled data problems given the inherent danger of overfitting already with little noise. While the approach covers various constraints that can be imposed on the underlying permutation, it is not exhaustive. For example, suppose we have information on the cycles of the permutation (numbers or lengths). Such information cannot be expressed in terms of index pairs, and hence requires a different paradigm. Kondor et al. (2007); Huang et al. (2009) use Fourier analysis on the permutation group (Diaconis, 1988) to facilitate learning of permutations, and it is an interesting direction of future research to study how that approach can be leveraged for the type of shuffled data problems considered in the present paper.
A Detailed derivations of the expected complete data negative log-likelihood for selected models
In this section, we provide detailed steps for deriving the expected complete data negative likelihood for least squares regression and precision matrix estimation for multivariate Normal data.
Recall that the expected complete data negative likelihood is given by
In this case, we take log p(x, y; , 2 ) = (y x > ) 2 /(2 2 ), and we have that
(ii) Precision matrix estimation for multivariate normal data. In this case, log p(x, y; ⌦) = log det ⌦ + tr(⌦zz > ), where z = [x > y > ] > ; zz > has diagonal blocks xx > , yy > and off-diagonal blocks
and analyzing each of the terms accordingly.
C Proof of Theorem 3.1
In light of relation ( 10), straightforward manipulations and omission of terms not depending on ⇧ show that the MAP estimator b ⇧ is the minimizer of the optimization problem
Since b ⇧ minimizes (S.2), the following basic inequality holds true:
Decomposing Y = µ + ⇠ with ⇠ = ⇤ ⇧ ⇤ ✏ and re-arranging terms in the above inequality yields that
where we have substituted d H (⇧ ⇤ , I n ) = k. By the Cauchy-Schwarz inequality, h⇧ ⇤ µ, b ⇧µi k⇧ ⇤ µk 2 2 , which implies that the first term in the previous inequality is non-negative. This in turn yields that h⇠, ( b
In the sequel, we will derive a probabilistic lower bound on the first term on the left hand side.
For any integer 1 s n, consider the event where for any integer 1 ` n, B 0 (`) here denotes the set of all unit vectors having at most `non-zero entries.
Denote by w(B(`)) = E g⇠N (0,In) [sup u2B(`) hu, gi] the corresponding Gaussian width (Vershynin, 2018, §7.5). Choosing t = w(B(m s )) + c 1 p 2 log n in (S.5) for c 1 1, standard tail bounds for the suprema of Gaussian processes (Boucheron et al., 2013, Theorem 5.8) yield
Using that w(B 0 (m s )) 4 p m s log(en/m s ) (e.g., Plan Vershynin, 2013, Lemma 2.3) and the fact that s log(en/s) log n for any n s 1, we have with c 1 = p 2
Combining this with the definition of v yields that
, and note that
Using the same argument as above, we choose t = ⇤ ⌧ s {w(B 0 (m s )) + c 2 p 2 log n} with c 2 = p 2. Putting together the pieces as above, we obtain
Now let 0 = 72 p SNR log(en/k) 72 p SNR log(en/m s ) and define the event
Observe that conditional on E s \ G s , the earlier inequality (S.4) implies that (recalling that m s = s + k)
Combination of the left and right hand sides and re-arranging terms implies the inequality 0 s k s + k Now for any s 2k, the right hand side is lower bounded by (1/3) . This in turn yields a contradiction whenever is chosen such that > 3 0 . In order to conclude that d H ( b ⇧, I n ) 2k with the stated probability in that case, it remains to provide a corresponding lower bound on the probability of the event S n s=1 (E s \ G s ), i.e., at least one of the events {E s \ G s } n s=1 occurs. Since the events inside the union are disjoint, we obtain that
(1 2/n 2 ) P(E s ) 1 2/n.
In order to prove the second part of Theorem 3.1, we first invoke the following basic inequality equivalent to (S.3)
Expanding the squares and re-arranging yields conditional on
where in the second inequality, we have used that if > 3 0 , conditional on on
After elementary manipulations, we obtain the inequality
with probability at least 1 2/n 1/n = 1 3/n, where the term sup u2B0(3k) is controlled similarly to (S.6).
The probability mass function of ( 8) is given by (Fligner and Verducci, 1986):
In the sequel, let us write {D(⇡) = d} as a shortcut for the event {d H (⇡, id) = d}. We then have
where !d denotes the number of derangements of {1, . . . , d}, i.e., the number of permutations ⌧ of d objects such that ⌧ (j) 6 = j for all 1 j d. Straightforward manipulations yield
For x 0 and integer m 1, define the (upper) incomplete Gamma function and its "normalized" counterpart by
where (m) = (m, 0) = (m 1)! denotes the Gamma function. It can be shown that (Abramowitz and Stegun, 1964
, we obtain the following for the right hand side of (S.8):
e (n + 1, exp( ) 1) .
At this point, we consider the upper bound on the probability of interest as stated in the proposition. We have
(1 + ) log n, which concludes the proof of the upper bound.
Regarding the lower bound, observe that in view of relation (S.9), the ratio of normalized incomplete Gamma functions can be expressed via the ratio of CDFs of two independent Poisson random variables, that is
where X 1 and X 2 are two independent Poisson random variables with parameters exp( ) and exp( ) 1, respectively. Setting = log(n k) yields that the right hand side is a function of the form c 1 (n, k) that is lower and upper bounded by 1/4 and 1, respectively, as
Similar to Eq. (S.3) in the proof of Proposition 3.1, we have the basic inequality
In the sequel, we will show that under the stated conditions, the left hand side must exceed the right hand side unless b ⇧ ij = 0 for all (i, j) such that |i j| > r. Expanding Y = ⇧ ⇤ µ + ⇤ ✏ and re-arranging terms yields the inequality
where we have used that k⇧ ⇤ µk 2 2 h b ⇧µ, ⇧ ⇤ µi 0. For the second term on the right hand side, the triangle inequality yields that for all indices i that are summed over, we have |b ⇡(i) ⇡ ⇤ (i)| 2r. Using the Cauchy-Schwarz inequality in combination with the Lipschitz property of the underlying function, we obtain that
By standard concentration results (e.g., Wainwright, 2019, § 2.3), the event E 1 = {k✏k 2 p 2n} holds with probability at least 1 exp(( p 2 1) 2 /2). We now turn to the first term on the right hand side of (S.11). We have the upper bound X ⇤ ensures that the left hand side exceeds the right hand side, which is a contradiction, and hence it must hold that |b ⇡(i) i| r, 1 i n.
Observe that conditional on the event {|b ⇡(i) i| r, 1 i n}, the basic inequality (S.10) reduces to hY, b ⇧µi hY, ⇧ ⇤ µi () kY b ⇧µk 2 2 kY ⇧ ⇤ µk 2 2 Substituting Y = ⇧ ⇤ µ + ✏ in the inequality of the right hand side and expanding squares, we obtain that conditional on {|b ⇡(i) i| r, 1 i n} and E 1
with the same arguments as used for (S.11) and (S.12). Dividing both sides by n yields the assertion. if r(e ⇡, ⇡ (k) ) > u: ⇡ (k+1) e ⇡.
else: ⇡ (k+1) ⇡ (k) .
end for return b E[⇡|D, ✓] as in ( 13) with m replaced by m invalid-mcmc-steps.
In this section, we provide references of each data set and regression model used in the real data analysis. A summary of each data set is shown in Table S.1 below. (Slawski et al., 2021) 2011 2 LR hamming 2k El Nino data(END) (Slawski et al., 2021) 93935 5 LR block 1.5k* CPS wages(CPS) (Slawski et al., 2021) 534 11 LR hamming 2k Bike sharing data(BSD) (Wang et al., 2020) 731 16 GLM block 1.5k* Flight Ticket Prices(FTP) (Slawski et al., 2020) 335 30 6 MVN hamming 2k Supply Chain Management(SCM) (Slawski et al., 2020) where Inv-2 (⌫, a 2 ) refers to the scaled inverse 2 -distribution with scale parameter a > 0 and ⌫ degrees of freedom (cf. §14 in Gelman et al. (2013) for more details on Bayesian inference for linear regression models).
For this illustration, we use m = 100, where each sequence {⇡ (j) } is generated by uniform thinning of Markov chains of length 4, 000 generated by Algorithm 2. The number of samples ( (k) , 2 (k) ) obtained via the above scheme is taken as 1,000. The sampling procedure is initialized from step (II) with the identity permutation. We compare both the unregularized case with the uniform prior for ⇡ as well as the regularized case with the Hamming prior (8) ( = log n in view of Proposition 3.2). Figure S.1 confirms that the proposed approach achieves visible improvements over the unregularized approach which suffers from serious amplification bias affecting the slope parameter 1 and serious underestimation of the error variance, as predicted by the brief analysis accompanying the first introductory example in §2.1.
In this section, we outline how the hyperparameter of the proposed prior on ⇡ can be selected based on Empirical and Hierarchical Bayes approaches. For simplicity, these approaches are presented for the linear regression model ( 14) in an exemplary fashion.
Consider the following hierarchical model specification: where ( ) denotes the terms in the normalization constant in the prior p(⇡| ) that depend on .