# A hybrid model-based and learning-based approach for classification with
limited number of data samples

Alireza Nooraiepour, Waheed U. Bajwa and Narayan B. Mandayam
The authors are with WINLAB, Department of Electrical and Computer
Engineering, Rutgers University, NJ, USA. Emails:
<EMAIL_ADDRESS><EMAIL_ADDRESS>This
work was supported by the National Science Foundation (NSF) under grants
ECCS-2028823 and OAC-1940074, and in part by ACI-1541069.

###### Abstract

The fundamental task of classification is considered for physical systems with
known parametric statistical models given a limited number of data samples.
The proposed solution, HyPhyLearn, is a hybrid classification method, which
exploits both the physics-based statistical models and the learning-based
classifiers. Notably, the standalone learning-based and statistical
classifiers face major challenges towards the fulfillment of the
classification task. Specifically, the physics-based statistical models
usually suffer from the inability to properly tuning the underlying
unobservable parameters, which would lead to a mismatched representation of
the system’s behaviors. Learning-based classifiers, on the other hand, rely on
a large number of data from the underlying physical process which might not be
accessible in most practical scenarios. In this vein, we conjecture that
HyPhyLearn would alleviate the difficulties associated with each approach by
fusing their individual strength. The proposed hybrid approach estimates the
the unobservable model parameters using the available (suboptimal) estimation
procedures, and subsequently use the physics-based statistical models to
generate synthetic data. Then, the data samples are incorporated with the
synthetic data in a learning-based classifier powered from adversarial
training of the neural networks. In particular, to address the mismatch
problem the classifier learns a mapping to a common feature space from the
data samples and the synthetic data. Simultaneously, the classifier is trained
to find discriminative features from this space in order to fulfill the
classification task. Two case studies revolving around two concrete
communication problems are presented in order to highlight the applicability
of HyPhyLearn. Numerical results demonstrate that the proposed approach leads
to major classification improvements in comparison to the existing standalone
or hybrid classification methods.

##  I Introduction

We revisit the problem of classification with limited number of training data
samples in this paper. The fundamental task of classification comes up in
various fields and is traditionally tackled within two frameworks: $1)$
statistical setting, and $2)$ fully data-driven setting. In the first case,
the main assumption is that data generation adheres to a known probabilistic
model of the underlying physical process. Subsequently, the classification
problem is usually dealt with within a hypothesis testing (HT) framework aimed
at testing between two (or more) hypotheses. Here, optimality in both the
Bayesian sense and the Neyman–Pearson sense relies on computation of the
likelihood-ratio terms, which requires clairvoyant knowledge of the
probabilistic models under different hypotheses [1]. However, accurate
modeling of the physical processes in increasingly complex engineered systems
is either not tractable or it relies on a large number of unobservable
parameters, estimation of which from limited number of data samples could be a
major hurdle [2, 3]. As a result, a mismatch between the physics-based
statistical models and the real physical processes is inevitable. This
precludes exact computation of the likelihood-ratio values, which deteriorates
the classification performance [4]. The fully data-driven (i.e., learning
based) setting, on the other hand, relies on a large number of data samples
for finding an optimal mapping from the data samples to the corresponding
labels. But availability of such data in many real-world problems, e.g.,
channel-based spoofing detection [5] and signal identification [6], is
generally limited, which might lead to learning of a suboptimal map. Moreover,
one should always expect mislabeled data in many applications, since the
employed labeling procedures might not be error free. Consequently,
classification performance of data-driven models can be seriously limited for
many real-world applications.

The overarching objective of this paper is to develop an algorithmic framework
for classification from limited number of training data samples in
applications in which neither model-based nor learning-based approaches alone
result in very good classification performance. To this end, note that
learning-based approaches traditionally tend to disregard the physics-based
models developed to describe the physical phenomena through tractable
mathematical analysis. For instance, in the context of wireless
communications, numerous theoretical models for channels and resource
management have been developed over the years [2, 5, 7]. Despite being
approximations in many cases, these models provide important prior information
about the corresponding physical systems that might be utilized to facilitate
the subsequent classification tasks. At the same time, physics-based models
consist of numerous unobservable parameters, the tuning of which is a major
hurdle for complex systems [3]. For example, physical channel models in the
multi-input multi-output (MIMO) and 5G communications scenarios rely on a
large number of multidimensional parameters that are defined over a mixed set
of discrete and continuous spaces [8, 9]. In such cases, the maximum
likelihood estimation (MLE) of the parameters could incur a formidable
computational cost [9, 10, 11]. Our goal in this context is to develop a
classification framework that can deal with these practical considerations
through a hybrid approach that consolidates physics-based and fully data-
driven classification approaches. The expectation is that the hybrid approach
would fuse the strengths of the two approaches towards achieving an overall
superior classification performance.

Our proposed hybrid approach first employs the (necessarily) suboptimal
parameter estimation methods to estimate the unobservable parameters. Then, it
utilizes them in the physics-based models to generate synthetic data, which
enables us to leverage learning-based classification approaches. The mismatch
between the physics-based models and the underlying physical process is
addressed in a learning setting. Specifically, a neural network is trained to
map the training and synthetic data to a common discriminative feature space,
which is often referred to as domain-invariant space in the domain adaptation
literature [12, 13]. Meanwhile, a neural network-based classifier is trained
on the mapped synthetic data to extract class-specific discriminative features
from them. The resulting classifier in this way is expected to perform well on
both synthetic and training data distributions.

###  I-A Relation to prior works

In the realm of statistical model-based classifiers, the difficulties
associated with estimating the parameters of the physics-based models are
recognized in various works [14, 4]. This is mainly attributed to the inherent
difficulties associated with determining probability distributions from only a
limited number of data samples. Along these lines, classification under the
assumption of mismatched models is considered in several works [15, 14, 16,
4]. Specifically, [14, 16] derive bounds on the probability of classification
error in the presence of mismatch via the $f$-divergence between the true and
mismatched distributions. In contrast to these bounds that are general in the
sense that no assumption is made regarding the underlying distributions, [4]
considers data that are contained in a linear subspace. This enables the
authors to derive an upper bound on the classification error of the mismatched
model that predicts the presence/absence of an error floor. The analyses in
these works, however, do not lead to a classification algorithm for the
mismatched setting as they merely analyze the mismatch problem itself.

The mismatch problem for the learning-based classifiers corresponds to the
cases where the distribution of the available training data is different from
that of the test data. Such mismatches are primarily studied in the transfer
learning (TL) and the data-shift literature [13]. In particular, covariate
shift [17], which is also studied under the name of transductive TL [18],
refers to the case where the underlying data distributions for the test and
training data are different. Concept shift [19], also known as inductive TL
[18], on the other hand, deals with situations in which the posterior
distribution of the labels given the data is not the same for the training and
the test data. A wide range of algorithms have been proposed in order to
alleviate the performance loss due to such shifts. For example, importance-
weighting technique [20, 21] is proposed for the covariate shift scenario to
remove the bias from the training data. Furthermore, algorithms based on
subspace mapping [22] and learning domain-invariant representations [12] have
also been proposed in the literature to address the mismatch problem. The
authors in [22] propose a transfer component analysis method aimed at finding
a transformation under which the maximum mean discrepancy between the true and
mismatched distributions is small. The work in [12] aims at finding a
representation that is invariant for the training and test distributions in
order to mitigate the effect of discrepancies in the subsequent learning
tasks. For the specific task of classification, the authors in [23] introduce
the domain-adversarial neural network (DANN) framework, which extracts domain-
invariant representations via (deep) neural networks that are discriminative
for the training data in order to devise a classifier on the test data.

Deep transfer learning (DTL) is another prime subject related to our work that
studies the transfer learning concept in the context of deep neural networks
(DNNs). DTL considers a DNN that has been pre-trained on the training data as
transferable knowledge useful for the test data. This knowledge can be
transferred based on different strategies. The pre-trained DNNs can either be
used directly for the test data, or serve as an intermediate feature
extracting step that could facilitate the subsequent learning process for the
test data. In another DTL strategy called fine-tuning, the pre-trained DNN or,
certain parts of it, is refined using the available test data to further
improve the effectiveness of transfer knowledge. We refer the reader to [24,
25] for a survey on DTL methods.

Model-based deep learning is another related line of work that aims at
designing systems whose operation combines physics-based models (domain
knowledge) and data. To this end, two main strategies are typically exploited
in such works, known as model-aided networks and DNN-aided inference. The
former results in specialized DNN architectures by identifying structures in a
model-based algorithm; e.g., an iterative structure for the case of deep
unfolding [26]. The latter primarily utilizes model-based methods for
inference, but replaces explicit domain-specific computations with dedicated
DNNs in order to facilitate operation in complex environments; e.g., using
generative models for compressed sensing applications [27]. We refer the
readers to [28] and references therein for the state-of-the-art strategies in
model-based deep learning methods.

There also have been previous attempts to incorporate physics-inferred
information in the fully data-driven setting. In the field of wireless
communications, for instance, the authors in [2] employ DTL to solve a
specific resource management problem. Similarly, the task of signal
classification is tackled via DTL under different practical assumptions, such
as real propagation effects [29], hardware impairments [30] and weak received
signal strength [31]. These works utilize abundant data from an approximate
model along with limited data from the real-world model in the DTL fine-tuning
approach. More closely to the idea of physics-guided machine learning (ML), a
recurrent neural network (RNN) is modified in [32] to incorporate information
from the physics-based model as an internal state of the RNN. Furthermore,
parameters of the physics-based models are combined with sensor readings and
used as input to a DNN to develop a hybrid prognostics model in [3].

We note that the aforementioned works in domain adaptation literature do not
employ any available physics-based statistical models and, consequently, rely
on large number of training data samples for dealing with the mismatch
problem. In addition, model-based deep learning strategies might not be
applicable to the statistical classification problem in general due to the
lack of algorithmic structure such as an iterative structure. Equally
importantly, DTL fine-tuning and physics-guided learning approaches do not
consider the difficulties associated with estimating the physics-based
parameters, which would indeed lead to inaccurate physics-based statistical
models. The resulting discrepancy between the model and the underlying
physical process necessitates a learning-based classifier that is capable of
leveraging the data in a way to alleviate this mismatch problem.

###  I-B Our contributions

The main contributions of this work are as follows.

  * •

We focus on the task of classification for a physical process assuming that a
limited number of training data samples, with possibly mislabeled instances,
is available. We consider the case where the physical process (or its
approximation) can be described by physics-based parametric statistical
models. As these models tend to be complex in general, estimation of the
unknown model parameters using the maximum likelihood estimation (MLE)
procedure could be computationally prohibitive.111As discussed later in
Section II, even using the MLE does not always provide any optimality
guarantees in general for the classification problem in a HT setting [33]. We
instead propose HyPhyLearn—a novel hybrid classification method—as a solution,
which exploits both physics-based statistical models and learning-based
classifiers. This approach makes use of (necessarily suboptimal) parameter
estimation algorithms/heuristics to obtain (approximate) parameter estimates.
Next, plugging in these estimates in the physics-based statistical models
enables us to generate synthetic data. HyPhyLearn then relies on neural
networks (NNs), which are powerful tools for finding a discriminative feature
space, towards obtaining a learning-based classifier. Specifically, the
learning process involves training a NN to map the training and synthetic data
to a common space under which they are not distinguishable. In the mean time,
a learning-based classifier is trained on the synthetic data mapped to the new
space to find discriminative class-level features. Indeed, learning the common
feature space addresses the distribution mismatch problem between the training
data samples and the generated synthetic data due to the errors in parameter
estimation. It is then expected that the classifier trained on the mapped
synthetic data will perform well on both data distributions. We repurpose
theories from the domain adaptation literature based on learning invariant
representations for our specific problem to justify the proposed hybrid
approach. A schematic of HyPhyLearn for a binary classification example is
illustrated in Fig. 1.

Figure 1: A schematic of our proposed hybrid classification approach
(HyPhyLearn) illustrated for a binary classification setting, which exploits
both physics-based statistical models and learning-based classifiers.

  * •

We also consider two prototypical problems from the wireless communications
literature to investigate the performance of our proposed approach and show
its superiority in comparison to the stand-alone statistical model-based
classifiers as well as the fine-tuning approach as the best existing hybrid
approach applicable to these problems. We first consider the problem of
channel spoofing in the wireless communications setting, where an adversary
(Eve) spoofs a legitimate transmitter (Alice) and sends a message to a
legitimate receiver (Bob) [5, 34, 35]. The spoofing detection at Bob involves
making a decision on whether an incoming message corresponds to Alice or Eve.
This can be cast as a binary classification problem at Bob. Second, we revisit
the problem of multi-user detection (MUD) in the uplink of a cellular network,
where different users are asynchronously sharing a channel with a base station
[11]. For a $K$-user system, MUD is basically a $2^{K}$-ary classification
problem in which the goal is to infer $K$ binary information bits from a given
observation. By obtaining likelihood ratio test (LRT) for each problem, we
show that statistical model-based classifiers rely heavily on the wireless
channel parameters in the above problems. However, estimation performance of
these parameters suffers from both the paucity of training data and complexity
of the physics-based statistical models. In fact, these models are complex in
the sense that MLEs of the corresponding parameters require an exhaustive
search over the space of the parameters, which is not feasible for many
communication scenarios including MIMO transmissions in a 5G setting [8]. For
both problems, numerical results show that HyPhyLearn provides major
improvements in terms of the classification accuracy in comparison to the best
existing approaches.

###  I-C Notation and organization

Throughout the paper, vectors are denoted with lowercase bold letters, while
uppercase bold letters are reserved for matrices. Furthermore, equality by
definition is expressed through the symbol $\overset{\bigtriangleup}{=}$. Non-
bold letters are used to denote scalar values and calligraphic letters denote
sets. Furthermore, the cardinality of a set $\mathcal{S}$ is denoted by
$\lvert\mathcal{S}\rvert$. The spaces of real and complex vectors of length
$d$ are denoted by $\mathbb{R}^{d}$ and $\mathbb{C}^{d}$, respectively. The
$m$th element of a vector $\mathbf{u}$ and the trace of a matrix $\mathbf{U}$
are shown by $\mathbf{u}[m]$ and $\operatorname{Tr}(\mathbf{U})$,
respectively. Also, real and imaginary parts of a complex number $a$ are
denoted by $\Re\\{a\\}$ and $\Im\\{a\\}$, respectively. The probability
density function and expectation of a random variable $w$ are denoted by
$p(w)$ and $\mathbb{E}_{p}(w)$, respectively, while $\mathbb{P}[\cdot]$ is
used to denote the probability of an event. The Gaussian and circularly-
symmetric complex Gaussian distributions are denoted by $\mathcal{N}$ and
$\mathcal{CN}$, respectively, while the uniform distribution supported between
two real numbers $a$ and $b$ is denoted by $\text{unif}(a,b)$. We denote the
$k$th standard basis vector of length $N$ in $\mathbb{R}^{N}$ by
$\mathbf{e}_{k}$, and use $\mathbf{\|\mathbf{u}\|}$ to refer to the Euclidean
norm of the vector $\mathbf{u}$. We refer to identity matrix of size $N$ and
the indicator function by $\mathbf{I}_{N}$ and
$\mathds{1}_{\mathcal{A}}(\mathbf{x})\overset{\bigtriangleup}{=}\begin{cases}1,\mathbf{x}\in\mathcal{A}\\\
0,\mathbf{x}\notin\mathcal{A}\end{cases}$, respectively. Transpose and
conjugate transpose of $\mathbf{u}$ are denoted by $\mathbf{u}^{T}$ and
$\mathbf{u}^{H}$, respectively. Furthermore, $\mathbf{e}_{n}(y)$ refers to a
one-hot encoded version of a non-negative integer $y$, which equals to an all-
zero vector of length $n$ except for the $y$th element which is set to $1$.
Also, $\circ$ and $\odot$ denote the Schur componentwise and the Khatri-Rao
product, respectively, while $\otimes$ is reserved for the Kronecker product.
Finally, given two vectors $\mathbf{a}$ and $\mathbf{b}$ of length $M$,
Toeplitz matrix of size $M\times M$ is defined as
$\text{toep}(\mathbf{a},\mathbf{b})\overset{\bigtriangleup}{=}\begin{bmatrix}\mathbf{a}[1]&\mathbf{b}[2]&\dots&\mathbf{b}[M]\\\
\mathbf{a}[2]&\ddots&\ddots&\mathbf{b}[M-1]\\\ \vdots&\ddots&\ddots&\vdots\\\
\mathbf{a}[M]&\mathbf{a}[M-1]&\dots&\mathbf{a}[1]\end{bmatrix}$.

The rest of the paper is organized as follows. The problem is formally posed
in Section II. Our proposed solution is described in Section III, which
discusses various pieces of HyPhyLearn approach. We introduce the first case
study involving the spoofing detection problem in Section IV. The second case
study, which concerns the multi-user detection problem, is presented in
Section V. We present numerical results concerning the application of our
proposed approach in the above two case studies in Section VI, and contrast it
with the existing methods. Finally, the paper is concluded in Section VII.

##  II Problem Formulation

Consider a physical process consisting of $C$ distinct behaviors where the
physics-based parametric statistical model for the $i$th behavior is available
in the form of a parametric probability density function (PDF) denoted by the
conditional prior $p_{i}(\mathbf{x};{\bm{\theta}_{i}})$ on observations
$\mathbf{x}$ that belong to an observation space $\mathcal{X}$. Assuming the
true underlying parameter for the $i$th behavior is $\bm{\theta}^{*}_{i}$, the
data for this behavior is generated by drawing independent and identically
distributed (i.i.d.) samples from $p_{i}(\mathbf{x};{\bm{\theta}^{*}_{i}})$.
Assuming further that the $i$th behavior is chosen with a prior probability
$\pi_{i}$, our goal is to devise a decision rule to determine a given sample
$\mathbf{x}=[x_{1},\dots,x_{n}]^{T}$ is generated under which behavior.
Clearly, this can be cast as a $C$-ary classification problem via
$H_{i}:\mathbf{x}\sim p_{i}(\mathbf{x};{\bm{\theta}^{*}_{i}}),\
i=0,\dots,C-1$. We consider the case where this decision is made by a
classifier $h_{\bm{\phi}}(\cdot)$ parameterized by
$\bm{\phi}\in\mathbb{R}^{d}$,
$h_{\bm{\phi}}(\mathbf{x}):\mathcal{X}\rightarrow\\{0,\dots,C-1\\}$, which
partitions $\mathcal{X}$ into $C$ disjoint sets, $\\{\mathcal{X}_{i}\\}$, and
decides in favor of $H_{i}$ if $\mathbf{x}\in\mathcal{X}_{i}$. Defining
${\bm{\theta}^{*}}{\overset{\bigtriangleup}{=}}[{\bm{\theta}^{*}_{0}},\dots,{\bm{\theta}^{*}_{C-1}}]$,
we denote the probability of error associated with $h_{\bm{\phi}}(\mathbf{x})$
by $\mathbb{P}_{\bm{\theta}^{*}}[{e}_{\mathbf{\phi}}]$, which can be computed
as

$\displaystyle\mathbb{P}_{\bm{\theta}^{*}}[{e}_{\bm{\phi}}]=\sum_{i=0}^{C-1}\pi_{i}\int_{\mathcal{X}}p_{i}(\mathbf{x};{\bm{\theta}^{*}_{i}})\mathds{1}_{\\{h_{\bm{\phi}}(\mathbf{x})\neq
i\\}}(\mathbf{x})d\mathbf{x},$ (1)

where ${e}_{\bm{\phi}}$ indicates the event that $h_{\bm{\phi}}(\mathbf{x})$
makes an erroneous decision. The optimal classifier
$h_{\bm{\phi}^{*}}(\mathbf{x})$ that minimizes the error probability is given
by the Bayes decision rule, i.e.,
$h_{\phi^{*}}(\mathbf{x})=\operatorname*{\arg\\!\max}_{i=0,\dots,C-1}\
\pi_{i}p_{i}(\mathbf{x};{\bm{\theta}^{*}_{i}})$ [1]. For the specific case of
$C=2$, this rule takes the famous form of the likelihood ratio test,
$\frac{p_{1}(\mathbf{x};{\bm{\theta}^{*}_{1}})}{p_{0}(\mathbf{x};{\bm{\theta}^{*}_{0}})}\underset{y=0}{\overset{y=1}{\gtreqless}}\frac{\pi_{0}}{\pi_{1}}$,
where $y=i$ implies making a decision in favor of the $i$th behavior.

We focus in this paper on the case where although the parametric model
$p_{i}(\mathbf{x};{\bm{\theta}_{i}})$ is known for the $i$th behavior, one
does not have access to the corresponding underlying true parameter
$\bm{\theta}^{*}_{i}$. Instead, only a small number of training data generated
in an i.i.d. manner from $p_{i}(\mathbf{x};{\bm{\theta}^{*}_{i}}),\forall i,$
are available. Specifically, we denote the available dataset by
$\mathcal{D}_{r}=\\{\mathbf{x}_{r,n}\\}_{n=1}^{N_{r}}$, where $N_{r}$ is the
total number of data samples. Also, the corresponding ground-truth label for
the $n$th sample is denoted by $y_{r,n}$ which is only given for $N_{r,l}$
number of data samples where $N_{r,l}\leq N_{r}$. Furthermore, we consider the
case where the model $p_{i}(\mathbf{x};{\bm{\theta}_{i}})$ under the $i$th
behavior is a non-trivial function of the underlying parameter for which
conventional estimation procedures such as maximum likelihood estimation (MLE)
are either not available or are computationally prohibitive to implement. The
implication of this aspect of the problem formulation is that the performance
of any suboptimal parameter estimation method is bound to be limited. As a
result, statistical model-based classifiers, which plug-in these estimates in
$p_{i}(\mathbf{x};{\bm{\theta}_{i}})$, would have a deteriorated performance
as well.

Unlike these classifiers that rely heavily on the knowledge of the parametric
statistical models and the estimated parameters, a purely data-driven approach
can result in a classifier that disregards the available parametric models.
However, as the data generation processes are governed by non-trivial models,
a large number of data is needed in this case to extract related patterns from
each behavior that would lead to a highly discriminative feature space. By
noting that the performance of the fully data-driven and the statistical
model-based classifiers is particularly curbed when they are used in a stand-
alone fashion, we conjecture that fusing the strengths of the two can lead to
a superior classification algorithm in our setting, as described in the next
section.

Before delving into the proposed solution for the described problem setting,
we discuss further two existing approaches towards obtaining a statistical
model-based classifier for the benefit of the reader. Recall that within the
framework of statistical model-based classification, one would first estimate
the unknown model parameters as $\widehat{\bm{\theta}}_{i}$’s, $i=1,\dots,C$,
and plug them in the available models to obtain
$p_{i}(\mathbf{x};{\widehat{\bm{\theta}}}_{i})$. The resulting plug-in models
are then used in practice in lieu of the true models within the optimal Bayes
decision rule. The parameters, $\bm{\phi}$, of the resulting _plug-in_
classifier consist solely of the parameters of physics-based statistical
models, i.e.,
$\bm{\phi}={\bm{\theta}}=[{\bm{\theta}}_{0},\dots,{\bm{\theta}}_{C-1}]$.222For
notational simplicity and without loss of generality, we have not included the
priors as part of the unknown parameters in the current discussion. Based on
this fact, we denote the _plug-in_ classifier by $h_{\bm{\theta}}(\mathbf{x})$
in the remainder of this section. The unknown model parameters can be
estimated using numerous approaches. In the following, we discuss two of the
most popular ways to estimate them as well as the shortcomings of these
approaches that warrant a new approach to classification.

Empirical error minimizer: Given a set of training data with their
corresponding labels, $\\{\mathbf{x}_{r,n},y_{r,n}\\}_{n=1}^{N_{r}}$, the most
natural approach for parameter estimation corresponds to the setting in which
the resulting plug-in classifier, $h_{{\bm{\theta}}}(\mathbf{x})$, minimizes
the empirical error probability defined by
$\widehat{\mathbb{P}}^{N_{r}}[{e}_{\bm{\theta}}]\overset{\bigtriangleup}{=}\frac{1}{N_{r}}\sum_{n=1}^{N_{r}}\mathds{1}_{\\{h_{\bm{\theta}}(\mathbf{x}_{r,n})\neq
y_{r,n}\\}}$. Specifically, for the case of $C=2$ consider the family of the
classifiers ${h_{{\bm{\theta}}}(\mathbf{x})}=\begin{cases}0,\quad\pi
p_{{\bm{\theta}}_{0}}(\mathbf{x})>(1-\pi)p_{{\bm{\theta}}_{1}}(\mathbf{x}),\\\
1,\quad\text{otherwise},\end{cases}$ for which the parameter values
${{\bm{\theta}}_{0}}$ and ${{\bm{\theta}}_{1}}$ are chosen from a space
$\bm{\Theta}$. The parameter estimates that minimize the empirical error are
obtained as
$\widehat{\bm{\theta}}=[\widehat{\bm{\theta}}_{0},\widehat{\bm{\theta}}_{1}]\in\operatorname*{\arg\\!\min}_{{\bm{\theta}}}\widehat{\mathbb{P}}^{N_{r}}[{e}_{{\bm{\theta}}}]$.
The following lemma, which is a direct result of Corollary $16.1$ in [36],
presents an upper bound on the performance of the Bayes decision rule in terms
of that of the plug-in classifier that is obtained using empirical error
minimization.

######  Lemma 1.

If $\bm{\theta}^{*}_{0},\ \bm{\theta}^{*}_{1}\in\bm{\Theta}$, then the error
probability of the Bayes decision rule, with the probability at least
$1-\delta$, is bounded by

$\displaystyle{\mathbb{P}}_{{\bm{\theta}^{*}}}[{e}_{{\bm{\theta}^{*}}}]\leq\widehat{\mathbb{P}}^{N_{r}}[{e}_{\widehat{\bm{\theta}}}]+8\sqrt{\frac{2}{N_{r}}\log\frac{8b}{\delta}},$
(2)

where $b$ denotes the Vapnik–Chervonenkis (VC) dimension [36] of the family of
classifiers, $h_{{\bm{\theta}}}(\mathbf{x})$, defined above.

The above lemma guarantees a $\mathcal{O}(\sqrt{\log{N_{r}}/{N_{r}}})$ rate of
convergence to the Bayes error for $h_{\hat{\bm{\theta}}}(\mathbf{x})$ when
$\widehat{\bm{\theta}}$ is chosen to minimize the empirical error. However,
obtaining such $\widehat{\bm{\theta}}$ is computationally expensive in general
as the empirical error probability might be a non-trivial function of the
parameters.

Maximum likelihood estimator: In practice, the unknown model parameters are
commonly replaced with their corresponding MLEs under each beahvior; the
resulting plug-in classifier gives rise to the well-known generalized
likelihood ratio test (GLRT) for the binary case ($C=2$) [1]. Specifically,
assuming the training data and their corresponding labels are available in the
form of $\\{\mathbf{x}_{r,n},y_{r,n}\\}_{n=1}^{N_{i}}$ for the $i$th
hypothesis, the MLE of $\bm{\theta}_{i}$ is obtained by
$\widehat{\bm{\theta}}^{MLE}_{i}=\operatorname*{\arg\\!\max}_{\bm{\theta}_{i}}\mathcal{L}(\mathcal{D}_{i}|\bm{\theta}_{i})$,
where $\mathcal{L}$ denotes the likelihood function. For the binary case where
$\frac{\pi
p_{1}(\mathbf{x};{{\bm{\theta}}_{1}})}{(1-\pi)p_{0}(\mathbf{x};{{\bm{\theta}}_{0}})+\pi
p_{1}(\mathbf{x};{{\bm{\theta}}_{1}})}$ is continuous in
$({\bm{\theta}}_{0},{\bm{\theta}}_{1},\pi)$, as the parameters’ estimates
converge to the true values, the error of the plug-in classifier also
converges to that of the Bayes decision rule. However, not only no optimality
condition can be stated in general for the plug-in classifier relying on MLEs
[33], obtaining such estimates might also be computationally prohibitive for
system with complex likelihood functions.

##  III Proposed Solution: HyPhyLearn

The main deciding factor in superiority of a solution for the problem setup
introduced in Section II is the extent to which it exploits the available
information, i.e., training data and the parametric statistical models. In
particular, the plug-in classifiers tend not to exploit this information in
the most optimal fashion as performance of the parameter estimation procedures
can be curbed due to the complexity of the underlying models and lack of the
corresponding ground-truth labels. We instead propose a novel hybrid
classification method to make use of the available information in learning-
based classifiers, which are powerful tools for finding discriminative feature
spaces. Specifically, our proposed solution relies on the parametric models to
generate synthetic data and incorporate them with the training data in a
classifier that makes use of adversarial training between NNs. Next, we
describe the various steps of the proposed solution that is termed HyPhyLearn
in detail.

Step $\mathbf{1}$—Imperfect labeling: As the available data are not assumed
completely labeled in our problem setup, the first step in our solution deals
with assigning labels to the unlabeled data samples in $\mathcal{D}_{r}$. This
involves a clustering step that partitions the dataset $\mathcal{D}_{r}$ into
$C$ distinct groups. Then, the groups are labeled using the available
$N_{r,l}$ labels. For example, a label can be assigned to a group based on the
number of labeled training data it includes from each behavior; If the
majority of such samples corresponds to the $i$th behavior, the group is
labeled as $i$. Subsequently, we refer to a group assigned with the label $i$
by $\mathcal{D}_{r,i}$ for $i=0,\dots,C-1$. Denoting this imperfect labeling
process by $g(\mathbf{x}):\mathcal{X}\rightarrow\\{0,\dots,C-1\\}$, a non-
trivial labeling error over $\mathcal{D}_{r}$ is associated with
$g(\mathbf{x})$ that can be computed via
$e_{r}=\frac{1}{N_{r}}\sum_{n=1}^{N_{r}}\mathds{1}_{\\{g(\mathbf{x}_{r,n})\neq
y_{r,n}\\}}$. In the remainder of this paper, we refer to the number of
samples in the cluster labeled as $i$ by $N_{r,i}$. The function
$g(\mathbf{x})$ may be obtained based on any one of the simple clustering
algorithms from the ML literature, such as the Gaussian mixture model [37], or
it may be a decision rule obtained based on the statistical analysis of the
parametric models. For instance, for the problem of channel spoofing
detection, a hypothesis test is proposed in [5] that assigns labels to
unlabeled samples based on their similarity, measured in terms of the
Euclidean distance, to a reference data sample.

Step $\mathbf{2}$—Parameter estimation: Based on the labels assigned in Step
$1$ to the unlabeled data samples, we estimate the parameters of the physics-
based statistical models under each behavior. To this end, we utilize
$\mathcal{D}_{r,i}$ to estimate the parameter vector $\bm{\theta}^{*}_{i}$
corresponding to the $i$th behavior. Furthermore, the priors are estimated as
$\widehat{\pi}_{i}=N_{r,i}/N_{r}$. We note that the procedure for estimating
$\bm{\theta}^{*}_{i}$ depends on the available parametric models corresponding
to the $i$th behavior, i.e., $p_{i}(\mathbf{x};\bm{\theta}_{i})$. We recall
from our problem setup that the MLE, which is usually utilized for parameter
estimation purposes, might not be employed here due to the formidable
complexity of optimizing $p_{i}(\mathbf{x};{\bm{\theta}_{i}})$ over
${\bm{\theta}_{i}}$. Instead, a (necessarily) suboptimal estimator,
$T(\cdot)$, built upon either heuristics or optimization techniques like
alternate maximization (see Sections IV-C and V-B) could be utilized to
estimate the parameters as $\widehat{\bm{\theta}}_{i}=T(\mathcal{D}_{r,i})$
for all the behaviors. The parameter estimation performance is therefore
limited here due to both the suboptimality of $T(\cdot)$ and presence of the
mislabeled samples in $\mathcal{D}_{r,i},\forall i$.

Step $\mathbf{3}$—Forming a synthetic dataset: The paucity of available data
in our problem formulation seems to preclude utilization of a learning-based
classifier as part of the solution. However, we note that the available
physics-based statistical models, in the form of parametric PDFs, enable us to
generate synthetic data to augment the available data, and make it possible to
exploit the discriminative power of learning-based classifiers. Having access
to the estimated parameter $\widehat{\bm{\theta}}_{i}$ obtained in Step $2$,
we plug it in the available physics-based statistical model to obtain a PDF
$p_{i}(\mathbf{x};{\widehat{\bm{\theta}}}_{i})$ for the $i$th behavior. In
order to generate a synthetic dataset, we first sample $w$ from a categorical
distribution parameterized by
$\widehat{\bm{\pi}}=[\widehat{\pi}_{0},\dots,\widehat{\pi}_{C-1}]$ over the
sample space of $\\{0,\dots,C-1\\}$. Then, we sample a data point
$\mathbf{x}_{s,i}$ according to $\mathbf{x}_{s,i}\sim
p_{w}(\mathbf{x};{\widehat{\bm{\theta}}_{w}})$ with the associated label
$y_{s,i}=w$. Repeating this process $N_{s}$ number of times, we obtain a
synthetic dataset
$\mathcal{D}_{s}=\\{\mathbf{x}_{s,i},y_{s,i}\\}_{i=1}^{N_{s}}$ in which the
data samples are generated in a statistically independent fashion.

Step $\mathbf{4}$—Incorporating synthetic and training data in a learning-
based classifier: The synthetic data generated in Step $3$, besides retaining
essential information about the underlying physics-based statistical models,
enables us to utilize the discriminative power of learning-based classifiers.
However, the errors introduced during the labeling and the parameter
estimation steps that precede the synthetic data generation process incur a
mismatch between the distributions corresponding to the training and synthetic
datasets. This mismatch is bound to deteriorate the performance of a
classifier trained on the synthetic data alone, when utilized in a real-world
setting. Then the question is how a learning-based classifier can be trained
to alleviate this problem. For example, in the fine-tuning approach [2], a NN-
based classifier will be trained on the synthetic data first, and then,
training data are used to refine the weights of the corresponding NN. However,
we conjecture that such learning strategies that utilize the training and
synthetic data in the separate stages of training are not the best solution
here; rather, synthetic and training data should jointly be incorporated in a
learning-based classifier. To this end, inspired by the works in the domain-
adaptation literature and specifically feature space mapping [12], we propose
to map the synthetic and training data through a data-driven function
$M_{\bm{\psi}}:\mathcal{X}\rightarrow\mathcal{Z}$, which is parameterized by a
real vector $\bm{\psi}$, into a common feature space $\mathcal{Z}$.
Consequently, a classifier $h_{\bm{\phi}_{1}}(\mathbf{z})$, parameterized by
$\bm{\phi}_{1}$, which is trained on the synthetic data within the space
$\mathcal{Z}$ is expected to perform well on both training and synthetic data.
To this end, we choose $M_{\bm{\psi}}$ and $h_{\bm{\phi}_{1}}$ to be NNs,
which are powerful tools for finding discriminative features from a given
dataset. We discuss this step in detail in the following subsection.

HyPhyLearn: We now present our final solution as an algorithmic framework
composed of the aforementioned four steps. In a nutshell, HyPhyLearn generates
synthetic data based on the physics-based parametric statistical models and
utilizes them along with the available data in a learning-based classifier
powered from the adversarial training of the NNs (see the following
subsection). In order to train the NNs based on their specific loss functions,
described in the following subsection, we utilize the stochastic gradient
descent method [37] along with mini-batches consisting of random samples from
the training and synthetic datasets in an iterative manner. The details of the
whole process is presented in Algorithm 1.

1Input: Parametric models $p_{i}(\mathbf{x};{\bm{\theta}_{i}})$
($i=0,\dots,C-1$); Training dataset
$\mathcal{D}_{r}=\\{\mathbf{x}_{r,n}\\}_{n=1}^{N_{r}}$; learning rates
$\mu_{r_{1}}$, $\mu_{r_{2}}$, $\mu_{r_{3}}$; Number of training steps
$N_{tr}$; Mini-batch size $N_{b}<N_{r}$; Number of synthetic data samples
$N_{s}$ to be generated

2 Output: The mapping $M_{\bm{\psi}}(\cdot)$ and the classifier
$h_{\bm{\phi}_{1}}(\cdot)$, parameterized by the real vectors ${\bm{\psi}}$
and ${\bm{\phi}_{1}}$, respectively

// Step $1$ \- Imperfect labeling

3 $\\{\mathcal{D}_{r,0},\dots,\mathcal{D}_{r,C-1}\\}\leftarrow$ Applying
$g(\mathbf{x})$ to unlabeled samples

// Step $2$ \- Parameter estimation

4 $\widehat{\bm{\theta}}_{i}\leftarrow T_{i}(\mathcal{D}_{r,i})$,
$\widehat{\pi}_{i}\leftarrow\frac{\lvert\mathcal{D}_{r,i}\rvert}{N_{r}}$ for
$i=0,\dots,C-1$

// Step $3$ \- Forming a synthetic dataset

5 $p_{i}(\mathbf{x};{\widehat{\bm{\theta}}_{i}})\leftarrow$Plug
$\widehat{\bm{\theta}}_{i}$ in $p_{i}(\mathbf{x};{{{\theta}}_{i}})$ for
$i=0,\dots,C-1$

6 for  _$n=1$ to $N_{s}$_ do

// Choosing a behavior

7 $r\sim\text{unif}(0,1)$,
$w=\operatorname*{\arg\\!\min}_{k}\sum_{i=0}^{k-1}\widehat{\pi}_{i}\geq r$

// Synthetic data generation

8 $\mathbf{x}_{s,n}{\sim}p_{w}(\mathbf{x};{\widehat{\bm{\theta}}_{w}})$,
$y_{s,n}=w$

9 Add $\\{\mathbf{x}_{s,n},y_{s,n}\\}$ to $\mathcal{D}_{s}$

10 end for

// Step $4$ \- Training the learning-based classifier

11 for  _$n_{tr}=1$ to $N_{tr}$_ do

12 $\mathcal{D}_{r,b}\leftarrow$ $N_{b}$ random samples from
$\mathcal{D}_{r}$, $\mathcal{D}_{s,b}\leftarrow$ $N_{b}$ random samples from
$\mathcal{D}_{s}$

// Forward propagation via (12), (14)

13 $L_{s}\leftarrow\mathcal{L}_{s}(\bm{\psi},\bm{\phi}_{1}|\mathcal{D}_{s,b})$

14
$L_{c}\leftarrow\mathcal{L}_{c}(\bm{\psi},\bm{\zeta}|\mathcal{D}_{r,b},\mathcal{D}_{s,b})$

// Backward propagation

15 Computing gradients:
$\mathcal{G}_{s,\bm{\phi}_{1}}\leftarrow\nabla_{\bm{\phi}_{1}}{L}_{s}$,
$\mathcal{G}_{s,\bm{\psi}}\leftarrow\nabla_{\bm{\psi}}{L}_{s}$

16 Computing gradients:
$\mathcal{G}_{c,\bm{\zeta}}\leftarrow\nabla_{\bm{\zeta}}{L}_{c}$,
$\mathcal{G}_{c,\bm{\psi}}\leftarrow\nabla_{\bm{\psi}}{L}_{c}$

// Update network parameters via (15)

17
$\bm{\psi}\leftarrow\bm{\psi}-\mu_{r_{1}}(\mathcal{G}_{s,\bm{\psi}}-\mathcal{G}_{c,\bm{\psi}})$,
$\bm{\phi}_{1}\leftarrow\bm{\phi}_{1}-\mu_{r_{2}}\mathcal{G}_{s,\bm{\phi}_{1}}$,
$\bm{\zeta}\leftarrow\bm{\zeta}-\mu_{r_{3}}\mathcal{G}_{c,\bm{\zeta}}$

18 end for

Algorithm 1 HyPhyLearn

###  III-A Incorporating synthetic and training data in a learning-based
classifier for HyPhyLearn

To elaborate further on Step $4$, we first denote the distributions
corresponding to the real and synthetic data as
$p_{{{\bm{\theta}^{*}}}}(\mathbf{x})=\sum_{i=0}^{C-1}\pi_{i}p_{i}(\mathbf{x};{{{\bm{\theta}}}^{*}_{i}})$
and
$p_{\widehat{\bm{\theta}}}(\mathbf{x})=\sum_{i=0}^{C-1}\widehat{\pi}_{i}p_{i}(\mathbf{x};{{\widehat{\bm{\theta}}}_{i}})$,
respectively. We refer to $p_{{{\bm{\theta}^{*}}}}(\mathbf{x})$ and
$p_{\widehat{\bm{\theta}}}(\mathbf{x})$ as the true and estimated
distributions, respectively. For each distribution, applying the mapping
$M_{\bm{\psi}}(\cdot)$ to the input space $\mathcal{X}$ would induce a
distribution over the feature space $\mathcal{Z}$. Specifically, we denote the
mapping of the true distribution $p_{{{\bm{\theta}^{*}}}}(\mathbf{x})$ to
$\mathcal{Z}$ by $p_{\bm{\psi},{\bm{\theta}^{*}}}(\mathbf{z})$, where
$\mathbf{z}=M_{\bm{\psi}}(\mathbf{x})$, $\mathbf{x}\sim
p_{{{\bm{\theta}^{*}}}}(\mathbf{x})$. Assuming that $\mathcal{X}$ and
$\mathcal{Z}$ are topological spaces, for any $\mathcal{A}\subset\mathcal{Z}$
the probability of $\mathcal{A}$ in space $\mathcal{Z}$ is

$\displaystyle\mathbb{P}_{\mathbf{z}}[{\mathcal{A}}]\overset{\bigtriangleup}{=}\mathbb{P}_{\mathbf{x}}\big{[}M^{-1}_{\bm{\psi}}({\mathcal{A}})\big{]}=\sum_{i=0}^{C-1}\pi_{i}\int_{M^{-1}_{\bm{\psi}}({\mathcal{A}})}p_{i}(\mathbf{x};{\bm{\theta}^{*}_{i}})d\mathbf{x},$
(3)

where the pre-image $M^{-1}_{\bm{\psi}}({\mathcal{A}})$ belongs to the Borel
$\sigma$-algebra over $\mathcal{X}$. Subsequently, the probability of error
corresponding to a classifier $h_{\bm{\phi}_{1}}(\mathbf{z})$, parameterized
by a real vector $\bm{\phi}_{1}$, with respect to the mapping of the true
distribution to the $\mathcal{Z}$ space is computed via

$\displaystyle\mathbb{P}_{\bm{\psi},\bm{\theta}^{*}}[{e}_{{\bm{\phi}}_{1}}]=\sum_{i=0}^{C-1}\pi_{i}\int_{\mathcal{Z}}p_{\bm{\psi},\bm{\theta}^{*}_{i}}(\mathbf{z})\mathds{1}_{\\{h_{{\bm{\phi}}_{1}}(\mathbf{z})\neq
i\\}}(\mathbf{z})d\mathbf{z},$ (4)

where the dependence of $\mathbb{P}$ on $\pi_{i}$’s is suppressed for
notational simplicity. Similarly, mapping of the estimated distribution to the
space $\mathcal{Z}$ is characterized by a distribution denoted by
$p_{\bm{\psi},{\widehat{\bm{\theta}}}}(\mathbf{z})$. Furthermore, the
probability of error for a classifier $h_{\bm{\phi}_{1}}(\mathbf{z})$ with
respect to $p_{\bm{\psi},{\widehat{\bm{\theta}}}}(\mathbf{z})$ can be computed
similar to (4), which we refer to as
$\mathbb{P}_{\bm{\psi},\widehat{\bm{\theta}}}[{e}_{\bm{\phi}}]$.

Our main goal is to learn a map $M_{\bm{\psi}}(\cdot)$ and a classifier
$h_{\bm{\phi}_{1}}(\mathbf{z})$ in a way that the probability of error of
$h_{\bm{\phi}_{1}}(\mathbf{z})$ with respect to the mapping of the true
distribution to $\mathcal{Z}$, i.e.,
$\mathbb{P}_{{\bm{\psi},{\bm{\theta}^{*}}}}[{e}_{\bm{\phi}_{1}}]$, is small.
To this end, we repurpose theories from the domain-adaptation literature in
the following to obtain an upper bound on
$\mathbb{P}_{{\bm{\psi},{\bm{\theta}^{*}}}}[{e}_{\bm{\phi}_{1}}]$, which leads
to explicit loss functions for the joint learning of $M_{\bm{\psi}}$ and
$h_{\bm{\phi}_{1}}(\mathbf{z})$ using both the training and synthetic
datasets. Specifically, it is desired for the mapping $M_{\bm{\psi}}(\cdot)$
from $\mathcal{X}$ to $\mathcal{Z}$ to transform the true and estimated
distributions in a way that $p_{{\bm{\psi},{\bm{\theta}^{*}}}}(\mathbf{z})$
and $p_{\bm{\psi},{\widehat{\bm{\theta}}}}(\mathbf{z})$, which are defined in
the feature space $\mathcal{Z}$, are similar. Mathematically, this similarity
should be measured in terms of a distance metric. However, as there are only a
limited number of samples available from
$p_{{\bm{\psi},{\bm{\theta}^{*}}}}(\mathbf{z})$, we need to be able to
approximate this distance from a finite number of samples. We expand further
on this idea by primarily focusing on binary classification in this section,
although the results are extendable to the classification task in general. We
begin with the following distance definitions.

######  Definition 1.

For a family of binary-valued functions
$\mathcal{H}_{\bm{\Phi}}=\\{h_{\bm{\phi}}:\mathcal{Z}\rightarrow\\{0,1\\}\\}$,
in which every member $h_{\bm{\phi}}\in\mathcal{H}_{\bm{\Phi}}$ is
parameterized by a real vector $\bm{\phi}\in\bm{\Phi}$, and the set
${A}_{\phi}=\\{\mathbf{z}|h_{\phi}(\mathbf{z})=1,\mathbf{z}\in\mathcal{Z}\\}$,
the $\mathcal{A}_{\Phi}$-distance between
$p_{{\bm{\psi},{\bm{\theta}^{*}}}}(\mathbf{z})$ and
$p_{\bm{\psi},{\widehat{\bm{\theta}}}}(\mathbf{z})$ is defined as

$d_{\mathcal{A}_{\Phi}}\big{(}p_{{\bm{\psi},{\bm{\theta}^{*}}}}(\mathbf{z}),p_{\bm{\psi},{\widehat{\bm{\theta}}}}(\mathbf{z})\big{)}\overset{\bigtriangleup}{=}2\sup_{h_{\phi}\in\mathcal{H}_{\Phi}}\bigg{|}\int_{A_{\phi}}(p_{{\bm{\psi},{\bm{\theta}^{*}}}}(\mathbf{z})-p_{\bm{\psi},{\widehat{\bm{\theta}}}}(\mathbf{z}))d\mathbf{z}\bigg{|}.$
(5)

Similarly, for
${B}_{\phi_{1},\phi_{2}}=\\{\mathbf{z}|h_{\phi_{1}}(\mathbf{z})\neq
h_{\phi_{2}}(\mathbf{z}),\mathbf{z}\in\mathcal{Z}\\}$, the
$\mathcal{B}_{\Phi}$-distance refers to 333Similar to the total variation
distance, it can be readily verified that $d_{\mathcal{A}_{\Phi}}$ and
$d_{\mathcal{B}_{\Phi}}$ are also distance metrics.

$\displaystyle
d_{\mathcal{B}_{\Phi}}\big{(}p_{{\bm{\psi},{\bm{\theta}^{*}}}}(\mathbf{z})$
$\displaystyle,p_{\bm{\psi},{\widehat{\bm{\theta}}}}(\mathbf{z})\big{)}\overset{\bigtriangleup}{=}$
$\displaystyle
2\sup_{h_{\phi_{1}},h_{\phi_{2}}\in\mathcal{H}_{\Phi}}\bigg{|}\int_{B_{\phi_{1},\phi_{2}}}(p_{{\bm{\psi},{\bm{\theta}^{*}}}}(\mathbf{z})-p_{\bm{\psi},{\widehat{\bm{\theta}}}}(\mathbf{z}))d\mathbf{z}\bigg{|}.$
(6)

The $\mathcal{A}_{\Phi}$-distance is also referred to via other names like
$A$-distance and $\mathcal{H}$-distance in [38, 23]. By looking at the
following extreme choices of $\mathcal{H}_{\bm{\Phi}}$, these distances are
clearly a function of richness of the class $\mathcal{H}_{\Phi}$. For a very
restrictive choice of only constant functions, i.e.,
$\mathcal{H}_{\Phi}=\\{h_{\bm{\phi}}|h_{\bm{\phi}}(\mathbf{z})=0,\forall\mathbf{z}\\}\bigcup\\{h_{\bm{\phi}}|h_{\bm{\phi}}(\mathbf{z})=1,\forall\mathbf{z}\\}$,
$d_{\mathcal{A}_{\Phi}}$ is always zero as the only possible choice for
$A_{\bm{\phi}}$ is either the empty set or $\mathcal{Z}$. On the other hand,
for $\mathcal{H}_{\Phi}=\\{h_{\bm{\phi}}|h_{\bm{\phi}}(\mathbf{z})=0\
\text{or}\ h_{\bm{\phi}}(\mathbf{z})=1,\forall\mathbf{z}\\}$, which represents
all the binary functions, $d_{\mathcal{A}_{\Phi}}$ is identical to definition
of the total variation distance [39] as the $\sup$ in (5) will effectively be
over the $\sigma$-algebra of subsets of the $\mathcal{Z}$ space. This
dependence of $d_{\mathcal{A}_{\Phi}}$ on the underlying family of functions
makes it possible to obtain an expression for the
${\mathcal{A}_{\Phi}}$-distance based on the finite set of samples from each
distribution. Specifically, consider two sets
$\mathcal{Z}_{{\bm{\psi},{\bm{\theta}^{*}}}}=\\{\mathbf{z}_{r,i}\\}_{i=1}^{N_{r}}$
and
$\mathcal{Z}_{\bm{\psi},{\widehat{\bm{\theta}}}}=\\{\mathbf{z}_{s,i}\\}_{i=1}^{N_{s}}$
sampled from the distributions $p_{{\bm{\psi},{\bm{\theta}^{*}}}}(\mathbf{z})$
and $p_{\bm{\psi},{\widehat{\bm{\theta}}}}(\mathbf{z})$ in an i.i.d. fashion,
respectively. In this case, for a family $\mathcal{H}_{\bm{\Phi}}$ that
satisfies the condition that if $h_{\phi}\in\mathcal{H}_{\Phi}$ then
$1-h_{\phi}\in\mathcal{H}_{\Phi}$, the $\mathcal{A}_{\Phi}$-distance can be
approximated from $\mathcal{Z}_{{\bm{\psi},{\bm{\theta}^{*}}}}$ and
$\mathcal{Z}_{\bm{\psi},{\widehat{\bm{\theta}}}}$ using [38]

$\displaystyle\widehat{d}_{\mathcal{A}_{\Phi}}(\mathcal{Z}_{{\bm{\psi},{\bm{\theta}^{*}}}},\mathcal{Z}_{\bm{\psi},{\widehat{\bm{\theta}}}})=$
$\displaystyle
2\bigg{(}1-\inf_{h_{\phi}\in\mathcal{H}_{\Phi}}\Big{(}\frac{1}{N_{r}}\sum_{i=1}^{N_{r}}\mathds{1}_{\\{h_{\phi}(\mathbf{z}_{r,i})=0\\}}+\frac{1}{N_{s}}\sum_{i=1}^{N_{s}}\mathds{1}_{\\{h_{\phi}({\mathbf{z}}_{s,i})=1\\}}\Big{)}\bigg{)}.$
(7)

As the bound on
$\mathbb{P}_{{\bm{\psi},{\bm{\theta}^{*}}}}[{e}_{\bm{\phi}_{1}}]$ should be
obtained based on a finite number of training and synthetic samples, it is
then of interest to see how far $\widehat{d}_{\mathcal{A}_{\Phi}}$ is from
$d_{\mathcal{A}_{\Phi}}$. To answer this question, one needs to rely on a
measure of complexity for a given class of functions such as the VC dimension
[36] and Rademacher complexity [40]. As we have chosen the mapping function
$M_{\bm{\psi}}(\mathbf{x})$ and the classifier $h_{\bm{\phi}_{1}}(\mathbf{z})$
to be NNs, we present the results based on the Rademacher complexity defined
as follows, which can be computed for certain classes of neural networks in a
closed-form fashion [40].

######  Definition 2.

Let $\mathcal{Z}_{1}=\\{\mathbf{z}_{i}\\}_{i=1}^{N}$ be a set of i.i.d.
samples drawn from a distribution $p(\mathbf{z})$ that is supported on
$\mathcal{Z}$. For $\mathcal{H}_{\Phi}$, a family of real-valued functions
over $\mathcal{Z}$, the empirical Rademacher complexity of
$\mathcal{H}_{\Phi}$, given a dataset $\mathcal{Z}_{1}$, is defined as

$R_{\mathcal{Z}_{1}}(\mathcal{H}_{\Phi})\overset{\bigtriangleup}{=}\underset{\underset{i=1,\dots,N}{\sigma_{i}\sim\\{-1,+1\\}}}{\mathbb{E}}\Bigg{[}\sup_{h_{\phi}\in\mathcal{H}_{\Phi}}\Bigg{(}\frac{1}{N}\sum_{i=1}^{N}\sigma_{i}h_{\phi}(\mathbf{z}_{i})\Bigg{)}\Bigg{]},$
(8)

where the expectation is over all the $\sigma_{i}$’s, each taking a binary
value with equal probability.

######  Lemma 2 ([40]).

Consider a family of functions
$\mathcal{H}_{\bm{\Phi}}=\\{h_{\bm{\phi}}:\mathcal{Z}\rightarrow\\{0,1\\}\\}$
and a distribution $p(\mathbf{z})$ over $\mathcal{Z}$. For a set
$\mathcal{Z}_{1}=\\{\mathbf{z}_{i}\\}_{i=1}^{N}$ of $N$ i.i.d. samples from
$p(\mathbf{z})$ and any $0<\delta<1$, the following holds $\forall
h_{\bm{\phi}}\in\mathcal{H}_{\Phi}$ with probability at least $1-\delta$:

$\displaystyle\mathbb{E}_{\mathbf{z}\sim
p(\mathbf{z})}[h_{\bm{\phi}}(\mathbf{z})]\leq\frac{1}{N}\sum_{i=1}^{N}h_{\bm{\phi}}(\mathbf{z}_{i})+2R_{\mathcal{Z}_{1}}(\mathcal{H}_{\Phi})+3\sqrt{\frac{\log(2/\delta)}{2N}}.$
(9)

Now, the difference between $d_{\mathcal{A}_{\Phi}}$ and
$\widehat{d}_{\mathcal{A}_{\Phi}}$ can be bounded in terms of the complexity
of the underlying family of functions and the number of available samples as
stated in the following lemma.

######  Lemma 3.

Let
$\mathcal{Z}_{{\bm{\psi},{\bm{\theta}^{*}}}}=\\{\mathbf{z}_{r,i}\\}_{i=1}^{N_{r}}$
and
$\mathcal{Z}_{\bm{\psi},{\widehat{\bm{\theta}}}}=\\{\mathbf{z}_{s,i}\\}_{i=1}^{N_{s}}$
be sets of i.i.d. samples corresponding to the distributions
$p_{{\bm{\psi},{\bm{\theta}^{*}}}}(\mathbf{z})$ and
$p_{\bm{\psi},{\widehat{\bm{\theta}}}}(\mathbf{z})$ on the space
$\mathcal{Z}$, respectively. Then, for any $0<\delta<1$ and a family of
functions
$\mathcal{H}_{\bm{\Phi}}=\\{h_{\bm{\phi}}:\mathcal{Z}\rightarrow\\{0,1\\}\\}$,
we have

$\displaystyle d_{\mathcal{A}_{\Phi}}$
$\displaystyle\big{(}p_{{\bm{\psi},{\bm{\theta}^{*}}}}(\mathbf{z}),p_{\bm{\psi},{\widehat{\bm{\theta}}}}(\mathbf{z})\big{)}\leq\widehat{d}_{\mathcal{A}_{\Phi}}(\mathcal{Z}_{{\bm{\psi},{\bm{\theta}^{*}}}},\mathcal{Z}_{\bm{\psi},{\widehat{\bm{\theta}}}})+2R_{\mathcal{Z}_{{\bm{\psi},{\bm{\theta}^{*}}}}}(\mathcal{H}_{\Phi})$
$\displaystyle+2R_{\mathcal{Z}_{\bm{\psi},{\widehat{\bm{\theta}}}}}(\mathcal{H}_{\Phi})+3\sqrt{(\log{2/\delta})/2N_{r}}+3\sqrt{(\log{2/\delta})/2N_{s}}$
(10)

with probability at least $1-\delta$.

###### Proof.

See Appendix A. ∎

The above lemma enables us to bound the $\mathcal{A}_{\bm{\Phi}}$ distance
between two distributions in terms of the collected samples from each.
Equipped with this result, we are able to bound the probability of error
$\mathbb{P}_{{\bm{\psi},{\bm{\theta}^{*}}}}[{e}_{\bm{\phi}_{1}}]$ via the
following theorem.

######  Theorem 1.

Assume that the training and synthetic datasets are mapped into the feature
space $\mathcal{Z}$ through the mapping function $M_{\bm{\psi}}(\mathbf{x})$,
with the resulting samples denoted by
$\mathcal{Z}_{{\bm{\psi},{\bm{\theta}^{*}}}}=\\{\mathbf{z}_{r,i}\\}_{i=1}^{N_{r}}$
and
$\mathcal{Z}_{\bm{\psi},{\widehat{\bm{\theta}}}}=\\{\mathbf{z}_{s,i}\\}_{i=1}^{N_{s}}$,
respectively. Then, for any $0<\delta<1$ and a family of functions
$\mathcal{H}_{\bm{\Phi}}:\mathcal{Z}\rightarrow\\{0,1\\}$,
$\mathbb{P}_{{\bm{\psi},{\bm{\theta}^{*}}}}[{e}_{\bm{\phi}_{1}}],\forall
h_{\bm{\phi}_{1}}\in\mathcal{H}_{\bm{\Phi}}$ is bounded by

$\displaystyle\mathbb{P}$
${}_{{\bm{\psi},{\bm{\theta}^{*}}}}[{e}_{\bm{\phi}_{1}}]\leq\mathbb{P}_{\bm{\psi},{\widehat{\bm{\theta}}}}[{e}_{\bm{\phi}_{1}}]+\frac{1}{2}\widehat{d}_{\mathcal{A}_{\Phi}}(\mathcal{Z}_{{\bm{\psi},{\bm{\theta}^{*}}}},\mathcal{Z}_{\bm{\psi},{\widehat{\bm{\theta}}}})+R_{\mathcal{Z}_{{\bm{\psi},{\bm{\theta}^{*}}}}}(\mathcal{H}_{\bm{\Phi}})+$
$\displaystyle
R_{\mathcal{Z}_{\bm{\psi},{\widehat{\bm{\theta}}}}}(\mathcal{H}_{\bm{\Phi}})+\frac{3}{2}\sqrt{(\log{2/\delta})/2N_{r}}+\frac{3}{2}\sqrt{(\log{2/\delta})/2N_{s})}.$
(11)

###### Proof.

See Appendix B. ∎

The above theorem bounds the probability of error with respect to
$p_{\bm{\psi},\bm{\theta}^{*}}(\mathbf{z})$ associated with a classifier
$h_{\bm{\phi}_{1}}(\mathbf{z})$ in terms of the quantities that do not depend
on the the unknown true parameters $\bm{\theta}^{*}$. As our primary goal is
to make $\mathbb{P}_{{\bm{\psi},{\bm{\theta}^{*}}}}[{e}_{\bm{\phi}_{1}}]$ as
small as possible, the mapping function $M_{\bm{\psi}}(\mathbf{x})$ and the
classifier $h_{\bm{\phi}_{1}}(\mathbf{z})$ should be chosen in a way to
minimize the above upper bound. We note that the complexity related terms in
the above bound are fixed for a chosen family of the functions and the bound
is primarily controlled by the first two terms. In other words,
$M_{\bm{\psi}}(\mathbf{x})$ and $h_{\bm{\phi}_{1}}(\mathbf{z})$ should be
chosen such that the probability of classification error with respect to the
mapping of the estimated distribution in the $\mathcal{Z}$ space, i.e.,
$\mathbb{P}_{\bm{\psi},{\widehat{\bm{\theta}}}}[{e}_{\bm{\phi}_{1}}]$, and the
approximated $\mathcal{A}_{\bm{\Phi}}$-distance between the synthetic and
training datasets are minimized simultaneously. To achieve this goal, we
restrict ourselves to $M_{\bm{\psi}}(\mathbf{x})$ and
$h_{\bm{\phi}_{1}}(\mathbf{z})$ that correspond to NNs that are trained to
minimize a loss function in accordance with the first two terms of the above
bound. One can efficiently solve the resulting optimization problem via the
stochastic gradient descent method as described in the following.

Joint learning of the feature map and the classifier: In terms of specifics,
we assume $M_{\bm{\psi}}(\mathbf{x})$ and $h_{\bm{\phi}_{1}}(\mathbf{z})$
belong to the class of feed-forward (deep) NNs whose parameters, i.e.,
$\bm{\psi}$ and $\bm{\phi}_{1}$, correspond to the weights and biases of each
network. The input and output layers of the NNs corresponding to
$M_{\bm{\psi}}$ have $n_{x}$ and $n_{z}$ number of neurons, respectively,
which denote the dimensions of the spaces $\mathcal{X}$ and $\mathcal{Z}$,
respectively. We note that $n_{x}$ is chosen according to the length of the
observation vector as part of the problem formulation, while $n_{z}$ can be
picked as a hyper-parameter to facilitate the training process. Subsequently,
the input layer of $h_{\bm{\phi}_{1}}$ has $n_{z}$ neurons while its output
layer contains $C$ neurons whose activation function is chosen to be the
softmax function $\bm{\sigma}(\mathbf{z})$ for which the $i$th element is
given by $\frac{e^{\mathbf{z}[i]}}{\sum_{i=1}^{n_{z}}e^{\mathbf{z}[i]}}$. In
this way, the $i$th component of the vector
${\mathbf{y}}_{\bm{\psi},\bm{\phi}_{1},\mathbf{x}}\overset{\bigtriangleup}{=}h_{\bm{\phi}_{1}}\big{(}M_{\bm{\psi}}(\mathbf{x})\big{)}$
denotes the probability that the classifier assigns to the input $\mathbf{x}$
that it belongs to the $i$th class for $i=0,\dots,C-1$. Consequently, the
averaged cross-entropy loss, minimizing of which leads to minimizing the
classification error associated with $h_{\bm{\phi}_{1}}$, over the synthetic
dataset $\mathcal{D}_{s}$ equals

$\displaystyle\mathcal{L}_{s}({\bm{\psi},\bm{\phi}_{1}}|\mathcal{D}_{s})=\frac{1}{n_{s}}\sum_{n=1}^{n_{s}}\sum_{i=1}^{C}\mathbf{l}_{s,n}[n]\log{{\mathbf{y}}_{\bm{\psi},\bm{\phi}_{1},\mathbf{x}_{s,n}}[n]},$
(12)

where $\mathbf{l}_{s,n}=\mathbf{e}_{C}({y}_{s,n})$ denotes the one-hot encoded
version of the label ${y}_{s,n}$ corresponding to the $n$th sample. Regrading
the computation of $\widehat{d}_{\mathcal{A}_{\Phi}}$ between the two sets
$\mathcal{Z}_{{\bm{\psi},{\bm{\theta}^{*}}}}$ and
$\mathcal{Z}_{\bm{\psi},{\widehat{\bm{\theta}}}}$, it is suggested by the
authors in [38, 23] that the classification accuracy corresponding to a
classifier trained to distinguish between the samples from the two sets can be
used as a surrogate for the $\inf$ part in (III-A) that can be readily
computed during the learning process. To train such classifier, we consider a
NN $d_{\bm{\zeta}}$ with $n_{z}$ input neurons and $2$ output neurons with
softmax activation function, which is trained to distinguish between
$\mathcal{Z}_{{\bm{\psi},{\bm{\theta}^{*}}}}$ and
$\mathcal{Z}_{\bm{\psi},{\widehat{\bm{\theta}}}}$ labeled as $0$ and $1$,
respectively. Consequently, by defining a two-dimensional vector
$d_{\bm{\psi},{\bm{\zeta}},\mathbf{x}}\overset{\bigtriangleup}{=}d_{\bm{\zeta}}\big{(}M_{\bm{\psi}}(\mathbf{x})\big{)}$,
the $\widehat{d}_{\mathcal{A}_{\Phi}}$ term can be approximated by the cross-
entropy loss associated with $d_{\bm{\zeta}}$ as follows:

$\displaystyle\mathcal{L}_{d}({\bm{\psi},\bm{\zeta}}|\mathcal{D}_{s},\mathcal{D}_{r})=2\big{(}1-2\mathcal{L}_{c}({\bm{\psi},\bm{\zeta}}|\mathcal{D}_{s},\mathcal{D}_{r})\big{)},$
(13)
$\displaystyle\mathcal{L}_{c}({\bm{\psi},\bm{\zeta}}|\mathcal{D}_{s},\mathcal{D}_{r})=\frac{1}{n_{r}}\sum_{i=1}^{n_{r}}\log{d_{\bm{\psi},{\bm{\zeta}},\mathbf{x}_{r,n}}[1]}+$
$\displaystyle\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\ \ \ \ \
\frac{1}{n_{s}}\sum_{n=1}^{n_{s}}\log{d_{\bm{\psi},{\bm{\zeta}},\mathbf{x}_{s,n}}[2]}.$
(14)

Now, using Theorem 1 the training goal for the constituent NNs is set to
simultaneously minimize the classification error corresponding to the
synthetic data and the distance between the real and synthetic data, both
measured in the mapped space $\mathcal{Z}$. Specifically, the NNs
$M_{\bm{\psi}}$ and $h_{\bm{\phi}_{1}}$ should be trained to minimize the sum
of the losses in (12) and (13), while the classifier $d_{\bm{\zeta}}$ is
trained to minimize (14). As $M_{\bm{\psi}}$ is trained to maximize
$\mathcal{L}_{c}({\bm{\psi},\bm{\zeta}})$ despite $d_{\bm{\zeta}}$’s goal to
minimize $\mathcal{L}_{c}({\bm{\psi},\bm{\zeta}})$, the learning process
involves adversarial training between these two NNs. Based on the approach
taken in [23] for adversarial training in the context of domain adaptation, we
train the above three NNs for finding the saddle points $\widehat{\bm{\psi}}$,
$\widehat{\bm{\phi}}_{1}$ and $\widehat{\bm{\zeta}}$, such that

$\displaystyle\widehat{\bm{\psi}},\widehat{\bm{\phi}}_{1}=\operatorname*{\arg\\!\min}_{{\bm{\psi}},{\bm{\phi}}_{1}}\mathcal{L}_{t}({\bm{\psi},\bm{\phi}_{1},\widehat{\bm{\zeta}}}|\mathcal{D}_{s},\mathcal{D}_{r}),$
(15)
$\displaystyle\widehat{\bm{\zeta}}=\operatorname*{\arg\\!\min}_{{\bm{\zeta}}}-\mathcal{L}_{t}({\widehat{\bm{\psi}},\widehat{\bm{\phi}}_{1},\bm{\zeta}}|\mathcal{D}_{s},\mathcal{D}_{r}),$
(16)
$\displaystyle\mathcal{L}_{t}({\bm{\psi},\bm{\phi}_{1},\bm{\zeta}}|\mathcal{D}_{s},\mathcal{D}_{r})=\mathcal{L}_{s}({\bm{\psi},\bm{\phi}_{1}}|\mathcal{D}_{s})+\mathcal{L}_{d}({\bm{\psi},\bm{\zeta}}|\mathcal{D}_{s},\mathcal{D}_{r}),$
(17)

which can be achieved by utilizing the stochastic gradient descent algorithm
for each minimization task. To this end, the minimization is performed over
the NN’s parameters, $\bm{\psi}$, $\bm{\phi}_{1}$ and $\bm{\zeta}$, that are
real vectors whose dimensions are determined by the architecture of each
network.

###  III-B An illustrative example: The case of two-dimensional Gaussian data

Next, we show how the learning-based classifier in Section III-A performs on
simple training and synthetic datasets in an illustrative manner. To this end,
we consider a toy example where the true and estimated distributions are a
mixture of two bivariate Gaussian distributions with full-rank covariance
matrix each. In particular, we focus on the problem of binary classification
where the distribution for the $i$th class is denoted by
$p_{i}(\mathbf{x};{{{\bm{\theta}}}^{*}_{i}})=\mathcal{N}(\bm{\mu}_{i},\bm{\Sigma})$
for $i=0,1$, $\bm{\mu}_{i}\in\mathbb{R}^{2\times 1}$,
$\bm{\Sigma}\in\mathbb{R}^{2\times 2}$, and equal priors. In order to
investigate the effect of mismatch between only mean parameters, the
corresponding estimated distributions are assumed to have the same covariance
but different means, i.e.,
$p_{i}(\mathbf{x};{\widehat{\bm{\theta}}}_{i})=\mathcal{N}(\widehat{\bm{\mu}}_{i},\bm{\Sigma})$
for $i=0,1$ and equal priors. For two multivariate Gaussian distributions, the
authors in [39] have proposed a bound for the corresponding total variation as
part of the following theorem.

(a) Samples corresponding to the true and estimated distributions in space
$\mathcal{X}$.

(b) Mapping of the samples via the function $M_{\bm{\psi}}$ to the space
$\mathcal{Z}.$

(c) Position of the means in the original space $\mathcal{X}$ and the space
$\mathcal{Z}$.

Figure 2: Visualization of the true and estimated distributions and their
mappings to the space $\mathcal{Z}$ for the case of $2$D Gaussian datasets.

######  Theorem 2 ([39]).

Consider two $d$-dimensional Gaussian distributions
$\mathcal{N}({{\bm{\mu}}}_{1},{\bm{\Sigma}}_{1})$ and
$\mathcal{N}({{\bm{\mu}}}_{2},{\bm{\Sigma}}_{2})$ where
${{\bm{\mu}}}_{1}\neq{{\bm{\mu}}}_{2}$ and $\bm{\Sigma}_{1}$ and
$\bm{\Sigma}_{2}$ are positive definite. Let
$\mathbf{v}={{\bm{\mu}}}_{1}-{{\bm{\mu}}}_{2}$ and $\mathbf{\Pi}$ be a
$d\times(d-1)$ matrix whose columns form a basis for the subspace orthogonal
to $\mathbf{v}$. Denote the eigenvalues of
$(\mathbf{\Pi}^{T}{\bm{\Sigma}}_{1}\mathbf{\Pi})^{-1}\mathbf{\Pi}^{T}{\bm{\Sigma}}_{2}\mathbf{\Pi}-\mathbf{I}_{d-1}$
by $\rho_{1},\dots,\rho_{d-1}$. Then, the total variation between the two
distribution can be bounded as

$\displaystyle\frac{1}{200}\leq\frac{TV(\mathcal{N}({{\bm{\mu}}}_{1},{\bm{\Sigma}}_{1}),\mathcal{N}({{\bm{\mu}}}_{2},{\bm{\Sigma}}_{2}))}{\min(1,V)}\leq\frac{9}{2}$
(18)

where
$V\stackrel{{\scriptstyle\text{def}}}{{=}}\max\Bigg{\\{}\frac{|\mathbf{v}^{T}({\bm{\Sigma}}_{1}-{\bm{\Sigma}}_{2})\mathbf{v}|}{\mathbf{v}^{T}{\bm{\Sigma}}_{1}\mathbf{v}},\frac{\mathbf{v}^{T}\mathbf{v}}{\sqrt{\mathbf{v}^{T}{\bm{\Sigma}}_{1}\mathbf{v}}},\sqrt{\sum_{i=1}^{d-1}\rho_{i}^{2}}\Bigg{\\}}.$

We note that a bound on total variation would also bound the
$\mathcal{A}_{\Phi}$ distance following the discussion after Definition 1.
Using the above result, we can bound the total variation distance between
$p_{i}(\mathbf{x};{{{\bm{\theta}}}^{*}_{i}})$ and
$p_{i}(\mathbf{x};{\widehat{\bm{\theta}}}_{i})$ as follows, which will provide
useful insights in the remainder of this section about the learning process
described in Section III-A.

######  Corollary 1.

For two Gaussian distributions $\mathcal{N}({{\bm{\mu}}}_{0},{\bm{\Sigma}})$
and $\mathcal{N}(\widehat{{{\bm{\mu}}}}_{0},{\bm{\Sigma}})$ with the same
positive definite covariance matrix ${\bm{\Sigma}}$, the corresponding total
variation is bounded from the above by
$\frac{9}{2}\min\Big{(}1,\frac{({{\bm{\mu}}}_{0}-\widehat{{{\bm{\mu}}}}_{0})^{T}({{\bm{\mu}}}_{0}-\widehat{{{\bm{\mu}}}}_{0})}{\sqrt{({{\bm{\mu}}}_{i}-\widehat{{{\bm{\mu}}}}_{0})^{T}{\bm{\Sigma}}({{\bm{\mu}}}_{0}-\widehat{{{\bm{\mu}}}}_{0})}}\Big{)}.$

Regarding the specific architecture for the NNs utilized in Section III-A, let
us now choose the mapping function $M_{\bm{\psi}}$ to be
$M_{\bm{\psi}}(\mathbf{x})=\mathbf{W}_{\bm{\psi},2}\mathbf{W}_{\bm{\psi},1}\mathbf{x}$
parameterized by $\bm{\psi}=\\{\mathbf{W}_{\bm{\psi},1}\in\mathbb{R}^{20\times
2},\mathbf{W}_{\bm{\psi},2}\in\mathbb{R}^{2\times 20}\\}$. In particular, we
have set the dimension of the space $\mathcal{Z}$ to $n_{z}=2$ in order to be
able to readily visualize it within the $2$D coordinate system. For each
$h_{\bm{\phi}_{1}}$ and $d_{\bm{\zeta}}$, we choose a two-layer NN with
softmax activation function. Specifically, for $h_{\bm{\phi}_{1}}$ we have
$h_{\bm{\phi}_{1}}\big{(}M_{\bm{\psi}}(\mathbf{x})\big{)}=\text{softmax}\big{(}\mathbf{V}_{\bm{\phi}_{1},2}(\mathbf{V}_{\bm{\phi}_{1},1}M_{\bm{\psi}}(\mathbf{x})+\mathbf{b}_{\bm{\phi}_{1},1})+\mathbf{b}_{\bm{\phi}_{1},2}\big{)}$
where
$\bm{\phi}_{1}=\\{\mathbf{V}_{\bm{\phi}_{1},1}\in\mathbb{R}^{20\times{2}},\mathbf{b}_{\bm{\phi}_{1},1}\in\mathbb{R}^{20},\mathbf{V}_{\bm{\phi}_{1},2}\in\mathbb{R}^{2\times
20},\mathbf{b}_{\bm{\phi}_{1},2}\in\mathbb{R}^{2}\\}$. Similarly,
$d_{\bm{\zeta}}$ is chosen to be
$d_{\bm{\zeta}}\big{(}M_{\bm{\psi}}(\mathbf{x})\big{)}=\text{softmax}\big{(}\mathbf{U}_{\bm{\phi}_{1},2}(\mathbf{U}_{\bm{\phi}_{1},1}M_{\bm{\psi}}(\mathbf{x})+\mathbf{b}_{\bm{\phi}_{1},1})+\mathbf{b}_{\bm{\phi}_{1},2}\big{)}$
for
$\bm{\zeta}=\\{\mathbf{U}_{\bm{\phi}_{1},1}\in\mathbb{R}^{20\times{2}},\mathbf{b}_{\bm{\phi}_{1},1}\in\mathbb{R}^{20},\mathbf{U}_{\bm{\phi}_{1},2}\in\mathbb{R}^{2\times
20},\mathbf{b}_{\bm{\phi}_{1},2}\in\mathbb{R}^{2}\\}$. Training of these NNs
involves finding the saddle points of (15) based on the available training and
synthetic datasets which would lead to the learning-based classifier
$h_{\bm{\phi}_{1}}$. We note that the above simple choice of the mapping
function maps $p_{i}(\mathbf{x};\mathbf{\theta}_{i}^{*}),\ i=0,\ 1$ to
Gaussian distributions in the $\mathcal{Z}$ space which allows us to utilize
Corollary 1 for analyzing the total variation distance between these mappings
in the following.

We now resort to numerical results for further illustration of this example.
To this end, we set ${{\bm{\mu}}}_{0}=[2.9,4.4]$, ${{\bm{\mu}}}_{1}=[5,6.4]$,
$\widehat{{{\bm{\mu}}}}_{0}=[2,3]$, $\widehat{{{\bm{\mu}}}}_{1}=[4,5]$ and
${\bm{\Sigma}}=\begin{bmatrix}0.15&0.11\\\ 0.11&0.15\end{bmatrix}$. Also, we
generate $n_{r}=40$ samples from the true distribution, while $n_{s}=2000$
samples are generated from the estimated distribution. The Figs. 2(a) and 2(b)
depict the samples from the true and estimated distributions and their mapping
through the function $M_{\bm{\psi}}$ into the $\mathcal{Z}$ space,
respectively. Furthermore, the positions of the means corresponding to the
samples from the real and estimated distributions in both space $\mathcal{X}$
and $\mathcal{Z}$ are illustrated in Fig. 2(c). An important observation in
relation to the Corollary 1 can be made by noting that the total variation
between $\mathcal{N}({{\bm{\mu}}}_{i},{\bm{\Sigma}})$ and
$\mathcal{N}(\widehat{{{\bm{\mu}}}}_{i},{\bm{\Sigma}})$ is bounded by the term
$||\mathbf{v}||\frac{\mathbf{e}_{v}^{T}\mathbf{e}_{v}}{\sqrt{\mathbf{e}_{v}^{T}{\bm{\Sigma}}\mathbf{e}_{v}}}$
where $\mathbf{v}=\bm{\mu}_{i}-\widehat{\bm{\mu}}_{i}$ and
$\mathbf{e}_{v}=\mathbf{v}/||\mathbf{v}||$. Assuming $\lambda_{1}$ and
$\lambda_{2}$ are eigenvalues of ${\bm{\Sigma}}$ with corresponding
eigenvectors $\mathbf{u}_{1}$ and $\mathbf{u}_{2}$ such that
$\lambda_{1}>\lambda_{2}$, it is straightforward to show that the maximal
value of
$\mathbf{e}_{v}^{T}{\bm{\Sigma}}\mathbf{e}_{v}=\lambda_{1}(\mathbf{u}_{1}^{T}\mathbf{e}_{v})+\lambda_{2}(\mathbf{u}_{2}^{T}\mathbf{e}_{v})$
is achieved when $\mathbf{e}_{v}\perp\mathbf{u}_{2}$. Therefore, for
$||\mathbf{v}||\frac{\mathbf{e}_{v}^{T}\mathbf{e}_{v}}{\sqrt{\mathbf{e}_{v}^{T}{\bm{\Sigma}}\mathbf{e}_{v}}}$
to be minimized $\mathbf{e}_{v}$ ought to be in the same direction of
$\mathbf{u}_{1}$ while $||\mathbf{v}||$ become minimum. Notably, Figs. 2(b)
and 2(c) highlight the fact that finding the saddle points in (15) in part
corresponds to mapping the datasets to a feature space $\mathcal{Z}$ that
satisfy both these two criteria.

##  IV Case study 1: Channel-based Spoofing Detection for Physical Layer
Security

We now present the first case study concerning channel spoofing detection,
which can be posed as a binary hypothesis testing problem. We obtain the
likelihood ratio test based on the time-variant channel models. Then, we
discuss parameter estimation procedures for the likelihood function
corresponding to channel frequency responses (CFRs), and show how our proposed
solution is applicable for achieving an enhanced spoofing detection system.

###  IV-A System Model

The problem of channel spoofing detection arises in a wireless communication
environment where a legitimate transmitter (Alice) is transmitting signals to
a legitimate receiver (Bob) in the presence of an adversary (Eve). Eve aims at
spoofing the Alice-Bob’s channel by using the Alice’s MAC address [5, 7].
Bob’s goal, in this setting, is to distinguish between the signals coming from
Alice and Eve based on the corresponding channel frequency responses (CFRs).

We envision the communication parties in a $5$G propagation setting relying on
MIMO-OFDM wideband communications where the number of antennas are set to
$N_{Tx}$ and $N_{Rx}$ at the transmitter (Tx) and the receiver (Rx),
respectively. We assume bob measures and stores channel frequency response
samples at $M$ tones, across an overall system bandwidth of $W$. We consider a
generalized time-variant channel model, where each measured frequency response
sample is made up of three components: $1)$ specular paths
($\overline{\mathbf{h}}$), $2)$ time-varying part ${\mathbf{d}}_{u}$, and $3)$
noise ${\mathbf{n}}$, which are complex vectors of size $M\times 1$. The
specular paths model the dominant portion of the channel which remains
unchanged within a coherence time. The time-varying part models the dense
multipath components which accounts for the diffuse scattering between two
transceivers. Finally, the noise part models the measurement noise. The
measured CFR at Bob at time $t=uT$, $u\in\mathbb{N}$, is denoted by
$\mathbf{h}_{k}$ which is a $M\times 1$ vector such that

$\displaystyle\mathbf{h}_{u}=\overline{\mathbf{h}}+{\mathbf{d}}_{u}+{\mathbf{n}}.$
(19)

We first introduce the dominant paths model suitable for MIMO-OFDM
communication, under a frequency-dependent array response [41]. For this
scenario, the $N_{Rx}\times N_{Tx}$ channel matrix associated with the $n$th
subcarrier ($n=1,\dots,N_{f}$) is expressed as

$\displaystyle\mathbf{H}[n]=\mathbf{A}_{R}[n]\mathbf{\Gamma}[n]\mathbf{A}^{H}_{T}[n].$
(20)

The antenna steering response vectors are defined as

$\displaystyle\mathbf{A}_{T}[n]=[\mathbf{a}_{T,n}(\psi_{T,0}),\dots,\mathbf{a}_{T,n}(\psi_{T,K-1})],$
(21)
$\displaystyle\mathbf{A}_{R}[n]=[\mathbf{a}_{Rx,n}(\psi_{R,0}),\dots,\mathbf{a}_{R,n}(\psi_{Rx,K-1})],$
(22)

where $K$ is the total number of dominant paths,
$\psi_{T,k}=\frac{2\pi}{\lambda_{n}}d\sin(\theta_{Tx,k})$,
$\lambda_{n}=c(NT_{s}+f_{c})/n$ denotes the signal bandwidth at the $n$th
subcarrier, and $d$ refers to the distance between two antenna elements. The
structure of the frequency-dependent antenna steering and response vectors
$\mathbf{a}_{T,n}(\psi_{T,K-1})$ and $\mathbf{a}_{R,n}(\psi_{R,K-1})$ depends
on the specific array structure. For the case of a uniform linear array (ULA)
which we consider in this work, we have

$\displaystyle\mathbf{a}_{T,n}(\psi_{T,k})=\frac{1}{N_{Tx}}[e^{-j\frac{N_{Tx}-1}{2}\psi_{T,k}},\dots,e^{j\frac{N_{Tx}-1}{2}\psi_{T,k}}],$
(23)

Similarly, $\mathbf{a}_{Rx,n}(\psi_{Rx,k})$ can be defined for the receiver’s
antennas. The path gain matrix is obtained by

$\displaystyle\mathbf{\Gamma}[n]$ $\displaystyle=\sqrt{N_{Rx}N_{Tx}}$ (24)
$\displaystyle\text{diag}\bigg{\\{}\rho_{0}e^{-j2\pi
n\tau_{0}/(NT_{s})},\dots,\rho_{K-1}e^{-j2\pi
n\tau_{K-1}/(NT_{s})}\bigg{\\}},$ (25)

where $h_{k}$ and $\tau_{k}$ denotes the complex channel gain and delays of
the $k$th path while $T_{s}$ is the sampling interval. Then,
$\bar{\mathbf{h}}$ is defined as concatenation of the vectorized version of
$\mathbf{H}[n]$ for all the subcarriers $n=1,\dots,N$, i.e.,

$\displaystyle\bar{\mathbf{h}}=\big{[}\text{vec}\\{\mathbf{H}[1]\\}^{T},\dots,\text{vec}\\{\mathbf{H}[N]\\}^{T}\big{]}^{T},$
(26)

where $\text{vec}\\{\cdot\\}$ denotes the column vector operator. We denote
the associated parameters with the specular paths contribution,
$\overline{\mathbf{h}},$ which remains constant during a coherence time
$T_{c}$, via a $4K\times 1$ vector $\bm{\theta}_{sp}$ defined as

$\displaystyle\bm{\theta}_{sp}=[\bm{\psi}_{T},\bm{\psi}_{R},\bm{\tau},\bm{\rho}]^{T},$
(27)

where $\bm{\psi}_{T}=[\psi_{T,0},\dots,\psi_{T,K-1}]$,
$\bm{\psi}_{R}=[\psi_{R,0},\dots,\psi_{R,K-1}]$,
$\bm{\tau}=[\tau_{0},\dots,\tau_{K-1}]$ and
$\bm{\rho}=[\rho_{0},\dots,\rho_{K-1}]$.

For modeling the variable part of the channel we first assume that the wide-
sense stationary uncorrelated scattering (WSSUS) assumption holds, and use a
multipath tapped delay line,
$h(t,\tau)=\sum_{l=0}^{L-1}A_{l}(t)\delta(\tau-l\Delta\tau),$ to model the
impulse response at time $t$ between two antennas. Here, $A_{l}(t)$ and
$\Delta\tau=1/W$ denote the (complex) amplitude of the $l$th path and the
delay between two consecutive paths, respectively. Sampling the impulse
response at time $t=uT$, followed by taking the Fourier transform w.r.t.
$\tau$ would result in a vector $\mathbf{q}_{{u}}$ whose $n$th element is
denoted by

$\displaystyle\mathbf{q}_{u}[n]=\mathcal{F}\\{h(kT,\tau)\\}|_{f=f_{0}-W/2+n\Delta
f}=$ $\displaystyle\sum_{l=0}^{L-1}A_{u,l}e^{-j2\pi(f_{0}-W/2+n\Delta
f)l/W},n=1,\dots,N_{f},$ (28)

where $\Delta f$ and $A_{u,l}$ denotes the subcarrier width and the $l$th
channel gain at time $u$, respectively. Following the exponential decay model
which holds for the power delay profile of $\mathbf{q}_{u}$ based on various
experimental observations [5], $A_{u,l}$ is modeled with zero-mean Gaussian
distribution with variance
$\text{Var}(A_{u,l})=\alpha^{2}(1-e^{-2\pi\beta})e^{-2\pi\beta}$. Here,
$\alpha^{2}$ and $\beta$ denotes the average power and the normalized
coherence bandwidth, respectively. The distribution of $\mathbf{q}_{u}$ is
given in the following lemma.

######  Lemma 4.

The vector $\mathbf{q}_{u}$ has a multivariate Gaussian distribution
$\mathcal{CN}(\mathbf{0},\mathbf{R}_{\mathbf{q}})$ with a Toeplitz covarinace
matrix
$\mathbf{R}_{\mathbf{q}}=\text{toep}(\bm{\nu}_{\mathbf{q}},\bm{\nu}_{\mathbf{q}}^{H})$
assuming

$\displaystyle\bm{\nu}_{\mathbf{q}}\overset{\bigtriangleup}{=}\big{[}\kappa(\bm{\theta}_{\mathbf{q}},0),\kappa(\bm{\theta}_{\mathbf{q}},\frac{1}{M}),\dots,\kappa(\bm{\theta}_{\mathbf{q}},1-\frac{1}{M})\big{]},$
(29)

where
$\kappa(\bm{\theta}_{\mathbf{q}},m)\overset{\bigtriangleup}{=}\frac{\alpha^{2}(1-e^{-2\pi\beta})(1-e^{-2\pi
L(\beta-mj)})}{(1-e^{-2\pi(\beta-mj)})}$ and
$\bm{\theta}_{\mathbf{q}}=[\alpha^{2},\beta,L]$.

###### Proof.

See Appendix C. ∎

Next, the contribution of measurement noise $\mathbf{n}$ is modelled with a
zero-mean complex multivariate Gaussian random variable as
$\mathbf{n}\sim\mathcal{CN}(\mathbf{0},\sigma^{2}\mathbf{I})$ where
$\sigma^{2}$ denotes the variance of the noise. This can be incorporated in
the above lemma by defining
$\bm{\nu}_{\mathbf{q},\mathbf{n}}\overset{\bigtriangleup}{=}\bm{\nu}_{\mathbf{q}}+[\sigma^{2},0,\dots,0]$
and
$\mathbf{R}_{\mathbf{q},\mathbf{n}}=toep(\bm{\nu}_{\mathbf{q},\mathbf{n}},\bm{\nu}_{\mathbf{q},\mathbf{n}}^{H})$.
We then follow the Kronecker model to obtain the covariance matrix of the CFR,
which holds when the diffuse spectrum contribution in the angular domains is
independent from that in the frequency domain [9, 42]. Under the Kronecker
model, the covariance matrix of the CFR can be decomposed as
$\mathbf{R}=\mathbf{I}_{N_{Rx}}\otimes\mathbf{I}_{N_{Tx}}\otimes\mathbf{R}_{\mathbf{q},\mathbf{n}}$.
Therefore, the distribution of the CFR in (19) within the above model can be
given as $\mathbf{h}_{u}\sim\mathcal{CN}(\overline{\mathbf{h}},\mathbf{R})$.
We denote the parameters associated with the covariance matrix by
$\bm{\theta}_{vn}=[\alpha,\beta,L,\sigma]$, which relates to the variable part
of the CFR and noise. As mentioned earlier, the mean $\overline{\mathbf{h}}$
solely depends on the specular paths parameters $\bm{\theta}_{sp}$.

###  IV-B Channel spoofing detection problem

Channel-based spoofing detection is generally studied [5, 7] in the “snapshot”
scenario where Bob receives a new message claiming to be sent by Alice, and he
has to check whether the claim is true. To this end, we assume that Bob is
able to measure and store a noisy version of the CFR corresponding to a
transmitting terminal. Based on the CFRs associated with the incoming
messages, the goal in this scenario is to determine whether a message at time
$t=(u+1)T$ belongs to Alice or Eve given a reference message from Alice444In
the remaining of this section, we use $A$ or $E$ in the superscript of vector
or scalar to indicate that it belongs to Alice or Eve, respectively.,
$\mathbf{h}^{A}_{k}$, at time $t=uT$. In this setup, we use the terms message
and CFR interchangeably. One can pose the spoofing detection as a binary
classification problem for which two hypotheses can be made:

$\displaystyle\mathcal{H}_{0}:\mathbf{h}_{u+1}=\mathbf{h}^{A}_{u+1},$ (30)
$\displaystyle\mathcal{H}_{1}:\mathbf{h}_{u+1}=\mathbf{h}^{E}_{u+1}.$ (31)

Under the null hypothesis, $\mathcal{H}_{0}$, the message at time $t=(u+1)T$
belongs to Alice, while under the alternative hypothesis $\mathcal{H}_{1}$ a
spoofing attack has occurred, i.e., the message belongs to Eve.

For the data acquisition phase, we consider a setting where Bob spots a finite
number of snapshots from a coherence time and stores the observed CFRs.
Furthermore, in order to label the incoming CFRs, we use a heuristic method
given by

$\displaystyle\lVert\mathbf{h}_{u+1}-\mathbf{h}^{A}_{u}\lVert^{2}\underset{\mathcal{H}_{0}}{\overset{\mathcal{H}_{1}}{\gtreqless}}\eta,$
(32)

in the lieu of the channel parameters $\bm{\theta}_{vn}$ and
$\bm{\theta}_{sp}$. This method can be viewed as an imperfect labeling
mechanism which decides in favor of $\mathcal{H}_{0}$ if the Euclidean
distance between an incoming CFR and the reference CFR is smaller than a
predefined threshold $\eta$.

From a statistical perspective, likelihood ratio test is the main approach for
deciding between the two hypothesis, which relies on the knowledge of unknown
channel parameters as obtained in the following. The likelihood ratio test for
the snapshot scenario at time $t=(u+1)T$ is defined as

$\displaystyle\mathds{L}\big{(}\mathbf{h}_{u+1}|\mathbf{h}^{A}_{u}\big{)}\overset{\bigtriangleup}{=}\frac{p(\mathbf{h}_{u+1}-\mathbf{h}^{A}_{u}|\mathcal{H}_{0})}{p(\mathbf{h}_{u+1}-\mathbf{h}^{A}_{u}|\mathcal{H}_{1})}\underset{\mathcal{H}_{0}}{\overset{\mathcal{H}_{1}}{\gtreqless}}\zeta,$
(33)

for a predefined threshold $\zeta$, where the conditional probability
distribution of $\mathbf{h}_{k+1}-\mathbf{h}^{A}_{k}$ serves as the likelihood
function under each behavior. In the following, we obtain closed-form
expressions for these likelihood functions assuming the statistical dependence
on the reference CFR for each hypothesis is specified via the conditional
distributions $p(\mathbf{q}_{u+1}^{A}|\mathbf{q}_{u}^{A})$ and
$p(\mathbf{q}_{u+1}^{E}|\mathbf{q}_{u}^{A})$ under $\mathcal{H}_{0}$ and
$\mathcal{H}_{1}$, respectively. Specifically, the dependence of
$\mathbf{q}_{u+1}^{A}$ on $\mathbf{q}_{u}^{A}$ is characterized through the
corresponding channel gains in terms of an order-$1$ auto-regressive model
(AR-$1$) [5], i.e.,

$\displaystyle
A^{A}_{u+1,l}=a^{A}A^{A}_{u,l}+\sqrt{(1-({a^{A}})^{2})\text{Var}(A^{A}_{u+1,l})}w_{u+1,l}$
(34)

where $a^{A}$ denotes the similarity parameter, and
$w_{u+1,l}\sim\mathcal{CN}(0,1)$ is independent of $A_{u,l}$. Similarly, the
$l$th path gain corresponding to $\mathbf{q}_{u+1}^{E}$ and
$\mathbf{q}_{u}^{A}$ are related with and AR-$1$ model with similarity
parameter $a^{E}$.

######  Lemma 5.

Under the null hypothesis,
$p(\mathbf{q}_{u+1}-\mathbf{q}^{A}_{u}|\mathcal{H}_{0})=\mathcal{CN}(\mathbf{0},\mathbf{R}_{q,\mathcal{H}_{0}})$
for

$\displaystyle\mathbf{R}_{q,\mathcal{H}_{0}}=\text{toep}(\bm{\nu}_{\mathcal{H}_{0}},\bm{\nu}_{\mathcal{H}_{0}}^{H}),$
(35)
$\displaystyle\bm{\nu}_{\mathcal{H}_{0}}\overset{\bigtriangleup}{=}\Big{[}2(1-a^{A})\kappa(\bm{\theta}_{\mathbf{q}}^{A},0),2(1-a^{A})\kappa\big{(}\bm{\theta}_{\mathbf{q}}^{A},\frac{1}{M}\big{)},\dots$
$\displaystyle\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
,2(1-a^{a})\kappa\big{(}\bm{\theta}_{\mathbf{q}}^{A},\frac{M-1}{M}\big{)}\Big{]},$
(36)

where
$\bm{\theta}_{\mathbf{q}}^{A}\overset{\bigtriangleup}{=}[\alpha^{A},\beta^{A},L^{A}]$
and the $\kappa$ function defined in Lemma 4.

###### Proof.

See Appendix D. ∎

######  Lemma 6.

Under the alternative hypothesis,
$p(\mathbf{q}_{u+1}-\mathbf{q}^{A}_{u}|\mathcal{H}_{1})=\mathcal{CN}(\mathbf{0},\mathbf{R}_{\mathbf{q},\mathcal{H}_{1}})$
for

$\displaystyle\mathbf{R}_{\mathbf{q},\mathcal{H}_{1}}=\text{toep}(\bm{\nu}_{\mathcal{H}_{1}},\bm{\nu}_{\mathcal{H}_{1}}^{H}),$
(37)
$\displaystyle\bm{\nu}_{\mathcal{H}_{1}}\overset{\bigtriangleup}{=}[\kappa^{{}^{\prime}}(a^{E},\theta_{\mathbf{q}}^{A},\theta_{\mathbf{q}}^{E},0),\kappa^{{}^{\prime}}(a^{E},\theta_{\mathbf{q}}^{A},\theta_{\mathbf{q}}^{E},\frac{1}{M}),\dots,$
$\displaystyle\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\ \ \
\kappa^{{}^{\prime}}(a^{E},\theta_{\mathbf{q}}^{A},\theta_{\mathbf{q}}^{E},\frac{M-1}{M})],$
(38)
$\displaystyle\kappa^{{}^{\prime}}(a^{E},\theta_{\mathbf{q}}^{A},\theta_{\mathbf{q}}^{E},m)\overset{\bigtriangleup}{=}\kappa(\theta_{\mathbf{q}}^{E},m)-2a^{E}\kappa(\theta_{\mathbf{q}}^{A},m)+\kappa(\theta_{\mathbf{q}}^{A},m),$
(39)

where
$\theta_{\mathbf{q}}^{A}\overset{\bigtriangleup}{=}[\alpha^{A},\beta^{A},L^{A}]$,
$\theta_{\mathbf{q}}^{E}\overset{\bigtriangleup}{=}[\alpha^{E},\beta^{E},L^{E}]$,
and the $\kappa$ function is defined in Lemma 4.

###### Proof.

See Appendix E. ∎

The above two lemmas enables us to obtain the likelihood functions for (33).
Regarding the null hypothesis $\mathcal{H}_{0}$, using the Kronecker model for
the covariance matrix [9] and the fact that the measurement noise is
independent from the other CFR’s components, we obtain the covariance matrix
of $\mathbf{h}_{u+1}-\mathbf{h}^{A}_{u}$ as
$\mathbf{R}_{\mathcal{H}_{0}}=\mathbf{I}_{N_{Rx}}\otimes\mathbf{I}_{N_{Tx}}\otimes\mathbf{R}_{\mathbf{q},\mathcal{H}_{0}}+2(\sigma^{A})^{2}\mathbf{I}_{M}$
where $\mathbf{R}_{\mathbf{q},\mathcal{H}_{0}}$ is given in Lemma 5.
Furthermore, as a snapshot is captured in one coherence time, the specular
paths contribution of the CFRs remains the same ($\overline{\mathbf{h}}^{A}$)
between two consecutive times in this case. Therefore, under $\mathcal{H}_{0}$
the likelihood function is
$\mathcal{CN}(\mathbf{0},\mathbf{R}_{\mathcal{H}_{0}})$. Similarly, for the
alternate hypothesis $\mathcal{H}_{1}$, the likelihood function can be
obtained as
$\mathcal{CN}(\overline{\mathbf{h}}^{E}-\overline{\mathbf{h}}^{A},\mathbf{R}_{\mathcal{H}_{1}})$
where
$\mathbf{R}_{\mathcal{H}_{1}}=\mathbf{I}_{N_{Rx}}\otimes\mathbf{I}_{N_{Tx}}\otimes\mathbf{R}_{\mathbf{q},\mathcal{H}_{1}}+(\sigma^{A})^{2}\mathbf{I}_{M}+(\sigma^{E})^{2}\mathbf{I}_{M}$
and $\mathbf{R}_{\mathbf{q},\mathcal{H}_{1}}$ is given in Lemma 6.

###  IV-C Parameter estimation

In order to utilize the likelihood ratio test in (33), Bob requires the
knowledge of the channel parameters $\bm{\theta}_{sp}$, $\bm{\theta}_{vn}$
corresponding to Alice-Bob and Eve-Bob channels along with the similarity
parameters. In practice, these parameters should be estimated based on the
data collected from the observed snapshots. Recall from Lemma 4 that a CFR
associated with a terminal has a Gaussian distribution of the from
$\mathcal{CN}\big{(}\overline{\mathbf{h}}_{\bm{\theta}_{sp}},\mathbf{R}_{\bm{\theta}_{vn}}\big{)}$.
The MLE estimates of the parameters are obtained via

$\displaystyle\hat{\bm{\theta}}_{sp},\hat{\bm{\theta}}_{vn}=\operatorname*{\arg\\!\max}_{\bm{\theta}_{sp},\bm{\theta}_{vn}}\mathcal{L}\big{(}\mathbf{h}|\bm{\theta}_{sp},\mathbf{R}_{\bm{\theta}_{vn}}\big{)},$
(40a)
$\displaystyle\mathcal{L}(\mathbf{h}|\bm{\theta}_{sp},\mathbf{R}_{\bm{\theta}_{vn}}\big{)}=$
(40b)
$\displaystyle-M\ln\pi-\ln{\det{\mathbf{R}_{\bm{\theta}_{vn}}}}-\big{(}\mathbf{h}-\overline{\mathbf{h}}_{\bm{\theta}_{sp}}\big{)}^{H}\mathbf{R}^{-1}_{\bm{\theta}_{vn}}\big{(}\mathbf{h}-\overline{\mathbf{h}}_{\bm{\theta}_{sp}}\big{)},$
(40c)

which amounts to jointly maximizing the arguments of some nonlinear objective
function. Besides, it can be proved that (40c) is not a convex function of
$\bm{\theta}_{sp}$, and as a result there is no unique solution set for the
optimization problem in (40a). In practice, solving such problem is far from
trivial, especially since the number of nonlinear parameters
($\bm{\theta}_{sp}$) is large and a multidimensional exhaustive search is not
feasible. As a workaround, the authors in [9, 42] propose a suboptimal
procedure to break the problem into two sub-problems and estimate
$\bm{\theta}_{sp}$ and $\bm{\theta}_{vn}$ in a separate manner. In our
problem, there is also similarity parameters $a^{A}$ and $a^{E}$ which similar
to $\bm{\theta}_{vn}$ appear in the covariance matrix of a Gaussian
distribution, as obtained in Lemmas 5 and 6. Therefore, based on the approach
taken in [9, 42], we break the parameter estimation problem into three sub-
problems. Each sub-problem involves numerically maximizing the objective
function of the form (40c) w.r.t. $\bm{\theta}_{sp}$ or $\bm{\theta}_{vn}$ or
similarity parameters via an iterative local optimization technique, such as
Gauss-Newton algorithm. In particular, the maximization processes are done
sequentially and in an alternating manner between the three sets of parameters
towards convergence. In the following, we elaborate on each sub-problem for
the specific channel model we described earlier.

####  IV-C1 Estimating the specular path parameters $\bm{\theta}_{sp}$

The main goal here is to obtain an estimate of $\bm{\theta}_{sp}$ which
maximizes (40c) for a given estimate of $\bm{\theta}_{vn}$. In the following,
we use the $N$-exponential basis function defined as

$\displaystyle\mathbf{U}_{N}^{\mathbf{v}}=\begin{bmatrix}e^{-j\big{(}-\frac{N-1}{2}\big{)}\mathbf{v}[1]}&\dots&e^{-j\big{(}-\frac{N-1}{2}\big{)}\mathbf{v}[n]}\\\
\vdots&\ddots&\vdots\\\
e^{-j\big{(}\frac{N-1}{2}\big{)}\mathbf{v}[1]}&\dots&e^{-j\big{(}\frac{N-1}{2}\big{)}\mathbf{v}[n]}\end{bmatrix},$
(41)

for a vector $\mathbf{v}$ of length $N$. Partial derivative of
$\mathbf{U}_{N}^{\mathbf{v}}$ w.r.t $\mathbf{v}$ is readily computed as
$\mathbf{D}_{N}^{\mathbf{v}}=\frac{\partial\mathbf{U}_{N}^{\mathbf{v}}}{\partial\mathbf{v}}=-j\Xi_{N}\mathbf{U}_{N}^{\mathbf{v}}$
where $\Xi_{N}=\text{diag}([-(N-1)/2,\dots,(N-1)/2])$. Furthermore, we recall
that for arbitrary matrices $\mathbf{A}\in\mathbb{C}^{N\times P},\
\mathbf{B}\in\mathbb{C}^{M\times P},\ \mathbf{Q}^{P\times
P}=\text{diag}(\mathbf{q})$ and a vector $\mathbf{q}\in\mathbb{C}^{P\times
1}$, one can write
$\text{vec}\\{\mathbf{B}\mathbf{Q}\mathbf{A}^{T}\\}=(\mathbf{A}\odot\mathbf{B})\mathbf{q}$.
Utilizing this result along with the exponential basis function we can rewrite
the specular path contribution introduced in (26) for the CFR model as

$\displaystyle\overline{\mathbf{h}}=\big{(}\mathbf{U}_{N_{Rx}}^{\bm{\psi}_{T}}\odot\mathbf{U}_{N_{Tx}}^{\bm{\psi}_{R}}\odot\mathbf{U}_{N_{f}}^{\bm{\tau}}\big{)}\bm{\rho},$
(42)

which greatly simplifies the calculation of the first and second derivatives
of $\overline{\mathbf{h}}$ w.r.t $\bm{\theta}_{sp}$. Specifically, the
Jacobian matrix for the above model is obtained via
$\mathbf{J}(\bm{\theta}_{sp})=\mathbf{J}_{{\bm{\psi}_{T}}}\odot\mathbf{J}_{\bm{\psi}_{R}}\odot\mathbf{J}_{\mathbf{\rho}}\odot\mathbf{J}_{\bm{\tau}}$
where the Jacobian matrix’s components are given by

$\displaystyle\mathbf{J}_{{\bm{\psi}_{T}}}$
$\displaystyle=\begin{bmatrix}\mathbf{D}_{N_{Tx}}^{\bm{\psi}_{T}}&\mathbf{U}_{N_{Tx}}^{\bm{\psi}_{T}}&\mathbf{U}_{N_{Tx}}^{\bm{\psi}_{T}}&\mathbf{U}_{N_{Tx}}^{\bm{\psi}_{T}}&\mathbf{U}_{N_{Tx}}^{\bm{\psi}_{T}}\end{bmatrix},$
(43a) $\displaystyle\mathbf{J}_{\bm{\psi}_{R}}$
$\displaystyle=\begin{bmatrix}\mathbf{U}_{N_{Rx}}^{\bm{\psi}_{R}}&\mathbf{D}_{N_{Rx}}^{\bm{\psi}_{R}}&\mathbf{U}_{N_{Rx}}^{\bm{\psi}_{R}}&\mathbf{U}_{N_{Rx}}^{\bm{\psi}_{R}}&\mathbf{U}_{N_{Rx}}^{\bm{\psi}_{R}}\end{bmatrix},$
(43b) $\displaystyle\mathbf{J}_{\bm{\tau}}$
$\displaystyle=\begin{bmatrix}\mathbf{U}_{N_{f}}^{\bm{\tau}}&\mathbf{U}_{N_{f}}^{\bm{\tau}}&\mathbf{D}_{N_{f}}^{\bm{\tau}}&\mathbf{U}_{N_{f}}^{\bm{\tau}}&\mathbf{U}_{N_{f}}^{\bm{\tau}}\end{bmatrix},$
(43c) $\displaystyle\mathbf{J}_{\bm{\rho}}$
$\displaystyle=\begin{bmatrix}{\bm{\rho}}^{T}&{\bm{\rho}}^{T}&{\bm{\rho}}^{T}&\mathbf{1}^{T}&\mathbf{1}^{T}j\end{bmatrix}.$
(43d)

The first order partial derivative of the log likelihood function (40c) given
an observation $\mathbf{h}$ with respect to the parameters $\bm{\theta}_{sp}$
is denoted by
$\mathbf{q}_{\bm{\theta}_{sp}}\big{(}\mathbf{h}|\mathbf{R}_{\bm{\theta}_{vn}}\big{)}$.
For a given Jacobian matrix $\mathbf{J}(\bm{\theta}_{s}p)$, one can compute
[42]

$\displaystyle\mathbf{q}_{\bm{\theta}_{sp}}\big{(}\mathbf{h}|\mathbf{R}_{\bm{\theta}_{vn}}\big{)}=2\Re{\big{\\{}\mathbf{J}^{H}(\bm{\theta}_{sp})\mathbf{R}^{-1}_{\bm{\theta}_{vn}}(\mathbf{h}-\overline{\mathbf{h}}_{\bm{\theta}_{sp}}\big{\\}}}.$
(44)

Furthermore, the negative covariance matrix of the above first order
derivative, i.e.,

$\displaystyle-\mathbb{E}\bigg{[}\frac{\partial\mathcal{L}(\mathbf{h}|\bm{\theta}_{sp},\mathbf{R}_{\bm{\theta}_{vn}}\big{)}}{\partial\bm{\theta}_{sp}}\Big{(}\frac{\partial\mathcal{L}(\mathbf{h}|\bm{\theta}_{sp},\mathbf{R}_{\bm{\theta}_{vn}}\big{)}}{\partial\bm{\theta}_{sp}}\Big{)}^{T}\bigg{]},$
(45)

is called the Fisher information matrix (FIM),
$\mathbf{F}\big{(}{\bm{\theta}_{sp}}|\mathbf{R}_{\bm{\theta}_{vn}}\big{)}$,
which can be expressed in terms of the Jacobian matrix
$\mathbf{J}(\bm{\theta}_{sp})$ as [42]

$\displaystyle\mathbf{F}\big{(}{\bm{\theta}_{sp}}|\mathbf{R}_{\bm{\theta}_{vn}}\big{)}=2\Re\\{{\mathbf{J}^{H}(\bm{\theta}_{sp})\mathbf{R}^{-1}_{\bm{\theta}_{vn}}\mathbf{J}(\bm{\theta}_{sp})}\\}.$
(46)

Having obtained (44) and (46), a local optimization technique can be utilized
to obtain an iterative rule for estimation of $\bm{\theta}_{sp}$. To this end,
we employ the Gauss-Newton algorithm as

$\displaystyle\hat{\bm{\theta}}^{i+1}_{sp}=\hat{\bm{\theta}}^{i}_{sp}+\zeta\
\mathbf{F}^{-1}\big{(}\hat{\bm{\theta}}^{i}_{sp}|\mathbf{R}_{\bm{\theta}_{vn}}\big{)}\mathbf{q}_{\hat{\bm{\theta}}^{i}_{sp}}\big{(}\mathbf{h}|\mathbf{R}_{\bm{\theta}_{vn}}\big{)}$
(47)

for a step length $\zeta$ which should be chosen such that
$\mathcal{L}(\mathbf{h}|\bm{\theta}^{i+1}_{sp},\mathbf{R}_{\bm{\theta}_{vn}}\big{)}>\mathcal{L}(\mathbf{h}|\bm{\theta}^{i}_{sp},\mathbf{R}_{\bm{\theta}_{vn}}\big{)}$.

####  IV-C2 Estimating the variable part and noise parameters
$\bm{\theta}_{vn}$

For a given estimate of $\bm{\theta}_{sp}$, the goal here is to estimate
${\bm{\theta}_{vn}}$. To this end, we assume $N$ number of CFRs denoted with
$\\{\mathbf{h}_{i}\\}_{i=1}^{N}$ are available. For each CFR, the approach
presented in IV-C1 is utilized to estimate the corresponding specular paths
parameters as $\\{\hat{\bm{\theta}}_{sp,i}\\}_{i=1}^{N}$. In order to remove
the contribution of the specular paths from the CFRs, we from an $M\times N$
matrix $\mathbf{H}$ whose $i$th column amounts to
$\mathbf{h}_{i}-\overline{\mathbf{h}}(\overline{\bm{\theta}}_{sp})$ where
$\overline{\bm{\theta}}_{sp}$ denotes the average value of
$\hat{\bm{\theta}}_{sp,i}$’s. In the following, we propose an estimation
procedure for $\bm{\theta}_{vn}$ using $\mathbf{H}$. We first note that all
the parameters in $\bm{\theta}_{vn}$ are continuous except for the number of
diffuse spectrum paths $L$ which is an integer value. As a result, the
objective function is not continuous in $L$ and the gradient of (40c) does not
exist w.r.t. $L$. In this vein, we take a sub-optimal approach towards
estimating $L$ by separating it from the rest of the parameters in
$\bm{\theta}_{vn}$. The authors in [43] have proposed an eigenvalue ratio
method for estimating the number of harmonics present in the signals from $N$
available observations. For our specific case, we denote the MLE of the
covariance of $\mathbf{H}$ by $\mathbf{C}_{\mathbf{H}}$, and the corresponding
eigenvalues by $\mathbf{e}_{i},\ i=1,\dots,M$. The propose heuristic approach
is to choose $\hat{L}$ such that
$\frac{\sum_{i=1}^{\hat{L}}\mathbf{e}_{i}}{\sum_{i=1}^{M}\mathbf{e}_{i}}\geq\eta,$
for a predefined value of$\eta$.

$\displaystyle\frac{\partial\bm{\nu}_{\mathbf{q},\mathbf{n}}}{\partial\sigma^{2}}=\begin{bmatrix}1,0,\dots,0\end{bmatrix},$
(48a)
$\displaystyle\frac{\partial\bm{\nu}_{\mathbf{q},\mathbf{n}}}{\partial\alpha^{2}}=\begin{bmatrix}1-e^{-2\pi
L\beta},\frac{(1-e^{-2\pi\beta})\big{(}1-f^{L}(\frac{1}{M})\big{)}}{1-f(\frac{1}{M})},\dots,\frac{(1-e^{-2\pi\beta})\big{(}1-f^{L}(1-\frac{1}{M})\big{)}}{1-f(1-\frac{1}{M})}\end{bmatrix},$
(48b)
$\displaystyle\frac{\partial\bm{\nu}_{\mathbf{q},\mathbf{n}}}{\partial\beta}=\Big{[}\begin{matrix}2\pi\alpha^{2}Le^{-2\pi\beta
L},\frac{2\pi
e^{-2\pi\beta}\big{(}f^{L}(\frac{1}{M})-1\big{)}}{f(\frac{1}{M})-1}+\frac{2L\pi
f^{L}(\frac{1}{M})(e^{-2\pi\beta}-1)}{f(\frac{1}{M})-1}-\frac{2\pi
f(\frac{1}{M})\big{(}f^{L}(\frac{1}{M})-1\big{)}(e^{-2\pi\beta}-1)}{\big{(}f(\frac{1}{M})-1\big{)}^{2}},\end{matrix}$
$\displaystyle\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \begin{matrix}\dots,\frac{2\pi
e^{-2\pi\beta}\big{(}f^{L}(1-\frac{1}{M})-1\big{)}}{f(1-\frac{1}{M})-1}+\frac{2L\pi
f^{L}(1-\frac{1}{M})(e^{-2\pi\beta}-1)}{f(1-\frac{1}{M})-1}-\frac{2\pi
f(1-\frac{1}{M})\big{(}f^{L}(1-\frac{1}{M})-1\big{)}(e^{-2\pi\beta}-1)}{\big{(}f(1-\frac{1}{M})-1\big{)}^{2}}\end{matrix}\Big{]}.$
(48c)

We plug the estimated value of $L$ in the parameter vector to obtain
$\bm{\theta}_{vn}=[\sigma,\alpha,\beta,\hat{L}]$. Then, the log-likelihood
function for the zero-mean CFRs can be written as

$\displaystyle\mathcal{L}(\mathbf{H}|{\bm{\theta}_{vn}})=-MN\ln\pi-N\ln{\det{\mathbf{R}_{\bm{\theta}_{vn}}}}-\operatorname{Tr}\big{(}\mathbf{H}^{H}\mathbf{R}^{-1}_{\bm{\theta}_{vn}}\mathbf{H}\big{)}.$
(49)

The first-order gradient of $\mathcal{L}(\mathbf{H}|{\bm{\theta}_{vn}})$
w.r.t. to each parameter can be computed as [42]

$\displaystyle\frac{\partial\mathcal{L}(\mathbf{H}|{\bm{\theta}_{vn}})}{\partial{\bm{\theta}_{vn}}[i]}=N\operatorname{Tr}\Big{(}\mathbf{R}^{-1}_{\bm{\theta}_{vn}}\frac{\partial\mathbf{R}_{\bm{\theta}_{vn}}}{\partial{\bm{\theta}_{vn}}[i]}\mathbf{R}^{-1}_{\bm{\theta}_{vn}}(\widehat{\mathbf{R}}-\mathbf{R}_{\bm{\theta}_{vn}})\Big{)}.$
(50)

for $i=1,2,3$. Subsequently, the $(i,j)$th element of the FIM corresponding to
$\mathcal{L}(\mathbf{H}|{\bm{\theta}_{vn}})$ equals [42]

$\displaystyle-\mathbb{E}\Big{[}\frac{\partial^{2}\mathcal{L}(\mathbf{H}|{\bm{\theta}_{vn}})}{\partial{\bm{\theta}_{vn}}[i]\partial{\bm{\theta}_{vn}}(j)}\Big{]}=N\operatorname{Tr}\Big{(}\mathbf{R}^{-1}_{\bm{\theta}_{vn}}\frac{\partial\mathbf{R}_{\bm{\theta}_{vn}}}{\partial{\bm{\theta}_{vn}}[i]}\mathbf{R}^{-1}_{\bm{\theta}_{vn}}\frac{\partial\mathbf{R}_{\bm{\theta}_{vn}}}{\partial{\bm{\theta}_{vn}}(j)}\Big{)}.$
(51)

To obtain explicit expressions for (50) and (51), one needs to compute partial
derivatives terms, i.e.,
$\frac{\partial\mathbf{R}_{\bm{\theta}_{vn}}}{\partial{\bm{\theta}_{vn}}[i]}$.
Considering the Toeplitz structure of the covariance model described in Lemma
4, we can write

$\displaystyle\frac{\partial\mathbf{R}_{\mathbf{q},\mathbf{n}}(\bm{\theta}_{vn})}{\partial{\bm{\theta}_{vn}}[i]}=toep\Big{(}\frac{\partial\bm{\nu}_{\mathbf{q},\mathbf{n}}}{\partial{\bm{\theta}_{vn}}[i]},\frac{\partial\bm{\nu}^{H}_{\mathbf{q},\mathbf{n}}}{\partial{\bm{\theta}_{vn}}[i]}\Big{)},$
(52)
$\displaystyle\frac{\partial\mathbf{R}_{\bm{\theta}_{vn}}}{\partial{\bm{\theta}_{vn}}[i]}=\mathbf{I}_{N_{Rx}}\otimes\mathbf{I}_{N_{Tx}}\otimes\frac{\partial\mathbf{R}_{\mathbf{q},\mathbf{n}}(\bm{\theta}_{vn})}{\partial{\bm{\theta}_{vn}}[i]},$
(53)

where the partial derivative for each parameter is obtained in (48a)-(48c) for
$f(m)=e^{-2\pi(\beta-jm)}$. Plugging this in (50) and (51) leads to
computation of first-order gradient and the FIM of the likelihood function.
Then, an iterative approach like the Gauss-Newton algorithm can be employed
for estimating ${\bm{\theta}_{vn}}$ in a similar fashion to the case of
$\bm{\theta}_{sp}$ in (47).

####  IV-C3 Estimating the similarity parameters $a^{A}$ and $a^{E}$

In order to estimate the parameters $a^{A}$ and $a^{E}$, one can use the
likelihood functions obtained in Lemmas 5 and 6, respectively, while using the
estimates obtained in Sections IV-C1 and IV-C2 for the values of all the other
parameters. Specifically, as these parameters appear in the covariance matrix
of a Gaussian distribution, a similar estimation procedure to that of
$\bm{\theta}_{sp}$ can be employed here as well. In fact, the expressions for
the first-order gradient and the FIM of the likelihood function in this case
is the same as those in (50) and (51), respectively, except for the fact that
there is only one parameter to estimate in this case. For example, for
estimation of $a^{A}$, by considering the Toeplitz structure of the covariance
model described in Lemma 5, we can write

$\displaystyle\frac{\partial\mathbf{R}_{\mathbf{q},\mathcal{H}_{0}}(a^{A})}{\partial
a^{A}}=toep\Big{(}\frac{\partial\bm{\nu}_{\mathcal{H}_{0}}}{\partial
a^{A}},\frac{\partial\bm{\nu}^{H}_{\mathcal{H}_{0}}}{\partial a^{A}}\Big{)},$
(54) $\displaystyle\frac{\partial\mathbf{R}_{\mathcal{H}_{0}}(a^{A})}{\partial
a^{A}}=\mathbf{I}_{N_{Rx}}\otimes\mathbf{I}_{N_{Tx}}\otimes\frac{\partial\mathbf{R}_{\mathbf{q},\mathcal{H}_{0}}(a^{A})}{\partial
a^{A}},$ (55)

where the partial derivative w.r.t. $a^{A}$ can be obtained as

$\displaystyle\frac{\partial\bm{\nu}_{\mathcal{H}_{0}}}{\partial
a^{A}}=-2\Big{[}$
$\displaystyle\frac{(\alpha^{A})^{2}(1-e^{-2\pi\beta^{A}})(1-f^{L^{A}}(0))}{1-f(0)},$
(56)
$\displaystyle\frac{(\alpha^{A})^{2}(1-e^{-2\pi\beta^{A}})(1-f^{L^{A}}(\frac{1}{M}))}{1-f(\frac{1}{M})},\dots,$
$\displaystyle\frac{(\alpha^{A})^{2}(1-e^{-2\pi\beta^{A}})(1-f^{L^{A}}(1-\frac{1}{M}))}{1-f(1-\frac{1}{M})}\Big{]},$

for $f(m)=e^{-2\pi(\beta-jm)}$. Subsequently, using the first-order gradient
and the FIM of the likelihood function, we utilize the Gauss-Newton algorithm
to estimate $a^{A}$. Similar approach can be taken for estimating $a^{E}$
using the covariance model described in Lemma 6 which we omit here for
brevity.

###  IV-D HyPhyLearn for channel spoofing detection

We propose to utilize HyPhyLearn described in Algorithm 1 for solving the
spoofing detection problem which can be seen as a binary instance of the
problem formulation described in Section II with two behaviors, as described
in (33). Besides, statistical parametric models are available for each
behavior, the high complexity of which makes one to resort to suboptimal
parameter estimation procedure. As mentioned in Section IV-B the data
corresponding to the Alice and Eve are collected in the snapshot setting, and
subsequently (imperfectly) labeled according to (32). Then, using these
collected CFRs, the underlying parameters of each likelihood function in (33)
is estimated using the sub-optimal parameter estimation procedure described in
Section IV-C. Next, the estimated parameters are plugged in the available
parametric models $\mathcal{CN}(\mathbf{0},\mathbf{R}_{\mathcal{H}_{0}})$ and
$\mathcal{CN}(\overline{\mathbf{h}}^{E}-\overline{\mathbf{h}}^{A},\mathbf{R}_{\mathcal{H}_{1}})$,
which subsequently are used to generate synthetic CFRs. Finally, the collected
and synthetic CFRs are incorporated in the Step $4$ of Algorithm 1 to train
the learning-based classifier which can be utilized as a spoofing detector.

##  V Case study 2: Multi-user detection

An important problem in multipoint-to-point digital communication networks
(e.g., radio networks, local-area networks, and uplink satellite channels) is
the optimum centralized demodulation of the information sent simultaneously by
several users through a Gaussian multiple-access channel. Even though the
users may not employ a protocol to coordinate their transmission epochs,
effective sharing of the channel is possible because each user modulates a
different signature signal waveform that is known by the intended receiver
(Code Division Multiple Access (CDMA)). In this section, we consider the
uplink of a cellular communication system where $K$ users are asynchronously
sharing a channel to communicate with a base station (BS). The problem of
multi-user detection in this setting amounts to inferring the information
associated with each user from this multiple access channel.

###  V-A Multi-user detection problem

Consider the uplink of an asynchronous direct-sequence (DS) CDMA system shared
by $K$ users, employing long spreading codes, bandlimited chip pulses and
operating over a frequency-selective fading channel. The baseband equivalent
of the received signal may be written as

$\displaystyle
r(t)=\sum_{p=0}^{P-1}\sum_{k=0}^{K-1}A_{k}b_{k}(p)s^{\prime}_{k,p}(t-\tau_{k}-pT_{b})^{*}c_{k}(t)+w(t),$
(57)

where $P$ is the number of transmitted packets and $s^{\prime}_{k,p}(t)$
denotes the $k$th user signature waveform. Furthermore, $T_{b}$ is the bit-
interval duration, $A_{k}$ and $\tau_{k}$ denote the respective complex
amplitude and timing offset of $k$th user, and $b_{k}(p)$ is the $k$th user’s
information bit in the $p$th signaling interval, whereas $w(t)$ is the complex
envelope of the additive noise term, which is assumed to be a zero-mean, wide-
sense stationary complex white Gaussian process. Moreover, $c_{k}(t)$ is the
impulse response modeling the channel effects between the BS and the $k$th
user. Note that the channel impulse responses $c_{k}(t)$ are assumed to be
time-invariant over each transmitted frame [44] under the assumption that the
channel coherence time exceeds the packet duration $BT_{b}$. Regarding the
$k$th user signature waveform, we have

$\displaystyle
s^{\prime}_{k,p}(t-\tau_{k}-pT_{b})=\sum_{n=0}^{N-1}\beta_{k,p}^{(n)}h_{SRRC}(t-nT_{c}),$
(58)

where $\\{\beta^{(n)}_{k,p}\\}_{n=0}^{N-1}$ is the pseudo-noise (PN) code
employed by user $k$ for spreading its data bit on the $p$th symbol interval,
$N$ is the processing gain, and $T_{c}=T_{b}/N$ is the chip interval.
Furthermore, $h_{SRRC}(t)$ denotes the square root raised-cosine waveform as
the bandlimited chip pulse which following [44] is time-limited to $[0,4Tc]$.

In the BS, chip-matched filtering and chip-rate sampling is done in order to
convert the received signal to discrete time domain. To this end, $r(t)$ is
convolved with chip-matched filter $h_{SRRC}(4T_{c}-t)$ followed by a sampler
at a rate $2/T_{c}$ (Nyquist rate). The convolution operation results in

$\displaystyle y(t)$ $\displaystyle=r(t)*h_{SRRC}(4T_{c}-t)$
$\displaystyle=\sum_{p=0}^{P-1}\sum_{k=0}^{K-1}b_{k}(p)h_{k,p}(t-pT_{b},\tau_{k})+n(t),$
(59)

where $h_{k,p}(t,\tau_{k})=A_{k}s_{k,p}(t-\tau_{k})^{*}c_{k}(t)$ is called the
effective signature waveform for
$s_{k,p}(t)=\sum_{n=0}^{N-1}\beta_{k,p}^{(n)}h_{RC}(t-nT_{c})$, and
$h_{RC}(t)$ represents a raised cosine chip waveform time-limited to
$[0,8Tc)$. As $h_{k,p}(t-pT_{b},\tau_{k})$ has a time domain support of
$[pT_{b},(p+2)T_{b}+7T_{c}]$, during the $p$th symbol interval
$\mathcal{I}_{p}=[pT_{b},(p+1)T_{b}]$, the contribution from at most three
bits for each user, i.e., the $p$th, the $p-1$th and the $p-2$th ones, is
observed assuming that $\tau_{k}+T_{m}<T_{b}$ where $T_{m}$ stands for the
maximum delay spread. Therefore, sampling the waveform $y(t)$ at rate
$M/T_{c}$, the MN-dimensional vector $y(p)$ collecting the data samples of the
interval $\mathcal{I}_{p}$ can be expressed as

$\displaystyle\mathbf{y}(p)=\sum_{k=0}^{K-1}[b_{k}(p-2)\mathbf{h}_{k,p-2}(p)+b_{k}(p-1)\mathbf{h}_{k,p-1}(p)$
$\displaystyle\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\ \ \ \ +b_{k}(p)\mathbf{h}_{k,p}(p)]+\mathbf{n}(p)$ (60)

where $\mathbf{h}_{k,p-i}(p)$ and $\mathbf{n}(p)$ comprise the $MN$ samples of
$h_{k,p-i}(t-(p-i)T_{b},\tau_{k})$, $i\in\\{0,1,2\\}$ and $n(t)$,
respectively, during $\mathcal{I}_{p}$. We set $M=2$ in the following
discussion. A compact representation of $\mathbf{y}(p)$ can be obtained by
relying on the notion of effective chip pulse defined as
$g_{k}(t,\tau_{k})=A_{k}h_{RC}(t-\tau_{k})^{*}c_{k}(t)$ which is supported on
the interval $[0,T_{b}+8T_{c}]$. Noting that
$h_{k,p}(t,\tau_{k})=\sum_{i=0}^{N-1}\beta_{k,p}^{n}g_{k}(t-nT_{c},\tau_{k})$,
and defining $\mathbf{g}_{k}\in\mathbb{C}^{MN+8M-1\times 1}$ as

$\displaystyle\mathbf{g}_{k}=\big{[}g_{k}(T_{c}/M,\tau_{k}),g_{k}(2T_{c}/M,\tau_{k}),\dots,$
$\displaystyle\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
g_{k}(T_{b}+(8M-1)T_{c}/M,\tau_{k})\big{]}^{T},$ (61)

one can write $\mathbf{h}_{k,p-i}(p)=\mathbf{C}_{k,p-i}(p)\mathbf{g}_{k}$,
where $\mathbf{C}_{k,p-i}(p)$ is a $MN\times(MN+8M-1)$ dimensional matrix as a
function of $\beta_{k,p}^{n}$, and obtained in details in ($9$)-($11$) of
[44]. Then, we have

$\displaystyle\mathbf{y}(p)=\sum_{k=0}^{K-1}\mathbf{A}_{k}(p)\mathbf{g}_{k}+\mathbf{n}(p)=\mathbf{A}(p)\mathbf{g}+\mathbf{n}(p),$
(62)

for
$\mathbf{A}_{k}(p)=b_{k}(p-2)\mathbf{C}_{k,p-2}(p)+b_{k}(p-1)\mathbf{C}_{k,p-1}(p)+b_{k}(p)\mathbf{C}_{k,p}(p)$,
$\mathbf{A}(p)=[\mathbf{A}_{0}(p),\dots,\mathbf{A}_{K-1}(p)]$, and
$\mathbf{g}=[\mathbf{g}_{0}^{T},\dots,\mathbf{g}_{K-1}^{T}]^{T}$.

The multiuser detection problem can be cast as $2^{K}$-ary classification
problem where the goal is to find the vector of information bits
$\mathbf{b}=[b_{0}(p),\dots,b_{K-1}(p)]$ given a observation vector
$\mathbf{y}^{\prime}(p)$. Assuming all the vectors
$\mathbf{b}\in\\{0,1\\}^{K}$ are a priori equiprobable the minimum distance
rule gives the maximum a posteriori decision [45]. Mathematically, the
multiuser detection is equivalent to solving the minimization problem
$\operatorname*{\arg\\!\min}_{\mathbf{b}\in\\{0,1\\}^{K}}\
\mathbf{y}^{\prime}(p)-\sum_{k=0}^{K-1}\mathbf{A}_{k}(p)\mathbf{g}_{k}$.
However, the complexity of such detector is exponential in the number of users
[45] and in practice sub-optimal methods like minimum mean square error (MMSE)
detector [45] is utilized for performance evaluation. We note that the
multiuser detection methods relies on the channel parameters and the spreading
codes corresponding to each user, and we assume true knowledge of both are not
available at the BS. Specifically, we consider a case where a mismatch exists
between the true spreading codes [46], utilized by the users, and the
corresponding ones available at the BS. Besides, we assume that BS has access
to $N_{T}$ number of training data from the $k$th user. The channel parameters
are then estimated based on the available training data using the procedure
described in the following section.

###  V-B Parameter estimation

The performance of multi-user detection relies heavily on the estimation of
the channel parameters. We assume the channel impulse response (CIR),
$c_{k}(t)$, takes the form of a time-invariant multipath channel with $L$
paths, i.e.,
$c_{k}(t)=\sum_{l=0}^{L-1}\alpha_{k,l}\delta(t-\tau^{\prime}_{k,l})$, which is
parameterized by the complex path gains $\alpha_{k,l}$ and the corresponding
path delays $\tau^{\prime}_{k,l}$. The joint ML estimate of these parameters
requires an exhaustive search over the continuous $K$-dimensional space
$[0,T_{b})^{K}$ which is computationally prohibitive. It is shown in [11] that
even using the conventional grid search based scheme to find a near-ML
solution NP-hard. As a workaround, alternative sub-optimal estimation methods
of low-complexity are proposed for practical settings. Notably, the authors in
[44] propose a two-step approach which first estimates the the channel impulse
response (CIR) using the Least Squares (LS) criterion, and then extracts the
underlying channel parameters. In particular, given the knowledge of the
spreading codes and information bits for all the users in the training
dataset, the overall CIR $\mathbf{g}$ may be directly estimated by invoking
the LS estimation procedure

$\displaystyle\hat{\mathbf{g}}$
$\displaystyle=\operatorname*{\arg\\!\min}_{\mathbf{x}}\sum_{i=1}^{N_{T}-1}||\mathbf{y}(p)-\mathbf{A}(p)\mathbf{x}||^{2}$
(63)
$\displaystyle=\bigg{[}\sum_{p=0}^{N_{T}-1}\mathbf{A}^{H}(p)\mathbf{A}(p)\bigg{]}^{-1}\bigg{[}\sum_{p=0}^{N_{T}-1}\mathbf{A}^{H}(p)\mathbf{y}(p)\bigg{]}.$
(64)

Relying on $\hat{\mathbf{g}}$ the authors in [44] propose an ad-hoc algorithm
to estimate the channel parameters. Specifically, the explicit parameters to
be estimated include delays $\tau_{k,l}=\tau^{\prime}_{k,l}+\tau_{k}$,
amplitudes $a_{k,l}=A_{k}|\alpha_{k,l}|$ and the phases
$\phi_{k,l}=\arg(a_{k,l})$ for $k=0,\dots,K-1$ and $l=0,\dots,L-1$. We provide
an overview of the above ad-hoc parameter estimation procedure in Appendix F
for completeness.

###  V-C HyPhyLearn for multiuser detection

We utilize HyPhyLearn to solve the problem of multiuser detection described in
Section V-A as a $2^{K}$-ary classification problem. In particular, due to the
available statistical parametric models for each class on one hand, and lack
of an estimation procedure for the underlying channel parameters which is
optimal in some sense on the other hand, the multiuser detection can be framed
within the problem formulation setting described in Section II. We use the
available data corresponding to the users in the suboptimal estimation method
described in Section V-B to obtain the estimates of the channel parameters for
$K$ users
$\hat{{\bm{\tau}}}=[\hat{\tau}_{0,0},\dots,\hat{\tau}_{0,L-1},\dots,\hat{\tau}_{K-1,L-1}]$,
$\hat{\mathbf{a}}=[\hat{a}_{0,0},\dots,\hat{a}_{0,L-1},\dots,\hat{a}_{K-1,K-1}]$
and
$\hat{\mathbf{\Phi}}=[\hat{\phi}_{0,0},\dots,\hat{\phi}_{0,L-1},\dots,\hat{\phi}_{K-1,L-1}]$.
Using these estimates along with the imperfect knowledge of spreading codes
for the training data, one can utilize the parametric model (62) to generate a
synthetic data example corresponding to sample information bits $\mathbf{b}$.
This synthetic data sample is subsequently added to the synthetic dataset
along with its corresponding label $\mathbf{b}$. Then, the synthetic dataset
is incorporated with the available training dataset according to the Step $4$
of the Algorithm 1 to find the learning-based classifier. Specifically, this
classifier has $2^{K}$ output neurons, each corresponding to a specific
information bits vector, which enables it to to serve as a data detection
method for the $K$-user system.

##  VI Numerical Results

In this section, we numerically evaluate the performance of our proposed
solution, HyPhyLearn, described in Algorithm 1 for the two case studies
described in Sections IV and V. This involves comparing the resulting
performance against that of the existing statistical classifiers and other
hybrid classification methods, and highlighting the superiority of our
proposed solution for the problems under study.

###  VI-A Spoofing detection problem

In the Alice–Eve–Bob setting, we begin with a scenario where the coherence
time of the Alice–Bob and the Eve–Bob channel are very large, and therefore
the corresponding channel parameters are fixed between the training and
testing stages. As mentioned in Section IV-B, the training data in this
problem are collected by observing finite number of snapshots by Bob. The
training CFRs from each snapshot are subsequently labeled using the heuristic
test (32). The number of received antennas and transmit antennas at Alice and
Bob is set to $2$. Also, following the discussion in [7] we assume Eve also
uses the same number of antennas to impersonate Alice. The number of
subcarriers is set to $N_{f}=20$, which makes the total number of samples
associated with each CFR equal $M=80$. We assume the Alice–Bob parameters are
$\sigma^{2}_{A}=20$, $\alpha^{2}_{A}=200$, $\beta_{A}=0.02$ and $a^{A}=0.85$,
while $\sigma^{2}_{E}=26$, $\alpha^{2}_{E}=250$, $\beta_{E}=0.08$ and
$a^{E}=0.65$ are used for the Eve–Bob channel. Furthermore, we set $L_{A}=20$
and $L_{E}=16$ as the number of diffuse spectrum virtual paths, while the
number of specular paths are set to $4$ for both channels in accordance with
the experimental measurements reported in [8].

Fig. 3 illustrates the spoofing detection performance of different methods for
the above scenario averaged over $10^{5}$ CFRs from each Alice-Bob and Eve-Bob
channel at the test stage, where the x-axis denotes the number of snapshots
observed during the training stage. In particular, we have evaluated the
performance of HyPhyLearn for this problem, as described in Section IV-D, and
compared it with other classifiers designed based on the likelihood ratio test
with plug-in estimates or existing ML algorithms. By looking at the resulting
spoofing detection accuracy, it can be seen that the performance of the ML
algorithms based on support vector machine (SVM) and Gaussian mixture model
(GMM) is limited in this case due to limited (and mislabeled) training data.
We note that the GMM is used as a classifier here by assigning labels to the
clusters using the available labels corresponding to the reference CFRs.
Specifically, we have used the radial basis function kernel [37] for the SVM
and two components for the GMM for these simulations. Furthermore, one can see
that the LRT method obtained in Section IV-B can improve upon the performance
of these ML algorithms by plugging the estimated parameters, as in Section
IV-C, in the statistical parametric models. In these experiments, we also use
the shrinkage method [47] which improves the covariance matrix estimation for
each likelihood function. For this method, a performance gain can be observed
for this approach in comparison to the no shrinkage case, assuming the
shrinkage parameter $\alpha$ is clarivoyantly chosen to maximize the spoofing
detection accuracy over the test dataset. This method is labeled as ‘LRT (best
shrinkage)’ in Fig. 3. However, in practice the parameter $\alpha$ has to be
estimated from the training data, which—as shown in the figure with label ‘LRT
(shrinkage)’—could deteriorate the LRT performance as the available data
includes mislabeled samples.

Furthermore, we evaluate the performance of an existing hybrid classification
approach known as fine tuning [2, 24] in DTL literature for this problem. In
this method, we first generate $5\times 10^{5}$ synthetic data samples using
the available likelihood parametric functions with plugged-in estimates. Then,
a neural network with $3$ hidden layers of $400$ neurons each is trained to
classify the synthetic data for this example. The training data are used
afterwards to refine the weights of this neural network. Notably, HyPhyLearn
is shown to outperform the aforementioned existing classification methods by
relying on both available and synthetic data and jointly using them in a
learning-based classifier.

For the sake of comparison, we have also considered a variation of HyPhyLearn
that relies on a generative adversarial network (GAN) for generating synthetic
data, i.e., it disregards the available physic-based models. We have observed
that the performance of this approach is impacted in the limited data regime
as GANs rely merely on the available training data for generating further
synthetic data of similar distribution. In fact, for this example, we have
verified that HyPhyLearn based on GAN needs to be trained on $20000$ data
samples in order to achieve the same level of spoofing detection accuracy as
HyPhyLearn based on physics-based models with $4000$ samples. Regrading the
specifics of GAN, we have used a DNN of two hidden layers with $200$ neurons
each as the generator, and a DNN with three hidden layers with $300$ neurons
each as the discriminator. In our implementation of HyPhyLearn, the number of
generated synthetic data samples is set to $4\times 10^{5}$. We have also used
NNs with $3$ hidden layers of $400$ neurons each for $\mathbf{M}_{\bm{\psi}}$
and $h_{\bm{\phi}_{1}}$, while a NN with one hidden layer of $40$ neurons each
is used for ${d}_{\mathbf{\zeta}}$. For all hidden layers, the ReLU activation
function is used. Furthermore, Adam optimizer [37] with a learning rate of
$0.0001$ is used for training in this example. We also note that the optimal
Bayes decision rule, which relies on the knowledge of the true parameters,
results in the spoofing detection accuracy of $0.996$.

Figure 3: Spoofing detection accuracy for different classification algorithms
as a number of available training data for the case when training and test
stage belong to the same coherence time.

Next, we consider a more realistic scenario where the channels’ variations
cause the training and test stage to not fall in the same coherence time. In
this case, Bob uses the heuristic test (32) for some time as it does not have
access to the channel parameters in this period. Afterwards, it uses the data
collected in the previous coherence times to estimate the channel parameters
for the current one. Fig. 4 depicts this setting where the training stage
consists of $n_{c}$ coherence times corresponding to the Alice–Bob channel.
Furthermore, in contrast to Alice, Eve’s transmissions are assumed to be
intermittent due to the uncertainty associated with Eve’s behaviour. During
each coherence time corresponding to the Alice-Bob channel, it is assumed that
Bob collects $100$ training data. Then, the estimation technique described in
Section IV-C is utilized to estimate the channel parameters under each
coherence time. Fig. 5 demonstrates the system performance as a function of
number of coherence times in the training stage. Regarding the physical setup,
we have used the same system parameters as those in Fig. 3, and assumed that
the coherence time of the Alice–Bob channel is $4$ times that of the Eve–Bob
channel for illustrative purpose. For DTL fine-tuning approach and HyPhyLearn,
the number of synthetic data generated for each behavior in a coherence time
is set to $20000$. For these two learning-based approaches, the training
specifications for are chosen to be the same as the ones used in Fig. 3. The
performance comparison again highlights the superiority of HyPhyLearn in
comparison to the existing statistical and data-driven methods.

Figure 4: Training and testing stages for the spoofing detection problem.
$T_{C}^{AB}$ and $T_{C}^{EB}$ denote the coherence time corresponding to the
Alice-Bob and Eve-Bob channels, respectively. The green bar indicates the time
interval within which a snapshot is observed by Bob. Figure 5: Spoofing
detection accuracy for the case where Bob collects training data during
certain number of coherence times before employing a classification algorithm.

###  VI-B Multi-user detection problem

In this section we present results of numerical simulations to investigate the
effectiveness of HyPhyLearn described in Section V-C for the MUD problem. We
choose the simulation parameters based on the setting described in [44] and
consider a system with processing gain of $N=32$ where the number of users is
either $K=3$ or $K=5$. Golden codes of length $32$ are used by the BS as the
pseudo-noise code in (58) and the users’ amplitudes ($A_{k}$’s) are set to
$2$. In addition, a chip interval of length $T_{c}=0.001$ and a sampling rate
of $2/T_{c}$ is employed. A near-far ratio (NFR) of $10$ dB is assumed, which
means the users’ amplitude are randomly unbalanced around $2$ with a variance
of $\pm 5$ dB. The fading channel between the users and the BS consists of $3$
paths, which makes the total number of unknown parameters in Section V-B to be
$9K$. We further consider a setting where the BS might not have access to the
perfect knowledge of the pseudo-noise sequences for all the users at the time
of detection, which would lead to a mismatched situation. To account for this
phenomenon, we introduce a parameter $\rho$ that in order to quantify the
averaged error in the pseudo-noise sequences at the BS while decoding.

As the performance metric, we consider the bit error rate (BER) at the BS
while decoding the users’ information bits, which is of major interest in
digital communication systems. As the MUD algorithm we employ the minimum mean
square error (MMSE) decoder introduced in [48], which is shown to outperform
other existing detection methods including matched filter receiver and box-
constrained maximum likelihood detector [44]. As mentioned in Section V, MUD
can be also solved by a classifier aiming at distinguishing between $2^{K}$
different classes each representing a unique decoded sequence of information
bits. In this case, BER is directly related to the classification accuracy of
the trained classifier. For the asynchronous system discussed in Section V,
the interval $\mathcal{I}_{2p}=[pT_{b},(p+2)T_{b}]$ contains most of the
energy content of the information symbol $b_{k}(p)$. Therefore, it is
sufficient for the MUD detector to process the data in the interval
$\mathcal{I}_{2p}$ in order to obtain estimates of the symbols $b_{k}(p)$,
$\forall k=0,\dots,K-1$.

$\displaystyle
2R_{\mathcal{Z}_{{\bm{\psi},{\bm{\theta}^{*}}}}}(\mathcal{H}_{\Phi})$
$\displaystyle+2R_{\mathcal{Z}_{\bm{\psi},{\hat{\bm{\theta}}}}}(\mathcal{H}_{\Phi})+3\sqrt{(\log{2\delta})/2N_{r}}+3\sqrt{(\log{2\delta})/2N_{s}}\geq$
(65a)
$\displaystyle\sup_{h_{\phi}\in\mathcal{H}_{\Phi}}\bigg{|}\int_{A_{\phi}}p_{{\bm{\psi},{\bm{\theta}^{*}}}}(\mathbf{z})d\mathbf{z}-\sum_{i=1}^{N_{r}}\mathds{1}_{\\{h_{\phi}(\mathbf{z}_{r,i})=1\\}}\bigg{|}+\sup_{h_{\phi}\in\mathcal{H}_{\Phi}}\bigg{|}\int_{A_{\phi}}p_{\bm{\psi},{\hat{\bm{\theta}}}}(\mathbf{z})d\mathbf{z}-\sum_{i=1}^{N_{s}}\mathds{1}_{\\{h_{\phi}(\mathbf{z}_{s,i})=1\\}}\bigg{|}\geq$
(65b)
$\displaystyle\sup_{h_{\phi}\in\mathcal{H}_{\Phi}}\bigg{|}\int_{A_{\phi}}p_{{\bm{\psi},{\bm{\theta}^{*}}}}(\mathbf{z})d\mathbf{z}-\sum_{i=1}^{N_{r}}\mathds{1}_{\\{h_{\phi}(\mathbf{z}_{r,i})=1\\}}-\bigg{(}\int_{A_{\phi}}p_{\bm{\psi},{\hat{\bm{\theta}}}}(\mathbf{z})d\mathbf{z}-\sum_{i=1}^{N_{s}}\mathds{1}_{\\{h_{\phi}(\mathbf{z}_{s,i})=1\\}}\bigg{)}\bigg{|}\geq$
(65c)
$\displaystyle\sup_{h_{\phi}\in\mathcal{H}_{\Phi}}\bigg{|}\int_{A_{\phi}}p_{{\bm{\psi},{\bm{\theta}^{*}}}}(\mathbf{z})d\mathbf{z}-\int_{A_{\phi}}p_{\bm{\psi},{\hat{\bm{\theta}}}}(\mathbf{z})d\mathbf{z}\bigg{|}-\sup_{h_{\phi}\in\mathcal{H}_{\Phi}}\bigg{|}\sum_{i=1}^{N_{r}}\mathds{1}_{\\{h_{\phi}(\mathbf{z}_{r,i})=1\\}}-\sum_{i=1}^{N_{s}}\mathds{1}_{\\{h_{\phi}(\mathbf{z}_{s,i})=1\\}}\bigg{|}=$
(65d) $\displaystyle
d_{\mathcal{A}_{\Phi}}\big{(}p_{{\bm{\psi},{\bm{\theta}^{*}}}}(\mathbf{z}),p_{\bm{\psi},{\hat{\bm{\theta}}}}(\mathbf{z})\big{)}-\hat{d}_{\mathcal{A}_{\Phi}}(\mathcal{Z}_{{\bm{\psi},{\bm{\theta}^{*}}}},\mathcal{Z}_{\bm{\psi},{\hat{\bm{\theta}}}}),$
(65e)

We present simulation results for the performance of the MMSE detector in the
above setting in Fig. 6, and compare it with our proposed approach in Section
V-C. Specifically, the parameter estimation procedure for HyPhyLearn is done
under two different levels of model mismatch, i.e., $\rho=0.2$ and
$\rho=0.25$. Furthermore, the number of training data available from each user
$N_{T}$ is set to $40$. As a general observation, Fig. 6 demonstrates that the
performance of all the detectors is deteriorated as the number of users and
the value of $\rho$ is increased. The perfect MMSE is referred to the case
where the true pseudo-noise sequences are assumed to be known as part of the
implementation of the decoder. In particular, huge performance gap between the
perfect MMSE and the MMSE decoder indicates the high sensitivity of the MMSE
detector to the mismatch. On the other hand, it is also highlighted that our
proposed approach can achieve a substantial gain over a wide range of SNRs by
dealing with the mismatch problem. For HyPhyLearn, the number of generated
synthetic data is set to $10^{6}$ for this example. We have also used NNs with
$4$ hidden layers of $300$ neurons each for $\mathbf{M}_{\bm{\psi}}$ and
$h_{\bm{\phi}_{1}}$ here. Also, a shallow NN with one hidden layer of $40$
neurons is used for ${d}_{\mathbf{\zeta}}$, while ReLU activation function is
used for all the hidden layers. During training, Adam optimizer with a
learning rate of $0.0001$ is utilized as the stochastic gradient descent
algorithm.

Figure 6: BER performance of the MMSE multi-user detector and HyPhyLearn as a
function of SNR. The results are provided for two different parameters, i.e.,
the number of users ($K$) and the mismatch parameter ($\rho$).

In Fig. 7, the BER performance of the multi-user detectors is investigated as
a function of number of available training data. For this example, SNR at the
BS is assumed to be fixed at the BS according to $8$ dB. It is demonstrated
that increasing the number of data samples does not lead to substantial
performance improvements in the case of MMSE method. This is attributed to the
aforementioned mismatch phenomenon in the pseudo-noise sequences which
prevents the MMSE detector from benefiting from the larger amount of data
considerably. Furthermore, it is further shown that the performance gap
between HyPhyLearn and the perfect MMSE shrinks as the number of data
increases. However, the degree to which this gap decreases is higher for the
case of $\rho=0.1$ in comparison to that of $\rho=0.25$. Indeed, HyPhyLearn
gets more benefit from the data at lower levels of mismatch where the
parameter estimates enjoy higher levels of accuracy.

Figure 7: BER performance of the MUD as a function of the number of training
data available at each user.

##  VII Conclusions

We have considered the problem of hypothesis testing in the context of
parametric classification where there is known model for each hypotheses but
the corresponding parameters are unknown. Towards designing a classifier in
this setting, we have taken into account several practical considerations
including the assumptions that available training data are limited and there
could be labeling errors associated with them. Furthermore, the model under
each hypothesis is assumed to be complex such that the MLEs of its parameters