Signal Amplitude Estimation and Detection
from Unlabeled Binary Quantized Samples
Guanyu Wang, Jiang Zhu, Rick S. Blum, Peter Willett,

arXiv:1706.01174v3 [] 8 Mar 2018

Stefano Marano, Vincenzo Matta, Paolo Braca
Abstract
Signal amplitude estimation and detection from unlabeled quantized binary samples are studied,
assuming that the order of the time indexes is completely unknown. First, maximum likelihood (ML)
estimators are utilized to estimate both the permutation matrix and unknown signal amplitude under
arbitrary, but known signal shape and quantizer thresholds. Sufficient conditions are provided under
which an ML estimator can be found in polynomial time and an alternating maximization algorithm is
proposed to solve the general problem via good initial estimates. In addition, the statistical identifiability
of the model is studied.
Furthermore, the generalized likelihood ratio test (GLRT) detector is adopted to detect the presence
of signal. In addition, an accurate approximation to the probability of successful permutation matrix
recovery is derived, and explicit expressions are provided to reveal the relationship between the number
of signal samples and the number of quantizers. Finally, numerical simulations are performed to verify
the theoretical results.
Index Terms
Estimation, detection, permutation, unlabeled sensing, quantization, identifiability, alternating
maximization.

I. I NTRODUCTION
In many systems, the data is transmitted with time information, which may sometimes be imprecise
[1], [2], [3], [4], [5], [6], [7]. One example is the global positioning system (GPS) spoofing attack which
can alter the time stamps on electric grid measurements [1] to make them useless so that the data must
be processed without time stamps. Since the exact form of civilian GPS signals is publicly known and
the elements needed are inexpensive, building a circuit to generate signals to spoof the GPS is easy.
In [2], a refined assessment of the spoofing threat is provided. In addition, the detailed information
of receiver-spoofer architecture, its implementation and performance, and spoofing countermeasures are

DRAFT

1

introduced. As a case study in [3], the impact of the GPS spoofing attack on wireless communication
networks, more specifically, the frequency hopping code division multiple access (FH-CDMA) based
ad hoc network, is investigated. A timing synchronization attack (TSA) is coined to the wide area
monitoring systems (WAMSs), and its effectiveness is demonstrated for three applications of a phasor
measurement unit (PMU) [4]. In [5], the out-of-sequence measurement (OOSM) problem where sensors
produce observations that are sent to a fusion center over communication networks with random delays are
studied, and a Bayesian solution is provided. The problem of random delay and packet loss in networked
control systems (NCS) is studied in [6]. In addition, a minimum error covariance estimator for the system
is derived and two alternative estimator architectures are presented for efficient computation. In [7], the
effect of an unknown timestamp delay in Automatic Identification System (AIS) is studied, and a method
based on adaptive filtering is proposed.
In the above examples, the relative order of the data is unknown, i.e., the samples are unlabeled.
Estimation and detection from unlabeled samples have drawn a great deal of attention recently [8], [9],
[10], [11], [12], [13], [14], [15], [16], [17], [18], [19], [20]. In [8], it is shown that the convex relaxation
based on a Birkhoff polytope approach does not recover the permutation matrix, and a global branch and
bound algorithm is proposed to estimate the permutation matrix. In the noiseless case with a random
linear sensing matrix, it is shown that the permutation matrix can be recovered correctly with probability
1, given that the number of measurements is twice the number of unknowns [9], [18]. In [10], [19],

the noise is taken into account and a condition under which the permutation matrix can be recovered
with high probability is provided. In addition, a polynomial time algorithm is proposed for a scalar
parameter case. Denoising linear regression model with shuffled data and additive Gaussian noise are
studied in [11], and the characterization of minimax error rate is provided. In addition, an algorithm for
the noiseless problem is also proposed, and its performance is demonstrated on an image point-cloud
matching task [11]. In [12], several estimators are compared in recovering the weights of a noisy linear
model from shuffled labels, and an estimator based on the self-moments of the input features and labels is
introduced. For unlabeled ordered sampling problems where the relative order of observations is known,
an alternating maximization algorithm combined with dynamic programming is proposed [13]. In [15],
a signal detection problem where the known signal is permuted in an unknown way is studied.
Compared to the location parameter estimation problem (xi = θ + wi ) in [17], the model in this paper
is a scale parameter estimation problem (xi = hi θ + wi ), in which hi , i = 1, · · · , K is the shape of a
signal, and θ is an amplitude of signal. As a result, the scale parameter estimation problem is much more
difficult than the location estimation in several aspects, and the scale parameters are especially relevant

DRAFT

2

in relation to the mislabeling/permutation issue. First, the model in [17] is always identifiable, while our
model may be unidentifiable, as shown later. Second, the problem in [17] can be solved efficiently via
simple sorting, while we can only prove that problem in this paper can be solved efficiently provided
certain conditions are satisfied. Third, good initial points are proposed to improve the performance of
alternating maximization algorithm. Furthermore, we provide an approximation to the probability of
successful permutation matrix recovery, which reveals the relationship between the length of signal and
the number of quantizers.
In this paper, we focus on the problems of scale estimation and signal detection from unlabeled
quantized samples. The main contribution of this work can be summarized as follows. First, a sufficient
condition for the existence of a polynomial time algorithm is provided for the unlabeled estimation
problem, and the model is shown to be unidentifiable in some special cases. Second, good initial
points are provided to improve the performance of an alternating maximization algorithm. And third,
we provide analytic approximations on probability of permutation matrix recovery in the case of known
signal amplitude, which can be used to predict when the permutation matrix can be correctly recovered.
The organization of this paper is as follows. In Section II, the problem is described. Background on
ML estimation and generalized likelihood ratio test (GLRT) detection from labeled data are presented in
Section III. In section IV, the model identifiability is studied, and the estimation problem from unlabeled
data is studied. Section V extends the detection work to unlabeled data, and derives an approximate
analytic formula for permutation matrix recovery probability. Finally, numerical results are presented in
Section VI, and conclusion follows in Section VII.
Notation: The K ×1 vector of ones is 1K . For an unknown deterministic parameter θ, θ0 denotes its true
value. For an unknown permutation matrix Π, Π0 denotes its true value. For a random vector y, p(y; θ)
denotes the probability density function (PDF) of y parameterized by θ, and Ey [·] denotes the expectation
taken with respect to y. Let N (µ, σ 2 ) denote a Gaussian distribution with mean µ and variance σ 2 . Let
Φ(·) and ϕ(·) denote the cumulative distribution function (CDF) and probability density function (PDF)

of a standard Gaussian random variable respectively. Let U(a, b) denote an uniform distribution, whose
minimum and maximum values are a and b. Let B(N, p) denote a binomial distribution, where N and p
denote the number of trials and the probability of event, respectively.
II. P ROBLEM S ETUP
Consider a signal amplitude estimation and detection problem where a collection of N binary quantizers
generate binary quantized samples which will be utilized to estimate the unknown scaling factor θ of a
K length signal and detect the presence of the signal, as shown in Fig. 1. The binary quantized samples

DRAFT

3

Ampulitude
Noiseless
Signal

hi θ

···
0 1

···
i

···

···

K

Time

Noise
Quantizers
Q1

···

Qj

···

QN
τi

τK

b11

bi1

bK1

···

b1j

···

···
···

bij

···

bKj

···

···

···

Binary
quantized
samples

τ1

b1N

biN

bKN

Binary channel
u11

···

u1j

···

···

···

uij

···

u1N

uiN

···

uKj

···

···

Flipped
samples

uK1

ui1

uKN

Permuted

···

uij

···

uKj

···

···

···
u1N

···

u1j

···

···

Unlabeled
samples

uK1

ui1

u11

uiN

uKN

Fig. 1: System diagram of unlabeled binary quantized samples generation.

bij are obtained via
bij = Qi (hi θ + wij ),

i = 1, · · · , K,

j = 1, · · · , N,

(1)

and the corresponding hypothesis problem can be formulated as


 H0 : bij = Qi (wij ), i = 1, · · · , K, j = 1, · · · , N,

 H1 : bij = Qi (hi θ + wij ), i = 1, · · · , K, j = 1, · · · , N,
where i and j respectively denote one of the K time indexes and one of the N quantizers, hi is the known
2 -variance distribution
coefficient characterizing the signal shape, wij is the i.i.d. noise drawn from the σw

whose PDF is fw (x/σw )/σw and CDF is Fw (x/σw ), where fw (x) and Fw (x) are the corresponding
unit-variance PDF and CDF, and Qi (·) implies a binary quantizer which produces 1 if the argument is

DRAFT

4

larger than a scalar threshold τi and 0 otherwise. The thresholds of N quantizers are identical given
any time index1 . We assume that the PDF fw (w) is log-concave, which is often met in practice such
as Gaussian distributions, and the thresholds of N quantizers are identical given any time index. We
assume that the PDF fw (w) is logconcave, which is often met in practice by, for example, the Gaussian
distribution.
The quantized data {bij } are transmitted over a binary channel with flipping probabilities q0 and q1
which are defined as Pr(uij = 1|bij = 0) = q0 and Pr(uij = 0|bij = 1) = q1 , where uij is the sample
received at the output of the channel, which we call the fusion center (FC) [21].
We assume that all the sets of data {uij }N
j=1 are transmitted to the FC with permuted time indexes.
Accordingly, the FC receives sets of data, say {ũij }N
j=1 , in which the time reference (represented by the
index i) is invalid. Specifically, the FC does not know which time index the data {ũij }N
j=1 belongs to,
but knows that {ũij }N
j=1 belongs to one of the K time indexes. Let us introduce the matrix U whose
e , as follows:
(i, j)−th entry is uij . Then, the unlabeled samples can be collected in a matrix U
e = ΠU,
U

(2)

where Π ∈ RK×K is an unknown permutation matrix; that is, a matrix of {0,1} entries in which each
row and each column sums to unity. We assume that θ is constrained to an interval [−∆, ∆], for algorithm
and theoretical reasons [22].

III. P RELIMINARIES
In this section, standard materials of parameter estimation and signal detection using labeled data are
presented.

A. Maximum Likelihood Estimation
The probability mass function (PMF) of uij can be calculated as

Pr(uij = 1) = q0 + (1 − q0 − q1 )Fw

hi θ − τi
σw


, pi ,
(3)

Pr(uij = 0) = 1 − pi .

1

Here we have thresholds fixed across quantizers and varying with time, with permutation across time. We could, equivalently,
have fixed thresholds of quantizers across time but varying across sensors and permuted across sensors. Mathematically, it is
the same problem and the formulation could as easily encompass it.

DRAFT

5

The PMF of U is
p(U; θ) =

K Y
N
Y

Pr(uij = 1)uij Pr(uij = 0)(1−uij ) .

(4)

i=1 j=1

Let ηi denote the fraction of uij = 1 in {uij }N
j=1 , i.e.,
ηi =

N
X

uij /N.

(5)

j=1

Consequently, the log-likelihood function l(η; θ) is
l(η; θ) = N

K
X
(ηi log pi + (1 − ηi ) log(1 − pi )),

(6)

i=1

where pi is given in (3). Note that in an error free binary symmetric channel scenario, i.e., q0 = q1 = 0 or
q0 = q1 = 1, the CDF Fw (x) is log-concave as it is the integral of a log-concave PDF fw (x). Therefore

maximizing the log-likelihood function is a convex optimization problem, which can be solved efficiently
via numerical algorithms [23], [24], [25]. For 0 < q0 + q1 < 2, it is difficult to determine the convexity
of the negative log-likelihood function. In this case all that can be guaranteed is a local optimum. As we
show in numerical experiments, we found that the ML estimator using gradient descent algorithm works
well and approaches the Cramér Rao lower bound (CRLB).
In addition, the Fisher Information (FI) I(θ) is the expectation of the negative second derivative of the
log-likelihood function l(η; θ) (6) taken with respect to θ, i.e., [28],
( 
  


 
)
K
hi θ−τi
∂ ηi 1−ηi
∂
hi θ−τi
ηi 1−ηi
N (1−q0 −q1) X
hi fw
Eη
−
+ fw
Eη
−
I(θ) = −
σw
σw
∂θ pi 1−pi
∂θ
σw
pi 1−pi
i=1


K h2 f 2 hi θ−τi
2
X
i w
σw
N (1−q0 −q1)
=
,
(7)
2
σw
pi (1 − pi )
i=1

where (7) follows due to Eη [ηi /pi − (1 − ηi )/(1 − pi )] = 0. Consequently, the CRLB is
CRLB(θ) = 1/I(θ),

(8)

which is later used as a benchmark performance for ML estimation from labeled data in Section VI.

B. GLRT detection
In the case of known θ, the optimal detector according to the NP criterion is the log-likelihood ratio
test [26]. For unknown θ, the GLRT is usually adopted. Although there is no optimality associated with

DRAFT

6

the GLRT, it appears to work well in many scenarios of practical interest [27]. The GLRT replaces the
unknown parameter by its MLE and decides H1 if
T1 (η) =

l(η; θ) − l(η; 0) > γ,

max

(9)

θ∈[−∆,∆]

where γ is a threshold determined by the given false alarm probability PF A .

IV. E STIMATION FROM UNLABELED DATA
In this section, we study the estimation problem from unlabeled data. First, we delineate the model.
The statistical identifiability is investigated, and the estimation problem is studied separately in the cases
of known and unknown θ.

A. Maximum likelihood estimation
Introduce the function π(·) such that m = π(i) if the permutation matrix Π in (2) maps the ith row
of U to the mth row of Ũ. The PMF of Ũ is
e θ, Π) =
p(U;

K Y
N
Y

Pr(ũmj = 1)ũmj Pr(ũmj = 0)(1−ũmj )

m=1 j=1

=

K Y
N
Y

(10)

Pr(ũπ(i)j = 1)ũπ(i)j Pr(ũπ(i)j = 0)(1−ũπ(i)j ) ,

i=1 j=1

where (Pr(ũij = 1), Pr(ũij = 0)) is the PMF of ũij . The corresponding log-likelihood function l(η̃; θ, Π)
is
l(η̃; θ, Π) = N

K
X


η̃π(i) log pi + (1 − η̃π(i) ) log(1 − pi ) ,

(11)

i=1

where η̃π(i) =

PN

j=1 ũπ(i)j /N

=

PN

j=1 ũmj /N .

The ML estimation problem can be formulated as

max

l(η̃; θ, Π),

θ∈[∆,∆],Π∈PK

(12)

where PK denotes the set of all possible K × K permutation matrices.

B. Estimation with permuted data and known θ
In this subsection, the permutation matrix recovery problem is studied in the case of known θ. It is
shown that the permutation matrix under ML estimation criterion can be recovered efficiently under the
ML criterion.

DRAFT

7

Proposition 1 Given the ML estimation problem in (12) with known θ, the ML estimate of the permutation
e , and equivalently the elements of η̃ , to have the same relative order
matrix Π will reorder the rows of U
as the elements of (1 − q0 − q1 )(hθ − τ ).
Proof: Note that the objective function l(η̃; θ, Π) (11) can be decomposed as
l(η̃; θ, Π) = K

N
X

η̃π(i) si + K

i=1

N
X

log(1 − pi ),

(13)

i=1

where si = log(pi /(1 − pi )). From (13), the ML estimate of the permutation matrix Π will reorder the
rows of U, and equivalently the elements of η̃ to have the same relative order as the elements of s
[15], [17]. Because si is monotonically increasing with respect to (1 − q0 − q1 )(hi θ − τi ), the elements
of η̃ should be reordered by the permutation matrix to have the same relative order as the elements of
(1 − q0 − q1 )(hθ − τ ) to maximize the likelihood.

If τi = c0 hi , then changing θ − c0 in −(θ − c0 ) would reverse the ordering. This might help to explain
why two solutions appear in the subsequent Proposition 3 when θ is unknown.

C. Estimation with permuted data and unknown θ
In general, θ may be unknown. Consequently, we should jointly estimate θ and permutation matrix Π.
However, finding the best permutation matrix is very challenging in most problems due to non-convexity.
One could try all the possible permutation matrices, at complexity cost O(N !). Given a permutation
matrix, one obtains the ML estimate of θ via numerical algorithms and achieves global optimum under
q0 = q1 = 0 or 1. Under 0 < q0 + q1 < 2, we do not know whether the negative log-likelihood function

is convex or not, and local optimum is guaranteed. Given θ, the computation complexity of finding the
optimal permutation matrix is just reordering, which costs O(N log N ), as we show in subsection IV-B.
1) Alternating maximization algorithm for general case: The problem structure induces us to optimize
the two unknowns alternately as shown in Algorithm 1. The alternating maximization in Algorithm 1 can
Algorithm 1 Alternating Maximization
Initialize t = 1 and θ̂t−1 ;
Fix θ = θ̂t−1 , reorder η̃ according to (1 − q0 − q1 )(hθ − τ ) and obtain the corresponding permutation
matrix Π̂t−1 ;
3: Solve max l(η̃; θ, Π̂t−1 ) and obtain θ̂t ;
θ
4: Set t = t + 1 and return to step 2 until a sufficient number of iterations has been performed or
|θ̂t − θ̂t−1 | ≤ , where  is a tolerance parameter.

1:
2:

be viewed as the alternating projection with respect to θ and Π. The objective function is l(η̃; θ, Π). In

DRAFT

8

step 2, given θ̂t−1 , we update the permutation matrix as Π̂t−1 , and the objective value is l(η̃; θ̂t−1 , Π̂t−1 ).
Given Π̂t−1 , we obtain ML estimation of θ as θ̂t , and the objective value is l(η̃; θ̂t , Π̂t−1 ) satisfying
l(η̃; θ̂t , Π̂t−1 ) ≥ l(η̃; θ̂t−1 , Π̂t−1 ). Given θ̂t , we update the permutation matrix as Π̂t , and the objective

value is l(η̃; θ̂t , Π̂t ) satisfying l(η̃; θ̂t , Π̂t ) ≥ l(η̃; θ̂t , Π̂t−1 ). Consequently, we have
l(η̃; θ̂t , Π̂t ) ≥ l(η̃; θ̂t−1 , Π̂t−1 ).

(14)

Provided that the maximum with respect to each θ and Π is unique, any accumulation point of the
sequence generated by Algorithm 1 is a stationary point [29].
2) Special cases for efficient recovery of Π under unknown θ:
Proposition 2 Given the ML estimation problem in (12) with unknown θ, if there exist constants c, d, e ∈
R such that cτ + dh = e1, the elements of η̃ should be reordered according to the order of the elements
of (q0 + q1 − 1)τ if c = 0, otherwise reordered according to h or −h.
Proof: We separately address the cases c = 0 and c 6= 0. In the case of c = 0, h must be a
constant vector. Reordering according to (1 − q0 − q1 )(hθ − τ ) is equivalent to reordering according
to (q0 + q1 − 1)τ . Since (q0 , q1 ) are known in this problem, η̃ should be reordered according to τ if
q0 + q1 > 1 or −τ if q0 + q1 < 1. In the case of c 6= 0, we have τ = (e/c)1 − (d/c)h. Consequently,
hθ − τ = (θ + d/c)h − (e/c)1, and η̃ is reordered according to h or −h.

The above proposition deals with four cases, i.e., h is a constant vector (c = 0), τ is a constant
vector (d = 0), h is a multiple of τ (e = 0) and each pair of components of h and τ lies in the same
line cτi + dhi = e (cde 6= 0). In [17] it is shown that reordering yields the optimal MLE given h = 1.
Proposition 2 extends the special case in [17] to more general cases. Consequently, we propose Algorithm
2, an efficient algorithm for parameter estimation.
Algorithm 2 Reordering algorithm
1:

If c = 0, reorder the elements of η̃ according to the elements of (q0 + q1 − 1)τ . The corresponding
permutation matrix is Π̂s0 . Solve the parameter estimation problem by numerical algorithm and
obtain θ̂ML = argmax l(η̃; θ, Π̂s0 );
θ

If c 6= 0, reorder the elements of η̃ according to the elements of h and −h. The corresponding
permutation matrices are Π̂s1 and Π̂s2 ;
3: Solve the single variable optimization problems and obtain θ̂s1 = argmax l(η̃; θ, Π̂s1 ) and θ̂s2 =

2:

θ

argmax l(η̃; θ, Π̂s2 ). Choose θ̂ML = θ̂s1 given that l(η̃; θ̂s1 , Π̂s1 ) ≥ l(η̃; θ̂s2 , Π̂s2 ), otherwise θ̂ML =
θ

θ̂s2 .

DRAFT

9

D. Statistical identifiability
Note that Algorithm 2 may generate two solutions (θ̂s1 , Π̂s1 ) and (θ̂s2 , Π̂s2 ). Given system parameters
h and τ , it is important to determine whether the two solutions (θ̂s1 , Π̂s1 ) and (θ̂s2 , Π̂s2 ) will yield

the same log-likelihood l(η̃; θ̂s1 , Π̂s1 ) = l(η̃; θ̂s2 , Π̂s2 ). If l(η̃; θ̂s1 , Π̂s1 ) = l(η̃; θ̂s2 , Π̂s2 ), two pairs
of parameter values lead to the same maximum likelihood. In this situation, (θ, Π) clearly cannot be
estimated consistently since η̃ provide no information as to whether it is (θ̂s1 , Π̂s1 ) or (θ̂s2 , Π̂s2 ).
This phenomenon motivates us delving into the identifiability of the model. Statistical identifiability
is a property of a statistical model which describes one-to-one correspondence between parameters and
probability distributions [32]. In this subsection, we provide the following proposition which justifies
that there exist cases in which the model is unidentifiable, i.e., there exist two different parameter values
((θs1 , Πs1 ) and (θs2 , Πs2 )) leading to the same distribution of the observations η̃ [32].
Proposition 3 Let ha and hd denote the ascending and descending ordered versions of h, and Πa h = ha
and Πd h = hd , where Πa and Πd are permutation matrices. Given τ = c0 h and ha = −hd , the
model is unidentifiable, i.e., l(η̃; θ, Π)|θ=θs1 ,Π=Πs1 = l(η̃; θ, Π)|θ=θs2 ,Π=Πs2 , where θs2 = 2c0 − θs1 and
Πs2 = Πs1 ΠT
a Πd .

Proof: Let Πs1 be a permutation matrix such that ΠT
s1 η̃ has the same relative order as h. Now we
T
prove that ΠT
s2 η̃ has the same relative order as −h. Utilizing ha = −hd = −Πd h and Πd Πd = I,
T
T
T
T
we obtain ΠT
d ha = −h. Note that Πs2 η̃ = Πd Πa Πs1 η̃ . Because Πs1 η̃ has the same relative order as
T
T
h, Πa ΠT
s1 η̃ has the same relative order as Πa h = ha , and Πd Πa Πs1 η̃ has the same relative order as

ΠT
d ha = −h.

Next we prove that l(η̃; θ, Π)|θ=θs1 ,Π=Πs1 = l(η̃; θ, Π)|θ=θs2 ,Π=Πs2 holds. Because θs2 = 2c0 − θs1 ,
we have
hi θs1 − τi = hi (θs1 − c0 ),

hi θs2 − τi = −hi (θs1 − c0 ).

(15)

By examining l(η̃; θ, Π) (13) and utilizing ha = −hd , the second term of l(η̃; θ, Π)|θ=θs1 ,Π=Πs1 is
equal to that of l(η̃; θ, Π)|θ=θs2 ,Π=Πs2 . For the first term, note that given θs1 and θs2 , the corresponding
s1 and s2 in (13) can be viewed as evaluating at h and −h according to (15), respectively. Because
ha = −hd , we can conclude that s1 is a permutated version of s2 . The first term of (13) can be
T
T
T
T
T
expressed as either (ΠT
s1 η̃) s1 or (Πs2 η̃) s2 . Because (Πs1 η̃) and s1 have the same relative order as
T
T
T
T
T
h, and (ΠT
s2 η̃) and s2 have the same relative order as −h, one has (Πs1 η̃) s1 = (Πs2 η̃) s2 . Thus

l(η̃; θ, Π)|θ=θs1 ,Π=Πs1 = l(η̃; θ, Π)|θ=θs2 ,Π=Πs2 .

Now an example is presented to substantiate the above proposition. Let c0 = 0.5, the true value

DRAFT

10

θ0 = 1, h = [2, −1, −2, 1]T , η = [η1 , η2 , η3 , η4 ]T and Π0 = [0 0 1 0; 0 1 0 0; 0 0 0 1; 1 0 0 0]. Then
η̃ = [η3 , η2 , η4 , η1 ]T , ha = [−2, −1, 1, 2]T , hd = [2, 1, −1, −2]T and ha = −hd . We can conclude that
l(η̃; θ, Π)|θ=1,Π=Π0 = l(η̃; θ, Π)|θ=0,Π=Π0 , where Π0 = [0 0 0 1; 1 0 0 0; 0 0 1 0; 0 1 0 0].

In addition, given |c0 | ≥ ∆, only one of {θs1 , θs2 } lies in the interval [−∆, ∆], and the model is
identifiable. In the following, a method to select good initial points for the alternating maximization
algorithm is provided for the general case.

E. Good initial points
For alternating maximization algorithms dealing with nonconvex optimization problems, an initial point
is important for the algorithm to converge to the global optimum. In the following text, we provide good
initial points for the alternating maximization algorithm. The key idea is to obtain a coarse estimate of
θ via matching the expected and actual number of ones in observations, and utilizing the orthogonal

property of permutation matrix.
Suppose that the number of measurements K is large. Consequently, as the number of measurements
tends to infinity, the law of large numbers (LLN) implies

p

ηi −→ q0 + (1 − q0 − q1 )Fw ((hi θ − τi )/σw ) ,

(16)

p

where −→ denotes convergence in probability. Given θ ∈ [−∆, ∆], −|hi |∆ − τi ≤ hi θ − τi ≤ |hi |∆ −
τi . In the following text, we only deal with q0 + q1 < 1 case. The case that q0 + q1 > 1 is very

similar and is omitted here. Define l =
max
i∈[1,··· ,N ]

min
i∈[1,··· ,N ]

(q0 + (1 − q0 − q1 )Fw ((−|hi |∆ − τi )/σw )) and u =

(q0 + (1 − q0 − q1 )Fw ((|hi |∆ − τi )/σw )). Then ηi should satisfy l ≤ ηi ≤ u. Let Il,u (η̃i )

denotes the projection of η̃i onto the interval [l, u]. Note that this projection operation is needed because
(16) is valid in the limit as K goes to infinity. From (16) one obtains
p

m , σw Fw−1 ((Il,u (η̃) − q0 1N )/(1 − q0 − q1 )) −→Π(hθ − τ ).

Utilizing ΠΠT = I yields
p

mT m−→hT hθ2 − 2τ T hθ + τ T τ ,

(17)

which is a quadratic equation in θ. Accordingly, using the asymptotic properties of mT m, one obtains
(18) via inverting (17)
θ1,2

τ Th
= T ±
h h

r

mT m − τ T τ
τ Th
+ ( T )2 .
T
h h
h h

(18)

DRAFT

11

The above two solutions can be used for the alternating maximization algorithm as initial points. Finally,
the optimum with larger likelihood is chosen as ML estimator. In Section VI, to provide a fair comparison
of the alternating maximization algorithm with good initial points, −∆ and ∆ are used as two initial
points, and we choose as ML estimator the solution whose likelihood is larger.
The result of (18) is consistent with that of Proposition 3. Given that the conditions in Proposition 3
are satisfied, and substituting τ = c0 h into (18), the solutions are θ1 = θ and θ2 = 2c0 − θ.
V. D ETECTION FROM UNLABELED DATA
In this section, we study the detection problem from unlabeled data. The GLRT detector is studied
separately in the cases of known and unknown θ. In addition, we investigate the permutation matrix
recovery probability.
A. Detection with permuted data and known θ
In the case of known θ, the GLRT can be formulated as
T2 (η̃) = max l(η̃; θ, Π) − max l(η̃; 0, Π) > γ.
Π∈PN

Π∈PN

(19)

As shown in Proposition 1, the ML estimate of the permutation matrix Π corresponding to the first term in
(19) will reorder the elements of η̃ to have the same relative order as the elements of (1−q0 −q1 )(hθ−τ ).
Similarly, for the ML estimation problem corresponding to the second term in (19), we reorder the
elements of η̃ to have the same order as that of −(1 − q0 − q1 )τ .
B. Detection with permuted data and unknown θ
For the unknown θ and unknown Π case, a GLRT is used to decide H1 if
T3 (η̃) = max l(η̃; θ, Π) − max l(η̃; 0, Π) > γ.
θ,Π∈PN

Π∈PN

(20)

Algorithm 1 for joint estimation of θ and Π has been described in section IV-C1. The performance of
the GLRT (20) will be evaluated in the Algorithm 1 for joint estimation of θ and Π is necessary for the
first term and has been described in section IV-C1; the second term is as in section V-A.
C. Approximations on permutation matrix recovery probability
In this subsection, we investigate the permutation matrix recovery probability problem. Since errors
in permutation matrix recovery are more likely to happen in the relatively indistinguishable cases, the
performances in terms of signal detection or estimation tasks may not be closely related to the recovery of

DRAFT

12

permutation matrix. However, it is meaningful to extract the accurate timestamp information or sensors’
identity information which corresponds to recovery of permutation matrix, as presented in the following.
It is difficult to obtain the permutation matrix recovery probability in the case of unknown θ. Instead,
we assume that θ is known, and analyze the permutation matrix recovery probability in terms of the
recovery algorithm provided in Proposition 1. Without loss of generality, we also assume that q0 + q1 < 1
in the following analysis. The case that q0 + q1 > 1 is similar and is omitted here.
First, let pi be ordered such that p(1) > p(2) > · · · > p(K) . From (3) we have (hi θ − τi )(1) >
(hi θ − τi )(2) > · · · > (hi θ − τi )(K) . Provided q0 + q1 < 1, the elements of η̃ should be reordered

according to the order of the elements of hθ − τ in Proposition 1. Therefore the permutation matrix will
be correctly recovered if and only if η(1) > η(2) > · · · > η(K) . Note that the subscripts of (hi θ − τi )(·)
and η(·) also correspond to the order of pi , instead of the order of hi θ − τi or ηi .
Define Ei as the event such that η(i) > η(i+1) and Ēi as the corresponding complement event of Ei ,
namely, η(i) ≤ η(i+1) . The probability that permutation matrix is recovered correctly can be written as
Pr(Π̂ML = Π0 ) = Pr(η(1) > · · · > η(K) ) = Pr

K−1
\

!
Ei

i=1

=1 − Pr

K−1
[

!
≥1−

Ēi

K−1
X

i=1

K−1
S

Pr(Ēi ),

i=1

K−1
P
where union bound Pr
Ēi ≤
Pr(Ēi ) is utilized in (21). From (3), we have uij ∼ B(1, pi )
i=1
i=1
PN
and N ηi = j=1 uij ∼ B(N, pi ). When N is large, the De Moivre-Laplace theorem [30] implies that



the distribution of ηi can be approximated by N (pi , pi (1 − pi )/N ). As a consequence, η(i) − η(i+1) is
approximately distributed as N (p(i) − p(i+1) , p(i) (1 − p(i) )/N + p(i+1) (1 − p(i+1) )/N ), and
Pr(Π̂ML = Π0 ) ≥ 1 −

K−1
X

Pr(Ēi ) = 1 −

i=1

≈1 −

K−1
X


Φ q

K−1
X

Pr(η(i) − η(i+1) ≤ 0)

i=1

√
−(p(i) − p(i+1) ) N

p(i) (1 − p(i) ) + p(i+1) (1 − p(i+1) )
i=1
 √ 
≥1 − (K − 1)Φ −t N
≈1 − (K − 1) √

min

i=1,··· ,K−1

√



1
2
√ e−t N/2
2πt N

1
t2
1
=1 − √ eln(K−1)−ln t− 2 ln N − 2 N , Pr(K, N ),
2π

where t =



vi
,
p(i) (1−p(i) )+p(i+1) (1−p(i+1) )

(21)

vi = p(i) − p(i+1) and the approximation Φ(−x) ≈

DRAFT

13

2

x
√ 1 e− 2
2πx

(x  0) is utilized.

Utilizing p(i) (1 − p(i) ) + p(i+1) (1 − p(i+1) ) ≤ 1/2, we define t̃ satisfying
√
2
t̃ =
min
vi ≤
t.
i=1,··· ,K−1
2

(22)

We conjecture that t̃ is on the order of K −α , i.e., t̃ = O(K −α ), which means that there exists constant
ct such that
t̃ ≈ ct K −α .

(23)

In the following text, we show that we can construct h such that t̃ = O(K −1 ) and t̃ = O(K −2 ).
According to (22) and (23), the approximation Pr(K, N ) (21) can be further simplified and relaxed as
1
1
2
f
Pr(K,
N ) = 1 − √ eln(K−1)−ln t̃− 2 ln N −t̃ N
2 π
c2
1
1
t
≈ 1 − √ e(1+α) ln K− 2 ln N − K 2α N .
2 πct

(24)
(25)

2

t
N of (25) must be far less than 0 for the recovery
From (25), the exponent (1 + α) ln K − 12 ln N − Kc2α

of permutation matrix. Given N is large, the term − 21 ln N is small compared to N . Thus (1 + α) ln K −
c2t
K 2α N

< 0 will ensure that the permutation matrix can be recovered in high probability. Simplifying

(1 + α) ln K −

c2t
K 2α N

< 0 yields
N>

(1 + α) 2α
K ln K.
c2t

(26)

The following cases are examples to illustrate t̃ = O(K −α ). For simplicity, we assume 1 − q0 − q1 > 0,
τ = ch(c < θ) and
a , (θ − c)/σw > 0.

1) t̃ = O(K −1 ): Let h be the shape of a ramp signal such that hi = u −

(27)
(u−l)(i−1)
(u
K−1

> |l|), and

2 ). Then the ordered sequence p
wij ∼ N (0, σw
(i) = pi , and t̃ can be approximated as

t̃ =

min

i=1,··· ,K−1

pi − pi+1

a(1 − q0 − q1 )(u − l)
min
fw (aξi )
i=1,··· ,K−1
K −1
a(1 − q0 − q1 )(u − l)fw (au)
≈
K −1

=

(28)

≈ct K −1 ,

DRAFT

14

where mean value theorem is utilized for ξi ∈ (hi+1 , hi ), ξ1 ≈ h1 = u is utilized when K is large, and
ct = a(1 − q0 − q1 )(u − l)fw (au).

(29)

Therefore t̃ can be reshaped in the form of (23).
2) t̃ = O(K −2 ): Let hi be independently drawn from the same distribution of wij /σw . The CDF of
pi is
Fpi (x) = Pr(pi ≤ x)
= Pr(q0 + (1 − q0 − q1 )Fw (ahi ) ≤ x)



1 −1
x − q0
= Pr hi ≤ Fw
a
1 − q0 − q1



1 −1
x − q0
= Fw
F
.
a w
1 − q0 − q1

(30)

In this case, we conjecture that t̃ = O(K g(a) ), where g(a) is a function of a, and the numerical results
under different a are shown in Fig. 2. In addition, the case in which h is the shape of a sinusoidal signal
is also presented in Fig. 3.

10 -2

10 -2

(a)

(b)

10 -6

10 -6

t~

10 -4

t~

10 -4

10 -8

10

10 -8

-10

10 -10

a=0.5
10 -12
10 1

10 2

a=1
10 3

10 4

10 -12
10 1

10 5

K

10 2

10 0

10 3

10 4

10 0

(d)

10 -5

10 -20

10 -10

10 -40

t~

t~

(c)

10 -15

10 -60

a=4

a=2
10 -20
10 1

10 5

K

t~, realization one
t~, realization two
t~, mean of 1000 realizations

10 2

10 3

K

10 4

10 5

10 -80
10 1

10 2

10 3

10 4

10 5

K

Fig. 2: The relationship of t̃ and K under different a. Note that q0 = q1 = 0, hi ∼ N (0, 1) and
2 ).
wij ∼ N (0, σw
Now we prove that t̃ = O(K −2 ) under certain conditions. Given that hi and wij /σw are i.i.d. random
variables and a = 1, the CDF Fpi (x) = (x − q0 )/(1 − q0 − q1 ), and the PDF of pi is


1

1−q0 −q1 , q0 ≤ x ≤ 1 − q1 ,
fpi (x) =

 0,
otherwise.

(31)

DRAFT

15

10 -2

10 -2

(a)

(b)

-4

10 -4

10 -6

10 -6

t~

t~

10

10

-8

10 -8

10 -10

10 -10

a=0.5
10

a=1

-12

10 1

10 2

10 3

10 4

10 -12
10 1

10 5

K

10 2

10 -2

10 3

10 4

10 0

(d)

(c)
10

10 5

K

t~, realization one
t~, realization two
t~, mean of 1000 realizations

-4

10 -5

t~

t~

10 -6
10 -8

10 -10

10 -10

a=4

a=2
10 -12
10 1

10 2

10 3

10 4

10 -15
10 1

10 5

10 2

K

10 3

10 4

10 5

K

Fig. 3: The relationship of t̃ and K under different a. Note that q0 = q1 = 0, hi = sin(2πxi ),
2 ).
xi ∼ U(0, 1) and wij ∼ N (0, σw

Then the variates p(1) , p(2) · · · , p(K) are distributed as K descending ordered statistics from an uniform
(q0 , 1 − q1 ) parent. For x ≤ (1 − q0 − q1 )/(K − 1), the CDF of t̃ can be derived as [31] (page 135,

equation (6.4.3))

Ft̃ (x) =Pr


min

i=1,··· ,K−1

vi ≤ x

=1 − Pr(v1 > x, v2 > x, · · · , vK−1 > x)


(K − 1)x K
.
=1 − 1 −
1 − q0 − q1

For x ≥ (1 − q0 − q1 )/(K − 1), Ft̃ (x) = 1. Then the PDF of t̃ is

h
iK−1

 K(K−1) 1 − (K−1)x
, 0≤x≤
1−q0 −q1
1−q0 −q1
ft̃ (x) =

 0,
otherwise.

1−q0 −q1
K−1 ,

(32)

(33)

The expectation of t̃ is
1

Z
Et̃ [t̃] =

0

Z
xft̃ (x) dx =

0

1−q0 −q1
K−1

t̃ft̃ (x) dx


Z 1−q0 −q1 
K−1
K(K − 1)
(K − 1)x K−1
=
x 1−
dx
1 − q0 − q1 0
1 − q0 − q1
1 − q0 − q1
=
.
K2 − 1

(34)

DRAFT

16

Hence the probability that t̃ falls into [c1 /K 2 , c2 /K 2 ] is
Pr(c1 /K 2 ≤ t̃ ≤ c2 /K 2 ) = Ft̃ (c2 /K 2 ) − Ft̃ (c1 /K 2 )
K 
K

c1 (K − 1)
c2 (K − 1)
= 1−
− 1−
.
(1 − q0 − q1 )K 2
(1 − q0 − q1 )K 2

(35)

0

When K is large, (K−1)/K ≈ 1 and (1−1/(c0 K))c K ≈ 1/e(c0 > 0). Equation (35) can be approximated
as
Pr(c1 /K 2 ≤ t̃ ≤ c2 /K 2 ) ≈ e

c

− 1−q 1−q
0

c

1

− 1−q 2−q

−e

0

1

.

(36)

Provided that q0 = q1 = 0, when c1 = 0.1 and c2 = 10, Pr(0.1/K 2 ≤ t̃ ≤ 10/K 2 ) ≈ 0.94; when
c1 = 0.01 and c2 = 100, Pr(0.01/K 2 ≤ t̃ ≤ 100/K 2 ) ≈ 0.99. It can be seen that t̃ falls near the order

of magnitude of K −2 with high probabilities. Thus it is reasonable that t̃ = O(K −2 ).
According to the definition of pi (3), equations (22) and (23), ct ∝ 1 − q0 − q1 . From (26), the number
of quantizers Nreq required for permutation matrix recovery probability is
Nreq ∝ 1/(1 − q0 − q1 )2 .

(37)

From (37), one can conclude that the number of quantizers for permutation matrix recovery with high
probability is 1/(1 − q0 − q1 )2 times that of unflipped case where q0 = q1 = 0.

VI.

NUMERICAL

S IMULATIONS

In this section, numerical experiments are conducted to evaluate the theoretical results. For simplicity,
2 ).
the distribution of noise wij is selected as the Gaussian distribution N (0, σw

A. Parameter estimation
For the first two experiments, we evaluate the performance of ML estimators proposed in section IV.
2 = 1, ∆ = 2, q = 0.05, q = 0.05 and the tolerance
Parameters are set as follows: K = 20, θ = 1, σw
0
1

parameter  in Algorithm 1 is 10−7 . The number of Monte Carlo trials is 5000.
For the first experiment, the MSE performance of Algorithm 2 is evaluated in Fig. 4. We let τ =
0.5h, which is a special case mentioned in Proposition 2. The coefficients h is equispaced with h =
[−1.50, −1.29, −1.08, · · · , 2.50]T , which corresponds to a ramp signal. It can be seen that h does not

satisfy the condition in Proposition 3, thus the model may be identifiable. It can be seen that the ML
estimator from labeled data always works well. Given limited number of quantizers, there is an obvious

DRAFT

17

gap between the MSEs of two estimators. As the number of quantizers increases, the performance of the
estimator from unlabeled data approaches that from labeled data.

9

#10 -3
Labeled data
Unlabeled data
CRLB

8
7

MSE

6
5
4
3
2
1
0
20

40

60

80

100

120

140

160

180

200

N

Fig. 4: MSE of θ vs. number of quantizers for the ML estimators from labeled and unlabeled data,
compared with the CRLB (8) for ramp signal.

For the second experiment, the MSE performance of Algorithm 1 (for the general case) is evaluated In
Fig. 5. The elements of the vector h describe the shape of a sinusoidal signal such that hi = sin(2πxi ),
where xi is drawn independently and randomly from the uniform distribution U(0, 1) and then sorted
in ascending order. The elements of the vector τ is drawn independently and randomly from the
uniform distribution U(−∆, ∆). It can be seen that when N < 80, good initial points improve the
MSE performance of the alternating maximization algorithm from unlabeled data. As N increases to 80,
the MSE performances of both unlabeled ML estimators approach a common level which is larger than
that achieved by the labeled data. Finally, the MSEs of both estimators from unlabeled data approach to
that from labeled data around N = 3 × 104 .

B. Signal detection
In Fig. 6, the relationship between PD and the number of quantizers N is employed. Parameters are
2 = 9 and P
consistent with the first experiment, except that σw
F A = 0.05.

In subgraph (a), h and τ are the same as those in the first experiment. It can be seen that the number
of quantizers has a significant effect on the detection probability. As N increases, the performance of
all the detectors improves, and the detection performance of the unlabeled GLRT approaches to that of
labeled GLRT. In subgraph (b), h and τ are the same as those in the second experiment, and similar
phenomena are observed. It seems that in this case little is gained by good initialization.

DRAFT

18

10 -1

0.1

MSE

10 -2

0.05

10 -3

MSE

0
20

40

60

80

N
10

-4

Labeled data
Unlabeled data
Unlabeled data, good initial points
CRLB

10 -5

10 -6

10 2

10 3

10 4

10 5

N

Fig. 5: MSE of θ vs. number of quantizers for the three ML estimators from labeled data, unlabeled
data via initial points ±∆ and unlabeled data via good initial points (18), compared with the CRLB (8)
for sinusoidal signal.

1

1

(b)

(a)

0.9

0.8

0.8
0.7

PD

PD

0.6
0.6
0.5
0.4
0.4
NP
labeled GLRT, 3 unknown
reorderd GLRT, permutation unknown
unlabeled GLRT

0.3
0.2
0.1

NP
labeled GLRT, 3 unknown
reorderd GLRT, permutation unknown
unlabeled GLRT
unlabeled GLRT,good initial points

0.2

0
20

40

60

80

N

100

120

140

20

40

60

80

100

120

N

Fig. 6: Pd vs. number of quantizers N for the ramp signal in subgraph (a) and the sinusoidal signal in
subgraph (b).

C. Permutation matrix recovery
In this subsection, the approximations for permutation matrix recovery are verified. Parameters are set
2 = 1. The number of Monte Carlo trials is
as follows: K = 20 , θ = 1.5, ∆ = 2, q0 = 0, q1 = 0 and σw

1000.

First, the relationship of t and t̃ (22) and the conjecture of t̃ (23) are illustrated in three cases.
√
From Fig. 7, one obtains that t can be approximated as 2t̃ in practice. For a ramp signal, h =
√
[−0.800, −0.705, −0.610, · · · , 1.000]T and τ = 0.5h. t ≈ 2ce /K where ct = ce = 0.4355 is evaluated
√
via (29). Because of the gap between t and 2ce /K , we use linear regression to fit t and obtain
cea = 0.6717, which is much more accurate than ce and will be utilized later to predict the number of

quantizers for permutation matrix recovery. For random generated h, h is drawn from standard normal
distribution and τ = 0.5h. It can be seen that t can be approximated by 1/K 2 . For a sinusoidal signal,
h and τ are drawn in the same way of the second experiment. We use linear regression and obtain

DRAFT

19

t ≈ 0.71/K 2.23 ≈

√

2t̃ =

√

2ct,s /K αt,s , ct,s = 0.5020 and αt,s = 2.23.

100

10-2

t

10-4

10-6

10-8

10-10

10-12
101

102

103

104

105

Fig. 7: The relationship of t and K , including equispaced, randomly generated and sinusoidal h cases.

Next, the empirical permutation recovery probability Pr(Π̂ML = Π0 ) versus N or K are presented in
f
Fig. 8, and the theoretical approximations Pr(K, N ) (21) and Pr(K,
N ) (24) are plotted for comparison.
In subgraph (a), (b) and (c), we set K = 20. While in subgraph (d), we set N = 104 . All h are drawn
in the same way as the second experiment. We also evaluate the empirical permutation matrix recovery
probability in the case of unknown θ, which has negligible difference compared to that in the known θ
case.
In subgraph(a), it can be seen that the permutation matrix of the ramp signal can be recovered with
high probability given N ≥ 5000. From N >
can conclude that N >

2
2
0.43552 K ln K|K=20

1+α 2α
ln K
c2t K

(26) where ct = ce = 0.4355 and α = 1, one

≈ 12636, which is more than twice of 5000. Utilizing the

fitted parameter cea , one obtain a more accurate result that N >

2
2
0.67172 K ln K|K=20

≈ 5312 ensures

permutation matrix recovery with high probability. For random h, N > 3K 4 ln K|K=20 ≈ 1.438 × 106
ensures recovery with high probability, which is not accurate enough, as subgraph(b) shows that N ≈ 105
is enough for recovery of permutation matrix. In subgraph(c), it is shown that N ≈ 106 is enough for
recovery of permutation matrix, which is also inaccurate compared to the fitted results of the sinusoidal
signal N >

3.23
4.46 ln K|
K=20
0.50202 K

≈ 2.437 × 107 . The numerical results show that the theoretical

bound Pr(K, N ) is accurate in predicting N with high probability in permutation matrix recovery, which
f
demonstrates that Pr(K,
N ) may be too conservative in predicting the number of quantizers ensuring
perfect permutation matrix recovery. In subgraph(d), 10000 = N >

2
2
0.67172 K ln K|K=26

≈ 9763, thus

K ≤ 26 will ensure permutation matrix recovery with high probability, which is consistent with the

numerical results.
In Fig. 9, the relationship of flipping probabilities (q0 , q1 ) and number of quantizers Nreq (37) required

DRAFT

20

1

Permutation Recovery Probability

Permutation Recovery Probability

1

(a)
0.8

0.6

0.4

0.2

0
10 3

0.2

10 4

10 5

10 6

N

1

Permutation Recovery Probability

Permutation Recovery Probability

0.4

experiment, 3 known, reordering
experiment, 3 unknown, AM
theorical lower bound, Pr(K; N )
f
theorical lower bound, Pr(K;
N)

1

(c)
0.8

0.6

0.4

0.2

0
10 3

0.6

0
10 3

10 4

N

(b)
0.8

(d)

0.8

0.6

0.4

0.2

0
10 4

10 5

10 6

20

40

60

80

100

K

N

Fig. 8: Pr(Π̂ML = Π0 ) vs. N or K for the ramp signal in subgraph (a)(d), random generated h in
f
subgraph (b) and the sinusoidal signal in subgraph (c). Pr(K, N ) and Pr(K,
N ) are evaluated via (21)
and (24), respectively.

for permutation matrix recovery with high probability is verified. Parameters are the same as those in
Fig. 8-(a) except for (q0 , q1 ). We use the result of the experiment in which q0 = q1 = 0 to predict those
in which q0 = q1 = 0.05, q0 = q1 = 0.1 and q0 = q1 = 0.15, and plot the experimental results for
comparison. It can be seen that the predictions are basically consistent with the experimental results,
which verifies (37).

1

Permutation Recovery Probability

0.9
0.8
0.7

experiment, q0 = q1 = 0
prediction, q0 = q1 = 0:05
experiment, q0 = q1 = 0:05
prediction, q0 = q1 = 0:1
experiment, q0 = q1 = 0:1
prediction, q0 = q1 = 0:15
experiment, q0 = q1 = 0:15

0.6
0.5
0.4
0.3
0.2
0.1
0
10 3

10 4

N

Fig. 9: Pr(Π̂ML = Π0 ) vs. number of quantizers N for ramp signal under different flipping
probabilities (q0 , q1 ).

DRAFT

21

VII.

CONCLUSION

We study a scale parameter estimation and signal detection problem from unlabeled quantized data for
a canonical (known signal shape) sensing model. A sufficient condition under which the signal amplitude
estimation problem can be solved efficiently is provided. It is also shown that in some settings the
model can even be unidentifiable. Given that the number of quantizers is limited, the performance of the
unlabeled estimator via reordering and alternating maximization algorithms is good, although there is a
gap between the performances of labeled and unlabeled ML estimators. In addition, good initial points
are provided to improve the performance of an alternating maximization algorithm for general estimation
problems. As the number of quantizers increases, the performance of the unlabeled estimator approaches
that of the labeled estimator due to the recovery of permutation matrix.
Furthermore, the performance of GLRT detector under unlabeled samples is evaluated, and numerical
results show that the performance degradation of the GLRT detector under unlabeled samples is significant
in noisy environments, compared to the GLRT detector with labeled samples given that the number
of quantizers is small. As the number of quantizers increases, the performance of the GLRT under
unlabeled samples approaches that of the GLRT detector under labeled samples. The explicit approximated
permutation matrix recovery probability predicts that in order to find the true label of K time indexes,
the number of quantizers N should be on the order of K 2α log K , where α is a constant depending on
the signal shape and the distribution of noise.

R EFERENCES
[1] P. Pradhan, K. Nagananda, P. Venkitasubramaniam, S. Kishore and R. S. Blum, “GPS spoofing attack characterization and
detection in smart grids,” Communications and Network Security (CNS), 2016 IEEE Conference on, pp. 391-395, 2016.
[2] T. E. Humphreys, B. M. Ledvina, M. L. Psiaki, B. W. O’hanlon and P. M. Kintner, “Assessing the spoofing threat:
development of a portable GPS civilian spoofer,” in Proc. Int. Tech. Meet. Satellite Div. The Ins. Navigation, pp. 2314-2325,
2008.
[3] Q. Zeng, H. Li and L. Qian,“GPS spoofing attack on time synchronization in wireless networks and detection scheme
design,” in MILCOM’12, pp. 1-5, 2012.
[4] Z. Zhang, S. Gong, A. D. Dimitrovski and H. Li, “Time synchronization attack in smart grid: Impact and analysis,” IEEE
Trans. Smart Grid, vol. 4, no. 1, pp. 87-98, 2013.
[5] S. Challa, R. J. Evans and X. Wang, “A Bayesian solution and its approximations to out-of-sequence measurement problems,”
Information Fusion, vol. 4, no. 3, pp. 185-199, 2003.
[6] L. Schenato, “Optimal estimation in networked control systems subject to random delay and packet drop,” IEEE Trans.
Autom. Control, vol. 53, no. 5, pp. 1311-1317, 2008.
[7] L. M. Millefiori, P. Braca, K. Bryan and P. Willett, “Adaptive filtering of imprecisely time-stamped measurements with
application to AIS networks,” in Proc. of the 18th Intern. Conf. on Inform. Fusion (FUSION), pp. 359-365, 2015.

DRAFT

22

[8] V. Emiya, A. Bonnefoy, L. Daudet and R. Gribonval, “Compressed sensing with unknown sensor permutation,” ICASSP,
pp. 1040-1044, 2014.
[9] J. Unnikrishnan, S. Haghighatshoar and M. Vetterli, “Unlabeled sensing: Solving a linear system with unordered
measurements,” Communication, Control, and Computing (Allerton), 2015 53rd Annual Allerton Conference on. IEEE,
pp. 786-793, 2015.
[10] A. Pananjady, M. J. Wainwright and T. A. Courtade, “Linear regression with an unknown permutation: statistical and
computational limits,”Communication, Control, and Computing (Allerton), 2015 53rd Annual Allerton Conference on. IEEE,
pp. 417-424, 2016.
[11] A. Pananjady, M. J. Wainwright and T. A. Courtade, “Denoising linear models with permutated data,”
http://arxiv.org/abs/1704.07461, 2017.
[12] A. Abid, A. Poon and J. Zou, “Linear regression with shuffled labels,” http://arxiv.org/abs/1705.01342, 2017.
[13] S. Haghighatshoar and G. Caire, “Signal recovery from unlabeled samples,” ISIT, pp. 451-455, 2017.
[14] L. Keller, M. J. Siavoshani, C. Fragouli and K. Argyraki, “Identity aware sensor networks,” Proceedings - IEEE INFOCOM,
pp. 2177-2185, 2009.
[15] S. Marano, V. Matta, P. Willett, P. Braca and R. S. Blum, “Hypothesis testing in the presence of Maxwell’s daemon: Signal
detection by unlabeled observations,” ICASSP, pp. 3286-3290, 2017.
[16] P. Braca, S. Marano, V. Matta, P. Willett, “Asymptotic efficiency of the PHD in multitarget/multisensor estimation,” IEEE
Journal of Selected Topics in Signal Processing, vol. 7, no. 3, pp. 553-564, 2013.
[17] J. Zhu, H. Cao, C. Song and Z. Xu, “Parameter estimation via unlabeled sensing using distributed sensors,” IEEE commun.
lett., vol. 21, no. 10, pp. 2130-2133, 2017.
[18] J. Unnikrishnan, S. Haghighatshoar, M. Vetterli, “Unlabeled sensing with random linear measurements,” avaliable at
https://arxiv.org/pdf/1512.00115.pdf, 2015.
[19] A. Pananjady, M. J. Wainwright and T. A. Courtade, “Linear regression with an unknown permutation: statistical and
computational limits,” available at https://arxiv.org/abs/1608.02902, 2016.
[20] S.

Haghighatshoar

and

G.

Caire,

“Signal

recovery

from

unlabeled

samples,”

avaliable

at

https://arxiv.org/pdf/1701.08701.pdf, 2017.
[21] O. Ozdemir and P. K. Varshney, “Channel aware target location with quantized data in wireless sensor networks,” IEEE
Trans. Signal Process., vol. 57, pp. 1190-1202, 2009.
[22] H. C. Papadopoulos, G. W. Wornell and A. V. Oppenheim, “Sequential signal encoding from noisy measurements using
quantizers with dynamic bias control,” IEEE Trans. Inf. Theory, vol. 47, no. 3, pp. 978-1002, Mar. 2001.
[23] A. Ribeiro and G. B. Giannakis, “Bandwidth-constrained distributed estimation for wireless sensor networks-part I: Gaussian
case,”IEEE Trans. Signal Process., vol. 54, no. 3, pp. 1131-1143, 2006.
[24] A. Ribeiro and G. B. Giannakis, “Bandwidth-constrained distributed estimation for wireless sensor networks-part II:
unknown probability density function,”IEEE Trans. Signal Process., vol. 54, no. 7, pp. 2784-2796, 2006.
[25] S. Boyd and L. Vandenberghe, Convex Optimization, Cambridge University Press, 2004.
[26] S. M. Kay, Fundamentals of Statistical Signal Processing, Volume I: Estimation Theory, Englewood Cliffs, NJ: Prentice
Hall, 1993.
[27] S. M. Kay, Fundamentals of Statistical Signal Processing, Volume II: Detection Theory, Englewood Cliffs, NJ: Prentice
Hall, 1993.
[28] S. A. Kassam, Signal Detection in Non-Gaussian Noise, World Publishing Corp., 1992.
[29] D. P. Bertsekas, Nonlinear Programming, 2nd ed., Athena Scientific, Belmont, MA, 1999.

DRAFT

23

[30] A. Papoulis and S. U. Pillai, Probability, Random Variables, and Stochastic Processes, Fourth Edition, 2002.
[31] H. A. David and H. N. Nagaraja, Order Statistics: Third Edition, John Wiley, pp. 133-135, 2003.
[32] E. L. Lehmann, Elements of large-sample theory, Springer-Verlag New York, inc., pp. 456-457, 1999.

DRAFT

