Linear Regression with Sparsely Permuted Data

arXiv:1710.06030v2 [math.ST] 15 Nov 2017

Martin Slawski
Department of Statistics
George Mason University
Fairfax, VA 22030, USA
mslawsk3@gmu.edu

Emanuel Ben-David
Center for Statistical Research Methodology
U.S. Census Bureau
Suitland, MD 20746, USA
emanuel.ben.david@census.gov

Abstract
In regression analysis of multivariate data, it is tacitly assumed that response and
predictor variables in each observed response-predictor pair correspond to the same
entity or unit. In this paper, we consider the situation of “permuted data” in which this
basic correspondence has been lost. Several recent papers have considered this situation
without further assumptions on the underlying permutation. In applications, the latter
is often to known to have additional structure that can be leveraged. Specifically, we
herein consider the common scenario of ”sparsely permuted data” in which only a
small fraction of the data is affected by a mismatch between response and predictors.
However, an adverse effect already observed for sparsely permuted data is that the
least squares estimator as well as other estimators not accounting for such partial
mismatch are inconsistent. One approach studied in detail herein is to treat permuted
data as outliers which motivates the use of robust regression formulations to estimate
the regression parameter. The resulting estimate can subsequently be used to recover
the permutation. A notable benefit of the proposed approach is its computational
simplicity given the general lack of procedures for the above problem that are both
statistically sound and computationally appealing.

1

Introduction

A largely unquestioned assumption in regression analysis with response-predictor pairs
{(yi , xi )}ni=1 is that each yi corresponds to the same statistical unit as xi . In this paper,
we consider the situation where the identifiers of the predictors (or equivalently those of
the responses) are subject to an unknown permutation so that correspondence between yi
and xi may not be be taken for granted. We refer to this situation as “permuted data”
respectively “sparsely permuted data” when the permutation only affects a small fraction
of all response-predictor pairs. Restoring the original correspondence between responses
and predictors may not be achievable in practice for both computational and statistical
reasons, but may also not be required for consistent estimation of regression parameters.
Conversely, if the primary interest concerns recovering the permutation itself, which is the
case in entity resolution, the regression model can facilitate that task.

1.1

Background and Motivation

Large organizations that own or have access to multiple data sources regularly rely on data
integration for conducting large-scale scientific projects. Available datasets are gathered
or produced at different points of time and independently of one another. The main
reason for combining datasets is that no single dataset contains all variables of interest
that pertain the research questions. Data collected from a single source are often limited
and do not have all the variables needed for statistical analysis. Limited budget, time
and resources prevent each particular agency from collecting a comprehensive dataset. In
sparse_permutation_arxiv-v1.tex, version: November 17, 2017

2

●
●
●
●
●
●●
● ●
●
●●
●
●
●
●
●
●
●
● ●●● ●
●
● ●
●
●
● ●
● ●●●●
● ●
●●
● ●●●
●●
●●
●●
●●
●
●
●
●● ●
●
●● ● ●
●
●●
●
●
●
●
●
●
●
●
●● ● ●
●
●
●●
●
●
● ● ●
● ●
● ●● ●
●
●
●
●
●
●
●
●
●
●
●
●● ●
●●●
●
●
●
●●●
●
●
●
●
● ●●
●● ●
● ●
●
● ●●●●
●●
●
●● ●
●
●
●
●
●
●
●
●
●●
●●●
●
●
●●
●●●
●
●
●
●
●
●
● ●●●● ●●
●
●
●●
●● ●●
● ●●●
●
●
●
● ●
●
●●
●
●
●
● ● ● ●●
●
●
●
●
●
●●
●
●
●
●●
●
●●●
●●
●
●
●
●
● ●
●●
●
●
●
●
● ●
●
●
●●
●
●●● ●
● ●
●
● ●
●
●
●
●● ● ● ●
●●
●
●
● ●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
● ●● ● ●●●
●● ● ●
●●
●
●
●
●
● ●
●●
●●● ● ● ●●
●●● ● ●
●
●
●
●
●● ●
●
●
●●● ●
●
●
●
●
●
●
●
●
● ●●
● ●●
● ●
●
●
● ● ●
●● ●
●
●
●
●●
●
●
●
●
●
●● ● ●
●●
●
●
●
●
●●●
●●●
●●
●
●
●
●
●
●●
●
●●
● ●
●
● ●●
●
●
●
●●
●
●
●● ●
●
●
●
●
●
●
●
●
●●●●● ●
●
●
● ●●
●
●
●
●●
●
●
●
●●
●
●
●
●
● ●
●
● ●
●
●
●
●
●
● ●
●
●
●
●●
●
●
● ● ●
●
●●
●
●
●
●
●
●

3

●

●

2

●

0

●
●
●

−4

0

2

● ●
● ● ●
●●
●
● ● ●●
●
●
●
● ● ●●●●
●● ●
● ● ● ●●
●●
●

●
●

●●
●
●
●
●●

●
●
●

●
● ●
●
● ●
● ●
●
●
● ●●
●
●
●●

●

●
●

●
●●
●
●
●●●●● ●
●
●●
●●●●●
●
●
●● ●● ●
●
●●
●● ●
●●
●●
●●
●
● ● ● ●● ● ●● ●● ● ●
●● ● ●
●
●
● ●
●
●
● ●
●
●●●●●
●
● ● ●●● ●
●
●
●●
●
●
●
●
●
●
● ● ● ●●● ●
● ●● ●
●●
● ●● ●
●
●● ●
●●
●
●
●
●
●
●●
●
●
●
●●
●●●●●
●
●
●
●
●●
●
●
●
●●●●●
●●
●
●●
●●
●●
●
●●●●● ●●
●●
●
●
●
●
●
●●
●●●●●●●●●●
● ●●● ●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●●
●
●
●
●
●● ●
●
●
●●● ●
●●
●
●●
●
●
●
●
●●
●
●● ● ●●
●
●
●●
●
●
●
●●
●
●● ●
●●
● ●
●●
●●
●
●● ●
●
● ●●
●●●
●
● ● ●●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●● ●● ●
●
●●
● ●
●
● ●
●
●●●
● ●● ●
●●
●
●●
● ●●
●●
● ●● ●●
●●
● ●●
●●
●
●●
●
●
● ● ●●
● ● ●●
●
●● ●● ●
●
●
●●●● ●
●
●
●●●
●
●
●
● ●●●●
●
●●
●● ●
● ●●
●
●●●
●
●
● ● ●
●
● ●
●●
●
● ●● ● ●
●●
●
●
● ● ●
●
●
●
●●
●
● ● ●
●
●

−3
−2

●

●

●
●
●

●

●
●

●
●●

●
●

●

●

●

●

−1

Residuals

1

●

●

−2

−2

0

●

−4

Permuted response variables in the merged data

y (original) y (merged)
●

original data
merged data
−4

−3

−2

Response variables as observed in the original data

Fitted

−1

0

1

..
.
−0.96
−3.45
0.28
−3.06
−0.33
−2.98
−0.53
−1.35
..
.

..
.
−3.45
−0.96
−1.86
−0.33
−3.02
−0.53
−2.98
−3.73
..
.

Figure 1: Illustration of the effect of mismatchings in the response variable on linear
regression with merged data. For this analysis, a subset of the El Niño data set (cf. §3.2)
was artificially broken into two parts which were subsequently merged based on a nonunique identifier. Left: Original response vs. response observed after data merging. Middle:
Residuals vs. fitted values in linear regression when using the original (dots) and merged
data set (triangles), respectively. The circles and triangles colored in red correspond to the
excerpt of the response variable depicted on the right hand side.
most situations, however, relevant datasets from other sources are already available. For
example, if a demographic survey does not have all the relevant variables for a particular
research project, other existing surveys or administrative data can be used to include those
missing variables.
Record linkage, or entity resolution, is an essential task in data integration. The task is
to identify which records in different datasets belong to the same entity. In this context, the
term ”entity” is to be understood in a broad sense, and may refer to customers, tax payers
or patients, etc. Record linkage has a long history of uses in large enterprises, government
agencies and health care institutions. A business register consisting of names, addresses,
and other identifying information such as total financial receipts might be constructed
from tax and employment data bases; a survey of retail establishments or agricultural
establishments might combine results from an area frame and a list frame. In this case,
units from the area frame would need to be identified in the list frame [52].
Generally, due to the lack of unique entity identifiers in data files, record matching is
based on methods of probabilistic record linkage [24] that score similarity between common
quasi-identifiers. For example, when the data files contain health information of individuals,
the quasi-identifiers can be first and last names, addresses, dates of birth. Spelling errors,
typographical variation and missing values in records inhibit exact matching and lead to
matching errors: mismatches and missed-matches. Even low matching error rates can
lead to selection bias and pervasive outliers, particularly when the record linkage process
includes or excludes particular types of entities and their attributes [10]. Selection bias
and outliers contaminate subsequent statistical analysis. To reduce bias, it is of interest to
develop statistical methods that can alleviate the adverse effects of matching errors.
In this paper, we focus on the problem of linear regression with sparsely permuted data
as it arises frequently when linking survey records to external data. The canonical example
is regression estimation of a study variable y on auxiliary variables x that reside in another
data source such as master data or administrative records. The matching error, determined
by the fraction of mismatches and missed-matches, is expected to be only a small fraction
of the sample size n [52]. The paper [47] considers two such situations. To model the
2




(y1 , x1 )


..


.
(yn , xn )




(y1 , xϕ(1) )


..


.
(yn , xϕ(n) )




(y1 , xϕ(1) )


..


.
(yn , xϕ(n) )

(y, Π∗ X)
(Π∗> y, X)
data with true correspondence

(y, X)
observed data

Figure 2: Regression with lost correspondence between response and predictors.
energy economy properly, an economist might need company-specific data on the fuel and
feedstocks used by companies that are only available from Agency A and corresponding
microdata on the goods produced for companies that is only available from Agency B. To
model the health of individuals in society, a demographer or health science policy worker
might need individual-specific information on those receiving social benefits from Agencies
B1, B2, and B3, corresponding income information from Agency I, and information on
health services from Agencies H1 and H2.
Regression with permuted data arises in other applications as well. The papers [3, 39]
provide extensive lists including multi-target tracking in radar systems [43] and pose and
correspondence estimation in computer vision [19].

1.2

Problem statement

To setup the problem more concretely, suppose the data (y1 , x1 ), . . . , (yn , xn ) is obtained
from matching two data files A and B, where for each i, yi ∈ R resides in file A and
xi ∈ Rd resides in file B. Since the process of record linkage is not error-free, some xi
may have been paired up with a non-correspondent yi . If the number of such mismatches
is known to be less than or equal to k ≤ n, then there is an unknown permutation ϕ on
[n] := {1, . . . , n} so that: ϕ moves at most k of the indices, and (y1 , xϕ(1) ), . . . , (yn , xϕ(n) )
are independent realizations from the classical linear regression model
y = x> β ∗ + , where  ∼ N (0, σ 2 ) and x ⊥⊥ .

(1)

Let Π∗ and Π∗> bethe matrix representations of ϕ and its inverse, respectively, and let
X > = x1 · · · xd . Eq. (1) readily implies that
y = Π∗ Xβ ∗ + , where y = (y1 , . . . , yn )> and  = (1 , . . . , n )>
⇔ Π∗> y = Xβ ∗ + Π∗> .

(2)

For simplicity, since the distribution of Π∗>  and  is the same, we do not distinguish
between Π∗>  and  in the sequel and write  in place of Π∗> .
The main question to be studied here is how the parameters Π∗ and β ∗ can be accurately
and efficiently estimated from the partially mismatched pairs (yi , xi ). Efficiency here refers
to computational complexity of the procedure that outputs the desired estimates. Let us
consider the number of mismatches which can be expressed as dH (Π∗ , In ) := |{i : Π∗ii = 0}|,
the Hamming distance between Π∗ and the identity matrix In . We can naturally estimate
Π∗ and β ∗ by a least-squares approach as:
b = argmin kΠXβ − yk2
b β)
(Π,
2
Π,β

(3)

subject to dH (Π, In ) ≤ k,
3

where Π and β run over Pn , the set of all permutation matrices in Rn×n , and Rd , respectively. A permutation respectively its matrix representation Π is said to be k-sparse if
dH (Π, In ) ≤ k. We note that [39] have shown that unless d = 1 the optimization problem
(3) is NP-hard for k = n. The same is true for any k that is defined as a fraction of n.
Notation. We here gather some notation frequently used in the present paper. For a
positive integer m, Im denotes the m × m identity matrix, and Sm−1 denotes the unit
sphere in Rm . We write |S| for the cardinality of a set S. The complement of S with
respect to a given base set depending on the context is denoted by S c . For a matrix
A, kAk2 = σmax (A) denotes its spectral norm respectively maximum singular value, and
range(A) denotes the column space of A. For an index set I ⊆ [m] := {1, . . . , m} and
v ∈ Rm , vI denotes the subvector of v corresponding to I. We write a ∨ b = max{a, b}.
Positive constants are denoted by C, c, c1 etc. We make use of the usual Big-O notation
in terms of O, o, Ω and Θ.

1.3

Prior work

Linear regression with linked data files is a common scenario in which traditional methods
of statistical data analyses are error-prone. Early work in [37] recognizes the adverse
effect of matching error on the regression analysis of linked datasets. They show that
as a consequence of matching error, the ordinary least squares estimator βbols is generally
not an unbiased estimator of β ∗ . In the paper [37], the process of matching is regarded
as random. Letting zi denote the response that is actually in correspondence to xi , and
letting yi denote the response that has been matched to xi , the model in [37] assumes that
for i = 1, . . . , n,
(
zi with probability qii
yi =
(4)
zj with probability qij , i 6= j.
Assuming that these probabilities can be estimated [34, 35, 37, 46, 47] focus on constructing
an estimator that is unbiased or less biased than βbols . In fact, in [35] the matrix Q = (qij ) is
used to construct an unbiased estimator of β ∗ . The papers [29, 30, 49] propose to estimate
β ∗ using a Bayesian procedure. The main shortcoming of these approaches are that they
rely on the assumption that the doubly stochastic matrix Q is known or can be accurately
estimated from the data. In addition, while achieving reduction in bias, the proposed
estimators may still be inconsistent and may have large mean squared errors.
The classical papers [20, 21, 22] and the later note [53] consider the situation of permuted data under the term “broken sample”. A broken sample is a sample of (yi , xi )-pairs
that up to some permutation of the {xi } (or equivalently the {yi }) are generated from their
joint distribution. In other words, each component of (yi , xi ) is observed separately (as if it
were generated from its marginal distribution), with possibly different orders for each component. Assuming that the {(yi , xi )} are generated from a bi-variate normal distribution
up to a permutation, these papers discuss recovering the permutation or estimating the correlation parameter of the underlying bivariate normal distribution. In [5, 15], the authors
discuss whether that parameter can be consistently estimated from a broken sample.
In recent years, computer science and engineering have witnessed a surge of interest
in regression analysis of permuted data resulting from problems in multi-target tracking,
statistical seriation [25], and unlabeled compressed sensing [50], to mention just a few.
The papers [39, 40, 50] are particularly important as they provide rigorous results on
fundamental questions associated with the problem. We shall refer to some of these results
in the subsequent sections in more detail. While the paper [17] is not concerned with a

4

regression setting, it provides a detailed analysis of the problem of finding correct matches
between two sets of objects in the presence of noise, which bears some relation to the
problem discussed in §2.4 below.

1.4

Summary of Contributions

The analysis in this paper concerns estimators of Π∗ and β ∗ given model (1) under the
i.i.d.

additional assumption xi ∼ N (0, Id ), i ∈ [n], as also assumed in the recent work [39]. It
i.i.d.

is not hard to extend the results to the case where xi ∼ N (0, Σ), i ∈ [n], where Σ is a
d × d symmetric positive definite matrix.
We first provide a negative result concerning the least squares approach (3) if k = n,
i.e., if no constraint is imposed on the permutation Π. It is shown that optimizing over
all possible choices of Π leads to overfitting in the sense that if β ∗ = 0, one still has
b 2 = Ω(1) with high probability. This result complements another one of a similar spirit
kβk
2
in the recent paper [3] who show that for d = 1, the least squares estimator (3) converges
almost surely to a limit different from β ∗ . Our result also aligns well with a minimax lower
bound in [40].
Altogether, these negative results provide additional justification to consider the regime
of sparse permutations with k  n. We demonstrate a bound on the estimation error
kβb − β ∗ k2 with βb as in (3). Specifically, the bound implies that the error vanishes as
both d/n and k log(n/k)/n → 0. In view of the fact that the optimization problem (3)
is not computationally tractable, we consider a convex relaxation that takes the form
of a robust regression estimation problem as it has been considered before in different
contexts [26, 38, 48], with the permuted observations here being in correspondence to gross
errors. The robust regression formulation can be reduced to a specific sparse regression
problem in an underdetermined setting with n − d samples and n parameters, one for each
observation in a given sample that is affected by partial mismatches between the {xi }ni=1
and the {yi }ni=1 . Our analysis of the robust regression problem applies generically beyond
the specific setting of sparsely permuted data. Prior works [26, 38] have considered a more
general version of the problem in which β ∗ is assumed to be sparse as well, thereby being
able to cover the regime n < d, but it is not clear whether the results in [26, 38] can be
specialized to match those of the present paper for the traditional n > d case.
While the robust regression formulation gives rise to an error bound for estimating
∗
β that is comparable to that of the computationally hard formulation (3), it does not
immediately yield an estimator for the permutation Π∗ . We address this issue by adopting
a two-stage approach that uses an accurate estimator θb of β ∗ to match the fitted values
n
bn
{x>
i θ}i=1 to the responses {yi }i=1 . This reduces to simple sorting operations, thereby
avoiding the computational challenges associated with problem (3). We show that our
approach recovers the underlying permutation under qualitatively the same condition as
in [39] which is considerably more stringent in terms of required signal-to-noise ratio than
what is required for accurate estimation of the regression coefficients β ∗ . We complement
our result with a comparable lower bound on the signal-to-noise ratio that is required for
permutation recovery even if β ∗ itself is known.
As already pointed out in [39], the Hamming ball constraint in (3) does not substantially
change the fundamental statistical limits of permutation recovery. However, that constraint
does help in that it gives rise to a computationally efficient estimator of Π∗ , whereas the
statistical achievability result in [39] is for the computationally hard estimator (3).

5

Π∗



Theorem 3

β∗

Π∗ Xβ ∗ − Xβ ∗



Theorem 1
Theorem 2

Lemma 2 (B)
Lemma 6 (D)



S ∗ = {i : Π∗ii 6= 1}
Corollary 1

Figure 3: Schematic overview on targets in the setting of linear regression with sparsely
permuted data and pointers to corresponding results in the paper. The directions of 
indicate that if the object on the left hand side is known or can be accurately estimated,
it immediately yields or simplifies estimation of the object on the right hand side.

2

Main results

This section is structured as follows. We start by showing that the estimator βb in (3) lacks
statistical consistency if k = n. This sets the stage for considering the scenario k  n. We
first focus on recovery of β ∗ , treating Π∗ as a nuisance parameter. We then present and
analyze a post-processing step for recovering Π∗ given an accurate estimate of β ∗ .
In the sequel, unless stated otherwise, we assume Gaussian design for the {xi }ni=1 :
(G) :

i.i.d.

xi ∼ N (0, Id ), i ∈ [n].
i.i.d.

Our results generalize to the case xi ∼ N (0, Σ), i ∈ [n], for a symmetric positive definite
matrix Σ ∈ Rd×d in a straightforward manner by defining a new regression parameter
Σ1/2 β ∗ . An estimator of that parameter (as discussed in the sequel) and an estimator of Σ
can then be combined to form an estimator of β ∗ . Note that estimation of Σ is not affected
by the presence of an unknown permutation.

2.1

Least squares estimation of (Π∗ , β ∗ ) without additional constraints:
a negative result

Let us consider problem (3) for k = n, i.e., no further constraints are imposed on the
b It turns out that in this case, we cannot hope for βb to be a good estimator
b β).
solution (Π,
∗
of β . As can be seen from the proof of the following proposition, complete freedom in
choosing Π to fit the data results into over-adaptation to noise even if β ∗ is low-dimensional
– in fact, the phenomenon already manifests itself for d = 1.
Proposition 1. Let β ∗ = 0 and Π∗ = In , i.e., y =  and consider the estimator (3) with
k = n. Then there exist constants c, C > 0 such that with probability at least 1−C exp(−cn)
b 2≥
kβk
2

n
σ2
.
2n + d 32π 2

b 2 will be bounded away from zero by a constant
In other words, in a pure noise setting, kβk
(assuming d = O(n)). A result of a similar flavor is shown in [3]. For d = 1, they show
that the estimator βb converges almost surely to a limit different from β ∗ as n → ∞, and
they derive an explicit expression for that limit. Our result here is more of a qualitative
nature, but it is non-asymptotic and not limited to d = 1.
Both the result in [3] as well as Proposition 1 raise the question whether there exist
alternative estimators that do significantly better in a regime where Π∗ can be an arbitrary
6

element of Pn . Theorem 1 in [40] indicates that the answer is negative: they show that for
e one has
e β),
any estimator (Π,
h
i
e βe − Π∗ Xβ ∗ k2 ≥ cσ 2 .
sup
E n1 kΠX
2
Π∗ ∈Pn , β ∗ ∈Rd

While the above lower bound concerns estimation of Π∗ Xβ ∗ rather than β ∗ , accurate insample prediction (or synonymously denoising), i.e., recovery of Π∗ Xβ ∗ , is typically easier
than recovery of β ∗ . Along this direction, another negative result is shown in [32]. They
demonstrate that a lower bound on the signal-to-noise ratio (SNR)
SNR = kβ ∗ k22 /σ 2

(5)

of the order Ω(d/ log log n) is required for any estimator of β ∗ to achieve a non-trivial
expected relative estimation error1 . As shown in [39], the condition SNR > nc for c large
b of (3) to recover Π∗ , in which case kβb − β ∗ k2 scales
enough is sufficient for the solution Π
2
as O(d/n) as in the usual regression setting in the absence of an unknown permutation.

2.2

Least squares estimation of β ∗ if Π∗ is k-sparse, k  n

In summary, the previous section points to the fact that we cannot hope for accurate
estimation of β ∗ without additional assumptions on Π∗ and/or the SNR of the problem.
As motivated in the introduction, we henceforth turn our attention to the case of Π∗ being
k-sparse, with k “significantly smaller” than n. The allowable range of k is addressed in
our analysis presented below.
We start by fixing notation. For 0 ≤ k ≤ n, let us introduce the shorthand
Pn,k = {Π ∈ Pn : dH (Π, In ) ≤ k}
for the constraint set of Π in problem (3). Moreover, for a compact and symmetric2 set
K ⊂ Rn , its Gaussian width is defined by


w(K) = Eg∼N (0,In ) sup | hg, xi | .
(6)
x∈K

While originating in geometric analysis, the Gaussian width is a measure of complexity that
has been increasingly adopted in the analysis of high-dimensional linear inverse problems
[13, 16, 42] in connection with Gordon’s ”Escape Through a Mesh Theorem” [27], which
is the key component in the proof of Theorem 1 below as well. In our setting, we use the
Gaussian width (6) in conjunction with the set
o
[ n
T =
range(Π − Π∗> ) ∩ Sn−1 ,
(7)
Π∈Pn,k

where we recall that Π∗> ∈ Pn,k is the inverse of Π∗ . Let k·k0 denote the `0 -”norm”, i.e.,
the number of non-zero entries of a vector. A simple observation is that for any Π ∈ Pn,k
and any v ∈ Rn , it holds that
kΠv − Π∗> vk0 = k(Π − In )v − (Π∗> − In )vk0 ≤ k(Π − In )vk0 + k(Π∗> − In )vk0 ≤ 2k.
1
2

Non-trivial here refers to a relative estimation error of lower order than that of the estimator θb ≡ 0.
A set S ⊂ Rn is called symmetric if x ∈ S implies that −x ∈ S.

7

As a result, T ⊆ B0 (2k; n) ∩ Sn−1 , where for 0 ≤ r ≤ n, the set B0 (r; n) = {v ∈ Rn :
kvk0 ≤ r} denotes the `0 -”ball” of radius r in Rn . By well-known results (cf. [41], Lemma
2.3), we hence have that
p
w(T ) ≤ w(B0 (2k; n) ∩ Sn−1 ) ≤ 3.5 2k log(en/2k).
(8)
Moreover, it is not hard to show that w(T ) ≥ w(B0 (n − k; k) ∩ Sn−1 ), i.e., there is a lower
bound on w(T ) of the same order. After these preparations, we are in position to state
the following result.
Theorem 1. Consider
optimization
h
i problem (3) and for positive integers m, denote νm =
m
1/2
Eg∼N (0,Im ) [kgk2 ] ∈ (m+1)1/2 , m
and let ε ∈ (0, 1/2) be a number such that
νn−d −

2
ε
νn ≥
w(T ).
1−ε
1−ε

(9)

If furthermore n > 9 ∨ 4d, then with probability ≥ 1 − 72 exp(−cε n) − 2 exp − 12 (d ∨ log n) ,
!
r
√ −2
σ
5(d
∨
log
n)
2(1
+
2)ε
{w(T
)
∨
log
n}
∗
q
√
kβb − β k2 ≤
+
.
(10)
n
n
4d∨log n
1−
n

We start our interpretation of Theorem 1 by inspecting condition (9) which imposes a limit
on how large d and k can be in relation to n. Roughly speaking, the requirements are n =
Ω(d) and n = Ω(w2 (T )). In light of (8), the latter condition becomes n = Ω(k log(n/k)).
Specifically, let us fix ε = 1/4 and d = αn for α ∈ (0, 14 ), then (9) essentially evaluates as
2

8
√
n≥
w2 (T ).
3( 1 − α − 1/3)
The error bound (10) consists of two parts. The first part equals the error one would
have if the permutation Π∗ were known in advance and is thus inevitable. The second
term is a bound on the excess error incurred for not knowing Π∗ . That term is well
controlled as long as long as the fraction of permuted observations is small relative to n. A
b 2 . Under
crucial intermediate step in the proof of Theorem 1 is a bound on kΠ∗> y − Πyk
b identifies
an additional condition on the SNR, we may deduce from that bound that Π
∗
∗
∗
the “support” S = {i : Πii 6= 1} of Π . We may then re-fit with the corresponding
observations left out, to achieve a smaller error in estimating β ∗ .
Corollary 1. For any δ ∈ (0, 1), under the conditions of Theorem 1, if it additionally
holds that
√
2(1 + 2)2 ε−4 k 2 {w(T ) ∨ log n}2
SNR >
,
n − k ≥ 9 ∨ 4d,
·
δ2
n
the following events hold with probability at least 1 − δ − 2 exp(− 21 (d ∨ log n)):
r
σ
5(d ∨ log n)
bb
∗
∗
b
b
q
S := {i : Πii 6= 1} = S ,
kβ − β k2 ≤
.
n−k
n
1 − 4d∨log
n−k
b
b
where βb denotes the ordinary least squares estimator based on data {(yi , xi ) : i ∈ [n] \ S}.
Under the conditions of the corollary, as long as n − k = Ω(n), the error rate in estimating
β ∗ is of the same order as if Π∗ were known in advance; the second term in (10) gets
eliminated. It is important to note that support recovery (i.e., {Sb = S ∗ }) is a weaker
b
result than permutation recovery (i.e., {Π∗ = Π}).
As discussed in more detail in §2.4
below, the latter requires a considerably more stringent condition on the SNR than what
is required in Corollary 1.
8

2.3

Convex relaxation

While the approach of the previous section has appealing statistical properties, its computational hardness asks for computationally efficient alternatives with similar guarantees.
As long as Π∗ is treated as a nuisance parameter, we may eliminate it in the following way.
Introducing f ∗ = Π∗ Xβ ∗ − Xβ ∗ = (Π∗ − In )Xβ ∗ , model (2) can be re-expressed as
y = Xβ ∗ + f ∗ + .
This prompts the optimization problem
ky − Xβ − f k22
[
subject to f ∈
range(Π − In ).
min

β∈Rd ,f ∈Rn

Π∈Pn,k

A first relaxation is given by
min

β∈Rd ,f ∈Rn

ky − Xβ − f k22
(11)

subject to kf k0 ≤ k.
P
We note in passing that one could additionally impose the constraint ni=1 fi = 1. However, it turns out that its addition does not yield significant statistical benefits, and it is
thus omitted. Relaxation (11) is still not convex, but following the standard approach of
replacing the `0 -norm by the `1 -norm, we end up with the convex optimization problem
min

β∈Rd ,f ∈Rn

ky − Xβ − f k22
(12)

subject to kf k1 ≤ b.
Since it tends to be difficult to choose b appropriately in practice, it is more convenient to
√
work with the Lagrangian form of (12). After re-parameterizing e = f / n3
min

β∈Rd ,e∈Rn

√
1
ky − Xβ − nek22 + λkek1 ,
n

λ > 0.

(13)

Formulation (13) and variants thereof have been used in robust regression and signal recovery with gross corruptions (e.g., [26, 36, 38, 48]). In fact, (13) is equivalent to employing
Huber’s loss [33] instead of squared loss, cf. [48] and the references therein. The connection
to robust regression comes naturally as observations with mismatch between x and y are
likely to induce severe errors in model fitting beyond the usual noise, and hence take the
role of outliers. Indeed, this reasoning could have been used to motivate (13) right away
instead of via the sequence of relaxations stated above. Formulation (13) is related to least
absolute deviation regression
min ky − Xβk1
(14)
β∈Rd

in that (14) is obtained from (13) in the limit λ → 0 and the additional constraint y =
√
Xβ + ne. In that sense, (14) can be seen as the counterpart to (13) in a noiseless setup.
Problem (14) has been analyzed under assumption (G) in the landmark papers of [14, 45]
on the classical error correcting problem in coding theory.
Theorem 2 below provides an upper bound on the `2 -error of the estimator of β ∗
resulting as the optimal β for (13).
3

This reparameterization is merely for technical reasons

9

e ee) be a minimizer of (13) with λ = 4(1 + M )σ
Theorem 2. Let (β,
M > 0. There exist constants c1 , c2 , ε > 0 so that if
k ≤ c1

p
2 log(n)/n for some

n−d
,
log(n/k)


2
the following holds with probability ≥ 1−2 exp(−c2 (n−d))−2n−M −2 exp − 21 (d ∨ log n) :
!
r
r
5(d
∨
log
n)
2k
log
n
σ
n
q
.
(15)
kβe − β ∗ k2 ≤
+ 48(1 + M )
ε−1
n
n
−
d
n
4d∨log n
1−
n
Comparing the upper bounds (10) and (15) of the original problem and its relaxation,
respectively, we observe a close agreement given that w(T ) = Θ(k log(n/k)). Apart from
the slight change of the order in the second term, only constants differ. Similarly to
Corollary 1, we may consider a two-stage procedure in order to further improve upon (15):
given ee and a suitable threshold t, we let Set = {i : |e
ei | ≥ t} and obtain a plain least squares
e
e With k assumed to be known, we may
fit based on data {(yi , xi ) : i ∈ [n] \ Set } yielding β.
take t as the k-th largest largest entry of ee in absolute magnitude. The formal analysis is
similar to Corollary 1 and is hence omitted.

2.4

Estimation of Π∗ given an estimator of β ∗

Having discussed estimation of β ∗ , we now come back to the problem of estimating Π∗ . As
indicated above, the latter problem turns out to be significantly more challenging than the
former problem. In the sequel, we study the question of how the availability of an accurate
estimate of β ∗ can be leveraged to construct an estimator of Π∗ that is computationally
feasible including the case d > 1. As mentioned earlier, except for d = 1, the optimization
problem in Eq. (3) is in general NP-hard. When d = 1, that is there is only one predictor
b is determined as a minimizer of the optimization problem
in the regression model, Π
max {hΠX, yi , hΠX, −yi}

subject to dH (Π, In ) ≤ k.

Π∈Pn

(16)

For k = n, maximizing each term inside the curly brackets is a specific instance of a linear
assignment problem [12], a class of problems that can be solved in polynomial time despite
their combinatorial nature. In fact, it easy to see that
max hΠX, yi =

Π∈Pn

n
X

X(i) y(i) ,

(17)

i=1

where X(i) and y(i) denote the i-th order statistic of X and y, respectively, i ∈ [n]. Hence,
for k = n, problem (16) reduces to two vanilla sorting operations. At this point, it is not
well understood yet whether the Hamming constraint dH (Π, In ) ≤ k for general k causes
problem (16) to be NP-hard again. So far, we are not aware of any computationally efficient
algorithm for general k.
Note that when β ∗ is known, computing the least squares estimator of Π∗ reduces to
solving precisely one of the two optimization problem encountered in (16) for d = 1:
max hΠXβ ∗ , yi

Π∈Pn

subject to dH (Π, In ) ≤ k.

(18)

At this point, a natural idea is to replace β ∗ by an accurate and computationally feasible
estimator like the one discussed in §2.3. From a computational point of view, this already
10

constitutes a simplification as (18) reduces to an integer linear program; while problems
in this class are still NP-hard in general, problem (18) can be considered as much more
benign than the original problem (3) which belongs to the class of quadratic assignment
problems notorious for their computational hardness. Due to recent advances in integer
programming that have meanwhile been taken advantage of for a series of other statistical
problems [7, 8], it turns out that problem (18) is practically feasible at least for n in the
order of a few thousands. For large n, we instead recommend estimating the support of Π∗
and then solve the unconstrained problem (17) restricted to observations in the estimated
support. Formally, denote by Se an estimator of the support of Π∗ and let ySe and XSe
e
be the sub-vector of y and the row submatrix of X corresponding to observations in S,
∗
e
respectively. We then estimate Π by Π defined by
!
ee e
Π
S,S
e =
ΠX
X,
(19)
In−|S|
e
where


e eX =
Π
S


XSe
,
XSec

e e e = argmax ΠX e, y e .
and Π
S,S
S
S
Π∈P|S|
e

We now turn to the statistical limits of permutation recovery, including a lower bound on
the SNR in the idealized case with known β ∗ . For simplicity, the theorem stated below is
for the case k = n, but it generalizes to k < n conditional on having {Se = S ∗ } in that all
expressions in n below would get replaced by k.
b
Theorem 3. For b ∈ Rd , consider Π(b)
= argmaxΠ∈Pn hΠXb, yi.
(a) Let δ, ∆, η > 0, and let θb be an estimator s.t. the event {kX θb − Xβ ∗ k∞ ≤ σ∆} holds
with probability at least 1 − η. If


kβ ∗ k22
n(n − 1) 2
n2 (n − 1)2
b 6= Π∗ ) ≤ 2δ + η.
b θ)
∆
+
2
log
>
, then P(Π(
σ2
4δ 2 π
δ
(b) Suppose that n ≥ 5. If
kβ ∗ k22
b ∗ ) 6= Π∗ ) ≥ .35.
< Cn2 , for some constant C > 0, then P(Π(β
σ2
e 4 ) is a
If β ∗ is known, specializing part (a) to the case ∆ = 0 asserts that SNR = Ω(n
sufficient condition for permutation recovery. This is to be compared to a result in [39]
b in (3) recovers Π∗ if SNR = Ω(nc ), where the constant c > 0
stating that the estimator Π
is not specified. Next, let us now turn to the case where β ∗ is substituted by an estimator
b Using standard concentration arguments, one shows that
θ.
p
√
kX θb − Xβ ∗ k∞ ≤ max kxi k2 kθb − β ∗ k2 ≤ C( d + log n)kθb − β ∗ k2
1≤i≤n

holds with high probability. Specifically, considering the estimator βe in §2.3, Theorem 2
√
then yields kX βe − Xβ ∗ k∞ ≤ Cσ n for k small enough, with high probability. Inserting
√
e 5 ). We stress again
∆ = C n into part (a) then results into the requirement SNR = Ω(n
e
b
e
that as opposed to Π, the estimator Π(β) is computationally appealing as it is obtained
from a quadratic program (13) and subsequent sorting.
Finally, part (b) provides evidence that the condition SNR = Ω(n2 ) is necessary for
permutation recovery. While the result in (b) concerns a specific estimator, there does not
seem to be much indication for the existence of a substantially better estimator. In [39], it
is shown that SNR = Ω(n8/9 ) is necessary for any estimator.
11

3

Numerical results

3.1

Simulated data

Below, we present the results of two synthetic data experiments that are meant to serve
as illustration of the developments in the previous section. In the first set of experiments, we generate n = 200 observations from model (2) under (G) with d = 10,
σ ∈ {.01, .02, .05, .1, .2, .5, 1}, and k/n ∈ {.01, .02, .05, .1, .15, .2, . . . , .5}, where the support
of Π∗ is selected uniformly at random. The parameter β ∗ is generated from the uniform
distribution on Sd−1 . For each pair (σ, k/n), we conduct 100 independent replications.
We compare the following estimators:
(i) The naive estimator (ordinary least squares estimator) of β ∗ that ignores the fact
that a fraction of the data is mismatched; this corresponds to the choice k = 0 in (3).
p
e ee) of §2.3 with the choice λ = 0.2σ log(n)/n.
(ii) The estimator (β,
(iii) A thresholded version of the estimator in (ii) that discards those observations corresponds to the largest k elements of {|e
ei |, i ∈ [n]}, and performs a least squares re-fit
using the remaining n − k observations.
e ee) is computed in CVX [28].
The estimator (β,
The estimator βb in §2.2 is compared to βe in a second set of simulations only for the
case d = 1. In that case, computation of βb reduces to (16) where each of the two inner
optimization problems can be cast as an integer linear program with n2 variables and 2n+1
linear constraints. It turns out that the general purpose routine cplexbilp in IBM CPLEX
[1] is able to solve such problems in only a few seconds for n = 200. Again, such reduction
is limited to the case d = 1.

-2
-4
-6
-8
0

0.1

0.2

0.3

0.4

(iii)
log2(estimation error)

(ii)
log2(estimation error)

log2(estimation error)

(i)
-2
-4
-6
-8
0

0.1

0.2

k/n

0.3

k/n

0.4

-2
-4
-6
-8
0

0.1

0.2

0.3

0.4

k/n

Figure 4: Comparison of the average `2 -errors (log2 -scale) in recovering β ∗ for the approaches (i) to (iii) described above. Each curve corresponds to a different value of the
noise level σ. The black dashed lines correspond to the oracle (ordinary least squares with
knowledge of Π∗ ).

Inspection of Figure 4 shows that the approach (ii) improves dramatically over the naive
estimator (i) as long as the SNR is not too small (for σ = 1 and σ = 0.5, there is no longer a
visible improvement). The results look promising in that the tolerable fraction of permuted
observations can be as much as 0.5 as long as the noise level is small; ”tolerable” here
refers to the fact that the estimation error increases gently with the fraction of permuted
observations as opposed to (i) with a sharp increase in error as k/n moves away from zero.
12

Approach (iii) appears to yield further improvements over (ii) for large SNR. For small
SNR, (ii) and (iii) are not distinguishable. This observation aligns well with Corollary 1.
Figure 5 reveals that for d = 1, the estimator βb performs substantially better than the

-2
-3
-4
-5
-6
-7
-8
-9
-10
-11
0

βe

log2(estimation error)

log2(estimation error)

βb

0.1

0.2

0.3

0.4

-2
-3
-4
-5
-6
-7
-8
-9
-10
-11
0

k/n

0.1

0.2

0.3

0.4

k/n

Figure 5: Comparison of the average `2 -errors (log2 -scale) in recovering β ∗ for the estimators βb and βe for d = 1. Each curve corresponds to a different value of the noise level σ.
The black dashed lines correspond to the oracle (least squares with knowledge of Π∗ ).

e We observe that βb is on par with the oracle for σ < 0.2,
convex formulation yielding β.
being hardly effected by an increase in k/n. That starts changing for σ = 0.2, while for
σ = 0.5, the adverse effect of an increase in k/n becomes clearly visible. At this point, it
is not clear whether similar improvements of βb over βe would be observed for d > 1. The
analysis in §2 only provides (comparable) upper bounds for both approaches, and lower
bounds would be needed to answer that question. The observation made from Figure 5
hence points towards an important open direction.

3.2

Real data

The El Niño data set in the UCI Machine Learning Repository [2] contains oceanographic
and surface meteorological readings taken from a series of buoys positioned throughout
the equatorial Pacific. The data set consists of n = 93, 935 records with the following attributes: buoy identifier, date, location (latitude and longitude), zonal and meridional
wind speeds (zon and mer), relative humidity (humidity), air temperature (air.temp),
sea surface temperature and subsurface temperatures down to a depth of 500 meters
(s.s.temp). The following linear regression model is considered:
air.temp = β0 + βzon · zon + βmer · mer+
+ βhumidity · humidity + βs.s.temp · s.s.temp + .

(20)

The results of least squares regression indicate an excellent fit of this model with an Rsquared of about 0.9. With the goal to mimic the situation of data merging based on nonunique identifiers, the data set is divided into two data sets A and B, with A containing the
response variable and B containing the predictors. The variables latitude and longitude
are maintained in both data sets, and used as quasi-identifiers for record linkage of A and
B with the R package fastLink [23]. Since (latitude, longitude)-pairs are generally
associated with multiple observations, the linkage process is not free of mismatches. The
merged data set in turn follows the format depicted in Figure 2 with a permutation affecting
13

interc

zon

mer

humidity

s.s.temp

RMSE

`2 -dist.

∅-hold-out error

βboracle
βbols

5.15

−.056

−.031

−.022

.844

.509

0

.260 (4.6·10−3 )

6.72

−.037

−.045

−.017

.774

.771

1.57

.276 (4.9·10−3 )

βe

5.74

−.044

−.037

−.016

.806

.773

.59

.267 (4.8·10−3 )

Table 1: Results of the real data analysis. Here, βboracle denotes the least squares estimator
based on the original data set, βbols denotes the naive least squares estimator, and βe denotes
the estimator from formulation (13). The columns labeled zon, . . ., s.s.temp contain the
regression parameter estimates for the respective variable; interc refers to the intercept.
The three rightmost columns contain the root mean square errors (RMSEs) on the original
e the `2 -distance to βboracle , and the average
(βboracle ) respectively the merged data set (βbols , β),
hold-out errors (standard errors in parentheses), respectively.
the correspondence between responses and predictors. As shown in Table 1, least squares
regression with the thus merged data set leads to an increase of the residual sum of squares
by 27%. In addition, the estimates of the regression parameters of the linear model (20)
change noticeably. As an alternative, we consider the approach (13) studied in detail
√
herein, where the parameter λ is chosen as λ = 2cb
σ / n with c = 1.345 and σ
b being a
robust estimator of the residual standard error. More specifically, we use the function
rlm (short for “robust linear model fitting”) in the R package MASS [44] with its default
arguments which is in exact correspondence to the above choice of λ, as follows from the
connection between linear regression with the Huber loss and formulation (13) (cf., e.g.,
Proposition 3.1 in [48]). The figures in Table 1 indicate that this approach provides some
remedy relative to the naive use of the least squares estimator based on the merged data
set. The estimates of the regression parameters from (13) reduce the gap to those obtained
with the original data set by about a factor of one half in terms of `2 -distance. Moreover,
approach (13) also yields a better fit in terms of the mean squared prediction error on a
collection of twenty hold-out sets obtained from leaving out different random subsets of
10% of the observations; in each of those twenty runs, the hold-in data set is split and
subsequently re-linked in the same manner as described above. The reduction in excess
risk (relative to the least squares estimator) achieved by (13) in comparison to the naive
least squares solution is again a factor of roughly one half.

4

Conclusion

Regression problems in which the correspondence between responses and predictors has
been lost can be traced back to classical work in statistics and naturally arises in the context of record linkage. Despite its long history, the problem has seen revived interest lately
and has been analyzed through the lens of modern non-asymptotic statistical theory. An
additional layer of complexity to be dealt with is the computational hardness of (quadratic)
assignment problems. In the present paper, both aspects are taken into account under the
assumption of sparsely permuted data. One central insight of our work is that the availability of an accurate estimator of β ∗ may be helpful in circumventing the computational
barrier associated with permutation recovery. On the other hand, preliminary numerical
results for d = 1 indicate a gap in statistical performance between a presumably hard
formulation capturing precisely the notion of a sparse permutation and a computationally convenient relaxation that introduces a loss of information. This observation deserves
further investigation.
14

Acknowledgments
The second author would like to thank William Winkler for his guidance and Anindya Roy
for helpful discussions related to this paper.

References
[1] IBM ILOG CPLEX Optimization
marketplace/ibm-ilog-cplex.

Studio.

http://www.ibm.com/us-en/

[2] UCI Machine Learning Repository. https://archive.ics.uci.edu/ml/datasets.
html.
[3] A. Abid, A. Poon, and J. Zou, Linear Regression with Shuffled Labels.
arXiv:1705.01342, 2017.
[4] R. Adler and J. Taylor, Random Fields and Geometry, Springer, 2007.
[5] Z. Bai and T. Hsing, The broken sample problem, Probability Theory and Related
Fields, 131 (2005), pp. 528–552.
[6] R. Baraniuk, M. Davenport, R. DeVore, and M. Wakin, A Simple Proof of
the Restricted Isometry Property for Random Matrices, Constructive Approximation,
28 (2006), pp. 253–263.
[7] D. Bertsimas, A. King, and R. Mazumder, Best subset selection via a modern
optimization lens, The Annals of Statistics, 44 (2016), pp. 813–852.
[8] D. Bertsimas and R. Mazumder, Least quantile regression via modern optimization, The Annals of Statistics, 42 (2014), pp. 2494–2525.
[9] P. Bickel, Y. Ritov, and A. Tsybakov, Simultaneous analysis of Lasso and
Dantzig selector, The Annals of Statistics, 37 (2009), pp. 1705–1732.
[10] M. Bohensky, Methodological Developments in Data Linkage, John Wiley & Sons,
Ltd, 2015, ch. Bias in data linkage studies.
[11] S. Boucheron, G. Lugosi, and P. Massart, Concentration Inequalities: A
Nonasymptotic Theory of Independence, Oxford University Press, 2013.
[12] R. Burkard, M. Dell’Amico, and S. Martello, Assignment Problems: Revised
Reprint, SIAM, 2009.
[13] T. Cai, T. Liang, and A. Rakhlin, Geometric inference for general highdimensional linear inverse problems, The Annals of Statistics, 44 (2016), pp. 1536–
1563.
[14] E. Candes and T. Tao, Decoding by Linear Programming, IEEE Transactions on
Information Theory, 51 (2005), pp. 4203–4215.
[15] H.-P. Chan and W.-L. Loh, A file linkage problem of DeGroot and Goel revisited,
Statistica Sinica, 11 (2001), pp. 1031–1045.
[16] V. Chandrasekaran, B. Recht, P. Parrilo, and A. Willsky, The convex
geometry of linear inverse problems, Foundations of Computational Mathematics, 12
(2012), pp. 805–849.
15

[17] O. Collier and A. Dalalyan, Minimax Rates in Permutation Estimation for Feature Matching, Journal of Machine Learning Research, 17 (2016), pp. 1–31.
[18] S. DasGupta, An elementary proof of a theorem of Johnson and Lindenstrauss, Random Structures and Algorithms, 22 (2003), pp. 60–65.
[19] P. David, D. Dementhon, R. Duraiswami, and H. Samet, Softposit: Simultaneous pose and correspondence determination, Int. J. Comput. Vision, 59 (2004),
pp. 259–284.
[20] M. DeGroot, P. Feder, and P. Goel, Matchmaking, The Annals of Mathematical
Statistics, 42 (1971), pp. 578–593.
[21] M. DeGroot and P. Goel, The Matching Problem for Multivariate Normal Data,
Sankhya, Series B, 38 (1976), pp. 14–29.
[22]

, Estimation of the correlation coefficient from a broken random sample, The
Annals of Statistics, 8 (1980), pp. 264–278.

[23] E. Enamorado, B. Eifield, and K. Imai, Fast Probabilistic Record Linkage with
Missing Data. R-package, Version 0.2.0.
[24] I. P. Fellegi and A. B. Sunter, A theory for record linkage, Journal of the American Statistical Association, 64 (1969), pp. 1183–1210.
[25] N. Flammarion, C. Mao, and P. Rigollet, Optimal rates of statistical seriation.
arXiv:1607.02435, 2016.
[26] R. Foygel and L. Mackey, Corrupted Sensing: Novel Guarantees for Separating
Structured Signals, IEEE Transactions on Information Theory, 60 (2014), pp. 1223–
1247.
[27] Y. Gordon, On Milman’s inequality and random subspaces which escape through a
mesh in Rn , Springer Berlin Heidelberg, Berlin, Heidelberg, 1988, pp. 84–106.
[28] M. Grant and S. Boyd, CVX: Matlab software for disciplined convex programming,
version 2.1. http://cvxr.com/cvx, Mar. 2014.
[29] R. Gutman, C. Afendulis, and A. Zaslavsky, A Bayesian Procedure for File
Linking to Analyze End-of-Life Medical Costs, Journal of the American Statistical
Association, 108 (2013), pp. 34–47.
[30] M. H. P. Hof and A. H. Zwinderman, Methods for analyzing data from probabilistic linkage strategies based on partially identifying variables, Statistics in Medicine, 31
(2012), pp. 4231–4242.
[31] D. Hsu, S. Kakade, and T. Zhang, A tail inequality for quadratic forms of subGaussian random vectors, Electronic Communications in Probability, 52 (2012), pp. 1–
6.
[32] D. Hsu, K. Shi, and X. Sun, Linear regression without correspondence.
arXiv:1705.07048, 2017.
[33] P. Huber, Robust Estimation of a Location Parameter, The Annals of Mathematical
Statistics, 53 (1964), pp. 73–101.
16

[34] G. Kim and R. Chambers, Regression analysis under probabilistic multi-linkage,
Statistica Neerlandica, 66 (2012), pp. 64–79.
[35] P. Lahiri and M. D. Larsen, Regression analysis with linked data, Journal of the
American Statistical Association, 100 (2005), pp. 222–230.
[36] J. N. Laska, M. A. Davenport, and R. G. Baraniuk, Exact Signal Recovery
from Sparsely Corrupted Measurements through the Pursuit of Justice, in Asilomar
Conference on Signals, Systems and Computers, 2009, pp. 1556–1560.
[37] J. Neter, S. Maynes, and R. Ramanathan, The effect of mismatching on the
measurement of response error, Journal of the American Statistical Association, 60
(1965), pp. 1005–1027.
[38] N. Nguyen and T. Tran, Robust Lasso with Missing and Grossly Corrupted Observations, IEEE Transactions on Information Theory, 59 (2013), pp. 2036–2058.
[39] A. Pananjady, M. Wainwright, and T. Cortade, Linear Regression with an Unknown Permutation: Statistical and Computational Limits, in Communication, Control and Computing (Allerton), 2016, pp. 417–424.
[40]

, Denoising Linear Models with Permuted Data. arXiv:1704.07461, 2017.

[41] Y. Plan and R. Vershynin, Robust 1-bit compressed sensing and sparse logistic regression: a convex programming approach, IEEE Transactions on Information Theory,
59 (2013), pp. 482–494.
[42]

, The generalized Lasso with non-linear observations, IEEE Transactions on Information Theory, 62 (2016), pp. 1528–1537.

[43] A. B. Poore and S. Gadaleta, Some assignment problems arising from multiple
target tracking, Math. Comput. Model., 43 (2006), pp. 1074–1091.
[44] B. Ripley, B. Venables, D. Bates, K. Hornik, A. Gebhard, and D. Firth,
MASS: Support Functions and Datasets for Venables and Ripley’s MASS. R-package
version 7.3.-47.
[45] M. Rudelson and R. Vershynin, Geometric approach to error correcting codes
and reconstruction of signals, International Mathematical Research Notices, 64 (2005),
pp. 4019–4041.
[46] F. Scheuren and W. Winkler, Regression analysis of data files that are computer
matched I, Survey Methodology, 19 (1993), pp. 39–58.
[47]

, Regression analysis of data files that are computer matched II, Survey Methodology, 23 (1997), pp. 157–165.

[48] Y. She and A. Owen, Outlier Detection Using Nonconvex Penalized Regression,
Journal of the American Statistical Association, 106 (2012), pp. 626–639.
[49] A. Tancredi and B. Liseo, Regression analysis with linked data: problems and
possible solutions, Statistica, 75 (2015), pp. 19–35.
[50] J. Unnikrishnan, S. Haghighatshoar, and M. Vetterli, Unlabeled sensing with
random linear measurements. arXiv:1512.00115, 2015.
17

[51] R. Vershynin, In: Compressed Sensing: Theory and Applications, Cambridge University Press, 2012, ch. ’Introduction to the non-asymptotic analysis of random matrices’.
[52] W. E. Winkler, Matching and record linkage, Wiley Interdisciplinary Reviews: Computational Statistics, 6 (2014), pp. 313–325.
[53] Y. N. Wu, A note on broken sample problem, tech. rep., Department of Statistics,
University of Michigan, 1998.
[54] S. Zhou, Restricted Eigenvalue Conditions on Subgaussian Random Matrices.
arXiv:0912.4045.

A

Proof of Proposition 1

Consider the least squares problem (3) with y = :
mink − ΠXβk22 .
Π,β

Omitting the term not depending on β, an equivalent problem is given by
min −2 hΠ, Xβi + kXβk22 .
Π,β

Re-parameterizing Xβ = u · r, with r ≥ 0 for u ∈ range(X) ∩ Sn−1 , we obtain the optimization problem
min
−2 hΠ, ui r + r2 .
Π,u∈range(X)∩Sn−1 ,r≥0

Observe that for any u and Π, the optimal value of r is found as
rb(u) = max{hΠ, ui , 0}.
b u
In particular, for a minimizer (Π,
b),
b u
b u
rb := rb(Π,
b) = hΠ,
bi
b 2 , and the optimal value of the objective is given by −b
so that rb = kX βk
r2 . We hence have
rb =

max

u∈range(X)∩Sn−1 ,Π

hΠ, ui

≥ max hΠ, Xj /kXj k2 i
1≤j≤d,Π

≥ maxhΠ, X1 i
Π

1
.
kX1 k2

We now bound
maxhΠ, X1 i ≥ hΠ± , X1 i,
Π

where the permutation Π± is defined as follows. Define
n,+ = {i : i > 0},

n,− = n − n,+ ,

nX1 ,+ = {i : Xi1 > 0},

nX1 ,− = n − nX1 ,+ .

Let further
n+ = min{n,+ , nX1 ,+ },

n− = min{n,− , nX1 ,− },
18

n± = n − n+ − n− .

(21)

Denote by
+
i+
1 < . . . < inX,+

the indices of the positive entries of X1 ,

j1+ < . . . < jn+,+

the indices of the positive entries of ,

−
i−
1 < . . . < inX,−
j1− < . . . < jn−,−

the indices of the negative entries of X1 ,
the indices of the negative entries of .

Furthermore, let
±
i±
1 < . . . < in±

j1± < . . . < jn±±
be the indices in
−
+
−
{1, . . . , n} \ ({i+
1 , . . . , in+ } ∪ {i1 , . . . , in− })

respectively {1, . . . , n} \ ({j1+ , . . . , jn++ } ∪ {j1− , . . . , jn−− })
in increasing order. We then construct Π± such that
hΠ± , X1 i =

n+
X

X1i+ j + +
k

n−
X

X1i− j − +
k

k

k

X1ik ± jk ± .

k=1

k=1

k=1

n±
X

All terms inside the first two sums are i.i.d. from the same distribution as |g||h|, where
g and h are independent N (0, 1)-random variables, and accordingly the terms inside the
last summand are i.i.d. from the same distribution as −|g||h|. Conditional on n± , hence
hΠ± , X1 i follows the same distribution as
n−n±

X

n
X

|gi ||hi | −

i=1

|gi ||hi |,

i=n−n± +1

where {gi }ni=1 and {hi }ni=1 are two independent sequences of n i.i.d. random variables from
the N (0, 1)-distribution. The expectation of the above expression (conditional on n± ) is
given by
2
(n − 2n± ) E[|g||h|] = (n − 2n± ) E[|g|]2 = (n − 2n± ) .
π
From Proposition 5.16 in [51], we have the following concentration inequality:


n−n±
n
X
X
1
(|gi ||hi | − E[|g||h|]) −
(|gi ||hi | − E[|g||h|]) ≥ (n − 2n± ) n±  ≤ exp(−nc).
P
π
i=1

i=n−n± +1

Hence, conditional on n± , with probability at least 1 − exp(−nc),
n−n±

X
i=1

|gi ||hi | −

n
X

|gi ||hi | ≥

i=n−n± +1

1
(n − 2n± ).
π

Consider the definition of n+ and n− Observe that n,+ , nX1 ,+ , n,− and nX1 ,− are each
binomial random variables with n trials and probability of success 1/2. By repeated application of Hoeffding’s inequality, we hence have
P(n+ ≤ 3n/8) ≤ 2 exp(−n/32),

P(n− ≤ 3n/8) ≤ 2 exp(−n/32).
19

Conditional on the events {n+ > 3n/8} and {n− > 3n/8}, we have n± ≤ 1/4 and thus
with probability at least 1 − C exp(−nc)
n−n±

X

|gi ||hi | −

i=1

n
X
i=n−n± +1

|gi ||hi | ≥

n
.
2π

n
Accordingly, the same applies to the event {hΠ± , X1 i > 2π
}. Lemma 15 in Appendix H
√
yields kX1 k2 ≤ 2 n with probability at least 1 − exp(−n/2). Putting together the pieces,
it follows that with probability at least 1 − C exp(−nc)
√
√
n
b 2 = rb ≥ n ⇒ kβk
b 2≥
kX βk
,
4π
σmax (X)4π

where σmax (X) denotes the
singular value of X. Finally, Lemma 16 in Appendix
√
√ maximum
H yields that σmax (X) ≤ 2n+ d with probability at least 1−exp(−cn), which concludes
the proof.

B

Proof of Theorem 1

b
b and β.
We start with a basic, yet important observation that allows us to decouple Π
⊥
Let PX and PX denote the orthogonal projection on range(X) respectively its orthogonal
complement.
Lemma 1. Consider optimization problem (3). Then, if n ≥ d, with probability one
(
)

> X −1 X > Π
b >y
X
>
⊥
2
>
⊥
2
b ∈ argminkP Πyk ,
b ∈ argminkP Πyk . (22)
Π
βb ∈
, Π
X
2
X
2
n
n
Π∈Pn,k
Π∈Pn,k
Proof. We have
min kΠXβ − yk22 = min kXβ − Πyk22
Π,β
Π,β
n
o
2
2
= min k P⊥
X Πyk2 + kXβ − PX Πyk2 .
Π,β

The second part can always be made equal to zero by choosing β such that Xβ = PX Πy.
It hence suffices to minimize the first part w.r.t. Π and to back-substitute the result into
the second part. Under model (G), X is non-singular with probability one as long as n ≥ d
which yields the second expression in (22).
b > y.
In the next lemma, we bound the `2 -distance between Π∗> y and Π
Lemma 2.
 Under the
conditions of Theorem 1, the following holds with probability at least
2
ψ(ε,d/n)νn
7
− exp(−{w(T ) ∨ log n}), where ψ(ε, d/n) is defined in (24) below.
1 − 2 exp −
8
√
2(1 + 2)σε−2 {w(T ) ∨ log n}
>
∗>
b
√
kΠ y − Π yk2 ≤
.
n
b > and the fact that Π∗> is a feasible solution,
Proof. In view of Lemma 1, by definition of Π
we have
⊥ ∗>
2
b> 2
kP⊥
X Π yk2 ≤ kPX Π yk2
∗>
⊥ ∗>
2
⊥ ∗>
2
b>
⇒ kP⊥
X (Π − Π )y + PX Π yk2 ≤ kPX Π yk2
∗>
2
⊥ ∗>
∗>
b>
b>
⇒ kP⊥
X (Π − Π )yk2 ≤ 2|hPX Π y, (Π − Π )yi|.

20

Since
∗>
⊥
∗
⊥
P⊥
X Π y = PX (Xβ + ) = PX ,

we have the implication that
∗>
2
⊥
∗>
b>
b>
kP⊥
X (Π − Π )yk2 ≤ 2|hPX , (Π − Π )yi|.
∗>
b>
Dividing both sides by kP⊥
X (Π − Π )yk2 , we obtain that
⊥
∗>
b > − Π∗> )yk22 kP⊥
b 2
b
b>
k(Π
X ∆k2 ≤ 2|hPX , ∆i|k(Π − Π )yk2 ,

⇒

b :=
∆

⊥
b > − Π∗> )yk2 kP⊥
b 2
b
k(Π
X ∆k2 ≤ 2|hPX , ∆i|.

b > − Π∗> )y
(Π
b > − Π∗> )yk2
k(Π
(23)

i) Lower bounding the left hand side of (23).
b ∈ T as defined in (7). Moreover, note that under (G), the random
Observe that ∆
subspace range(X) ⊆ Rn follows the uniform distribution on the Grassmannian G(n, d).
Hence, applying Lemma 13 with K = T , V = range(X) and thus m = n and l = n − d,
we have for any ε ∈ (0, 1/2) in compliance with (9)


⊥
P(dist(T , range(X)) > ε) = P min kPX ∆k2 > ε
∆∈T

 !
7
1 (1 − ε)νn−d − ενn − w(T ) 2
≥ 1 − exp −
2
2
3 + ε + ενn /νn−d
!2


ν
(1 − ε) n−d
7
ψ(ε, d/n)νn2
νn − ε
≥ 1 − exp −
, ψ(ε, d/n) :=
,
ε2
2
8
3 + ε + 1−ε
(24)
where the last inequality follows from condition (9). Conditional on the event {kP⊥
X ∆k2 ≥
ε}, which holds with the probability as stated, the left hand side of (23) is lower bounded
b > − Π∗> )yk2 .
by ε2 k(Π
ii) Upper bounding the right hand side of (23).
Consider t > 0 arbitrary, to be chosen later. We have


⊥
⊥
b
P(|hPX , ∆i| ≥ t) ≤ P sup |hPX , ∆i| > t .
∆∈T

⊥
Now note that  and P⊥
X are independent. As a result, when conditioning on PX , we have
from Lemma 14 that






⊥
⊥
⊥
⊥
⊥
⊥
P sup |hPX , ∆i| > t PX ≤ P sup |hPX , ∆i| > E sup |hPX , ∆i| PX + τ PX
∆∈T
∆∈T
∆∈T


⊥
⊥
≤ P sup |hPX , ∆i| > σw(T ) + τ PX
∆∈T


2
2
⊥
⊥
≤ exp(−τ /(2σ )), τ := t − E sup |hPX , ∆i| PX
∆∈T

≥ t − σw(T ).
The second equality is a consequence of the Sudakov-Fernique comparison
√ inequality (Theorem 2.2.3 in [4]) and the fact that kP⊥
k
≤
1.
Choosing
t
=
(1
+
2)σ{w(T ) ∨ log n},
X 2
combining this result with i), and putting together the pieces in (23), the assertion of the
lemma follows.
21

In order to complete the proof of Theorem 1, we need one last lemma whose proof is
deferred to Appendix G.
Lemma 3. Let σmin (·) denote the minimum singular value of a matrix. If n ≥ 9 ∨ 4d, it
holds that

(
)  
r
p
 X > X −1 X > 

√
5(d
∨
log
n)}
4d
∨
log
n

√
q
P  σmin (X/ n) < 1 −
∪
≥σ

n
n
n
4d∨log n 
1
−
2
n


1
≤ 2 exp − (d ∨ log n) .
2
Equipped with Lemmas 2 and 3, we are in position to conclude Theorem 1. Expanding
the expression for βb in (22), we obtain that

βb =

X >X
n

−1

−1 > ∗>
 > −1 > b >
X >X
X Π y
X X
X (Π − Π∗> )y
+
n
n
n
n
 > −1 >
>
∗>
b
X X
X ( + (Π − Π )y)
= β∗ +
,
n
n

b >y
X >Π
=
n



and thus

kβb − β ∗ k2 ≤

X>X
n

−1
√



X>X
n

≤
1−


2

n

−1
√

≤

>
X
√ 
n

>
X
√ 
n

2

n

σ
q

4d∨log n
n

r

X>X
n

−1

>
X
√
n

+

+

b > − Π∗> )yk2
k(Π
√2

1 √
b>
k(Π
σmin (X/ n)

√

n

− Π∗> )yk2

n

5(d ∨ log n) 2(1 +
+
n

√

2)σε−2 {w(T ) ∨ log n}
√
n

with the stated probability by combining Lemmas 2 and 3.

C

Proof of Corollary 1

Observe that the event {Sb = S ∗ } is implied by the event


∗>
>
b
min∗ |yi − yϕ−1 (i) | > kΠ y − Π yk∞ ,
i∈S

which is in turn implied by


b > yk2 =: γ ,
min∗ |yi − yϕ−1 (i) | > kΠ∗> y − Π
i∈S

We now have




P min∗ |yi − yϕ−1 (i) | ≤ γ
i∈S

≤

X

P |yi − yϕ−1 (i) | ≤ γ

i∈S ∗

≤ k max P (|yi − yj | ≤ γ) .
i6=j

22



!

Note that for any i 6= j, yi − yj ∼ N (0, 2(kβ ∗ k22 + σ 2 )). Invoking Lemma 19, we thus obtain
that for any i 6= j
γ
.
P (|yi − yj | ≤ γ) ≤ p
∗
2(kβ k22 + σ 2 )
It hence follows from the previous two displays that for any δ > 0,


P min∗ |yi − yϕ−1 (i) | ≤ γ ≤ δ
i∈S

if

k2 γ 2
.
2δ 2
Substituting γ by the bound from Lemma 2, the above condition becomes
√
2(1 + 2)2 ε−4 k 2 {w(T ) ∨ log n}2
SNR >
·
,
δ2
n
kβ ∗ k22 >

and the first part of the corollary follows from Lemma 2. Turning to the second part,
consider the least squares estimator using the reduced data set {(xi , yi ), i ∈ [n] \ S ∗ } and
b ∗ ) the corresponding least squares estimator. We then have for any δ > 0
denote by β(S
b ∗ ) − β ∗ k2 ≤ δ) = P(kβ(S
b ∗ ) − β ∗ k2 ≤ δ, Sb = S ∗ ) + P(kβ(S
b ∗ ) − β ∗ k2 ≤ δ, Sb 6= S ∗ ),
P(kβ(S
which implies
b
b ∗ ) − β ∗ k2 ≤ δ) − P(Sb 6= S ∗ ).
P(kβb − β ∗ k2 ≤ δ) ≥ P(kβ(S
It thus remains to control the first probability on the left hand side, which can be done
with arguments similar to those used in the last part of the proof of Theorem 1.

D

Proof of Theorem 2

We start with a Lemma that parallels Lemma 1. The reasoning is analogous, and the proof
is hence omitted.
e ee). Then, if n ≥ d, with
Lemma 4. Consider optimization problem (13) with solution (β,
probability one
(

> X −1 X > (y − √ne
√
1 ⊥
X
e)
2
ee ∈ argmin kPX (y − ne)k2 + λkek1 , βe ∈
,
n
n
n
e
o
√
1
2
ee ∈ argmin kP⊥
(y
−
ne)k
+
λkek
1 .
2
n X
e
For our analysis of the optimization problem of Lemma 4, we need the following crucial
lemma whose proof is provided in Appendix F.
n−d
Lemma 5. There exists constants c1 , c2 , ε > 0 so that if k ≤ c1 log(n/k)
, the following event
holds with probability at least 1 − 2 exp(−c2 (n − d)).


1 ⊥ 2
n−d
2
∗
kP vk ≥ ε
kvk2 , v ∈ C(S , 3) ,
n X 2
n

with S ∗ = {j : Π∗jj 6= 1} and C(S ∗ , 3) defined in (31).
23

√
In correspondence with Lemma 2, we first bound ke
e −e∗ k2 , where e∗ = (Π∗ Xβ ∗ −Xβ ∗ )/ n.
Lemma 6. Under the conditions of Theorem 2, with probability at least 1 − 2 exp(−cε (n −
2
d)) − 2n−M , we have
r
n −1 k log n
∗
ε
.
ke
e − e k2 ≤ 48(1 + M )σ
n−d
n
Proof. We first decompose
√ ∗
⊥
∗
⊥ √ ∗
P⊥
ne + e, e := P⊥
X y = PX (Xβ + ne + ) = PX
X .
According to Lemma 4, since ee is a minimizer, we have that
1 ⊥√ ∗
1
ek1 ≤ ke
kPX n(e − ee) + ek22 + λke
k22 + λke∗ k1 .
n
n
Expanding the square on the left-hand side and re-arranging terms, we obtain that
√
1 ⊥√ ∗
kPX n(e − ee)k22 + λke
ek1 ≤ 2| e∗ − ee, e/ n | + λke∗ k1 ,
n
which implies
√
λke
ek1 ≤ 2ke∗ − eek1 ke
/ nk∞ + λke∗ k1 ,

(25)

and in turn
√
λke
e(S ∗ )c − e∗(S ∗ )c k1 = λke
e(S ∗ )c k1 ≤ 2ke∗ − eek1 ke
/ nk∞ + λ(ke∗ k1 − ke
eS ∗ k1 )
√
∗
∗
≤ 2ke − eek1 ke
/ nk +λkeS ∗ − eeS ∗ k1 .
| {z ∞}
=:λ0

Re-arranging, we arrive at
(λ − 2λ0 )ke
e(S ∗ )c − e∗(S ∗ )c k1 ≤ (λ + 2λ0 )ke∗S ∗ − eeS ∗ k1 .
Conditional on the event Λ = {λ ≥ 4λ0 }, we hence have that
ke
e(S ∗ )c − e∗(S ∗ )c k1 ≤ 3ke
eS ∗ − e∗S ∗ k1 .

(26)

Equipped with this is intermediate result, we go back to (25). We first obtain by rearranging that
1 ⊥√ ∗
n(e − ee)k22 ≤ (2λ0 + λ)ke
e − e∗ k1
kP
n X
≤ (2λ0 + λ)4ke
eS ∗ − e∗S ∗ k1
√
≤ (2λ0 + λ)4 kke
eS ∗ − e∗S ∗ k2 ,
where we have used (26) in the second inequality. In order to lower bound the left hand
side, note that (26) can be written as ee − e∗ ∈ C(S ∗ , 3) according to (31). Conditional on
the event in Lemma 5,
n−d ∗
1 ⊥√ ∗
kP
n(e − ee)k22 ≥ ε
ke − eek22 .
n X
n
Combining this with the previous display, we find that
√
n
ke∗ − eek2 ≤ 12ε−1
λ k.
n−d
p
In view of Lemma 18, the choice λ = 4(1 + M )σ 2 log(n)/n for M ≥ 0 guarantees that
2
the event Λ defined above (26) occurs with probability at least 1 − 2n−M . This completes
the proof of the Lemma.
Using Lemmas 2 and 4, the proof of the theorem is completed along the lines of the last
part of the proof of Theorem 1. For brevity, we omit those steps here.
24

E

Proof of Theorem 3

The proof of part (a) of the theorem relies on two lemmas, which are stated and proved
first.
Lemma 7. Let γ, δ > 0. Consider the event


> ∗
> ∗ 2
2
Aγ = min(xi β − xj β ) ≤ γ .
i<j

(27)

Then P (Aγ ) < δ if
n2 (n − 1)2 2
kβ ∗ k22 >
γ .
4δ 2 π
o
S n
∗ − x> β ∗ )2 ≤ γ 2 . Using the fact that x> β ∗ −
Proof. It is clear that Aγ ⊆ i<j (x>
β
i
j
i
∗
∗ 2
x>
j β ∼ N (0, 2kβ k2 ), we have



∗
[
γ
β
γ

≤ (xi − xj )> √
≤√
−√
P (Aγ ) ≤ P 
∗
∗
∗k
2kβ
k
2kβ
k
2kβ
2
2
2
i<j

X 
γ
γ
≤
P −√
≤g≤ √
, g ∼ N (0, 1)
2kβ ∗ k2
2kβ ∗ k2
i<j

≤

n(n − 1)
γ
√
∗
2
kβ k2 π

Thus, it suffices to require that

by Lemma 19 in Appendix H.

n(n − 1)
γ
√ < δ. Resolving for kβ ∗ k2 , this condition
∗
2
kβ k2 π

becomes
kβ ∗ k22 >

n2 (n − 1)2 2
γ .
4δ 2 π

The following lemma is of interest in its own when the design matrix is fixed and both β ∗
and σ 2 are known or accurate estimates of these quantities are available. Condition (28)
can then be evaluated explicitly, at least after substituting β ∗ and σ 2 by their respective
estimates. The result is close in spirit to those in the framework in [17].


b ∗ ) 6= Π∗ | X < δ if
Lemma 8. Conditional on X, we have P Π(β

2
n(n − 1)
∗
> ∗
min x>
β
−
x
β
> 4σ 2 log
.
i
j
i<j
δ

(28)

b ∗ ),
Proof. Without loss of generality, we may assume that Π∗ = In . By construction of Π(β
∗
∗
>
∗
>
∗
b
it is clear that Π(β ) 6= Π whenever there exists i 6= j such that (yi −yj )(xi β −xj β ) < 0,
∗
> ∗
that is, the order between x>
i β and xj β is different than that of yi , yj . First, conditioning
on the random design matrix X, we have


o
[n
∗
> ∗
b ∗ ) 6= Π∗ | X) = 2 P 

P(Π(β
yi − yj > 0 | x>
i β − xj β < 0
i<j

≤2

X



∗
> ∗
P yi − yj > 0 | x>
β
−
x
β
<
0
i
j

i<j

25

∗
> ∗
> ∗
> ∗
2
Since yi − yj | x>
i β − xj β ∼ N (xi β − xj β , 2σ ), using the usual tail bound for the
Gaussian distribution for each term in the above sum, we obtain
 
2 
> β ∗ − x> β ∗
x
X
i
j


b ∗ ) 6= Π∗ | X) ≤ 2
P(Π(β
exp −

2
4σ
i<j



 mini<j
≤ n(n − 1) exp −



∗
x>
i β

−

4σ 2

∗
x>
j β

2 

.

b ∗ ) 6= Π∗ | X) < δ if
Hence, for any δ > 0, P(Π(β

2
n(n − 1)
∗
> ∗
> 4σ 2 log
min x>
.
i β − xj β
i<j
δ

b ∗ ) from the
Combining Lemmas 7 and 8, part (a) in Theorem 3 q
follows immediately for Π(β
b ∗ ) 6=
previous lemma by invoking Lemma 7 with γ = 2σ log n(n−1) . We then have P(Π(β
δ

Π∗ ) ≤ P(Aγ ) + (1 − P(Aγ )) δ < 2δ under the conditions of the theorem. The general
version of the theorem results from
>b
> ∗
> ∗
> ∗
b
b
min |x>
i θ − xj θ| ≥ min |xi β − xj β | − 2 max |xi (β − θ)|
i<j

i<j

≥

1≤i≤n

∗
min |x>
i β
i<j

−

∗
x>
j β |

− 2σ∆.

∗
b
conditional on the event {kXβ
The assertion follows from an application
q − X θk∞ ≤ σ∆}.

of Lemma 7 with γ = 2σ
log n(n−1)
+∆ .
δ

The proof of part (b) in Theorem 3 requires another lemma.
Lemma 9. Consider the event Aγ in (27). Then, for all n ≥ 5, P (Aγ ) ≥ 0.85 if
kβ ∗ k22 ≤

1 2 2
n γ .
25π

∗
∗
Proof. Since x>
i β /kβ k2 ∼ N (0, 1), we have for g ∼ N (0, 1)




γ
γ
γ
γ
γ
> ∗
√ ,
p := P − √ < x β < √
= P −√
<g< √
≥
2
2
2kβ ∗ k2
2kβ ∗ k2
kβ ∗ k2 2π

(29)

(30)

whenever γ < √32 kβ ∗ k2 , otherwise p > .86 as can be verified by evaluating the Gaussian
cumulative
 distribution functionnumerically.
γ
γ
Let B := − √ < x> β ∗ < √ . For i.i.d. observations x1 , . . . , xn , let
2
2
u :=

n
X

I(xi ∈ B) ∼ binom(n, p),

i=1

26

with p as defined in (30). The event Aγ in (27) occurs whenever at least two xi , xj ∈ B,
i.e., the event {u ≥ 2} ⊆ Aγ . Let us suppose for the moment that p ≥ 3.5/n for n ≥ 5.
Then
P (Aγ ) ≥ P (u ≥ 2)



 
3.5 n−1
3.5 n
3.5
1−
≥1− 1−
−n
n
n
n
> .85 for n ≥ 5.
Using Eq.(30), we now conclude that P (Aγ ) ≥ .85 whenever
holds if
kβ ∗ k22 ≤

γ
kβ ∗ k2

√

2π

≥ 3.5/n. This

1 2 2
n γ .
25π

We finally turn to the proof of part (b) of the theorem. Define indices i0 , j0 ⊂ [n] by the
∗
> ∗
> ∗
> ∗
relation x>
i0 β − xj0 β = mini<j |xi β − xj β |. For given γ > 0, we have

 

b ∗ ) 6= Π∗ ) ≥ P yi − yj < 0 | 0 < x> β ∗ − x> β ∗ < γ P 0 < x> β ∗ − x> β ∗ < γ
P(Π(β
i0
j0
i0
j0
0
0


∗
> ∗
> ∗
> ∗
= P x>
β
+
σ
−
x
β
−
σ
<
0
|
0
<
x
β
−
x
β
<
γ
×
i0
j0
i0
j0
i0
j0


∗
> ∗
× P 0 < x>
i0 β − xj0 β < γ


∗
> ∗
≥ P (σi0 − σj0 < −γ) P 0 < x>
β
−
x
β
<
γ
i0
j0
 


γ
∗
> ∗
P 0 < x>
β
−
x
β
<
γ
≥ P g < −√
i0
j0
2σ




√
√
4 πkβ ∗ k2
4 πkβ ∗ k2
> ∗
> ∗
c ,
c P 0 < xi0 β − xj0 β <
≥P g<− √
n
2σn
√
4 πkβ ∗ k2
by choosing γ =
c, for 0 < c ≤ 1.
n
√
4 πkβ ∗ k2
≥ Φ(−.1)(.85), by Lemma 9, for a choice of c s.t. √
c ≤ .1.
2σn
≥ .35,
Rx
where in the penultimate line, Φ(x) = −∞ √12π exp(−t2 /2) dt denotes the Gaussian cub ∗ ) 6= Π∗ ) ≥ .35 whenever
mulative distribution function. Therefore, P(Π(β
some C > 0. This concludes the proof of part (b) of the theorem.

F

kβ ∗ k22
σ2

<

C
n2

for

Proof of Lemma 5

The proof of Lemma 5 is given at the end of this section. The proof depends on several
lemmas collected from the literature which are stated first.
For 1 ≤ s ≤ p, let J (s) denote the set of all index subsets of [p] of size s.
Definition 1. For J ∈ J (s), and α ∈ [1, ∞), define
C(J, α) = {v ∈ Rp : kvJ c k1 ≤ αkvJ k1 }.
27

(31)

Definition 2. [14] Let A ∈ RN ×p and B0 (s; p) = {v ∈ Rp : kvk0 ≤ s}. Define the s-sparse
minimum and maximum eigenvalues associated with A by
λmin (s) =

kAvk22
,
v∈B0 (s;p)\{0} kvk2
2
min

λmax (s) =

kAvk22
.
v∈B0 (s;p)\{0} kvk2
2
max

Let δ ∈ (0, 1). We say that A is an (s, δ)-restricted isometry if
1 − δ ≤ λmin (s) ≤ λmax (s) ≤ 1 + δ.
Definition 3. [9] Let A ∈ RN ×p . We define the (s, α)-restricted eigenvalue by
φ(s, α) = min

min

J∈J (s) v∈C(J,α)\0

kAvk22
.
kvJ k22

(32)

Lemma 10. [9] Suppose that A is a (3s, δ ∗ /4)-restricted isometry for some δ ∗ ∈ (0, 1).
The (s, 3)-restricted eigenvalue then satisfies φ(s, 3) ≥ 1 − δ ∗ .
Lemma 11. [6] Let A be a random of matrix of size N × p satisfying
P((1 − ε)kvk22 ≤ kAvk22 ≤ (1 + ε)kvk22 ) ≤ 2 exp(−N c0 (ε))

(33)

for any ε ∈ (0, 1) and any v ∈ Rp , where c0 (ε) > 0 is a constant depending only on
ε. Then for any δ ∈ (0, 1), there exist constants c1 , c2 > 0 depending only on δ such
that A is a (q, δ)-restricted isometry for any q ≤ c1 N/ log(p/q) with probability at least
1 − 2 exp(−c2 N ).
Lemma 12. (from Lemma 2.2 in [18]) Let P denote the orthogonal projection on an N dimensional subspace of Rp chosen uniformly at random from the Grassmannian G(p, N ).
Then, for any v ∈ Rp and any ε ∈ (0, 1),


p
P (1 − ε)kvk22 ≤
kP vk22 ≤ (1 + ε)kvk22 ≤ 2 exp(−N (ε2 /4 − ε3 /6)})
N
Proof of Lemma 5: the lemma is essentially implied by the fact that under the stated
conditions, the matrix P⊥
X has its (k, 3)-restricted eigenvalue bounded away from zero by
a constant ε as to be shown below. ”Essentially” refers to the fact that it can be shown
(cf. [54]) that a matrix A with restricted eigenvalue φ(k, 3) obeys
[
1
kAvk22 ≥ φ(k, 3)kvk22
∀v ∈
C(J, 3).
32
J∈J (k)

As sketched above, a lower bound on φ(k, 3) follows if A is a (3k, δ ∗ /4) restricted isometry
for some δ ∗ ∈ (0, 1). That property is in turn satisfied with high probability if A satisfies
(33). Applying Lemma 12 with P = P⊥
X , p = n and N = n − d yields the conclusion.

G

Proof of Lemma 3

 > −1 >
X
√
We invoke Lemma 17 conditional on X, with X n X
in the role of A, so that
n
 > −2 >
X
√ . Observe that
Γ = √Xn X n X
n
 > −1 !
X X
d
1
√ , (ii) kΓk2 = 2
√ ,
(i) tr(Γ) = tr
≤ 2
n
σmin (X/ n)
σmin (X/ n)
v
(34)
u
√
 > −2 !
u
p
X
X
d
√ .
(iii)
tr(Γ2 ) = ttr
≤ 2
n
σmin (X/ n)
28

In view of the relations in (34), application of Lemma 17 with t = d ∨ log n and comparison
of terms yields that
!
p
 > −1 >
5{d ∨ log n}
X X
X 
√
√
P
>σ
X ≤ exp(−d ∨ log n).
n
n
σmin (X/ n)
2

Using Lemma 16 for σmin (X) with the choice t =

H

√

d ∨ log n, the assertion readily follows.

Additional lemmas

We here collect various auxiliary results acting as lemmas in the proofs of our main results.
Lemma 13. (Gordon’s Escape through a Mesh theorem [27]) Let K be a closed
subset of the unit sphere in Rm , let νr = Eg∼N (0,Ir ) [kgk2 ], and let ε > 0. If the Gaussian
width (6) of K obeys w(K) < (1 − ε)νl − ενm , then an (m − l)-dimensional subspace V
drawn uniformly from the Grassmannian G(m, m − l) satisfies

 !
1 (1 − ε)νl − ενm − w(K) 2
7
P(dist(K, V ) > ε) ≥ 1 − exp −
.
2
2
3 + ε + ενm /νl
Lemma 14. (Concentration of Gaussian processes [11]) Let K be a closed subset of
the unit sphere in Rm with Gaussian width w(K), and let g ∼ N (0, Im ). Then for all t > 0


P sup | hg, xi | ≥ w(K) + t ≤ exp(−t2 /2).
x∈K

Lemma 15. (Concentration of the norm of Gaussian random vectors [51])
√
Let g ∼ N (0, In ). Then for all t > 0, P(kgk2 ≥ (1 + t) n) ≤ exp(−t2 n/2).
Lemma 16.(Concentration of extreme singular values of Gaussian matrices [51])
Let X be an n × d Gaussian matrix with i.i.d. N (0, 1)-entries. Denote by σmin (X) and
σmax (X) the minimum respectively maximum singular value of X. Then for any t > 0
√
√
√
√
P(σmin (X) ≥ n − d − t) ≤ exp(−t2 /2),
P(σmax (X) ≥ n + d + t) ≤ exp(−t2 /2).
Lemma 17. (Concentration of quadratic forms [31])
Let A be a an m × n matrix, Γ = A> A, and g ∼ N (0, In ). Then for all t > 0
p
P(kAgk22 > tr(Γ) + 2 tr(Γ2 )t + 2kΓk2 t) ≤ exp(−t).
Lemma 18. (Maximal inequality for Gaussian random variables)
Let g1 , . . . , gn be zero-mean Gaussian random variables such that max1≤i≤n E[gi2 ] ≤ σ 2 .
√
2
Then for all M > 0, P(max1≤i≤n |gi | > (1 + M )σ 2 log n) ≤ 2n−M .
Lemma 19. (Small ball probability)
Let g ∼ N (0, 1). Then for all t > 0, P(|g| < t) ≤

29

q

2
π

t < t.

