Signal Amplitude Estimation and Detection from Unlabeled Binary Quantized Samples

TL;DR

This paper develops methods for estimating signal amplitude and detecting signals from unlabeled binary quantized samples, addressing unknown ordering and providing algorithms with theoretical guarantees.

Contribution

It introduces polynomial-time ML estimators for permutation and amplitude, analyzes model identifiability, and derives probability bounds for permutation recovery.

Findings

ML estimators effectively recover signal amplitude and permutation.

Explicit probability expressions relate sample size and quantizer count.

Numerical results confirm theoretical accuracy and robustness.

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.

Figures

Equations

$$b_{ij}=Q_i(h_i\theta +w_{ij}),i=1,\cdots ,K,j=1,\cdots ,N,$$

$$\{\begin{array}{ll}& {\mathcal{H}}_0:b_{ij}=Q_i(w_{ij}),i=1,\cdots ,K,j=1,\cdots ,N,\\ & {\mathcal{H}}_1:b_{ij}=Q_i(h_i\theta +w_{ij}),i=1,\cdots ,K,j=1,\cdots ,N,\end{array}$$

$$\tilde{\mathbf{U}}=\mathbf{\Pi }\mathbf{U},$$

$$\mathrm{Pr}(u_{ij}=1)$$

$$\mathrm{Pr}(u_{ij}=0)$$

$$p(\mathbf{U};\theta )=\prod _{i=1}^K\prod _{j=1}^N\mathrm{Pr}(u_{ij}=1{)}^{u_{ij}}\mathrm{Pr}(u_{ij}=0{)}^{(1-u_{ij})}.$$

$${\eta }_i=\sum _{j=1}^N{u}_{ij}/N.$$

$$l(\mathbf{\eta };\theta )=N\sum _{i=1}^K({\eta }_i\log p_i+(1-{\eta }_i)\log (1-p_i)),$$

$$I(\theta )$$

$$=\frac{N(1-q_0-q_1{)}^2}{{\sigma }_w^2}\sum _{i=1}^K\frac{h_i^2{f}_w^2\left(\frac{h_i\theta -{\tau }_i}{{\sigma }_w}\right)}{p_i(1-p_i)},$$

$$\mathrm{CRLB}(\theta )=1/I(\theta ),$$

$$T_1(\mathbf{\eta })=\underset{\theta \in [-\Delta ,\Delta ]}{\max }l(\mathbf{\eta };\theta )-l(\mathbf{\eta };0)>\gamma ,$$

$$\begin{array}{rl}& p(\tilde{\mathbf{U}};\theta ,\mathbf{\Pi })=\prod _{m=1}^K\prod _{j=1}^N\mathrm{Pr}({\tilde{u}}_{mj}=1{)}^{{\tilde{u}}_{mj}}\mathrm{Pr}({\tilde{u}}_{mj}=0{)}^{(1-{\tilde{u}}_{mj})}\\ & =\prod _{i=1}^K\prod _{j=1}^N\mathrm{Pr}({\tilde{u}}_{\pi (i)j}=1{)}^{{\tilde{u}}_{\pi (i)j}}\mathrm{Pr}({\tilde{u}}_{\pi (i)j}=0{)}^{(1-{\tilde{u}}_{\pi (i)j})},\end{array}$$

$$l(\tilde{\mathbf{\eta }};\theta ,\mathbf{\Pi })=N\sum _{i=1}^K\left({\tilde{\eta }}_{\pi (i)}\log p_i+(1-{\tilde{\eta }}_{\pi (i)})\log (1-p_i)\right),$$

$$\underset{\theta \in [\Delta ,\Delta ],\mathbf{\Pi }\in {\mathcal{P}}_K}{\max }l(\tilde{\mathbf{\eta }};\theta ,\mathbf{\Pi }),$$

$$l(\tilde{\mathbf{\eta }};\theta ,\mathbf{\Pi })=K\sum _{i=1}^N{\tilde{\eta }}_{\pi (i)}{s}_i+K\sum _{i=1}^N\log (1-p_i),$$

$$l(\tilde{\mathbf{\eta }};{\hat{\theta }}_t,{\hat{\mathbf{\Pi }}}_t)\ge l(\tilde{\mathbf{\eta }};{\hat{\theta }}_{t-1},{\hat{\mathbf{\Pi }}}_{t-1}).$$

$$h_i{\theta }_{s1}-{\tau }_i=h_i({\theta }_{s1}-c_0),h_i{\theta }_{s2}-{\tau }_i=-h_i({\theta }_{s1}-c_0).$$

$${\eta }_i\stackrel{\mathrm{p}}{⟶}{q}_0+(1-q_0-q_1)F_w\left((h_i\theta -{\tau }_i)/{\sigma }_w\right),$$

$$\mathbf{m}\triangleq {\sigma }_w{F}_w^{-1}\left(({\mathcal{I}}_{l,u}(\tilde{\mathbf{\eta }})-q_0{1}_N)/(1-q_0-q_1)\right)\stackrel{\mathrm{p}}{⟶}\mathbf{\Pi }(\mathbf{h}\theta -\mathbf{\tau }).$$

$${\mathbf{m}}^{\mathrm{T}}\mathbf{m}\stackrel{\mathrm{p}}{⟶}{\mathbf{h}}^{\mathrm{T}}\mathbf{h}{\theta }^2-2{\mathbf{\tau }}^{\mathrm{T}}\mathbf{h}\theta +{\mathbf{\tau }}^{\mathrm{T}}\mathbf{\tau },$$

$${\theta }_{1,2}=\frac{{\mathbf{\tau }}^{\mathrm{T}}\mathbf{h}}{{\mathbf{h}}^{\mathrm{T}}\mathbf{h}}\pm \sqrt{\frac{{\mathbf{m}}^{\mathrm{T}}\mathbf{m}-{\mathbf{\tau }}^{\mathrm{T}}\mathbf{\tau }}{{\mathbf{h}}^{\mathrm{T}}\mathbf{h}}+(\frac{{\mathbf{\tau }}^{\mathrm{T}}\mathbf{h}}{{\mathbf{h}}^{\mathrm{T}}\mathbf{h}}{)}^2}.$$

$$T_2(\tilde{\mathbf{\eta }})=\underset{\mathbf{\Pi }\in {\mathcal{P}}_N}{\max }l(\tilde{\mathbf{\eta }};\theta ,\mathbf{\Pi })-\underset{\mathbf{\Pi }\in {\mathcal{P}}_N}{\max }l(\tilde{\mathbf{\eta }};0,\mathbf{\Pi })>\gamma .$$

$$T_3(\tilde{\mathbf{\eta }})=\underset{\theta ,\mathbf{\Pi }\in {\mathcal{P}}_N}{\max }l(\tilde{\mathbf{\eta }};\theta ,\mathbf{\Pi })-\underset{\mathbf{\Pi }\in {\mathcal{P}}_N}{\max }l(\tilde{\mathbf{\eta }};0,\mathbf{\Pi })>\gamma .$$

$$\begin{array}{rl}& \mathrm{Pr}({\hat{\mathbf{\Pi }}}_{\mathrm{ML}}={\mathbf{\Pi }}_0)=\mathrm{Pr}({\eta }_{(1)}>\cdots >{\eta }_{(K)})=\mathrm{Pr}\left(\bigcap _{i=1}^{K-1}{E}_i\right)\\ =& 1-\mathrm{Pr}\left(\bigcup _{i=1}^{K-1}{\bar{E}}_i\right)\ge 1-\sum _{i=1}^{K-1}\mathrm{Pr}({\bar{E}}_i),\end{array}$$

$$\mathrm{Pr}({\hat{\mathbf{\Pi }}}_{\mathrm{ML}}={\mathbf{\Pi }}_0)\ge 1-\sum _{i=1}^{K-1}\mathrm{Pr}({\bar{E}}_i)=1-\sum _{i=1}^{K-1}\mathrm{Pr}({\eta }_{(i)}-{\eta }_{(i+1)}\le 0)$$

$$\approx$$

$$\ge$$

$$\approx$$

$$=$$

$$\tilde{t}=\underset{i=1,\cdots ,K-1}{\min }{v}_i\le \frac{\sqrt{2}}{2}t.$$

$$\tilde{t}\approx c_t{K}^{-\alpha }.$$

$$\widetilde{\mathrm{Pr}}(K,N)$$

Full text

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.

I Introduction

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 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 ($x_{i}=\theta+w_{i}$) in [17], the model in this paper is a scale parameter estimation problem ($x_{i}=h_{i}\theta+w_{i}$), in which $h_{i},~{}i=1,\cdots,K$ is the shape of a signal, and $\theta$ 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 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\times 1$ vector of ones is ${\mathbf{1}}_{K}$. For an unknown deterministic parameter $\theta$, $\theta_{0}$ denotes its true value. For an unknown permutation matrix $\boldsymbol{\Pi}$, $\boldsymbol{\Pi}_{0}$ denotes its true value. For a random vector $\mathbf{y}$, ${p}({\mathbf{y}};\theta)$ denotes the probability density function (PDF) of $\mathbf{y}$ parameterized by $\theta$, and ${\rm E}_{\mathbf{y}}[\cdot]$ denotes the expectation taken with respect to $\mathbf{y}$. Let ${\mathcal{N}}(\mu,\sigma^{2})$ denote a Gaussian distribution with mean $\mu$ and variance $\sigma^{2}$. Let $\Phi(\cdot)$ and $\varphi(\cdot)$ denote the cumulative distribution function (CDF) and probability density function (PDF) of a standard Gaussian random variable respectively. Let ${\mathcal{U}}(a,b)$ denote an uniform distribution, whose minimum and maximum values are $a$ and $b$. Let ${\mathcal{B}}(N,p)$ denote a binomial distribution, where $N$ and $p$ denote the number of trials and the probability of event, respectively.

II Problem Setup

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 $\theta$ of a $K$ length signal and detect the presence of the signal, as shown in Fig. 1.

The binary quantized samples $b_{ij}$ are obtained via

[TABLE]

and the corresponding hypothesis problem can be formulated as

[TABLE]

where $i$ and $j$ respectively denote one of the $K$ time indexes and one of the $N$ quantizers, $h_{i}$ is the known coefficient characterizing the signal shape, $w_{ij}$ is the i.i.d. noise drawn from the $\sigma_{w}^{2}$-variance distribution whose PDF is $f_{w}(x/\sigma_{w})/\sigma_{w}$ and CDF is $F_{w}(x/\sigma_{w})$, where $f_{w}(x)$ and $F_{w}(x)$ are the corresponding unit-variance PDF and CDF, and $Q_{i}(\cdot)$ implies a binary quantizer which produces $1$ if the argument is larger than a scalar threshold $\tau_{i}$ and [math] otherwise. The thresholds of $N$ quantizers are identical given any time index111Here 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.. We assume that the PDF $f_{w}(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 $f_{w}(w)$ is logconcave, which is often met in practice by, for example, the Gaussian distribution.

The quantized data $\{b_{ij}\}$ are transmitted over a binary channel with flipping probabilities $q_{0}$ and $q_{1}$ which are defined as ${\rm Pr}(u_{ij}=1|b_{ij}=0)=q_{0}$ and ${\rm Pr}(u_{ij}=0|b_{ij}=1)=q_{1}$, where $u_{ij}$ 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 $\{{u}_{ij}\}_{j=1}^{N}$ are transmitted to the FC with permuted time indexes. Accordingly, the FC receives sets of data, say $\{\tilde{u}_{ij}\}_{j=1}^{N}$, in which the time reference (represented by the index $i$) is invalid. Specifically, the FC does not know which time index the data $\{\tilde{u}_{ij}\}_{j=1}^{N}$ belongs to, but knows that $\{\tilde{u}_{ij}\}_{j=1}^{N}$ belongs to one of the $K$ time indexes. Let us introduce the matrix $\mathbf{U}$ whose $(i,j)-$th entry is $u_{ij}$. Then, the unlabeled samples can be collected in a matrix $\widetilde{\mathbf{U}}$, as follows:

[TABLE]

where ${\boldsymbol{\Pi}}\in{\mathbb{R}}^{K\times 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 $\theta$ is constrained to an interval $[-\Delta,\Delta]$, for algorithm and theoretical reasons [22].

III Preliminaries

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

III-A Maximum Likelihood Estimation

The probability mass function (PMF) of $u_{ij}$ can be calculated as

[TABLE]

The PMF of $\mathbf{U}$ is

[TABLE]

Let $\eta_{i}$ denote the fraction of $u_{ij}=1$ in $\{u_{ij}\}_{j=1}^{N}$, i.e.,

[TABLE]

Consequently, the log-likelihood function $l({\boldsymbol{\eta}};{\theta})$ is

[TABLE]

where $p_{i}$ is given in (3). Note that in an error free binary symmetric channel scenario, i.e., $q_{0}=q_{1}=0$ or $q_{0}=q_{1}=1$, the CDF $F_{w}(x)$ is log-concave as it is the integral of a log-concave PDF $f_{w}(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<q_{0}+q_{1}<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 Cramér Rao lower bound (CRLB).

In addition, the Fisher Information (FI) $I(\theta)$ is the expectation of the negative second derivative of the log-likelihood function $l({\boldsymbol{\eta}};{\theta})$ (6) taken with respect to $\theta$, i.e., [28],

[TABLE]

where (7) follows due to ${\rm E}_{\boldsymbol{\eta}}[\eta_{i}/p_{i}-(1-\eta_{i})/(1-p_{i})]=0$. Consequently, the CRLB is

[TABLE]

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

III-B GLRT detection

In the case of known $\theta$, the optimal detector according to the NP criterion is the log-likelihood ratio test [26]. For unknown $\theta$, the GLRT is usually adopted. Although there is no optimality associated with 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 ${\mathcal{H}}_{1}$ if

[TABLE]

where $\gamma$ is a threshold determined by the given false alarm probability $P_{FA}$.

IV Estimation 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 $\theta$.

IV-A Maximum likelihood estimation

Introduce the function $\pi(\cdot)$ such that $m=\pi(i)$ if the permutation matrix ${\boldsymbol{\Pi}}$ in (2) maps the $i$th row of ${\mathbf{U}}$ to the $m$th row of $\tilde{\mathbf{U}}$. The PMF of $\tilde{\mathbf{U}}$ is

[TABLE]

where $({\rm Pr}(\tilde{u}_{ij}=1),{\rm Pr}(\tilde{u}_{ij}=0))$ is the PMF of $\tilde{u}_{ij}$. The corresponding log-likelihood function $l({\tilde{\boldsymbol{\eta}}};{\theta},{\boldsymbol{\Pi}})$ is

[TABLE]

where $\tilde{\eta}_{\pi(i)}=\sum_{j=1}^{N}{\tilde{u}}_{\pi(i)j}/N=\sum_{j=1}^{N}{\tilde{u}}_{mj}/N$. The ML estimation problem can be formulated as

[TABLE]

where ${\mathcal{P}}_{K}$ denotes the set of all possible $K\times K$ permutation matrices.

IV-B Estimation with permuted data and known $\theta$

In this subsection, the permutation matrix recovery problem is studied in the case of known $\theta$. It is shown that the permutation matrix under ML estimation criterion can be recovered efficiently under the ML criterion.

Proposition 1

*Given the ML estimation problem in (12) with known $\theta$, the ML estimate of the permutation matrix ${\boldsymbol{\Pi}}$ will reorder the rows of $\widetilde{\mathbf{U}}$, and equivalently the elements of $\tilde{\boldsymbol{\eta}}$, to have the same relative order as the elements of $(1-q_{0}-q_{1})({\mathbf{h}}\theta-{\boldsymbol{\tau}})$.  *

Proof:

Note that the objective function $l({\tilde{\boldsymbol{\eta}}};{\theta},{\boldsymbol{\Pi}})$ (11) can be decomposed as

[TABLE]

where $s_{i}=\log(p_{i}/(1-p_{i}))$. From (13), the ML estimate of the permutation matrix ${\boldsymbol{\Pi}}$ will reorder the rows of $\mathbf{U}$, and equivalently the elements of $\tilde{\boldsymbol{\eta}}$ to have the same relative order as the elements of $\mathbf{s}$ [15, 17]. Because $s_{i}$ is monotonically increasing with respect to $(1-q_{0}-q_{1})(h_{i}\theta-\tau_{i})$, the elements of $\tilde{\boldsymbol{\eta}}$ should be reordered by the permutation matrix to have the same relative order as the elements of $(1-q_{0}-q_{1})({\mathbf{h}}\theta-{\boldsymbol{\tau}})$ to maximize the likelihood.

If $\tau_{i}=c_{0}h_{i}$, then changing $\theta-c_{0}$ in $-(\theta-c_{0})$ would reverse the ordering. This might help to explain why two solutions appear in the subsequent Proposition 3 when $\theta$ is unknown.

IV-C Estimation with permuted data and unknown $\theta$

In general, $\theta$ may be unknown. Consequently, we should jointly estimate $\theta$ and permutation matrix $\boldsymbol{\Pi}$. 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 $\theta$ via numerical algorithms and achieves global optimum under $q_{0}=q_{1}=0$ or $1$. Under $0<q_{0}+q_{1}<2$, we do not know whether the negative log-likelihood function is convex or not, and local optimum is guaranteed. Given $\theta$, 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.

IV-C1 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 be viewed as the alternating projection with respect to $\theta$ and $\boldsymbol{\Pi}$. The objective function is $l({\tilde{\boldsymbol{\eta}}};{\theta},{\boldsymbol{\Pi}})$. In step $2$, given $\hat{\theta}_{t-1}$, we update the permutation matrix as $\hat{\boldsymbol{\Pi}}_{t-1}$, and the objective value is $l({\tilde{\boldsymbol{\eta}}};{\hat{\theta}_{t-1}},\hat{\boldsymbol{\Pi}}_{t-1})$. Given $\hat{\boldsymbol{\Pi}}_{t-1}$, we obtain ML estimation of $\theta$ as $\hat{\theta}_{t}$, and the objective value is $l({\tilde{\boldsymbol{\eta}}};{\hat{\theta}_{t}},\hat{\boldsymbol{\Pi}}_{t-1})$ satisfying $l({\tilde{\boldsymbol{\eta}}};{\hat{\theta}_{t}},\hat{\boldsymbol{\Pi}}_{t-1})\geq l({\tilde{\boldsymbol{\eta}}};{\hat{\theta}_{t-1}},\hat{\boldsymbol{\Pi}}_{t-1})$. Given $\hat{\theta}_{t}$, we update the permutation matrix as $\hat{\boldsymbol{\Pi}}_{t}$, and the objective value is $l({\tilde{\boldsymbol{\eta}}};{\hat{\theta}_{t}},\hat{\boldsymbol{\Pi}}_{t})$ satisfying $l({\tilde{\boldsymbol{\eta}}};{\hat{\theta}_{t}},\hat{\boldsymbol{\Pi}}_{t})\geq l({\tilde{\boldsymbol{\eta}}};{\hat{\theta}_{t}},\hat{\boldsymbol{\Pi}}_{t-1})$. Consequently, we have

[TABLE]

Provided that the maximum with respect to each $\theta$ and $\boldsymbol{\Pi}$ is unique, any accumulation point of the sequence generated by Algorithm 1 is a stationary point [29].

IV-C2 Special cases for efficient recovery of $\boldsymbol{\Pi}$ under unknown $\theta$

Proposition 2

*Given the ML estimation problem in (12) with unknown $\theta$, if there exist constants $c,d,e\in{\mathbb{R}}$ such that $c{\boldsymbol{\tau}}+d{\mathbf{h}}=e{\mathbf{1}}$, the elements of $\tilde{\boldsymbol{\eta}}$ should be reordered according to the order of the elements of $(q_{0}+q_{1}-1){\boldsymbol{\tau}}$ if $c=0$, otherwise reordered according to ${\mathbf{h}}$ or $-{\mathbf{h}}$.  *

Proof:

We separately address the cases $c=0$ and $c\neq 0$. In the case of $c=0$, $\mathbf{h}$ must be a constant vector. Reordering according to $(1-q_{0}-q_{1})({\mathbf{h}}\theta-{\boldsymbol{\tau}})$ is equivalent to reordering according to $(q_{0}+q_{1}-1){\boldsymbol{\tau}}$. Since $(q_{0},q_{1})$ are known in this problem, $\tilde{\boldsymbol{\eta}}$ should be reordered according to ${\boldsymbol{\tau}}$ if $q_{0}+q_{1}>1$ or $-{\boldsymbol{\tau}}$ if $q_{0}+q_{1}<1$. In the case of $c\neq 0$, we have ${\boldsymbol{\tau}}=(e/c){\mathbf{1}}-(d/c){\mathbf{h}}$. Consequently, ${\mathbf{h}}\theta-{\boldsymbol{\tau}}=(\theta+d/c){\mathbf{h}}-(e/c){\mathbf{1}}$, and $\tilde{\boldsymbol{\eta}}$ is reordered according to ${\mathbf{h}}$ or $-{\mathbf{h}}$.

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

IV-D Statistical identifiability

Note that Algorithm 2 may generate two solutions $(\hat{\theta}_{s1},\hat{\boldsymbol{\Pi}}_{s1})$ and $(\hat{\theta}_{s2},\hat{\boldsymbol{\Pi}}_{s2})$. Given system parameters $\mathbf{h}$ and $\boldsymbol{\tau}$, it is important to determine whether the two solutions $(\hat{\theta}_{s1},\hat{\boldsymbol{\Pi}}_{s1})$ and $(\hat{\theta}_{s2},\hat{\boldsymbol{\Pi}}_{s2})$ will yield the same log-likelihood $l({\tilde{\boldsymbol{\eta}}};\hat{\theta}_{s1},\hat{\boldsymbol{\Pi}}_{s1})=l({\tilde{\boldsymbol{\eta}}};\hat{\theta}_{s2},\hat{\boldsymbol{\Pi}}_{s2})$. If $l({\tilde{\boldsymbol{\eta}}};\hat{\theta}_{s1},\hat{\boldsymbol{\Pi}}_{s1})=l({\tilde{\boldsymbol{\eta}}};\hat{\theta}_{s2},\hat{\boldsymbol{\Pi}}_{s2})$, two pairs of parameter values lead to the same maximum likelihood. In this situation, $(\theta,{\boldsymbol{\Pi}})$ clearly cannot be estimated consistently since $\tilde{\boldsymbol{\eta}}$ provide no information as to whether it is $(\hat{\theta}_{s1},\hat{\boldsymbol{\Pi}}_{s1})$ or $(\hat{\theta}_{s2},\hat{\boldsymbol{\Pi}}_{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 ($({\theta}_{s1},{\boldsymbol{\Pi}}_{s1})$ and $({\theta}_{s2},{\boldsymbol{\Pi}}_{s2})$) leading to the same distribution of the observations ${\tilde{\boldsymbol{\eta}}}$ [32].

Proposition 3

*Let ${\mathbf{h}}_{a}$ and ${\mathbf{h}}_{d}$ denote the ascending and descending ordered versions of $\mathbf{h}$, and ${\boldsymbol{\Pi}}_{a}{\mathbf{h}}={\mathbf{h}}_{a}$ and ${\boldsymbol{\Pi}}_{d}{\mathbf{h}}={\mathbf{h}}_{d}$, where ${\boldsymbol{\Pi}}_{a}$ and ${\boldsymbol{\Pi}}_{d}$ are permutation matrices. Given ${\boldsymbol{\tau}}=c_{0}{\mathbf{h}}$ and ${\mathbf{h}}_{a}=-{\mathbf{h}}_{d}$, the model is unidentifiable, i.e., $l({\tilde{\boldsymbol{\eta}}};{\theta},{\boldsymbol{\Pi}})|_{\theta={\theta}_{s1},{\boldsymbol{\Pi}}={\boldsymbol{\Pi}}_{s1}}=l({\tilde{\boldsymbol{\eta}}};{\theta},{\boldsymbol{\Pi}})|_{\theta={\theta}_{s2},{\boldsymbol{\Pi}}={\boldsymbol{\Pi}}_{s2}}$, where ${\theta}_{s2}=2c_{0}-{\theta}_{s1}$ and ${\boldsymbol{\Pi}}_{s2}={\boldsymbol{\Pi}}_{s1}{\boldsymbol{\Pi}}_{a}^{\rm T}{\boldsymbol{\Pi}}_{d}$.  *

Proof:

Let ${\boldsymbol{\Pi}}_{s1}$ be a permutation matrix such that ${\boldsymbol{\Pi}}_{s1}^{\rm T}{\tilde{\boldsymbol{\eta}}}$ has the same relative order as $\mathbf{h}$. Now we prove that ${\boldsymbol{\Pi}}_{s2}^{\rm T}{\tilde{\boldsymbol{\eta}}}$ has the same relative order as $-\mathbf{h}$. Utilizing ${\mathbf{h}}_{a}=-{\mathbf{h}}_{d}=-{\boldsymbol{\Pi}}_{d}{\mathbf{h}}$ and ${\boldsymbol{\Pi}}_{d}{\boldsymbol{\Pi}}_{d}^{\rm T}={\mathbf{I}}$, we obtain ${\boldsymbol{\Pi}}_{d}^{\rm T}{\mathbf{h}}_{a}=-{\mathbf{h}}$. Note that ${\boldsymbol{\Pi}}_{s2}^{\rm T}{\tilde{\boldsymbol{\eta}}}={\boldsymbol{\Pi}}_{d}^{\rm T}{\boldsymbol{\Pi}}_{a}{\boldsymbol{\Pi}}_{s1}^{\rm T}{\tilde{\boldsymbol{\eta}}}$. Because ${\boldsymbol{\Pi}}_{s1}^{\rm T}{\tilde{\boldsymbol{\eta}}}$ has the same relative order as $\mathbf{h}$, ${\boldsymbol{\Pi}}_{a}{\boldsymbol{\Pi}}_{s1}^{\rm T}{\tilde{\boldsymbol{\eta}}}$ has the same relative order as ${\boldsymbol{\Pi}}_{a}\mathbf{h}={\mathbf{h}}_{a}$, and ${\boldsymbol{\Pi}}_{d}^{\rm T}{\boldsymbol{\Pi}}_{a}{\boldsymbol{\Pi}}_{s1}^{\rm T}{\tilde{\boldsymbol{\eta}}}$ has the same relative order as ${\boldsymbol{\Pi}}_{d}^{\rm T}{\mathbf{h}}_{a}=-{\mathbf{h}}$.

Next we prove that $l({\tilde{\boldsymbol{\eta}}};{\theta},{\boldsymbol{\Pi}})|_{\theta={\theta}_{s1},{\boldsymbol{\Pi}}={\boldsymbol{\Pi}}_{s1}}=l({\tilde{\boldsymbol{\eta}}};{\theta},{\boldsymbol{\Pi}})|_{\theta={\theta}_{s2},{\boldsymbol{\Pi}}={\boldsymbol{\Pi}}_{s2}}$ holds. Because ${\theta}_{s2}=2c_{0}-{\theta}_{s1}$, we have

[TABLE]

By examining $l({\tilde{\boldsymbol{\eta}}};{\theta},{\boldsymbol{\Pi}})$ (13) and utilizing ${\mathbf{h}}_{a}=-{\mathbf{h}}_{d}$, the second term of $l({\tilde{\boldsymbol{\eta}}};{\theta},{\boldsymbol{\Pi}})|_{\theta={\theta}_{s1},{\boldsymbol{\Pi}}={\boldsymbol{\Pi}}_{s1}}$ is equal to that of $l({\tilde{\boldsymbol{\eta}}};{\theta},{\boldsymbol{\Pi}})|_{\theta={\theta}_{s2},{\boldsymbol{\Pi}}={\boldsymbol{\Pi}}_{s2}}$. For the first term, note that given ${\theta}_{s1}$ and ${\theta}_{s2}$, the corresponding ${\mathbf{s}}_{1}$ and ${\mathbf{s}}_{2}$ in (13) can be viewed as evaluating at $\mathbf{h}$ and $-\mathbf{h}$ according to (15), respectively. Because ${\mathbf{h}}_{a}=-{\mathbf{h}}_{d}$, we can conclude that ${\mathbf{s}}_{1}$ is a permutated version of ${\mathbf{s}}_{2}$. The first term of (13) can be expressed as either $({\boldsymbol{\Pi}}_{s1}^{\rm T}{\tilde{\boldsymbol{\eta}}})^{\rm T}{\mathbf{s}}_{1}$ or $({\boldsymbol{\Pi}}_{s2}^{\rm T}{\tilde{\boldsymbol{\eta}}})^{\rm T}{\mathbf{s}}_{2}$. Because $({\boldsymbol{\Pi}}_{s1}^{\rm T}{\tilde{\boldsymbol{\eta}}})^{\rm T}$ and ${\mathbf{s}}_{1}$ have the same relative order as $\mathbf{h}$, and $({\boldsymbol{\Pi}}_{s2}^{\rm T}{\tilde{\boldsymbol{\eta}}})^{\rm T}$ and ${\mathbf{s}}_{2}$ have the same relative order as $-\mathbf{h}$, one has $({\boldsymbol{\Pi}}_{s1}^{\rm T}{\tilde{\boldsymbol{\eta}}})^{\rm T}{\mathbf{s}}_{1}=({\boldsymbol{\Pi}}_{s2}^{\rm T}{\tilde{\boldsymbol{\eta}}})^{\rm T}{\mathbf{s}}_{2}$. Thus $l({\tilde{\boldsymbol{\eta}}};{\theta},{\boldsymbol{\Pi}})|_{\theta={\theta}_{s1},{\boldsymbol{\Pi}}={\boldsymbol{\Pi}}_{s1}}=l({\tilde{\boldsymbol{\eta}}};{\theta},{\boldsymbol{\Pi}})|_{\theta={\theta}_{s2},{\boldsymbol{\Pi}}={\boldsymbol{\Pi}}_{s2}}$.

Now an example is presented to substantiate the above proposition. Let $c_{0}=0.5$, the true value $\theta_{0}=1$, ${\mathbf{h}}=[2,-1,-2,1]^{\rm T}$, ${\boldsymbol{\eta}}=[\eta_{1},\eta_{2},\eta_{3},\eta_{4}]^{\rm T}$ and ${\boldsymbol{\Pi}}_{0}=[0~{}0~{}1~{}0;0~{}1~{}0~{}0;0~{}0~{}0~{}1;1~{}0~{}0~{}0]$. Then $\tilde{\boldsymbol{\eta}}=[\eta_{3},\eta_{2},\eta_{4},\eta_{1}]^{\rm T}$, ${\mathbf{h}}_{a}=[-2,-1,1,2]^{\rm T}$, ${\mathbf{h}}_{d}=[2,1,-1,-2]^{\rm T}$ and ${\mathbf{h}}_{a}=-{\mathbf{h}}_{d}$. We can conclude that $l({\tilde{\boldsymbol{\eta}}};{\theta},{\boldsymbol{\Pi}})|_{\theta=1,{\boldsymbol{\Pi}}={\boldsymbol{\Pi}}_{0}}=l({\tilde{\boldsymbol{\eta}}};{\theta},{\boldsymbol{\Pi}})|_{\theta=0,{\boldsymbol{\Pi}}={\boldsymbol{\Pi}}^{\prime}}$, where ${\boldsymbol{\Pi}}^{\prime}=[0~{}0~{}0~{}1;1~{}0~{}0~{}0;0~{}0~{}1~{}0;0~{}1~{}0~{}0]$.

In addition, given $|c_{0}|\geq\Delta$, only one of $\{\theta_{s1},\theta_{s2}\}$ lies in the interval $[-\Delta,\Delta]$, 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.

IV-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 $\theta$ 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

[TABLE]

where $\overset{\rm p}{\longrightarrow}$ denotes convergence in probability. Given $\theta\in[-\Delta,\Delta]$, $-|h_{i}|\Delta-\tau_{i}\leq h_{i}\theta-\tau_{i}\leq|h_{i}|\Delta-\tau_{i}$. In the following text, we only deal with $q_{0}+q_{1}<1$ case. The case that $q_{0}+q_{1}>1$ is very similar and is omitted here. Define $l=\underset{i\in[1,\cdots,N]}{\operatorname{min}}~{}\left(q_{0}+(1-q_{0}-q_{1})F_{w}(({-|h_{i}|\Delta-\tau_{i}})/{\sigma_{w}})\right)$ and $u=\underset{i\in[1,\cdots,N]}{\operatorname{max}}~{}\left(q_{0}+(1-q_{0}-q_{1})F_{w}(({|h_{i}|\Delta-\tau_{i}})/{\sigma_{w}})\right)$. Then ${\eta}_{i}$ should satisfy $l\leq{\eta}_{i}\leq u$. Let ${\mathcal{I}}_{l,u}(\tilde{\eta}_{i})$ denotes the projection of $\tilde{\eta}_{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

[TABLE]

Utilizing ${\boldsymbol{\Pi}}{\boldsymbol{\Pi}}^{\rm T}=\mathbf{I}$ yields

[TABLE]

which is a quadratic equation in $\theta$. Accordingly, using the asymptotic properties of ${\mathbf{m}}^{\rm T}{\mathbf{m}}$, one obtains (18) via inverting (17)

[TABLE]

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, $-\Delta$ and $\Delta$ 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 ${\boldsymbol{\tau}}=c_{0}{\mathbf{h}}$ into (18), the solutions are $\theta_{1}=\theta$ and $\theta_{2}=2c_{0}-\theta$.

V Detection 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 $\theta$. In addition, we investigate the permutation matrix recovery probability.

V-A Detection with permuted data and known $\theta$

In the case of known $\theta$, the GLRT can be formulated as

[TABLE]

As shown in Proposition 1, the ML estimate of the permutation matrix ${\boldsymbol{\Pi}}$ corresponding to the first term in (19) will reorder the elements of $\tilde{\boldsymbol{\eta}}$ to have the same relative order as the elements of $(1-q_{0}-q_{1})({\mathbf{h}}\theta-{\boldsymbol{\tau}})$. Similarly, for the ML estimation problem corresponding to the second term in (19), we reorder the elements of $\tilde{\boldsymbol{\eta}}$ to have the same order as that of $-(1-q_{0}-q_{1}){\boldsymbol{\tau}}$.

V-B Detection with permuted data and unknown $\theta$

For the unknown $\theta$ and unknown $\boldsymbol{\Pi}$ case, a GLRT is used to decide ${\mathcal{H}}_{1}$ if

[TABLE]

Algorithm 1 for joint estimation of $\theta$ and $\boldsymbol{\Pi}$ has been described in section IV-C1. The performance of the GLRT (20) will be evaluated in the Algorithm 1 for joint estimation of $\theta$ and $\boldsymbol{\Pi}$ is necessary for the first term and has been described in section IV-C1; the second term is as in section V-A.

V-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 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 $\theta$. Instead, we assume that $\theta$ 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 $q_{0}+q_{1}<1$ in the following analysis. The case that $q_{0}+q_{1}>1$ is similar and is omitted here.

First, let $p_{i}$ be ordered such that $p_{(1)}>p_{(2)}>\cdots>p_{(K)}$. From (3) we have $(h_{i}\theta-\tau_{i})_{(1)}>(h_{i}\theta-\tau_{i})_{(2)}>\cdots>(h_{i}\theta-\tau_{i})_{(K)}$. Provided $q_{0}+q_{1}<1$, the elements of $\tilde{\boldsymbol{\eta}}$ should be reordered according to the order of the elements of $\mathbf{h}{\theta}-\boldsymbol{\tau}$ in Proposition 1. Therefore the permutation matrix will be correctly recovered if and only if $\eta_{(1)}>\eta_{(2)}>\cdots>\eta_{(K)}$. Note that the subscripts of $(h_{i}\theta-\tau_{i})_{(\cdot)}$ and $\eta_{(\cdot)}$ also correspond to the order of $p_{i}$, instead of the order of $h_{i}\theta-\tau_{i}$ or $\eta_{i}$.

Define $E_{i}$ as the event such that $\eta_{(i)}>\eta_{(i+1)}$ and $\bar{E_{i}}$ as the corresponding complement event of $E_{i}$, namely, $\eta_{(i)}\leq\eta_{(i+1)}$. The probability that permutation matrix is recovered correctly can be written as

[TABLE]

where union bound ${\rm Pr}\left(\bigcup\limits_{i=1}^{K-1}{\bar{E}_{i}}\right)\leq\sum\limits_{i=1}^{K-1}{\rm Pr}(\bar{E}_{i})$ is utilized in (LABEL:unionbound). From (3), we have $u_{ij}\sim\mathcal{B}(1,p_{i})$ and $N\eta_{i}=\sum_{j=1}^{N}u_{ij}\sim\mathcal{B}(N,p_{i})$. When $N$ is large, the De Moivre-Laplace theorem [30] implies that the distribution of $\eta_{i}$ can be approximated by ${\mathcal{N}}(p_{i},{p_{i}(1-p_{i})/N})$. As a consequence, $\eta_{(i)}-\eta_{(i+1)}$ is approximately distributed as ${\mathcal{N}}(p_{(i)}-p_{(i+1)},p_{(i)}(1-p_{(i)})/N+p_{(i+1)}(1-p_{(i+1)})/N)$, and

[TABLE]

where $t=\underset{i=1,\cdots,K-1}{\operatorname{min}}~{}\frac{v_{i}}{\sqrt{p_{(i)}(1-p_{(i)})+p_{(i+1)}(1-p_{(i+1)})}}$, $v_{i}=p_{(i)}-p_{(i+1)}$ and the approximation $\Phi(-x)\approx\frac{1}{\sqrt{2\pi}x}{\rm e}^{-\frac{x^{2}}{2}}(x\gg 0)$ is utilized.

Utilizing $p_{(i)}(1-p_{(i)})+p_{(i+1)}(1-p_{(i+1)})\leq 1/2$, we define $\tilde{t}$ satisfying

[TABLE]

We conjecture that $\tilde{t}$ is on the order of $K^{-\alpha}$, i.e., $\tilde{t}=O(K^{-\alpha})$, which means that there exists constant $c_{t}$ such that

[TABLE]

In the following text, we show that we can construct $\mathbf{h}$ such that $\tilde{t}=O(K^{-1})$ and $\tilde{t}=O(K^{-2})$. According to (22) and (23), the approximation ${\rm Pr}(K,N)$ (21) can be further simplified and relaxed as

[TABLE]

From (25), the exponent $(1+\alpha)\ln K-\frac{1}{2}\ln N-\frac{c_{t}^{2}}{K^{2\alpha}}N$ of (25) must be far less than [math] for the recovery of permutation matrix. Given $N$ is large, the term $-\frac{1}{2}\ln N$ is small compared to $N$. Thus $(1+\alpha)\ln K-\frac{c_{t}^{2}}{K^{2\alpha}}N<0$ will ensure that the permutation matrix can be recovered in high probability. Simplifying $(1+\alpha)\ln K-\frac{c_{t}^{2}}{K^{2\alpha}}N<0$ yields

[TABLE]

The following cases are examples to illustrate $\tilde{t}=O(K^{-\alpha})$. For simplicity, we assume $1-q_{0}-q_{1}>0$, $\boldsymbol{\tau}=c\mathbf{h}(c<\theta)$ and

[TABLE]

V-C1 $\tilde{t}=O(K^{-1})$

Let $\mathbf{h}$ be the shape of a ramp signal such that $h_{i}=u-\frac{(u-l)(i-1)}{K-1}(u>|l|)$, and $w_{ij}\sim{\mathcal{N}}(0,\sigma_{w}^{2})$. Then the ordered sequence $p_{(i)}=p_{i}$, and $\tilde{t}$ can be approximated as

[TABLE]

where mean value theorem is utilized for $\xi_{i}\in(h_{i+1},h_{i})$, $\xi_{1}\approx h_{1}=u$ is utilized when $K$ is large, and

[TABLE]

Therefore $\tilde{t}$ can be reshaped in the form of (23).

V-C2 $\tilde{t}=O(K^{-2})$

Let $h_{i}$ be independently drawn from the same distribution of $w_{ij}/\sigma_{w}$. The CDF of $p_{i}$ is

[TABLE]

In this case, we conjecture that $\tilde{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 $\mathbf{h}$ is the shape of a sinusoidal signal is also presented in Fig. 3.

Now we prove that $\tilde{t}=O(K^{-2})$ under certain conditions. Given that $h_{i}$ and $w_{ij}/\sigma_{w}$ are i.i.d. random variables and $a=1$, the CDF $F_{p_{i}}(x)=(x-q_{0})/(1-q_{0}-q_{1})$, and the PDF of $p_{i}$ is

[TABLE]

Then the variates $p_{(1)},p_{(2)}\cdots,p_{(K)}$ are distributed as $K$ descending ordered statistics from an uniform $(q_{0},1-q_{1})$ parent. For $x\leq(1-q_{0}-q_{1})/(K-1)$, the CDF of $\tilde{t}$ can be derived as [31] (page 135, equation (6.4.3))

[TABLE]

For $x\geq(1-q_{0}-q_{1})/(K-1)$, $F_{\tilde{t}}(x)=1$. Then the PDF of $\tilde{t}$ is

[TABLE]

The expectation of $\tilde{t}$ is

[TABLE]

Hence the probability that $\tilde{t}$ falls into $[{c_{1}}/{K^{2}},{c_{2}}/{K^{2}}]$ is

[TABLE]

When $K$ is large, $(K-1)/{K}\approx 1$ and $(1-1/(c^{\prime}K))^{c^{\prime}K}\approx 1/e(c^{\prime}>0)$. Equation (V-C2) can be approximated as

[TABLE]

Provided that $q_{0}=q_{1}=0$, when $c_{1}=0.1$ and $c_{2}=10$, ${\rm Pr}({0.1}/{K^{2}}\leq\tilde{t}\leq{10}/{K^{2}})\approx 0.94$; when $c_{1}=0.01$ and $c_{2}=100$, ${\rm Pr}({0.01}/{K^{2}}\leq\tilde{t}\leq{100}/{K^{2}})\approx 0.99$. It can be seen that $\tilde{t}$ falls near the order of magnitude of $K^{-2}$ with high probabilities. Thus it is reasonable that $\tilde{t}=O(K^{-2})$.

According to the definition of $p_{i}$ (3), equations (22) and (23), $c_{t}\propto 1-q_{0}-q_{1}$. From (26), the number of quantizers $N_{\rm req}$ required for permutation matrix recovery probability is

[TABLE]

From (37), one can conclude that the number of quantizers for permutation matrix recovery with high probability is ${1}/{(1-q_{0}-q_{1})^{2}}$ times that of unflipped case where $q_{0}=q_{1}=0$.

VI numerical Simulations

In this section, numerical experiments are conducted to evaluate the theoretical results. For simplicity, the distribution of noise $w_{ij}$ is selected as the Gaussian distribution $\mathcal{N}(0,\sigma_{w}^{2})$.

VI-A Parameter estimation

For the first two experiments, we evaluate the performance of ML estimators proposed in section IV. Parameters are set as follows: $K=20$, $\theta=1$, $\sigma_{w}^{2}=1$, $\Delta=2$, $q_{0}=0.05$, $q_{1}=0.05$ and the tolerance parameter $\epsilon$ 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 $\boldsymbol{\tau}=0.5\mathbf{h}$, which is a special case mentioned in Proposition 2. The coefficients $\mathbf{h}$ is equispaced with ${\mathbf{h}}=[-1.50,-1.29,-1.08,\cdots,2.50]^{\rm T}$, which corresponds to a ramp signal. It can be seen that $\mathbf{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 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.

For the second experiment, the MSE performance of Algorithm 1 (for the general case) is evaluated In Fig. 5. The elements of the vector $\mathbf{h}$ describe the shape of a sinusoidal signal such that $h_{i}=\sin(2\pi x_{i})$, where $x_{i}$ is drawn independently and randomly from the uniform distribution $\mathcal{U}(0,1)$ and then sorted in ascending order. The elements of the vector ${\boldsymbol{\tau}}$ is drawn independently and randomly from the uniform distribution $\mathcal{U}(-\Delta,\Delta)$. 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\times 10^{4}$.

VI-B Signal detection

In Fig. 6, the relationship between $P_{D}$ and the number of quantizers $N$ is employed. Parameters are consistent with the first experiment, except that $\sigma_{w}^{2}=9$ and $P_{FA}=0.05$.

In subgraph ($a$), $\mathbf{h}$ and $\boldsymbol{\tau}$ 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$), $\mathbf{h}$ and $\boldsymbol{\tau}$ 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.

VI-C Permutation matrix recovery

In this subsection, the approximations for permutation matrix recovery are verified. Parameters are set as follows: $K=20$ , $\theta=1.5$, $\Delta=2$, $q_{0}=0$, $q_{1}=0$ and $\sigma_{w}^{2}=1$. The number of Monte Carlo trials is $1000$.

First, the relationship of $t$ and $\tilde{t}$ (22) and the conjecture of $\tilde{t}$ (23) are illustrated in three cases. From Fig. 7, one obtains that $t$ can be approximated as $\sqrt{2}{\tilde{t}}$ in practice. For a ramp signal, $\mathbf{h}=[-0.800,-0.705,-0.610,\cdots,1.000]^{\rm T}$ and $\boldsymbol{\tau}=0.5\mathbf{h}$. $t\approx\sqrt{2}c_{e}/K$ where $c_{t}=c_{e}=0.4355$ is evaluated via (29). Because of the gap between $t$ and $\sqrt{2}c_{e}/K$, we use linear regression to fit $t$ and obtain $c_{ea}=0.6717$, which is much more accurate than $c_{e}$ and will be utilized later to predict the number of quantizers for permutation matrix recovery. For random generated $\mathbf{h}$, $\mathbf{h}$ is drawn from standard normal distribution and $\boldsymbol{\tau}=0.5\mathbf{h}$. It can be seen that $t$ can be approximated by $1/K^{2}$. For a sinusoidal signal, $\mathbf{h}$ and $\boldsymbol{\tau}$ are drawn in the same way of the second experiment. We use linear regression and obtain $t\approx 0.71/K^{2.23}\approx{\sqrt{2}\tilde{t}}=\sqrt{2}c_{t,s}/K^{\alpha_{t,s}}$, $c_{t,s}=0.5020$ and $\alpha_{t,s}=2.23$.

Next, the empirical permutation recovery probability ${\rm Pr}(\hat{\boldsymbol{\Pi}}_{\rm ML}={\boldsymbol{\Pi}}_{0})$ versus $N$ or $K$ are presented in Fig. 8, and the theoretical approximations ${\rm Pr}(K,N)$ (21) and $\widetilde{\rm 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=10^{4}$. All $\mathbf{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 $\theta$, which has negligible difference compared to that in the known $\theta$ case.

In subgraph($a$), it can be seen that the permutation matrix of the ramp signal can be recovered with high probability given $N\geq 5000$. From $N>\frac{1+\alpha}{c_{t}^{2}}K^{2\alpha}\ln K$ (26) where $c_{t}=c_{e}=0.4355$ and $\alpha=1$, one can conclude that $N>\frac{2}{0.4355^{2}}K^{2}\ln K|_{K=20}\approx 12636$, which is more than twice of $5000$. Utilizing the fitted parameter $c_{ea}$, one obtain a more accurate result that $N>\frac{2}{0.6717^{2}}K^{2}\ln K|_{K=20}\approx 5312$ ensures permutation matrix recovery with high probability. For random $\mathbf{h}$, $N>3K^{4}\ln K|_{K=20}\approx 1.438\times 10^{6}$ ensures recovery with high probability, which is not accurate enough, as subgraph($b$) shows that $N\approx 10^{5}$ is enough for recovery of permutation matrix. In subgraph($c$), it is shown that $N\approx 10^{6}$ is enough for recovery of permutation matrix, which is also inaccurate compared to the fitted results of the sinusoidal signal $N>\frac{3.23}{0.5020^{2}}K^{4.46}\ln K|_{K=20}\approx 2.437\times 10^{7}$. The numerical results show that the theoretical bound ${\rm Pr}(K,N)$ is accurate in predicting $N$ with high probability in permutation matrix recovery, which demonstrates that $\widetilde{\rm Pr}(K,N)$ may be too conservative in predicting the number of quantizers ensuring perfect permutation matrix recovery. In subgraph($d$), $10000=N>\frac{2}{0.6717^{2}}K^{2}\ln K|_{K=26}\approx 9763$, thus $K\leq 26$ will ensure permutation matrix recovery with high probability, which is consistent with the numerical results.

In Fig. 9, the relationship of flipping probabilities $(q_{0},q_{1})$ and number of quantizers $N_{\rm req}$ (37) required for permutation matrix recovery with high probability is verified. Parameters are the same as those in Fig. 8-($a$) except for $(q_{0},q_{1})$. We use the result of the experiment in which $q_{0}=q_{1}=0$ to predict those in which $q_{0}=q_{1}=0.05$, $q_{0}=q_{1}=0.1$ and $q_{0}=q_{1}=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).

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\alpha}\log K$, where $\alpha$ is a constant depending on the signal shape and the distribution of noise.