# Information Leakage Detection through Approximate Bayes-optimal Prediction

Pritha Gupta<EMAIL_ADDRESS>
Department of Computer Science, Paderborn University
Paderborn, Germany Marcel Wever<EMAIL_ADDRESS>
MCML, LMU Munich
Munich, Germany Eyke Hüllermeier<EMAIL_ADDRESS>
MCML, LMU Munich
Munich, Germany

###### Abstract

In today’s data-driven world, the proliferation of publicly available
information intensifies the challenge of information leakage (IL), raising
security concerns. IL involves unintentionally exposing secret (sensitive)
information to unauthorized parties via systems observable information.
Conventional statistical approaches rely on estimating mutual information (MI)
between observable and secret information for detecting IL, facing challenges
such as the curse of dimensionality, convergence, computational complexity,
and MI misestimation. Furthermore, though proven effective, emerging
supervised machine learning (ML) methods are limited to binary system
sensitive information and lack a comprehensive theoretical framework. To
address these limitations, we establish a theoretical framework using
statistical learning theory and information theory to quantify and detect IL
accurately. We demonstrate that MI can be accurately estimated by
approximating the log-loss and accuracy of the Bayes predictor. As in
practice, the Bayes predictor is typically unknown, and we propose to
approximate it with the help of automated machine learning (AutoML). First, we
compare our MI estimation approaches against current baselines, using
synthetic data sets generated using the multivariate normal (MVN) distribution
with known MI. Second, we introduce a cut-off technique using one-sided
statistical tests to detect IL and employ the Holm-Bonferroni correction to
increase confidence in detection decisions. Our study evaluates IL detection
performance on real-world data sets, highlighting the effectiveness of Bayes
predictor’s log-loss’s estimation, and finds our proposed method to
effectively estimate MI on synthetic data sets and thus detect ILs accurately.

Keywords: Information Leakage Detection, Mutual Information, Bayes-optimal
Predictor, AutoML, Statistical Tests, Holm-Bonferroni, Privacy

##  1 Introduction

IL has become a significant and complex challenge in today’s data-driven world
(Shabtai et al., 2012). The rapid proliferation of publicly available data,
coupled with the increasing use of Internet of Things (IoT) technologies, has
magnified the threat of IL, posing a substantial risk to the security and
confidentiality of systems (Shabtai et al., 2012; Kelsey, 2002). This
heightened threat of IL is amplified by the increasing use of Internet of
Things (IoT) technologies, posing substantial risks to the security and
confidentiality of systems.

IL occurs when sensitive or confidential information is inadvertently exposed
to unauthorized individuals through observable information of a system
(Hettwer et al., 2019). This unintended disclosure of sensitive data can lead
to severe consequences, ranging from potential electrical blackouts to the
theft of critical information such as medical records and military secrets
(Hettwer et al., 2019; Shabtai et al., 2012). Consequently, the efficient
detection and quantification of IL are of paramount importance. The
theoretical foundations of quantifying IL in systems stem from the field of
information theory, proposes to estimate the MI between observable and secret
information (Chatzikokolakis et al., 2010). MI is shown to be difficult to
compute specifically for high-dimensional data, despite being a pivotal
measure (Gao et al., 2015; Eirola et al., 2014). Traditional solutions often
resort to using statistical estimation to calculate MI, but often encounter
challenges when dealing with high-dimensional data, a problem known as the
_curse of dimensionality_ (Gao et al., 2015; Kraskov et al., 2004). Recent
works offer more robust non-parametric approaches with improved convergence
rates for MI estimation in complex datasets (Singh and Póczos, 2016; Kandasamy
et al., 2015; Moon et al., 2021). However, their implementation in high-
dimensional scenarios remains challenging due to the sophisticated statistical
techniques involved (Kandasamy et al., 2015).

In recent years, ML techniques have become very popular in the field of
information leakage detection (ILD), particularly in the sub-field of
performing side-channel attacks on cryptographic systems (Picek et al., 2023).
These systems unintentionally release the _observable information_ via many
modes, such as network messages, CPU caches, power consumption, or
electromagnetic radiation called the side-channels (Picek et al., 2023). These
modes are exploited by the SCAs to reveal the secret inputs (secret keys,
plaintexts) to the adversaries, potentially rendering all implemented
cryptographic protections irrelevant (Moos et al., 2021; Hettwer et al.,
2019). A system is considered to contain a side-channel if it is susceptible
to SCAs, thereby making it _vulnerable_. Therefore, detecting the existence of
a side-channel in a cybersecurity system is equivalent to uncovering IL
(Hettwer et al., 2019).

In this field, the most relevant literature uses ML to perform SCAs, not
preventing side-channels through early detection of ILs (Hettwer et al.,
2019). Notably, current ML-based methods in this realm detect side-channels,
preventing SCAs and protecting the system on both algorithmic and hardware
levels (Mushtaq et al., 2018; Moos et al., 2021). These approaches leverage
the observable information of a system to classify it as vulnerable (with IL)
or non-vulnerable (without IL) (Perianin et al., 2020; Mushtaq et al., 2018).
These techniques extract observable information from secure systems used as
input, categorizing the systems as non-vulnerable (labeled as $0$).
Subsequently, they intentionally introduce known ILs in these systems,
categorizing them as vulnerable (labeled as $1$). Finally, the extracted
observable information with corresponding class labels is used to produce the
binary classification data set for the learning model. However, this approach
confines these techniques to domain-specific scenarios and cannot be easily
transferred to detect other unknown leakage patterns (Perianin et al., 2020).

ML-based techniques have shown promising potential in estimating MI within
classification system data sets (Belghazi et al., 2018; Qin and Kim, 2019;
Cristiani et al., 2020). Despite their promise, these methods grapple with
challenges related to convergence issues and high computational complexity
(Belghazi et al., 2018; Cristiani et al., 2020). For instance, the approach
proposed by Qin and Kim (2019) introduces a lower bound on MI using the
Kullback-Leibler (KL) divergence. However, this method tends to underestimate
MI and might not capture certain subclasses of IL. Recent advancements have
demonstrated the effectiveness of ML-based techniques in directly detecting IL
by analyzing the accuracy of the supervised learning models on extracted
system data (Moos et al., 2021; Drees et al., 2021). Yet, these methods
exhibit limitations in handling imbalanced and noisy real-world data sets,
commonly encountered in practical scenarios, and tend to miss ILs by producing
false negatives (Zhang et al., 2020; Picek et al., 2018).

To address these limitations, in our prior work, we propose utilizing weighted
binary classifiers integrated with statistical tests as Fisher’s exact test
(FET) and paired t-test (PTT) to account for imbalance (Gupta et al., 2022).
To mitigate noise, an ensemble of binary classifiers, including a deep multi-
layer perceptron, along with their derived results ($p$-values) from the
statistical tests, are aggregated using Holm-Bonferroni correction to enhance
ILD accuracy and increase confidence in systems. Despite its merits, this
approach remains confined to binary classification tasks and lacks a
comprehensive theoretical framework.

###  1.1 Our Contributions

We offer the following contributions in this work to address the existing
literature’s shortcomings.

  * •

In Section 2, we introduce a novel theoretical framework for addressing the
ILD problem in systems that generate classification data. In this context, we
draw upon the core concepts from the information theory and statistical
learning theory to leverage the connection between MI and the performance of
Bayes predictor, as explained in Sections 2.2.1 and 2.2.2, respectively. We
also quantify the IL and formalize the conditions under which it exists in a
system from both the statistical learning theory and information theory
perspectives in Section 2.3.

  * •

In Section 3, we propose two distinct approaches for quantifying IL using MI
within a given system based on approximating the Bayes predictor’s accuracy
and the log-loss. We employ probability calibration techniques to enhance the
precision of Log-Loss estimation, described in Section 3.2.4. We describe the
three state-of-the-art baselines in Section 3.3 (Eirola et al., 2014; Belghazi
et al., 2018; Qin and Kim, 2019). We propose to approximate the Bayes
predictor’s accuracy and the log-loss using the competitive AutoML tools,
described in Section 4. Section 5 describes the rigorous empirical evaluation
of our methods against three state-of-the-art approaches, which involves
applying them to synthetic data sets with known MI and utilizing the
normalized mean absolute error (NMAE) performance metric for assessment.

  * •

In Section 6.2, we formalize a technique for ILD using a cut-off on estimated
MI through the one-sample t-test (OTT) statistical test. We discuss
classification-based ILD approaches using FET and PTT statistical tests in
Section 6.3. Our prior work demonstrates the competence of these approaches in
accurately detecting padding side-channel leaks within OpenSSL TLS servers
(Gupta et al., 2022). To enhance the confidence of IL detection decisions, we
propose employing the Holm-Bonferroni correction on multiple models’
(pipelines’) estimates, as described in Section 6.4. This approach enhances
robustness against noise in estimating MI (e.g., accuracy, confusion matrix)
and variations in the quality of AutoML pipelines.

  * •

Finally, in Section 7, we conduct an extensive benchmarking study, comparing
our ILD methods against state-of-the-art approaches, including FET and PTT,
for detecting timing side-channels (ILs) to counter Bleichenbacher attacks.

##  2 Information Leakage Problem

This section will briefly explain the task of ILD using a mapping that
categorizes a given system as vulnerable and non-vulnerable. Additionally, we
will define a measure to quantify IL and formalize the conditions under which
IL occurs in a system using the foundational concepts of information theory
and statistical learning theory described in Section 2.2. The notations used
throughout the paper are listed in Appendix A.

###  2.1 Information Leakage Detection

ILD stands as a pivotal task in modern cybersecurity, aiming to identify any
unintended disclosure of confidential data or sensitive information through a
system’s emitted _observable information_. The ILD process involves a leakage
detector (ILD algorithm) analyzing the extracted system data set, denoted as
$\mathcal{D}=\\{(\boldsymbol{x}_{i},y_{i})\\}_{i=1}^{N}\subset\mathcal{X}\times\mathcal{Y}$
(size $N\in\mathbb{N}$), where $\mathcal{X}=\mathbb{R}^{d}$ denotes the space
of observable information and $\mathcal{Y}=[M]=\\{1,2,\ldots,M\\}$ represents
the secret information in terms of categorical class labels. The goal of ILD
boils down to associating this data set $\mathcal{D}$ with a binary label,
where $0$ indicates an absence of leakage (no IL) and $1$ suggests the
presence of leakage (IL) in the system. This objective is represented by the
mapping $L$ as

$L:\bigcup_{N\in\mathbb{N}}(\mathcal{X}\times\mathcal{Y})^{N}\rightarrow\\{0,1\\}\,,$

which takes a data set $\mathcal{D}$ (extracted from the system) of any size
as input and provides an assessment of potential IL in the system as an
output. This process is encapsulated in the predicted IL mapping $\hat{L}$
produced by an ILD approach.

Let $\mathcal{L}=\\{(\mathcal{D}_{i},z_{i})\\}_{i=1}^{N_{L}}$ be an IL-
Dataset, such that ${N_{L}}\in\mathbb{N},z_{i}\in\\{0,1\\}$. Further, let
$\boldsymbol{z}=(z_{1},\ldots,z_{N_{L}})$ be the ground-truth vector,
generated by $L$, such that $z_{i}=L(\mathcal{D}_{i}),\forall i\in[N_{L}]$.
Let $\boldsymbol{\hat{z}}=(\hat{z}_{1},\ldots,\hat{z}_{N_{L}})$ be the
corresponding prediction vector, such that
$\hat{z}_{i}=\hat{L}(\mathcal{D}_{i}),\forall i\in[N_{L}]$. Since the ILD task
produces binary decisions, the performance of predicted IL mapping $\hat{L}$
is measured using metrics commonly used in binary classification, e.g.,
accuracy. To avoid confusion, we will refer to $\mathcal{L}$ as IL-Dataset and
each $\mathcal{D}$ as the data set.

###  2.2 Fundamentals of Information and Statistical Learning Theory

This section recalls some fundamental concepts from information theory and
statistical learning theory required to formalize the conditions that signify
a system’s occurrence of IL. We provide a formal description of MI rooted in
information theory and how it is estimated for a given system data set.
Additionally, we formally describe the classification problem with the
relevant evaluation metrics from statistical learning theory. Building upon
this, we define the concept of Bayes predictor and the classification problem
handled by estimating class (conditional) class probabilities.

####  2.2.1 Mutual Information

MI, a fundamental concept in information theory, measures the extent to which
knowledge of one random variable informs about another, thereby quantifying
their degree of dependence (Cover and Thomas, 2005). Consider any pair of
random variables $X$ and $Y$ with joint distribution $P_{(X,Y)}(\cdot)$ on
$\mathcal{X}\times\mathcal{Y}$. Although these random variables can be of any
type, we subsequently assume that $X$ is a continuous $d$-dimensional real-
valued random variable ($\mathcal{X}=\mathbb{R}^{d}$) and $Y$ is a discrete
random variable with $M$ possible values ($\mathcal{Y}=\\{1,\ldots,M\\}$) — in
the information leakage scenario, this will be the relevant case for us.

In the following, for the sake of simplicity, we will use the same notation
for a probability measure and its probability mass or density function. Often,
the symbol for the former is capitalized, whereas small letters are used for
the latter. Assuming that the meaning is always clear from the context, we
will use capital letters throughout. In particular, we will denote the
marginals of the joint distribution $P_{(X,Y)}(\cdot)$ by $P_{X}(\cdot)$ and
$P_{Y}(\cdot)$, and the conditional distributions by $P_{{Y}|{X}}(\cdot)$ and
$P_{{X}|{Y}}(\cdot)$, respectively.

Recall that the entropy of a discrete random variable $Y$ is defined as

$H(Y)=-\sum_{y\in\mathcal{Y}}P_{Y}(Y=y)\cdot\lg(P_{Y}(Y=y))\,,$ (1)

where $0\cdot\lg(0)=0$ by definition; the continuous version — differential
entropy — is essentially obtained by replacing summation with integration. The
maximum entropy is $\lg(M)$ and occurs when all outcomes are equally likely,
indicating complete uncertainty about the outcome of $Y$, i.e.,
$P_{Y}(y)=1/M,\,\forall y\in[M]$. The minimum entropy is $0$, indicating
complete certainty about $Y$, and occurs for the case of a Dirac measure (one
outcome $y=m$ has a probability $P_{Y}(m)=1$, while all others have a
probability of $0$, i.e., $P_{Y}(i)=0\,,\forall i\in[M]\setminus\\{m\\}$). The
conditional entropy of $Y$ given $X$ is defined as

$H({Y}\,|\,{X})=-\int\limits_{\boldsymbol{x}\in\mathcal{X}}P_{X}(\boldsymbol{x})\sum_{y\in\mathcal{Y}}P_{{Y}|{X}}({Y=y}\,|\,{X=\boldsymbol{x}})\cdot\lg(P_{{Y}|{X}}({Y=y}\,|\,{X=\boldsymbol{x}}))\,\mathrm{d}\boldsymbol{x}\,.$

According to this definition, conditional entropy measures the residual
uncertainty in one random variable given knowledge of the other — more
specifically, it measures the _expected_ residual uncertainty, with the
expectation taken with respect to the marginal distribution of $Y$.
$H({Y}\,|\,{X})$ is bounded by $H(Y)$, reaching its upper bound when
$H({Y}\,|\,{X=\boldsymbol{x}})$ follows the marginal distribution
$P_{Y}(\cdot)$ for all $\boldsymbol{x}\in\mathcal{X}$, signifying that $X$
does not provide any information about $Y$. The minimum entropy is $0$,
indicating that $X$ completely determines $Y$, and again occurs in the case
where all conditionals $P_{{Y}|{X}}({\cdot}\,|\,{\boldsymbol{x}})$ are Dirac
distributions.

Now, MI between a pair of random variables $(X,Y)$ measures the reduction of
uncertainty about variable $Y$ by observing variable $X$:

$I(X;Y)=H(Y)-H({Y}\,|\,{X})\,.$ (2)

Thus defined, MI is a symmetric measure:

$\displaystyle I(X;Y)$ $\displaystyle=H(Y)-H({Y}\,|\,{X})$
$\displaystyle=H(Y)+H(X)-H(X,Y)$ $\displaystyle=H(X)+H(Y)-H(Y,X)$
$\displaystyle=H(X)-H({X}\,|\,{Y})=I(Y;X)\,.$

Consequently, its value ranges between $0$ and $\min(\\{H(X),H(Y)\\})$. A
value of $0$ indicates that $X$ and $Y$ are completely independent, while the
maximal value $\min(\\{H(X),H(Y)\\})$ signifies full dependence (Cover and
Thomas, 2005). By plugging in the expressions for (conditional) entropy and
rearranging terms, one easily shows that MI equals the KL divergence of the
joint distribution $P_{(X,Y)}(\cdot)$ from the product of the marginals (i.e.,
the joint distribution under the assumption of independence):

$\displaystyle I(X;Y)$ $\displaystyle=H(Y)-H({Y}\,|\,{X})$
$\displaystyle=\int\limits_{\boldsymbol{x}\in\mathcal{X}}\sum_{y\in\mathcal{Y}}{P_{(X,Y)}(X=\boldsymbol{x},Y=y)\cdot\lg{\left({\frac{P_{(X,Y)}(X=\boldsymbol{x},Y=y)}{P_{X}(X=\boldsymbol{x})\cdot
P_{Y}(Y=y)}}\right)}}\,d\boldsymbol{x}\,.$ (3)

The unit of MI varies based on the base of the logarithm: natural logarithm
$\ln$ uses _nats_ , decimal logarithm $\log$ employs _Hartleys_ , and binary
logarithm $\lg$ corresponds to _bits_. We calculate MI using the binary
logarithm $\lg$, thus using _bits_ as the unit of information.

####  2.2.2 Classification Problem

In the realm of classification, the learning algorithm (learner) is provided
with a training data set
$\mathcal{D}=\\{(\boldsymbol{x}_{i},y_{i})\\}_{i=1}^{N}\subset\mathcal{X}\times\mathcal{Y}$
of size $N\in\mathbb{N}$, where $\mathcal{X}=\mathbb{R}^{d}$ is the input
(instance) space and $\mathcal{Y}=\\{1,2,\ldots,M\\}=[M],\,M\in\mathbb{N}$ the
output (categorical class labels) space (Vapnik, 1991), and the
$(\boldsymbol{x}_{i},y_{i})$ are assumed to be independent and identically
distributed (i.i.d.) according to $P_{(X,Y)}(\cdot)$. According to statistical
learning theory, the primary goal of the learner in standard classification is
to induce a hypothesis
$h\colon\mathcal{X}\rightarrow\mathcal{Y},\,h\in\mathcal{H}$, with low
generalization error (risk)

$R(h)=\mathbb{E}[\ell(y,h(\boldsymbol{x}))]=\int\limits_{\mathcal{X}\times\mathcal{Y}}\ell(y,h(\boldsymbol{x}))\,d\,P_{(X,Y)}(X=\boldsymbol{x},Y=y)\,,$
(4)

where $\mathcal{H}$ is the underlying hypothesis space (set of candidate
functions the learner can choose from),
$\ell\colon\mathcal{Y}\times\mathcal{Y}\rightarrow\mathbb{R}$ a loss function,
and $P_{(X,Y)}(\cdot)$ is the joint probability measure modeling the
underlying data-generating process. A loss function commonly used in standard
classification is the $0\text{-}1$ loss
$\ell_{01}(y,\hat{y})\coloneq\llbracket\hat{y}\neq y\rrbracket$, where
$\hat{y}=h(\boldsymbol{x})$ and $\llbracket\cdot\rrbracket$ is the indicator
function as defined in Appendix A. The risk minimizer (denoted by $h^{*}$)
achieves the minimum expected loss $\mathbb{E}[\ell(\cdot)]$ in terms of the
loss function $\ell$ across the entire joint distribution $P_{(X,Y)}(\cdot)$
and is defined as $h^{*}=\operatorname*{\arg\,\min}_{h\in\mathcal{H}}R(h)$.
The measure $P_{(X,Y)}(\cdot)$ in (4) induces marginal probability (density or
mass) functions on $\mathcal{X}$ and $\mathcal{Y}$ as well as a conditional
probability of the class $Y$ given an input instance $\boldsymbol{x}$, i.e.,
$P_{(X,Y)}(\boldsymbol{x},y)=P_{{Y}|{X}}({y}\,|\,{\boldsymbol{x}})\times
P_{X}(\boldsymbol{x})$. In practice, these probabilities are not observed by
the learner, so directly minimizing the risk (4) is not feasible. Instead,
learning in a standard classification setting is commonly accomplished by
minimizing (a regularized version of) _empirical risk_ for $h$:

$R_{\text{emp}}(h)=\dfrac{1}{N}\sum_{i=1}^{N}\ell(y_{i},h(\boldsymbol{x}_{i}))\,.$
(5)

In the subsequent discussions, we denote by
$g=\operatorname*{\arg\,\min}_{h\in\mathcal{H}}R_{\text{emp}}(h)$ the
(learned) hypothesis that minimizes (5), i.e., the empirical risk minimizer
(Vapnik, 1991). In practice, for the available finite sampled data set
$\mathcal{D}$, $g$ is the best possible approximation (empirical estimation)
of the true risk minimizer $h^{*}$ and serves as an appropriation thereof.

2.2.2.1 Bayes-optimal Predictor In statistical learning theory, the _Bayes
predictor_ is the optimal (pointwise) classification function
$g^{b}\colon\mathcal{X}\rightarrow\mathcal{Y}$ producing a minimum expected
risk, i.e., minimizing (4) for a given loss function $\ell(\cdot)$:

$\displaystyle g^{b}(\boldsymbol{x})$
$\displaystyle=\operatorname*{\arg\,\min}_{\hat{y}\in\mathcal{Y}}\sum_{y\in\mathcal{Y}}\ell(y,\hat{y})\,\cdot
P_{{Y}|{X}}({Y=y}\,|\,{X=\boldsymbol{x}})=\operatorname*{\arg\,\min}_{\hat{y}\in\mathcal{Y}}\mathbb{E}_{y}[{\ell(y,\hat{y})}\,|\,{\boldsymbol{x}}]\,,$

where $\mathbb{E}_{y}[\ell]$ is the expected loss of the prediction $\hat{y}$
with respect to $y\in\mathcal{Y}$, and $P_{{Y}|{X}}({y}\,|\,{\boldsymbol{x}})$
is the conditional probability of the class $y$ given an input instance
$\boldsymbol{x}$ (Devroye et al., 1996). In the case where the loss is
$0\text{-}1$, the Bayes predictor simplifies to

$g^{bc}(\boldsymbol{x})=\operatorname*{\arg\,\max}_{y\in\mathcal{Y}}P_{{Y}|{X}}({Y=y}\,|\,{X=\boldsymbol{x}})\,.$
(6)

The expected loss of this predictor is known as the _Bayes error-rate_
$\operatorname{m}_{\textsc{{Err}}}{(\textsc{{Bc}})}$. To determine the Bayes
predictor, the conditional probability distribution $P_{{Y}|{X}}(\cdot)$ must
be known for the underlying data set $\mathcal{D}$. In the case where
$\mathcal{X}$ and $\mathcal{Y}$ are independent of each other, i.e.,
$P_{{Y}|{X}}({y}\,|\,{\boldsymbol{x}})$ corresponds to the marginal
distribution $P_{Y}(y)$ as input features are completely uninformative, we
have

$g^{bc}(\boldsymbol{x})\equiv
g^{mc}(\boldsymbol{x})\equiv\operatorname*{\arg\,\max}_{y\in\mathcal{Y}}P_{Y}(Y=y)\,.$
(7)

We refer to this predictor as the _marginal Bayes predictor_ on $\mathcal{X}$.
The marginal Bayes predictor assigns to each input $\boldsymbol{x}$ a class
label from the set of labels with the highest marginal probability. Strictly
speaking, as the maximum in (7) is not necessarily unique, the marginal Bayes
predictor may pick any label with the highest probability — in that case, we
nevertheless assume that it picks the same label for every $\boldsymbol{x}$,
so that it is a constant function.

2.2.2.2 Probabilistic Classification Distinct from standard classifiers,
probabilistic classifiers focus on estimating the (conditional) class
probabilities $P_{{Y}|{X}}({y}\,|\,{\boldsymbol{x}})$ for each class $y$ in
$\mathcal{Y}$, given an input instance $\boldsymbol{x}\in\mathcal{X}$. We
denote predictions of that kind by
$\hat{P}_{{Y}|{X}}({y}\,|\,{\boldsymbol{x}})$. As before, training data comes
in the form
$\mathcal{D}=\\{(\boldsymbol{x}_{i},y_{i})\\}_{i=1}^{N}\subset\mathcal{X}\times\mathcal{Y}$,
where the $(\boldsymbol{x}_{i},y_{i})$ are i.i.d. according to
$P_{(X,Y)}(\cdot)$, and the goal to induce a hypothesis
$h_{p}\colon\mathcal{X}\rightarrow\mathbb{P}(\mathcal{Y})$ with low
generalization error (risk)

$R_{p}(h_{p})=\mathbb{E}[\ell_{p}(y,h_{p}(\boldsymbol{x}))]=\int\limits_{\mathcal{X}\times\mathcal{Y}}\ell_{p}(y,h_{p}(\boldsymbol{x}))\,d\,P_{(X,Y)}(X=\boldsymbol{x},Y=y)\,.$

Now, however, instead of comparing a predicted class $\hat{y}$ with a true
class $y$, the loss $\ell_{p}$ compares a predicted probability distribution
with $y$ — thus, the loss is a mapping
$\ell_{p}\colon\mathcal{Y}\times\mathbb{P}(\mathcal{Y})\rightarrow\overline{\mathbb{R}}$,
where $\mathbb{P}(\mathcal{Y})$ denotes the set of probability mass functions
on $\mathcal{Y}$.

Since $\mathcal{Y}$ is finite and consists of $M$ classes,
$\mathbb{P}(\mathcal{Y})$ can be represented by the $(M-1)$-simplex, i.e.,
predictions can be represented as probability vectors
$h_{p}(\boldsymbol{x})=\boldsymbol{\hat{p}}=(\hat{p}_{1},\ldots,\hat{p}_{M})$,
where
$\hat{p}_{m}=\boldsymbol{\hat{p}}[m]=\hat{P}_{{Y}|{X}}({m}\,|\,{\boldsymbol{x}})$
is the probability assigned to class $m$. Typically, learning predictors of
that kind involve minimizing the _empirical risk_

$R_{{\text{emp}}|{p}}(h_{p})=\dfrac{1}{N}\sum_{i=1}^{N}\ell_{p}(y_{i},h_{p}(\boldsymbol{x}_{i}))\,.$
(8)

Subsequently, we denote the empirical risk minimizer by
$g_{p}=\operatorname*{\arg\,\min}_{h_{p}\in\mathcal{H}_{p}}R_{{\text{emp}}|{p}}(h_{p})$
(Vapnik, 1991; Bishop, 2006, chap. 4). In the case where $Y$ is independent of
$X$, the marginal Bayes predictor or marginal classifier (denoted by
$g_{p}^{mc}$) is again the best constant probability predictor, i.e., the one
with the lowest risk (8) among all constant predictors.

_Proper Scoring Rules_ So-called (strictly) _proper scoring rules_ constitute
an important class of loss functions is probabilistic classification (Gneiting
and Raftery, 2007). Roughly speaking, such losses incentivice the learner to
predict accurate probabilities. Formally, a loss $\ell_{p}(\cdot)$ is a proper
scoring rule if

$\mathbb{E}_{Y\sim\boldsymbol{p}}\,\ell_{p}(Y,\boldsymbol{p})\leq\mathbb{E}_{Y\sim\boldsymbol{p}}\,\ell_{p}(Y,\hat{\boldsymbol{p}})$

for all distributions $\boldsymbol{p},\hat{\boldsymbol{p}}$, and a _strictly_
proper scoring rule is the above inequality is strict whenever
$\boldsymbol{p}\neq\hat{\boldsymbol{p}}$. Thus, by predicting the true
distribution $\boldsymbol{p}$ of the (categorical) random variable $Y$, the
learner minimizes its loss in expectation. Important examples of strictly
proper scoring rules include the Brier score and the categorical cross-entropy
(CCE). The latter is defined as
$\ell_{\textsc{{CCE}}}(y,\boldsymbol{\hat{p}})\coloneq-\ln(\hat{p}_{y})$
(Bishop, 2006, chap. 4). The CCE loss (logarithmic scoring rule) is widely
recognized for its information-theoretic interpretations and practical
effectiveness (Good, 1992; Roulston and Smith, 2002; Daley and Vere-Jones,
2004). It can also be motivated in the realm of Neyman–Pearson theory
(Feuerverger and Rahman, 1992).

_Probability Calibration_ Despite its sound theoretical grounding, the
(empirical) minimization of a proper scoring rule is not guaranteed to produce
accurate probability estimates in practice. Instead, predictions are often
noisy or systematically biased — neural networks, for instance, tend to be
over-confident and predict probabilities that are biased toward the extremes
(Szegedy et al., 2016). So-called _calibration_ methods seek to improve
probability estimates through suitable post-hoc adjustments. Broadly speaking,
such methods leverage validation data to learn calibration functions in the
form of (monotone) mappings from predicted probabilities to improved
probabilities (Silva Filho et al., 2023). In this study, five prevalent
techniques for calibrating multi-class probabilities will be used and briefly
described in Section 3.2.4.

_Deterministic Prediction_ Obviously, a probabilistic prediction may also
serve as a basis for a deterministic prediction, if needed. Typically, the
class label with the highest (predicted) probability is adopted in that case,
i.e.,
$\hat{y}=\operatorname*{\arg\,\max}_{y\in\mathcal{Y}}\hat{P}_{{Y}|{X}}({y}\,|\,{\boldsymbol{x}})$.
As readily seen, this prediction minimizes the standard $0\text{-}1$ loss in
expectation. The other way around, this also shows that to perform well in
terms of standard $0\text{-}1$ loss, the learner merely needs to identify the
most probable class, or, stated differently, strong performance can even be
achieved with relatively inaccurate estimates, provided the highest predicted
probability is still assigned to the truly most probable class label (Domingos
and Pazzani, 1997; Cheng and Hüllermeier, 2012).

2.2.2.3 Performance Evaluation Metrics Koyejo et al. (2015) and Powers (2011)
define evaluation measures used for evaluating the performance of classifiers,
using the ground-truth labels denoted by $\boldsymbol{y}=(y_{1},\ldots,y_{N})$
for a given $\mathcal{D}=\\{(\boldsymbol{x}_{i},y_{i})\\}_{i=1}^{N}$ and the
predictions denoted by the vector
$\boldsymbol{\hat{y}}=(\hat{y}_{1},\ldots,\hat{y}_{N})$. The predictions could
be obtained by the empirical risk minimizer $g$ derived using (5) for a
standard classification, such that $\hat{y}_{i}=g(\boldsymbol{x_{i}}),\forall
i\in[N]$ or by the deterministic predictions produced by a probabilistic
classifier $g_{p}$ derived using (8), such that
$\hat{y}_{i}=\operatorname*{\arg\,\max}_{m\in
M}g_{p}(\boldsymbol{x_{i}})[m],\forall i\in[N]$.

_Accuracy_ is defined as the proportion of correct predictions

$\operatorname{m}_{\textsc{{Acc}}}(\boldsymbol{y},\boldsymbol{\hat{y}})\coloneq\dfrac{1}{N}\sum_{i=1}^{N}\llbracket
y_{i}=\hat{y}_{i}\rrbracket\,.$

_Error-rate_ is defined as the proportion of incorrect predictions

$\operatorname{m}_{\textsc{{Err}}}(\boldsymbol{y},\boldsymbol{\hat{y}})\coloneq\dfrac{1}{N}\sum_{i=1}^{N}\llbracket
y_{i}\neq\hat{y}_{i}\rrbracket\,.$

_CM_ Many evaluation metrics are based on the concepts, _true positive_
($\operatorname{m}_{\textsc{{tp}}}$), _true negative_
($\operatorname{m}_{\textsc{{tn}}}$), _false positive_
($\operatorname{m}_{\textsc{{fp}}}$), and _false negative_
($\operatorname{m}_{\textsc{{fn}}}$), which are defined as

$\displaystyle\operatorname{m}_{\textsc{{tn}}}(\boldsymbol{y},\boldsymbol{\hat{y}})=\sum_{i=1}^{N}\llbracket
y_{i}=0,\hat{y}_{i}=0\rrbracket,\,\operatorname{m}_{\textsc{{tp}}}(\boldsymbol{y},\boldsymbol{\hat{y}})=\sum_{i=1}^{N}\llbracket
y_{i}=1,\hat{y}_{i}=1\rrbracket,$
$\displaystyle\operatorname{m}_{\textsc{{fp}}}(\boldsymbol{y},\boldsymbol{\hat{y}})=\sum_{i=1}^{N}\llbracket
y_{i}=0,\hat{y}_{i}=1\rrbracket,\,\operatorname{m}_{\textsc{{fn}}}(\boldsymbol{y},\boldsymbol{\hat{y}})=\sum_{i=1}^{N}\llbracket
y_{i}=1,\hat{y}_{i}=0\rrbracket\,.$

The confusion matrix is defined using these metrics as

$\operatorname{m}_{\textsc{{CM}}}(\boldsymbol{y},\boldsymbol{\hat{y}})=\begin{pmatrix}\operatorname{m}_{\textsc{{tn}}}(\boldsymbol{y},\boldsymbol{\hat{y}})&\operatorname{m}_{\textsc{{fp}}}(\boldsymbol{y},\boldsymbol{\hat{y}})\\\
\operatorname{m}_{\textsc{{fn}}}(\boldsymbol{y},\boldsymbol{\hat{y}})&\operatorname{m}_{\textsc{{tp}}}(\boldsymbol{y},\boldsymbol{\hat{y}})\end{pmatrix}\,.$

_FNR_ is defined as the ratio of _false negatives_ to the number of positive
instances

$\operatorname{m}_{\textsc{{FNR}}}(\boldsymbol{y},\boldsymbol{\hat{y}})=\tfrac{\operatorname{m}_{\textsc{{fn}}}(\boldsymbol{y},\boldsymbol{\hat{y}})}{\operatorname{m}_{\textsc{{fn}}}(\boldsymbol{y},\boldsymbol{\hat{y}})+\operatorname{m}_{\textsc{{tp}}}(\boldsymbol{y},\boldsymbol{\hat{y}})}\,.$

_FPR_ is defined as the ratio of _false positives_ to the number of negative
instances

$\operatorname{m}_{\textsc{{FPR}}}(\boldsymbol{y},\boldsymbol{\hat{y}})=\tfrac{\operatorname{m}_{\textsc{{fp}}}(\boldsymbol{y},\boldsymbol{\hat{y}})}{\operatorname{m}_{\textsc{{fp}}}(\boldsymbol{y},\boldsymbol{\hat{y}})+\operatorname{m}_{\textsc{{tn}}}(\boldsymbol{y},\boldsymbol{\hat{y}})}\,.$

_MCC_ is another balanced measure that considers true and false positives,
including true and false negatives. It produces a value between $-1$ and $+1$,
where $+1$ represents a perfect prediction, $0$ is random, and $-1$ indicates
total disagreement between the predictions and ground truths. It is formally
defined as

$\displaystyle\operatorname{m}_{\textsc{{MCC}}}(\boldsymbol{y},\boldsymbol{\hat{y}})$
$\displaystyle=\left(\tfrac{(\operatorname{m}_{\textsc{{tp}}}(\boldsymbol{y},\boldsymbol{\hat{y}})\times\operatorname{m}_{\textsc{{tn}}}(\boldsymbol{y},\boldsymbol{\hat{y}})-\operatorname{m}_{\textsc{{fp}}}(\boldsymbol{y},\boldsymbol{\hat{y}})\times\operatorname{m}_{\textsc{{fn}}}(\boldsymbol{y},\boldsymbol{\hat{y}}))}{\sqrt{(\operatorname{m}_{\textsc{{tp}}}(\boldsymbol{y},\boldsymbol{\hat{y}})+\operatorname{m}_{\textsc{{fp}}}(\boldsymbol{y},\boldsymbol{\hat{y}}))(\operatorname{m}_{\textsc{{tp}}}(\boldsymbol{y},\boldsymbol{\hat{y}})+\operatorname{m}_{\textsc{{fn}}}(\boldsymbol{y},\boldsymbol{\hat{y}}))}}\right)$
$\displaystyle\quad\times\left(\tfrac{1}{(\operatorname{m}_{\textsc{{tn}}}(\boldsymbol{y},\boldsymbol{\hat{y}})+\operatorname{m}_{\textsc{{fp}}}(\boldsymbol{y},\boldsymbol{\hat{y}}))(\operatorname{m}_{\textsc{{tn}}}(\boldsymbol{y},\boldsymbol{\hat{y}})+\operatorname{m}_{\textsc{{fn}}}(\boldsymbol{y},\boldsymbol{\hat{y}}))}\right)\,.$

_BER_ is the average of both false positive rate (FPR) and false negative rate
(FNR), offering a balanced evaluation metric for a given classification model
as

$\displaystyle\operatorname{m}_{\textsc{{BER}}}(\boldsymbol{y},\boldsymbol{\hat{y}})$
$\displaystyle=\tfrac{1}{2}\left(\tfrac{\operatorname{m}_{\textsc{{fp}}}(\boldsymbol{y},\boldsymbol{\hat{y}})}{\operatorname{m}_{\textsc{{tn}}}(\boldsymbol{y},\boldsymbol{\hat{y}})+\operatorname{m}_{\textsc{{fp}}}(\boldsymbol{y},\boldsymbol{\hat{y}})}+\tfrac{\operatorname{m}_{\textsc{{fn}}}(\boldsymbol{y},\boldsymbol{\hat{y}})}{\operatorname{m}_{\textsc{{tp}}}(\boldsymbol{y},\boldsymbol{\hat{y}})+\operatorname{m}_{\textsc{{fn}}}(\boldsymbol{y},\boldsymbol{\hat{y}})}\right)=\tfrac{\operatorname{m}_{\textsc{{fnr}}}(\boldsymbol{y},\boldsymbol{\hat{y}})+\operatorname{m}_{\textsc{{fpr}}}(\boldsymbol{y},\boldsymbol{\hat{y}})}{2}\,.$

###  2.3 Quantifying and Detecting Information Leakage

This section describes a measure for quantifying IL in systems by
approximating the Bayes predictor, enabling IL detection. IL occurs in a
system when _observable information_ ($\boldsymbol{x}\in\mathcal{X}$,
represented by $X$) is directly or indirectly correlated with its secret
information (e.g., secret keys, plaintexts) ($y\in\mathcal{Y}$, represented by
$Y$), i.e., some information in $\boldsymbol{x}\in\mathcal{X}$ can be
exploited to infer the secret $y\in\mathcal{Y}$ (Kelsey, 2002; Chatzikokolakis
et al., 2010). The presence of IL in a system is determined by analyzing the
system’s classification data set
$\mathcal{D}=\\{(\boldsymbol{x}_{i},y_{i})\\}_{i=1}^{N}\subset\mathcal{X}\times\mathcal{Y}$
if the Bayes predictor (or a suitable proxy $g_{p}$ or $g$) significantly
outperforms the marginal Bayes predictor (or a suitable proxy $g_{p}^{mc}$).
The performance difference quantifies IL within a system and serves as a
foundational condition for its detection. Thus, the ILD involves analyzing the
learnability of $g_{p}$ minimizing (8) (proxy of Bayes predictor) on the given
system data set $\mathcal{D}$.

_IL Quantification_ IL measures the extent of information in the inputs or
observable information $\boldsymbol{x}\in\mathcal{X}$ that can be exploited to
deduce the outputs or secret class labels $y\in\mathcal{Y}$, quantified
directly through MI (Chatzikokolakis et al., 2010). Another method for
quantifying IL in a system uses Bayes predictor $g^{bc}$ to explore
dependencies (correlation) between $\mathcal{X}$ and $\mathcal{Y}$. As shown
in (7), for cases where $\mathcal{X}$ and $\mathcal{Y}$ are independent, the
Bayes predictor reduces to the marginal Bayes predictor. IL is quantified as
the difference in average penalties between the marginal Bayes predictor
$g^{mc}$ and the Bayes predictor $g^{bc}$ using loss functions or metrics like
$0\text{-}1$, CCE, log-loss etc as
$\delta_{\textsc{{Err}}}=\operatorname{m}_{\textsc{{Err}}}{(\textsc{{Mc}})}-\operatorname{m}_{\textsc{{Err}}}{(\textsc{{Bc}})}$.

In practice, since the ground-truth joint probability density function (PDF)
$P_{(X,Y)}(.)$ are seldom observed, we approximate the Bayes predictor and the
marginal Bayes predictor, using the empirical risk minimizers, $g_{p}$ (or
$g$) and $g_{p}^{mc}$ respectively as suitable proxies. When penalty or loss
is evaluated solely based on the classification decisions of the learners, we
may use the empirical risk minimizer $g$ or the decision rule
$\operatorname*{\arg\,\max}_{y\in\mathcal{Y}}g_{p}(\boldsymbol{x})[y]$ of
$g_{p}$ as a proxy for Bayes predictor. However, when using error-rate as the
loss function, $g$ serves as a better proxy as it directly targets and
minimizes the average $0\text{-}1$ loss. In scenarios requiring assessment
based on (conditional) class probabilities, the empirical risk minimizer
$g_{p}$ should be employed.

IL is then approximately quantified as the difference in the penalty for the
give loss function between $g_{p}^{mc}$ and the empirical risk minimizers,
i.e.,
$\delta_{\textsc{{Err}}}\approxeq\operatorname{m}_{\textsc{{Err}}}{(g_{p}^{mc})}-\operatorname{m}_{\textsc{{Err}}}{(g)}\approxeq\operatorname{m}_{\textsc{{Acc}}}{(g)}-\operatorname{m}_{\textsc{{Acc}}}{(g_{p}^{mc})}$
using error-rate (average $0\text{-}1$ loss). Notably, using log-loss as the
loss function estimates MI between $X$ and $Y$ to quantify IL, i.e.,
$\delta_{ll}=I(X;Y)=\mathbb{E}[\ell_{ll}(\textsc{{Mc}})]-\mathbb{E}[\ell_{ll}(\textsc{{Bc}})]\approxeq\mathbb{E}[\ell_{ll}(g_{p}^{mc})]-\mathbb{E}[\ell_{ll}(g_{p})]$
as per (11). We also propose estimating MI by calculating the midpoint between
its upper and lower bounds, defined through Bayes error-rate (proxy
$\operatorname{m}_{\textsc{{Err}}}{(g_{p})}$), as detailed in Section 3.1. Our
Log-Loss and Mid-Point approaches outperform state-of-the-art techniques in
accurately estimating MI on synthetic data set with known MI, as detailed in
Section 5.

_IL Detection_ The presence of IL in a system generating data set
$\mathcal{D}$ is indicated if the approximated MI or $\delta_{(\cdot)}$ is
significantly greater than $0$. Our prior work (Gupta et al., 2022) proposes
classification-based ILD approaches that are built on using average
$0\text{-}1$ loss function or error-rate to quantify and detect IL through
$\delta_{\textsc{{Err}}}$, serving as baselines for this study. We propose MI-
based ILD approaches which employ the OTT on multiple MI estimates to
determine the presence of IL, as explained in Section 6.2. Our ILD approach,
using log-loss for MI estimation, effectively detects leakage in OpenSSL TLS
server implementations as discussed in Section 7.

##  3 Mutual Information Estimation Techniques

In this section, we will introduce the two techniques for estimating the MI in
a system generating classification data set using statistical learning theory
rather than relying solely on standard statistical techniques. We briefly
describe Mid-Point and Log-Loss in Sections 3.1 and 3.2, respectively.
Additionally, Section 3.3 outlines three competitive state-of-the-art MI
estimation approaches (Eirola et al., 2014; Belghazi et al., 2018; Qin and
Kim, 2019). All these techniques are eventually used to perform ILD using OTT
in a given system, as explained in Section 6.2.

###  3.1 Mid-point Estimation

We propose the Mid-Point approach, which leverages the relationship between
Bayes error-rate and MI for a given classification task to estimate the MI and
eventually perform ILD in the given system. To estimate the MI, we use the
_empirical risk minimizer_ $g$ of (5) to approximate Bayes error-rate of the
given system data set. This approach uses lower and upper bounds on the
conditional entropy $H({Y}\,|\,{X})$ which are derived based on the Bayes
error-rate $\operatorname{m}_{\textsc{{Err}}}{(\textsc{{Bc}})}$ and number of
class labels $M$ (Tebbe and Dwyer, 1968).

Derived by Hellman and Raviv (1970) and Tebbe and Dwyer (1968), the lower
bound on the conditional entropy of $Y$ given $X$ is expressed as

$H_{l}(\operatorname{m}_{\textsc{{Err}}}{(\textsc{{Bc}})},M)=\lg(m)+m(m+1)\left(\lg\left(\frac{m+1}{m}\right)\right)\left(\operatorname{m}_{\textsc{{Err}}}{(\textsc{{Bc}})}-\frac{m-1}{m}\right)\,,$

where $H_{l}$ is valid for $\frac{1}{m+1}\leq
1-\operatorname{m}_{\textsc{{Err}}}{(\textsc{{Bc}})}\leq\frac{1}{m}\,,m=1,\ldots,M-1$.

Derived by Fano (1961) and Tebbe and Dwyer (1968), the upper bound on the
conditional entropy of $Y$ given $X$ is expressed as

$H_{u}(\operatorname{m}_{\textsc{{Err}}}{(\textsc{{Bc}})},M)=H_{2}(\operatorname{m}_{\textsc{{Err}}}{(\textsc{{Bc}})})+\operatorname{m}_{\textsc{{Err}}}{(\textsc{{Bc}})}\lg(M-1)\,,$

where $H_{2}(a)=-a\lg(a)-(1-a)\lg(1-a)$ is the binary cross-entropy function.

Using $H_{u}(\cdot)$ and $H_{l}(\cdot)$, the conditional entropy
$H({Y}\,|\,{X})$ is bounded as

$\displaystyle H_{l}(\operatorname{m}_{\textsc{{Err}}}{(\textsc{{Bc}})},M)\leq
H({Y}\,|\,{X})\leq
H_{u}(\operatorname{m}_{\textsc{{Err}}}{(\textsc{{Bc}})},M)$ $\displaystyle
2\operatorname{m}_{\textsc{{Err}}}{(\textsc{{Bc}})}\leq H({Y}\,|\,{X})\leq
H_{2}(\operatorname{m}_{\textsc{{Err}}}{(\textsc{{Bc}})})\,,\quad M=2\,.$

These bounds are substituted in (2) to derive the upper $F_{u}(\cdot)$ and
lower $F_{l}(\cdot)$ bounds on MI, in terms of the Bayes error-rate
$\operatorname{m}_{\textsc{{Err}}}{(\textsc{{Bc}})}$ and number of classes $M$

$\displaystyle
H(Y)-H_{u}(\operatorname{m}_{\textsc{{Err}}}{(\textsc{{Bc}})},M)\leq
I(X;Y)\leq H(Y)-H_{l}(\operatorname{m}_{\textsc{{Err}}}{(\textsc{{Bc}})},M)$
$\displaystyle F_{l}(\operatorname{m}_{\textsc{{Err}}}{(\textsc{{Bc}})},M)\leq
I(X;Y)\leq F_{u}(\operatorname{m}_{\textsc{{Err}}}{(\textsc{{Bc}})},M)$
$\displaystyle
H(Y)-H_{2}(\operatorname{m}_{\textsc{{Err}}}{(\textsc{{Bc}})})\leq I(X;Y)\leq
H(Y)-2\operatorname{m}_{\textsc{{Err}}}{(\textsc{{Bc}})}\,,\quad M=2\,.$

We propose employing the mid-point of the upper $F_{u}(\cdot)$ and the lower
$F_{l}(\cdot)$ bound to estimate the MI, by replacing
$\operatorname{m}_{\textsc{{Err}}}{(\textsc{{Bc}})}$ with the error-rate of
$g$ as

$\widehat{I}(X;Y)\approxeq\frac{F_{u}(\operatorname{m}_{\textsc{{Err}}}{(g)},M)+F_{l}(\operatorname{m}_{\textsc{{Err}}}{(g)},M)}{2}\,.$
(9)

Here $g$ is the _empirical risk minimizer_ of (5) and serves as an appropriate
approximation (proxy) of the Bayes predictor, as
$\operatorname{m}_{\textsc{{Err}}}{(\textsc{{Bc}})}\approxeq\operatorname{m}_{\textsc{{Err}}}{(g)}$.

_Overestimated MI in Imbalanced Data Sets_ As per Figure 1, we can see that in
the case of balanced binary classification, when the mid-point MI is $0$, the
Bayes error-rate is $0.5$, i.e.,
$\operatorname{m}_{\textsc{{Err}}}{(\textsc{{Bc}})}=\operatorname{m}_{\textsc{{Err}}}{(\textsc{{Mc}})}=0.5$
which implies that there exists no IL in the system. But for the imbalanced
case $r=0.3$, when the
$\operatorname{m}_{\textsc{{Err}}}{(\textsc{{Bc}})}=\operatorname{m}_{\textsc{{Err}}}{(\textsc{{Mc}})}=0.3$,
the mid-point MI is greater than $0$, while in this case there exists no IL in
the system according the statistical learning theory condition illustrated in
Section 2.3. This highlights the problem of overestimation associated with the
Mid-Point approach because the ground-truth MI might indeed be equivalent to
the lower bound
($F_{l}(\operatorname{m}_{\textsc{{Err}}}{(\textsc{{Bc}})},M)$).

Figure 1: Bayes error-rate versus MI estimated using Mid-Point approach

To illustrate this issue, let us consider an example where the class
imbalanced joint distribution $P_{(X,Y)}(\cdot)$, along with the calculated
marginals $P_{Y}(\cdot)$, $P_{X}(\cdot)$, and the conditional
$P_{{Y}|{X}}(\cdot)$, is defined as

$\begin{array}[]{|c|c|c|c|}\hline\cr\vrule\lx@intercol\hfil\text{Joint
Distribution}\hfil\lx@intercol\vrule\lx@intercol\\\ \hline\cr
P_{(X,Y)}(x,y)&y=0&y=1&P_{X}(x)\\\ \hline\cr x=0&0.5&0&0.5\\\ \hline\cr
x=1&0.45&0.05&0.5\\\ \hline\cr P_{Y}(y)&0.95&0.05&-\\\
\hline\cr\end{array}\quad\begin{array}[]{|c|c|c|c|}\hline\cr\vrule\lx@intercol\hfil\text{Conditionals
on $Y$ given $X$}\hfil\lx@intercol\vrule\lx@intercol\\\ \hline\cr
P_{{Y}|{X}}({y}\,|\,{x})&y=0&y=1&P_{X}(x)\\\ \hline\cr x=0&1&0&0.5\\\
\hline\cr x=1&0.9&0.1&0.5\\\ \hline\cr P_{Y}(y)&0.95&0.05&-\\\
\hline\cr\end{array}$ (10)

In the given scenario, we calculate the ground-truth MI by substituting the
joint and marginal PMFs from (2) in (3) as
$I(X;Y)=\sum_{x\in\\{0,1\\}}\sum_{y\in\\{0,1\\}}P_{(X,Y)}(x,y)\cdot\lg\left(\frac{P_{(X,Y)}(x,y)}{P_{X}(x)\cdot
P_{Y}(y)}\right)=$\mathbf{0.0519}\text{\,}\text{bits}$$. The Bayes predictor
producing optimal decisions is defined using (6) as $g^{bc}(x)=0$ with Bayes
error-rate of $0.05$, as
$\operatorname{m}_{\textsc{{Err}}}{(\textsc{{Bc}})}=\sum_{x\in\\{0,1\\}}P_{{Y}|{X}}({y=1}\,|\,{x})\cdot
P_{X}(x)=0.05$. Consequently, the binary entropy of Bayes error-rate and $Y$
are calculated using (1) as
$H(Y)=H_{2}(\operatorname{m}_{\textsc{{Err}}}{(\textsc{{Bc}})})\approxeq$\mathbf{0.286}\text{\,}\text{bits}$$.
Inserting them into (9), the Mid-Point MI estimation is
$\mathbf{0.0932}\text{\,}\text{bits}$ approximately.

In this example, the ground-truth MI is equal to
$\mathbf{0.0519}\text{\,}\text{bits}$, while the Mid-Point MI estimation is
$\mathbf{0.0932}\text{\,}\text{bits}$, indicating the problem of
overestimation. This arises due to the deterministic nature of error-rates in
capturing (quantifying) the probabilistic relationship (MI) between two random
variables. Furthermore, error-rates might not serve as an appropriate metric
for assessing the performance of imbalanced data sets, leading to detecting
non-existent ILs in the system and resulting in false positives. To enhance
estimation accuracy using these bounds, it becomes crucial to define a Bayes
predictor that minimizes metrics as balanced error-rate (BER) or maximizes
mathews correlation coefficient (MCC), F$1$-score, etc., as defined in Section
2.2.2. Zhao et al. (2013) derives a relationship between BER and MI for the
binary classification setting. Unfortunately, such a relationship doesn’t
extend seamlessly to handling multiple classes.

_Baseline ILDs Oversights_ One common approach to identifying IL is the PTT-
Majority baseline, which relies on uncovering significant differences in
error-rates between the Bayes predictor (or a suitable approximation $g$ or
$g_{p}$) and the marginal Bayes predictor (Gupta et al., 2022). However,
considering the example in (10), the Bayes predictor ($g^{bc}(x)=0$) coincides
with marginal Bayes predictor and always predicts class $0$ for any given
input $x$. This alignment between Bayes predictor and marginal Bayes predictor
leads to a situation where the IL present in a system goes undetected by PTT-
Majority. Furthermore, the FET based ILD approaches, including FET-Mean and
FET-Median, also fall short in identifying this IL. The failure to identify
the IL stems from the Bayes predictor’s confusion matrix, which in this case
is evaluated as $\begin{pmatrix}95&00\\\ 05&00\\\ \end{pmatrix}$, yielding an
exact $p$-value of $1$ when tested using the FET (Gupta et al., 2022).
Consequently, in cases when Bayes error-rate is very close to marginal Bayes
predictor’s error-rate, these baseline ILD approaches fail to detect the IL in
the system and produce false negatives (Gupta et al., 2022).

###  3.2 Log-Loss Estimation

To address the challenges posed by the Mid-Point approach’s overestimation,
false positives, and the false negatives produced by baseline ILD approaches,
we propose the Log-Loss approach for MI estimation. This approach utilizes the
empirical risk minimizer $g_{p}$ defined in (8) to approximate the _log-loss_
of the Bayes predictor to estimate MI. The MI between $X$ and $Y$ can be
calculated by substituting the log-loss of Bayes predictor and marginal Bayes
predictor in (2) as

$\displaystyle I(X;Y)$ $\displaystyle=-H({Y}\,|\,{X})+H(Y)$
$\displaystyle=\int\limits_{\boldsymbol{x}\in\mathcal{X}}P_{X}(X=\boldsymbol{x})\sum_{y\in\mathcal{Y}}P_{{Y}|{X}}({Y=y}\,|\,{X=\boldsymbol{x}})\cdot\lg\left(P_{{Y}|{X}}({Y=y}\,|\,{X=\boldsymbol{x}})\right)d\boldsymbol{x}$
$\displaystyle\quad-\sum_{y\in\mathcal{Y}}P_{Y}(Y=y)\cdot\lg{\left(P_{Y}(Y=y)\right)}$
$\displaystyle=-\mathbb{E}[\ell_{ll}(\text{Bayes
predictor})]+\mathbb{E}[\ell_{ll}(\text{marginal Bayes predictor})]$
$\displaystyle=-\mathbb{E}[\ell_{ll}(\textsc{{Bc}})]+\mathbb{E}[\ell_{ll}(\textsc{{Mc}})]\,,$

where $\mathbb{E}[\ell_{ll}(\cdot)]$ represent the expected log-loss. Note,
that the conditional entropy ($H({Y}\,|\,{X})$) is equal to the expected log-
loss of _Bayes predictor_ , i.e
$H({Y}\,|\,{X})=\mathbb{E}[\ell_{ll}(\textsc{{Bc}})]$), while the entropy of
$Y$ is the expected log-loss of _marginal Bayes predictor_ , i.e,
$H(Y)=\mathbb{E}[\ell_{ll}(\textsc{{Mc}})]$).

####  3.2.1 Log-Loss

As discussed in Section 2.2.2, the most common loss function, CCE, is used to
obtain a probabilistic classifier $g_{p}$ for a given classification task
($\mathcal{D}$). This classifier provides the estimated conditional _class
probabilities_ for an instance $\boldsymbol{x}$, forming a vector
$g_{p}(\boldsymbol{x})=\boldsymbol{\hat{p}}=(\hat{p}_{1},\ldots,\hat{p}_{M})$,
where $\hat{p}_{m}$ is the probability assigned to class $m$, and
$\sum_{m=1}^{M}\hat{p}_{m}=1$. (Bishop, 2006, chap. 4) defines the log-loss
for a probabilistic classifier $g_{p}$ as

$\ell_{ll}(y,\boldsymbol{\hat{p}})=-\sum_{m=1}^{M}\hat{p}_{m}\cdot\lg(\hat{p}_{m})=-\sum_{m=1}^{M}g_{p}(\boldsymbol{x})[m]\cdot\lg(g_{p}(\boldsymbol{x})[m])\,.$

The log-loss reaches its minimum value of $0$ if the probability is high,
i.e., $\hat{p}_{m}\approxeq 1\,,m\in[M]$ for one specific class $m$ and low
probabilities for all other classes, i.e., $\hat{p}_{j}\approxeq 0\,,\forall
j\in[M]\setminus\\{m\\}$. It’s akin to conditional entropy, where uncertainty
in the classification increases log-loss, i.e., if $g_{p}$ produces uniformly
distributed (conditional) class probabilities results in a higher penalty,
i.e., $\hat{p}_{m}=\tfrac{1}{M}\,,\forall
m\in[M]\implies\ell_{ll}(y,\boldsymbol{\hat{p}})=\lg(M)$ (Bishop, 2006, chap.
4). Therefore, expected log-loss of Bayes predictor and marginal Bayes
predictor serves as a direct substitution estimate for conditional entropy and
entropy of $Y$, receptively, i.e.,
$H({Y}\,|\,{X})\approxeq\ell_{ll}(y,\boldsymbol{\hat{p}})$. Note that CCE also
behaves similarly, but it is only minimized when the estimated probability of
the ground-truth class label $y$ is high $\approxeq 1$ for a given input
$\boldsymbol{x}\in\mathcal{X}$ and therefore makes it an appropriate loss
function to estimate conditional class-probabilities by learning $g_{p}$
(approximation of Bayes predictor) for a given data set $\mathcal{D}$ using
(8).

####  3.2.2 Marginal Bayes Predictor Estimation

The estimation of marginal Bayes predictor is referred to as marginal
classifier and denoted by $g_{p}^{mc}$, relies on estimating the class-
distribution from the provided data set
$\mathcal{D}=\\{(\boldsymbol{x}_{i},y_{i})\\}_{i=1}^{N}$ to make predictions
using the (8) and only the class labels $y$ of the given data set
$\mathcal{D}$. For each input $\boldsymbol{x}\in\mathcal{X}$ from the data set
$\mathcal{D}$, the estimated class probabilities
$g_{p}^{mc}(\boldsymbol{x})=\boldsymbol{\hat{p}}_{mc}=(\hat{p}_{1},\ldots,\hat{p}_{M})$
are determined, where
$\hat{p}_{m}=\tfrac{\lvert\\{(\boldsymbol{x}_{i},y_{i})\in\mathcal{D}\;\lvert\;y_{i}=m\\}\rvert}{\lvert\mathcal{D}\rvert}$
represents the estimated prior probability for class $M$. The log-loss of
marginal Bayes predictor is used to estimate the entropy of $Y$, i.e.,
$H(Y)\approxeq\ell_{ll}(y,\boldsymbol{\hat{p}}_{mc})$ The maximum value of
$\ell_{ll}(\textsc{{Mc}})$ being $\lg(M)$ for a balanced data set and any
value between $0$ and $\lg(M)$, indicates an imbalance in the data set. It
takes a minimum of $0$ for a highly imbalanced data set where one class
dominates the other classes.

####  3.2.3 MI Estimation

In the Log-Loss approach, we propose to induce the Bayes predictor using
$g_{p}$ and marginal Bayes predictor using $g_{p}^{mc}$. The expected log-loss
of Bayes predictor is approximated by evaluating the log-loss of $g_{p}$
(_empirical risk minimizer_ of (8)) for the given data set $\mathcal{D}$,
i.e., $\mathbb{E}[\ell_{ll}(\text{Bayes
predictor})\approxeq\mathbb{E}[\ell_{ll}(g_{p})]$

$\mathbb{E}[\ell_{ll}(g_{p})]=\dfrac{1}{N}\sum_{i=1}^{N}\ell_{ll}(y,g_{p}(\boldsymbol{x}_{i}))=-\dfrac{1}{N}\sum_{i=1}^{N}\sum_{m=1}^{M}g_{p}(\boldsymbol{x}_{i})[m]\cdot\lg(g_{p}(\boldsymbol{x}_{i})[m])\,.$

Similarly, the expected log-loss of marginal Bayes predictor is approximated
by evaluating the log-loss of $g_{p}^{mc}$, i.e.,
$\mathbb{E}[\ell_{ll}(\text{marginal Bayes
predictor})\approxeq\mathbb{E}[\ell_{ll}(g_{p}^{mc})]$

$\mathbb{E}[\ell_{ll}(g_{p}^{mc})]=\dfrac{1}{N}\sum_{i=1}^{N}\ell_{ll}(y,g_{p}^{mc}(\boldsymbol{x}_{i}))=-\sum_{m=1}^{M}\boldsymbol{\hat{p}}_{mc}[m]\cdot\lg(\boldsymbol{\hat{p}}_{mc}[m])\,.$

The MI is then estimated using approximations for the Log-Loss approach as

$\displaystyle\widehat{I}(X;Y)$
$\displaystyle\approxeq\mathbb{E}[\ell_{ll}(g_{p}^{mc})]-\mathbb{E}[\ell_{ll}(g_{p})]$
$\displaystyle\approxeq\sum_{m=1}^{M}\left(\dfrac{1}{N}\sum_{i=1}^{N}g_{p}(\boldsymbol{x}_{i})[m]\cdot\lg(g_{p}(\boldsymbol{x}_{i})[m])-\boldsymbol{\hat{p}}_{mc}[m]\cdot\lg(\boldsymbol{\hat{p}}_{mc}[m])\right)\,.$
(11)

_Precision_ Consider the example in (10), to illustrate the MI estimation
precision of the Log-Loss approach. While this example doesn’t conclusively
demonstrate suitability, it does highlight the potential precision of the Log-
Loss approach in estimating MI and detecting ILs. Since we propose to employ
the consistent learners using the AutoML tools to induce Bayes predictor, the
learned predictors $g_{p}$ and $g_{p}^{mc}$ closely approximate the
conditional and true marginal probabilities, respectively. Therefore, the
predicted marginal probabilities by marginal Bayes predictor ($g_{p}^{mc}$)
are approximately equal to $\boldsymbol{\hat{p}}_{mc}\approxeq(0.95,0.05)$,
(conditional) class probabilities by $g_{p}$ are approximately equal to
$g_{p}(x=0)\approxeq(P_{{Y}|{X}}({y=0}\,|\,{x=0}),P_{{Y}|{X}}({y=1}\,|\,{x=0}))=(1.0,0.0)$
for $x=0$ and
$g_{p}(x=1)\approxeq(P_{{Y}|{X}}({y=0}\,|\,{x=1}),P_{{Y}|{X}}({y=1}\,|\,{x=1}))=(0.9,0.1)$
for $x=1$. The MI estimated via (11) as $I(X;Y)\approxeq
0.286-0.234\approxeq$\mathbf{0.052}\text{\,}\text{bits}$$, will be close to
the actual MI. This underlines how the accurate approximation of learned
predictors contributes to the close estimation of MI and reinforces the Log-
Loss approach’s potential for precision in detecting ILs within a system.

####  3.2.4 Classifier Calibration

Ensuring the precision of probability estimates from the probabilistic risk
minimizer $g_{p}$ is crucial for the effectiveness of the Log-Loss approach.
However, even with the CCE loss function, these probabilities might not
accurately reflect the ground-truth conditionals $P_{{Y}|{X}}$, especially in
imbalanced class scenarios (Pakdaman Naeini et al., 2015; Domingos and
Pazzani, 1997). Famous classification methods like Naive Bayes, random forest
classifier (RF) produce inaccurate probability estimates despite their success
in accurate classification (Cheng and Hüllermeier, 2012; Domingos and Pazzani,
1997). Consequently, the utilization of classifier calibration techniques
becomes essential to produce predicted probabilities close to actual
frequencies, thereby enhancing the predictive value of probabilistic models
(Silva Filho et al., 2023). In the realm of classifier calibration, post-
processing methods are employed to map predicted class probabilities to more
reliable and calibrated probabilities, often achieved through distinct models
or transformation methods (Silva Filho et al., 2023). This study focuses on
five prevalent multi-class calibration approaches to enhance the accuracy of
MI estimation using the Log-Loss approach:

  * •

Isotonic Regression: This approach calibrates predicted probabilities by
fitting a monotonic function that transforms the original scores by minimizing
the mean squared error (MSE), ensuring improved calibration without
overfitting the training data (Barlow and Brunk, 1972; Zadrozny and Elkan,
2002). When combining this method with log-loss to estimate MI, we refer to it
IR Cal Log-Loss.

  * •

Platt’s Scaling: It involves training a logistic regression model on the
original classifier’s outputs and true labels to act as a calibrated version
of the original classifier to produce calibrated class probabilities (Platt,
1999, 2000). When combining this method with log-loss to estimate MI, we refer
to it PS Cal Log-Loss.

  * •

Beta Calibration: This approach employs a parametric transformation using a
beta distribution with adjustable hyperparameters to flexibly calibrate
probabilities for different data sets (Kull et al., 2017). When combining this
method with log-loss to estimate MI, we refer to it Beta Cal Log-Loss.

  * •

Temperature Scaling: By scaling the logits (log-odds) of predicted
probabilities using a “temperature” parameter, this method calibrates the
entire distribution while maintaining the rank order of class labels according
to original predicted class probabilities (Guo et al., 2017). When combining
this method with log-loss to estimate MI, we refer to it TS Cal Log-Loss.

  * •

Histogram Binning: Predicted probabilities are discretized into bins, and the
true class frequency in each bin is estimated to adjust predicted class
probabilities and achieve calibration (Zadrozny and Elkan, 2001). When
combining this method with log-loss to estimate MI, we refer to it HB Cal Log-
Loss.

These calibration techniques are pivotal in refining the predicted class
probabilities and ensuring accurate MI estimation through the Log-Loss
approach. When using these techniques to estimate MI and perform ILD using
Log-Loss, we refer to them collectively as Cal Log-Loss.

###  3.3 Baselines

In this section, we will briefly describe three baseline approaches, namely
Gaussian mixture model (GMM), mutual information neural estimation (MINE) and
PC-Softmax, recently introduced to estimate the MI (Eirola et al., 2014;
Belghazi et al., 2018; Qin and Kim, 2019).

####  3.3.1 Gaussian Mixture Model

GMM provide a robust solution to approximating unknown PDFs ($P_{(X,Y)}$,
$P_{X}$, $P_{Y}$) required to estimate MI via (3) effectively for
classification data sets $\mathcal{D}$ (Maia Polo and Vicente, 2022). The
process begins by fitting a GMM to the joint space
$\mathcal{X}\times\mathcal{Y}$, to estimate the joint PDF
$\hat{P}_{(X,Y)}(\cdot)$ between $X$ and $Y$. To deduce the marginal PDFs on
$X$ ($\hat{P}_{X}(\cdot)$) and $Y$ ($\hat{P}_{Y}(\cdot)$), fitted GMM is used
with appropriate variables. The MI is approximated using the estimated PDFs
($\hat{p}(\cdot)$) and the provided data set
$\mathcal{D}=\\{\boldsymbol{x}_{i},y_{i}\\}_{i=1}^{N}$ as

$\widehat{I}(X;Y)=\frac{1}{N}\sum_{i=1}^{N}\left(\lg(\hat{P}_{(X,Y)}(X=\boldsymbol{x}_{i},Y=y_{i}))-\lg(\hat{P}_{X}(\boldsymbol{x}_{i}))-\lg(\hat{P}_{Y}(y_{i}))\right)\,.$

While GMMs are robust, they face challenges of overfitting with high-
dimensional data sets, where model parameters grow quadratically with data set
dimensionality ($d$) (Maia Polo and Vicente, 2022; Eirola et al., 2014). The
high-dimensional data clustering (HDDC) variant restricts covariance matrices
(spherical, diagonal) to reduce parameters, allowing for efficient estimation
and compatibility with the expectation-maximization (EM) algorithm to address
the problem of curse of dimensionality (Eirola et al., 2014). GMM performs
poorly on high-dimensional data sets and is particularly suited for data that
follows a MVN distribution.

_Estimation Process_ Considering a GMM with $K$ components for the space
$\mathcal{X}\times\mathcal{Y}$, each component $k$ is associated with its
unique parameters: mixing coefficients ($\pi_{k}$), means ($\mu_{k}$), and
covariance matrices ($\Sigma_{k}$) (Eirola et al., 2014). These mixing
coefficients must be $0<\pi_{k}<1$ for each component, and the sum of all
coefficients over the components equals one, i.e., $\sum_{k=1}^{K}\pi_{k}=1$,
which are written as

$\mu_{k}=\begin{bmatrix}\mu_{X_{k}}\\\
\mu_{Y_{k}}\end{bmatrix},\quad\Sigma_{k}=\begin{bmatrix}\Sigma_{XX_{k}}&\Sigma_{XY_{k}}\\\
\Sigma_{YX_{k}}&\Sigma_{YY_{k}}\end{bmatrix}\,.$

The marginal model for $X$ ($\hat{P}_{X}(\cdot)$) is a GMM of $K$ components
with the same mixing coefficients ($\pi_{k}$), but only consider the means
($\mu_{X_{k}}$) and covariance matrices ($\Sigma_{XX_{k}}$) associated with
$X$ (Eirola et al., 2014). The marginals GMM for $Y$ ($\hat{P}_{Y}(\cdot)$) is
determined in a similar way. EM algorithm is employed to fit a GMM and
determine $K$ (number of components). The value of $K$ is essential since an
insufficient number of components can overlook data properties, while too many
may lead to overfitting. The Akaike information criterion (AIC) assists in
this decision, considering both the log-likelihood $\ell(\theta)$ of the
fitted model and the number of free parameters $F$, defined as
$\text{AIC}=-2\log\ell(\theta)+2F$ (Eirola et al., 2014). In order to make the
approach robust and comparison fair with our approaches, we perform the
hyperparameter optimization (HPO) using Akaike information criterion (AIC) as
the objective function.

####  3.3.2 Mutual Information Neural Estimation

MINE aims to estimate the MI between two variables ($X,\,Y$) by training a
neural network to optimize a specific objective function for a given data set
$\mathcal{D}=\\{\boldsymbol{x}_{i},y\\}_{i=1}^{N}$

$\ell_{\textsc{{MINE}}}(\boldsymbol{x},y)=\sup_{\theta}\left(\mathbb{E}_{(\boldsymbol{x},y)\sim
P_{(X,Y)}(\boldsymbol{x},y)}[T_{\theta}(\boldsymbol{x},y)]-\log(\mathbb{E}_{(\boldsymbol{x},y)\sim
P_{X}(\boldsymbol{x})\cdot
P_{Y}(y)}[\mathbf{e}^{(T_{\theta}(\boldsymbol{x},y))}])\right)\,,$

where $T_{\theta}(\boldsymbol{x},y):\mathbb{R}^{d+M}\rightarrow\mathbb{R}$
represents the neural network (learnable parameters $\theta$) function, with
input layer of size $d+M$ (units) (Belghazi et al., 2018). The $d$-dimensional
input $\boldsymbol{x}$ concatenated with one-hot-encoded $M$ length ground-
truth class label $y$ serves as an input to the network. The network’s
architecture encompasses an input layer (Linear units) and hidden layers using
ReLU activation, evaluating a real-valued output (Linear unit) quantifying the
dependency between $\boldsymbol{x}$ and $y$. MINE maximizes the discrepancy
between the network’s expectations over the joint ($P_{(X,Y)}(\cdot)$) and
marginal ($P_{X}(\cdot),P_{Y}(\cdot)$) distributions. Essentially, MINE serves
as an _upper-bound_ estimator for MI by capturing the variable’s statistical
dependency without having to explicitly determine their unobserved
(underlying) probability distributions (Belghazi et al., 2018). Through
iterative training, the process outlined in Algorithm 1, MINE efficiently
estimates MI, especially when data dimensionality is high with unobserved
PDFs.

Algorithm 1 Training algorithm: MINE

1:Dataset $\mathcal{D}=\\{(\boldsymbol{x}_{i},y_{i})\\}_{i=1}^{N}$, Neural
Network $T_{\theta}$

2:Estimated MI: $\hat{I}_{\theta}$

3:Initialize network parameters $\theta$ and hyperparameters

4:Set $N_{epochs}=10000$, Set $B$ such that $N\mod B=0$

5:for $n\in[N_{epochs}]$ do

6: Define a random permutation $\pi$ on set of indices $[N]$

7: for $b\in[B]$ do

8: Define batch start and end: $b_{s}=(b-1)\frac{N}{B}+1$,
$b_{e}=b\cdot\frac{N}{B}$

9: Evaluate the lower-bound:

10:
$V(\theta)\leftarrow\frac{1}{B}\sum_{i=b_{s}}^{b_{e}}T_{\theta}(\boldsymbol{x}_{i},y_{i})-\log\left(\frac{1}{B}\sum_{i=b_{s}}^{b_{e}}\left(\text{exp}{\left(T_{\theta}(\boldsymbol{x}_{i},y_{\pi(i)})\right)}\right)\right)$

11: Evaluate bias-corrected gradients (using exponential moving average) as:

12: $G_{b}(\theta)\leftarrow\nabla_{\theta}V(\theta)$

13: Update the statistics network parameters:
$\theta\leftarrow\theta+G_{b}(\theta)$

14: end for

15:end for

16:Compute
$\hat{I}_{\theta}=\dfrac{1}{N}\sum_{i=1}^{N}T_{\theta}(\boldsymbol{x}_{i},y_{i})-\log\left(\dfrac{1}{N}\sum_{i=1}^{N}\left(\text{exp}{\left(T_{\theta}(\boldsymbol{x}_{i},y_{\pi(i)})\right)}\right)\right)$

17:return Estimated MI: $\lg(\mathrm{e})\cdot\hat{I}_{\theta}$

_Overfitting Solution_ The main challenge for MINE lies in determining the
suitable neural network architecture for estimating MI. The HPO approach is
not straightforward due to the absence of a suitable objective function.
Discrepancy-based methods like MINE and KL divergence, which compare network
outputs between shuffled and original data sets, can potentially lead to
architectures that overfit and produce biased MI estimates. To tackle these
challenges, we propose an innovative solution that combines HPO with ensemble
methods. Instead of relying on a single architecture, we leverage the power of
diversity by aggregating the outputs of top-$5$ performing architectures
identified through HPO, with MSE as the objective function. This enhances MI
estimation precision and yields more reliable and robust MI estimates.
Moreover, it demonstrates strong generalization across diverse data sets,
making it a promising solution.

####  3.3.3 PC-Softmax

PC-Softmax offers an alternative approach to estimating MI in classification
data sets $\mathcal{D}$ using a modified version of _softmax_ function,
referred to as probabilitiy-corrected softmax (PC-softmax) (Qin and Kim,
2019). This approach combines PC-softmax with the CCE loss function to acquire
an empirical risk minimizer $g_{p}$ for classification data sets
$\mathcal{D}$, using (8). The (conditional) class probabilities predicted by
$g_{p}$ are then used to estimate MI using the PC-softmax function (Qin and
Kim, 2019).

This approach provides several advantages over MINE — it’s straightforward and
easy to implement, and the network architecture can be easily optimized using
classification evaluation metrics or loss functions as objective functions. It
approximates the lower bound of KL divergence between the joint
($P_{(X,Y)}(\cdot)$) and the product of marginals
($P_{X}(\cdot),P_{Y}(\cdot)$) for variables $X$ and $Y$ to estimate a lower
bound on the MI (Qin and Kim, 2019).

_Estimation Process_ The accurate estimation of (conditional) class
probabilities by the empirical risk minimizer $g_{p}$ is quite challenging,
particularly compared to standard classification. In practice, typically
$g_{p}$ is obtained by learning a scoring function
$h_{{s}|{\boldsymbol{\theta}}}(\boldsymbol{x})\colon\mathcal{X}\rightarrow\mathbb{R}^{M},\,h_{{s}|{\boldsymbol{\theta}}}\in\mathcal{H}_{s}$,
which assigns an un-normalized real-valued score to each class label,
parameterized by $\boldsymbol{\theta}\in\mathbb{R}^{n}$, where $n$ is the
number of trainable parameters (Bishop, 2006, chap. 4). Neural networks or
multi-layer perceptron (MLP) with $M$ linear nodes in the second last layer
are commonly used to estimate these parameters $\boldsymbol{\theta}$ of the
scoring function $h_{{s}|{\boldsymbol{\theta}}}(\cdot)$. MLP outputs un-
normalized scores
$h_{{s}|{\boldsymbol{\theta}}}(\boldsymbol{x})=\boldsymbol{s}=(s_{1},\ldots,s_{M})$
for a given instance $\boldsymbol{x}$, where
$s_{m}=h_{{s}|{\boldsymbol{\theta}}}(\boldsymbol{x})[m]=\boldsymbol{s}[m]$ is
the score assigned to class $m$. The _softmax_ function
$S(\boldsymbol{x})\colon\mathbb{R}^{M}\rightarrow[0,1]^{M}$ is then commonly
used to convert these scores into probabilities
$\boldsymbol{\hat{p}}=(\hat{p}_{1},\ldots,\hat{p}_{M})$, typically implemented
in the output layer with $M$ nodes of the MLP. For a given ground-truth class
$y\in[M]$, the score $\boldsymbol{s}[y]$ is transformed to the corresponding
(conditional) class probability $\boldsymbol{\hat{p}}[y]$ by the _softmax_
function

$S(\boldsymbol{x})[y]=\boldsymbol{\hat{p}}[y]=\frac{\text{exp}{\left(\boldsymbol{s}[y]\right)}}{\sum_{m=1}^{M}\text{exp}{\left(\boldsymbol{s}[m]\right)}}\,.$

The empirical risk minimizer $g_{p}$ is then obtained using the CCE loss,
computed for real-valued scores using the _softmax_ function, defined as
$\ell_{CCE}(\boldsymbol{x},y)\coloneq-\ln(S(\boldsymbol{x})[y])$.

Qin and Kim (2019) shows that the expected CCE is equivalent to the MI between
input $\mathcal{X}$ and output $Y$ random variable up to a constant $\lg(M)$
under a uniform label distribution, i.e.,
$I(X;Y)\geq\widehat{I}(X;Y)\approxeq-\sum_{(\boldsymbol{x},y)\in\mathcal{D}}\log(S(\boldsymbol{x})[y])$.
The authors introduce PC-softmax as an extension of traditional _softmax_ to
enhance MI estimation, providing mathematical justifications for its
effectiveness as an estimator

$S_{pc}(\boldsymbol{x})[y]=\frac{\text{exp}{\left(\boldsymbol{s}[y]\right)}}{\sum_{m=1}^{M}\text{exp}{\left(P_{Y}(m)\cdot\boldsymbol{s}[m]\right)}}\,.$

The PC-softmax function reduces to the traditional _softmax_ when
$P_{Y}(y)=\frac{1}{M},\,\forall y\in[M]$. The corresponding CCE loss is
defined as
$\ell_{\textsc{{Pcce}}}(\boldsymbol{x},y)\coloneq-\ln(S_{pc}(\boldsymbol{x})[y])$,
which is especially effective for enhancing MI estimation accuracy in
unbalanced data sets. The loss is evaluated by estimating the marginal using
the data set, such that
$\hat{P}_{Y}(y=m)=\tfrac{\lvert\\{(\boldsymbol{x}_{i},y_{i})\in\mathcal{D}\;\lvert\;y_{i}=m\\}\rvert}{\lvert\mathcal{D}\rvert}$,
since $P_{Y}(\cdot)$ is seldom observed. This loss function maximizes the
lower bound on the MI, making it a simple yet effective technique for
estimating MI in classification tasks. The MI for a given data set
$\mathcal{D}=\\{\boldsymbol{x}_{i},y\\}_{i=1}^{N}$ is estimated as

$I(X;Y)\geq\widehat{I}(X;Y)\approxeq\frac{1}{N}\sum_{(\boldsymbol{x},y)\in\mathcal{D}}\ell_{\textsc{{Pcce}}}(\boldsymbol{x},y)\times\lg(\mathrm{e})=-\frac{1}{N}\sum_{(\boldsymbol{x},y)\in\mathcal{D}}\lg(S_{pc}(\boldsymbol{x})[y])\,.$

This approach provides a _lower bound_ on MI, similar to Log-Loss, but Log-
Loss offers a simpler and potentially more accurate method for estimating MI
in classification tasks. Contrary to MINE, the MLP architecture for this
approach can be optimized for any given data set $\mathcal{D}$ using standard
objective functions used for optimizing MLP architectures in classification
tasks, like accuracy, validation CCE loss, etc.

##  4 AutoML tools

In the last decade, the field of AutoML has emerged as a response to the high
and unmet demand for engineering machine learning applications experts. While
AutoML envisions the automation of the entire data science process, the
arguably most studied problem is that of the combined algorithm selection and
hyperparameter optimization (CASH), which was first formally specified by
Thornton et al. (2013).

Since then, various systems have been devised (Thornton et al., 2013; Feurer
et al., 2019; Olson et al., 2016; Mohr et al., 2018), demonstrating promising
performances for tailoring the choice of ML algorithms and the setting of
their hyperparameters to a given task, typically comprising a data set and a
loss function. A more thorough overview of AutoML methods can be found in
(Zöller and Huber, 2021). While various AutoML systems with complementary
strengths have been proposed in the literature, a recent benchmark study
(Gijsbers et al., 2022) suggests AutoGluon (Erickson et al., 2020, 2022) as
the AutoML system with the best performance across data sets and different
tasks.

In (Hollmann et al., 2023), instead of tackling the CASH problem, a general
predictor called TabPFN is fitted across various different data sets and is
able to immediately return highly accurate predictions. Yielding competitive
performance to AutoGluon, TabPFN suggests itself as a state-of-the-art AutoML
tool that returns results substantially faster than other AutoML systems.

Due to their strong performance, we selected AutoGluon and TabPFN to
approximate the Bayes predictor in our method. In Sections 4.1 and 4.2, we
briefly describe the concepts AutoGluon and TabPFN are based on, respectively.

###  4.1 AutoGluon

In contrast to other AutoML systems that aim to select feature preprocessing
algorithms and a learning algorithm and tune their respective hyperparameters,
AutoGluon (Erickson et al., 2020, 2022) follows a different approach. More
specifically, AutoGluon focuses on building a stacking ensemble, greedily
adding more models as long as they are beneficial for the overall performance
of the ensemble. To this end, AutoGluon searches over a wide range of
different ML algorithms and pre-fitted models, including deep learning,
gradient boosting, and linear models. Simultaneously, AutoGluon tunes
hyperparameters for each model to optimize their performance. For performing
HPO, different standard techniques can be used such as ASHA (Li et al., 2020),
Hyperband (HB) (Li et al., 2017), Bayesian optimization (BO) (Frazier, 2012),
or Bayesian optimization and Hyperband (BOHB) (Falkner et al., 2018).

A recent benchmark study found AutoGluon to perform superior to other AutoML
systems in terms of predictive performance and robustness. The only drawback
noted by Gijsbers et al. (2022) is the prediction time, as it will grow with
the number of models included in the returned stacking ensemble.

Moreover, AutoGluon allows the customization of the set of considered ML
algorithms and models, including gradient boosting machine (GBM) tree-based
models and deep neural networks. For our study, we limited the search space of
AutoGluon to consistent classifiers (Mielniczuk and Tyrcha, 1993; Biau et al.,
2008). In Table 4, we present the models and learning algorithms with
corresponding hyperparameters, including their range values available to
AutoGluon in our experimental evaluation in Sections 5 and 7.

###  4.2 TabPFN

TabPFN (Hollmann et al., 2023) is an approach to AutoML that solves
classification tasks in a transductive way. Instead of fitting a model
specifically to a data set and returning this model to the user, TabPFN is a
transformer-based neural network architecture trained on vast number of
different probability distributions of tabular data. Instead of searching for
a proper learning algorithm and its hyperparameter setting in the first place,
TabPFN can be used to instantaneously make predictions as it is already fitted
on prior data.

Due to this, predictions can be obtained faster by orders of magnitude
compared to classical AutoML estimation baseline tools. Moreover, the
predictions have proven to be highly accurate and excel, particularly on small
tabular data sets. However, TabPFN was only considered for data sets
consisting of up to $1000$ training data points, $100$ numeric features, and
$10$ classes. Hence, in practical applications with more data points, TabPFN
is currently not reliably applicable.

Due to its highly accurate predictions, which furthermore appear to come along
with well-calibrated probability estimates for the classes, TabPFN is an
interesting pick for our approach to approximate the Bayes predictor. Since
TabPFN is limited to handling data sets with up to $100$ numeric features, our
real IL-Datasets obtained from OpenSSL TLS server implementations contain more
than $100$ features. To address this, we recommend applying dimensionality
reduction techniques, specifically RF and extra trees classifier (XT) (only
when $d>100$). We perform fine-tuning on TabPFN, adjusting parameters such as
the number of reduced features, the reduction model, and the number of prior-
fitted models using HPO. Table 4 lists the specific parameter ranges used for
conducting our experiments detailed in Sections 5 and 7.

##  5 Mutual Information Estimation: Empirical Evaluation

This section extensively evaluates our methods and compares them to baseline
approaches for estimating MI using data sets generated by synthetic systems,
subsequently referred to as _synthetic data sets_. We begin by outlining the
procedure for generating these synthetic data sets via MVN distribution, as
detailed in Section 5.1. Section 5.2 provides a comprehensive overview of the
empirical setup, utilizing the NMAE metric for performance assessment. To
summarize the results, we assess the overall performance of MI estimation
techniques across various systems in Section 5.3.1. Our evaluation also
explores generalization capabilities across factors, including the number of
classes ($M$), input dimensions ($d$), class imbalance ($r$), and noise level
($\epsilon$) present within the given synthetic data set. Concise results on
generalization capabilities for selected approaches are discussed in Section
5.3.2, with analysis detailed in Section C.2. Our results highlight our
proposed approaches’ superior performance and better generalization
capabilities compared to state-of-the-art methods.

###  5.1 Simulating Synthetic Systems

We employ the MVN distribution to simulate a real-world scenario involving
systems containing ILs and generating synthetic classification data set. MVN
is a natural choice due to its ability to mimic real-world data patterns
(Amancio et al., 2014; Bishop, 2006, chap. 2). Some illustrations of two-
dimensional data points generated using the two noise introduction procedures
are shown in Figure 3.

####  5.1.1 Introducing Class Imbalance ($r$)

In order to generate the synthetic data sets using the MVN distribution, we
define the vector $\boldsymbol{n}_{M}=(n_{1},\ldots,n_{M})$, where
$n_{m}=\lvert\\{(\boldsymbol{x}_{i},y_{i})\in\mathcal{D}\;\lvert\;y_{i}=m\\}\rvert$
denotes the sample counts for each class label $m$ in the data set
$\mathcal{D}$.

Figure 2: Class frequencies in the imbalanced synthetic data sets

The parameter $r$ indicates the minimum proportion of instances associated
with a particular class label in $\mathcal{D}$, offering insights into
frequency disparities in the class distribution of the data set.
Mathematically, it is defined as:
$r=\min_{m\in[M]}\frac{n_{m}}{\lvert\mathcal{D}\rvert}$, i.e.,
$r\in(0,\tfrac{1}{M}]$. So, when the underlying data set is balanced,
$r=\tfrac{1}{M}$, whereas in the case of imbalanced one, $r<\tfrac{1}{M}$.

In order to generate a multi-class imbalanced data set ($M>2$), we introduce
two generation methods termed _Minority_ and _Majority_. With the _Minority_
approach, we create data sets with a singular minority class, and the other
classes have an equal but larger number of data points ($>\frac{N}{M}$) in
comparison to the minority class ($\mathsf{g}_{r}=\text{Minority}$).
Conversely, with the _Majority_ method, we create data sets with just one
majority class, while the other classes have an equal but smaller number of
data points ($<\frac{N}{M}$) relative to the majority class
($\mathsf{g}_{r}=\text{Majority}$). The resulting sample counts vector
$\boldsymbol{n}_{M}^{\mathsf{g}_{r}}$ for both generation methods
$\mathsf{g}_{r}\in\\{\text{Minority},\text{Majority}\\}$ is defined as

$\displaystyle\boldsymbol{n}_{M}^{\text{Majority}}$ $\displaystyle=(\lceil
N\cdot r\rceil,\ldots,\lceil N\cdot r\rceil,N-(M-1)\lceil N\cdot r\rceil)$
$\displaystyle\boldsymbol{n}_{M}^{\text{Minority}}$
$\displaystyle=(\lceil\frac{N(1-r)}{(M-1)}\rceil,\ldots,\lceil\frac{N(1-r)}{(M-1)}\rceil,N-(M-1)\cdot\lceil\frac{N(1-r)}{(M-1)}\rceil)\,,$
(12)

where $\lvert\mathcal{D}\rvert=N=1000\cdot M$ is the number of instances or
the data set size, $M$ is the number of classes and $r$ is the class imbalance
parameter. The class frequencies for each generation method are depicted in
Figure 2.

####  5.1.2 Data Generation using MVN Distribution

We use the MVN distribution to generate the realistic synthetic classification
data sets, widely employed for benchmarking classifiers (Amancio et al.,
2014). The MVN distribution is a natural choice due to its simplicity, well-
established statistical properties, and ability to effectively model complex
inter-variable correlations in real-world data (Bishop, 2006, chap. 2).
Moreover, the central limit theorem provides additional justification for
preferring the MVN distribution over other multi-variate distributions, as it
states that the cumulative sum of a large number of independent random
variables often conforms to a normal distribution (Bishop, 2006, chap. 4).

_Formal Model_ In order to generate the synthetic system’s data set while
maintaining a clear means of obtaining the ground-truth MI, we define the
conditional probability on $X$ given $Y$ as

$P_{{X}|{Y}}({X=\boldsymbol{x}}\,|\,{Y=m})=\textsc{MVN}(\boldsymbol{\mu}_{m},\Sigma)\,,$
(13)

where $\boldsymbol{\mu}_{m}$ denotes the $d$-dimensional mean vector
associated with class label $m$ and $\Sigma$ is the covariance matrix of size
$d\times d$, expressing correlations among input features $\boldsymbol{x}$.
The marginal on $Y$ is estimated using the $\boldsymbol{n}_{M}$ as
$P_{Y}(Y=m)=\frac{n_{m}}{M}$, which boils down to $P_{Y}(Y=m)=\tfrac{1}{M}$
for balanced data sets, for which the entropy of $Y$ is maximal, i.e.
$H(Y)=\lg(M)$. The marginal on $Y$ for data sets generated using either the
“Minority” or “Majority” generation method are

$\mathsf{g}_{r}=\begin{cases}\text{Majority}&\begin{cases}P_{Y}(Y=m)=r,&\text{if
}m\in[M]\setminus M\\\ P_{Y}(Y=m)=1-r\cdot(M-1),&\text{if }m=M\end{cases}\\\
\hline\cr\text{Minority}&\begin{cases}P_{Y}(Y=m)=\frac{1-r}{M-1},&\text{if
}m\in[M]\setminus M\\\ P_{Y}(Y=m)=r,&\text{if }m=M\end{cases}\end{cases}$ (14)

The joint distribution between $X$ and $Y$ is
$P_{(X,Y)}(\boldsymbol{x},m)=P_{Y}(m)\times
P_{{X}|{Y}}({\boldsymbol{x}}\,|\,{m})$ and the marginal on $X$ is calculated
as

$\displaystyle P_{X}(\boldsymbol{x})$
$\displaystyle=\sum_{m=1}^{M}P_{Y}(m)\cdot
P_{{X}|{Y}}({X=\boldsymbol{x}}\,|\,{Y=m})=\sum_{m=1}^{M}P_{Y}(m)\cdot\textsc{MVN}(\mu_{m},\Sigma)\,.$
(15)

Using this, we can compute the exact class probabilities for a given instance
$\boldsymbol{x}$, with the conditional on $Y$ given $X$ is calculated as

$\displaystyle P_{{Y}|{X}}({Y=m}\,|\,{X=\boldsymbol{x}})$
$\displaystyle=\frac{P_{Y}(Y=m)\cdot
P_{{X}|{Y}}({X=\boldsymbol{x}}\,|\,{Y=m})}{\sum_{m=1}^{M}P_{Y}(Y=m)\cdot
P_{{X}|{Y}}({X=\boldsymbol{x}}\,|\,{Y=m})}$
$\displaystyle=\frac{P_{(X,Y)}(X=\boldsymbol{x},Y=m)}{P_{X}(X=\boldsymbol{x})}\,.$
(16)

####  5.1.3 Introducing Noise ($\epsilon$) in the System

To simulate systems generating classification data sets, we sample the
instances $(\boldsymbol{x}_{i},y_{i})\sim
P_{{Y}|{X}}({\boldsymbol{x}_{i}}\,|\,{y_{i}}),\,\forall i\in[N]$ (using
Equation 13) and the class-distribution follows $P_{Y}(y)$ (using Equation
14). However, we need to introduce noise into the data set to simulate real-
world IL scenarios in cryptographic systems. Specifically, we aim to replicate
situations where the certainty about the true output $y_{i}\in\mathcal{Y}$
diminishes as we observe more inputs $\boldsymbol{x}_{i}\in\mathcal{X}$, i.e,
$H({Y}\,|\,{X})>0,\,I(X;Y)<H(Y)$.

We propose two methods: the MVN perturbation and proximity techniques. The
perturbation technique introduces noise by flipping a certain percentage of
outputs or class labels ($y_{i},\,\forall i\in[N]$) in the generated synthetic
data set. On the other hand, the proximity technique introduces noise by
reducing the distance between the mean vectors ($\boldsymbol{\mu}_{m}$)
associated with each class. As a result, this leads to a higher degree of
overlap between the Gaussians $\textsc{MVN}(\boldsymbol{\mu}_{m},\Sigma)$
across all classes $m$. These methods enable us to simulate scenarios where
cryptographic systems leak no information (non-vulnerable), i.e., the case
when $H({Y}\,|\,{X})=H(Y)$ and $I(X;Y)=0$. Figure 3 illustrates an example of
two-dimensional data points generated by both vulnerable systems with noise
levels of $0$ and $0.5$ and non-vulnerable systems with noise level of $1.0$,
simulated using the outlined techniques to introduce noise.

Figure 3: Generated synthetic data sets using MVN perturbation and proximity
techniques

5.1.3.1 MVN Perturbation Technique Our initial approach involves introducing
noise by flipping a certain percentage of class labels in the simulated
system’s data set. The procedure for simulating perturbed systems to generate
classification data sets, for a unique configuration of the generation method
($\mathsf{g}_{r}$), number of classes ($M$), input dimensions ($d$), noise
($\epsilon$), and class imbalance ($r$), is outlined in Algorithm 2.

_Modified PDFs_: The introduction of noise by flipping $\epsilon$ percentage
of class labels ($y\sim\mathcal{Y}$) modifies the conditional probability
measures $P_{{Y}|{X}}(\cdot)$ and $P_{{X}|{Y}}(\cdot)$ defined in Section
5.1.2. For fixed $\epsilon$, let random variables
$B\sim\text{Bernoulli}(p=\epsilon,q=1-\epsilon)$,
$\tilde{Y}\sim\text{Categorical}(P_{Y}(1),\ldots,P_{Y}(M))$ such that
$\\{B,\tilde{Y}\\}$ are independent of $\\{X,Y\\}$. The modified output random
variable $Y_{\epsilon}$ is defined as

$Y_{\epsilon}=Y\cdot\llbracket B=0\rrbracket+\llbracket
B=1\rrbracket\cdot\tilde{Y}\,.$

The conditional on the modified output variable $Y_{\epsilon}$ given $X$
($P_{{Y_{\epsilon}}|{X}}(\cdot)$) is derived as

$\displaystyle
P_{{Y_{\epsilon}}|{X}}({Y_{\epsilon}=m}\,|\,{X=\boldsymbol{x}})$
$\displaystyle=P_{B}(B=0)\cdot
P_{{Y}|{X}}({m}\,|\,{\boldsymbol{x}})+P_{B}(B=1)\cdot
P_{\tilde{Y}}(\tilde{Y}=m)$ $\displaystyle=(1-\epsilon)\cdot
P_{{Y}|{X}}({Y=m}\,|\,{X=\boldsymbol{x}})+\epsilon\cdot P_{Y}(Y=m)\,.$

The marginal distributions on $X$ remain unchanged because only the class
labels are flipped. Accordingly, the modified marginal $Y_{\epsilon}$ is
derived as: $P_{Y_{\epsilon}}(Y_{\epsilon}=y)=\epsilon\cdot
P_{Y}(Y=y)+(1-\epsilon)\cdot P_{Y}(Y=y)=P_{Y}(Y=y)$. So, the marginal on
$Y_{\epsilon}$ is the same as that of $Y$. Consequently, we can evaluate the
ground-truth MI using the marginal $P_{Y}(Y=m)$ and modified conditional PDFs,
i.e., $P_{{Y_{\epsilon}}|{X}}({Y_{\epsilon}=m}\,|\,{X=\boldsymbol{x}})$ using
(17). The conditional on $X$ given $Y_{\epsilon}$
($P_{{X}|{Y_{\epsilon}}}(\cdot)$) is derived as

$P_{{X}|{Y_{\epsilon}}}({X=\boldsymbol{x}}\,|\,{Y_{\epsilon}=m})=\frac{P_{{Y_{\epsilon}}|{X}}({Y_{\epsilon}=m}\,|\,{X=\boldsymbol{x}})\cdot
P_{X}(X=\boldsymbol{x})}{P_{Y}(Y=m)}\,.$

The altered conditional distribution
$P_{{X}|{Y_{\epsilon}}}({X=\boldsymbol{x}}\,|\,{Y_{\epsilon}=m})$ appears to
be a combination of multiple MVN distributions, which introduces complexity
into the data set. This increased complexity poses a greater challenge for
models to learn from the data, making the estimation of MI more challenging.

Algorithm 2 Generate Perturbation MVN $\mathcal{D}$ for given
$\mathsf{g}_{r}$, $M$, $d$, $\epsilon$ and $r$

1:Define $\mathcal{D}=\\{\\},\,N=1000\cdot M$

2:Calculate $\boldsymbol{n}_{M}^{\mathsf{g}_{r}}=(n_{1},\ldots,n_{M})$ using
(12)

3:Sample positive semi-definite matrix for the MVN distribution
$\displaystyle\boldsymbol{Q}$
$\displaystyle\sim\text{OrthogonalMatrix}_{d\times
d}=\begin{pmatrix}q_{1,1}&\ldots&q_{1,d}\\\ \ldots&\ldots&\ldots\\\
q_{d,1}&\ldots&q_{d,d}\\\
\end{pmatrix},\,q_{i,j}\in[-1,1),\,\forall(i,j)\in[d]\times[d]$
$\displaystyle\boldsymbol{S}$
$\displaystyle\sim\text{Diagonal}(\text{RandomMatrix}_{d\times
d})=\begin{pmatrix}s_{1,1}&\ldots&0\\\ 0&\ldots&0\\\ 0&\ldots&s_{d,d}\\\
\end{pmatrix},\,s_{i,i}\in[0,1),\,\forall i\in[d]$ $\displaystyle\Sigma$
$\displaystyle=(\boldsymbol{Q}\cdot\boldsymbol{S})\cdot\boldsymbol{Q}^{T}$

4:for $m\in[M]$ do

5: Define $\boldsymbol{\mu_{m}}=(1.5m,\ldots,1.5m)\in\mathbb{R}^{d}$

6: for $i\in[n_{m}]$ do

7: Sample $\boldsymbol{x}_{i}\sim\text{MVN}(\boldsymbol{\mu_{m}},\Sigma)$

8: Assign label $y_{i}$ using $B\sim\text{Bernoulli}(p=\epsilon,q=1-\epsilon)$
$\triangleright$ Flip $\epsilon$ $\text{\,}\mathrm{\char 37\relax}$ labels

9: $\begin{cases}y_{i}=m,\quad B=0\\\
y_{i}\sim\text{Categorical}(P_{Y}(1),\ldots,P_{Y}(M)),\quad B=1\end{cases}$

10:

11: $\mathcal{D}=\mathcal{D}\cup\\{(\boldsymbol{x}_{i},y_{i})\\}$
$\triangleright$ Add instances

12: end for

13:end for

14:return $\mathcal{D}$

Algorithm 3 Generate Proximity MVN $\mathcal{D}$ for given $\mathsf{g}_{r}$,
$M$, $d$, $\epsilon$ and $r$

1:Define $\mathcal{D}=\\{\\},\,N=1000\cdot M$

2:Calculate $\boldsymbol{n}_{M}^{\mathsf{g}_{r}}=(n_{1},\ldots,n_{M})$ using
(12)

3:Sample positive semi-definite matrix for the MVN distribution
$\displaystyle\boldsymbol{Q}$
$\displaystyle\sim\text{OrthogonalMatrix}_{d\times
d}=\begin{pmatrix}q_{1,1}&\ldots&q_{1,d}\\\ \ldots&\ldots&\ldots\\\
q_{d,1}&\ldots&q_{d,d}\\\
\end{pmatrix},\,q_{i,j}\in[-1,1),\,\forall(i,j)\in[d]\times[d]$
$\displaystyle\boldsymbol{S}$
$\displaystyle\sim\text{Diagonal}(\text{RandomMatrix}_{d\times
d})=\begin{pmatrix}s_{1,1}&\ldots&0\\\ 0&\ldots&0\\\ 0&\ldots&s_{d,d}\\\
\end{pmatrix},\,s_{i,i}\in[0,1),\,\forall i\in[d]$ $\displaystyle\Sigma$
$\displaystyle=(\boldsymbol{Q}\cdot\boldsymbol{S})\cdot\boldsymbol{Q}^{T}$

4:for $m\in[M]$ do

5: Define $\boldsymbol{\mu_{m}}=(1.5\cdot m(1-\epsilon),\ldots,1.5\cdot
m(1-\epsilon))\in\mathbb{R}^{d}$ $\triangleright$ Distance $1.5(1-\epsilon)$

6: for $i\in[n_{m}]$ do

7: Sample $\boldsymbol{x}_{i}\sim\text{MVN}(\boldsymbol{\mu_{m}},\Sigma)$

8: $\mathcal{D}=\mathcal{D}\cup\\{(\boldsymbol{x}_{i},m)\\}$ $\triangleright$
Add instances

9: end for

10:end for

11:return $\mathcal{D}$

5.1.3.2 MVN Proximity Technique Our second approach introduces noise in the
simulated systems by reducing the distance between the Gaussians generated by
MVN distributions. This is achieved by moving the mean vectors
$\boldsymbol{\mu}^{\prime}_{m},\forall m\in[M]_{0}$ corresponding to MVN
distribution representing each class closer to each other. The updated MVNs
for each class $m$ is defined by
$\textsc{MVN}(\boldsymbol{\mu}^{\prime}_{m},\Sigma)$, which in turn represents
the conditional PDF $P_{{X}|{Y}}(\cdot)$ for the underlying generated system
data set. The process for simulating systems using the proximity technique to
generate classification data sets, for a unique configuration of the
generation method ($\mathsf{g}_{r}$), number of classes ($M$), input
dimensions ($d$), noise ($\epsilon$), and class imbalance ($r$), is outlined
in Algorithm 3. Notably, when introducing proximity, only the original MVN is
modified, which only modifies the conditional on $X$ given $Y$, i.e.,
$P_{{X}|{Y}}({X=\boldsymbol{x}}\,|\,{Y=m})=\textsc{MVN}(\boldsymbol{\mu}^{\prime}_{m},\Sigma)$.
Consequently, the conditional probability $P_{{Y}|{X}}(\cdot)$ and the
marginal on $X$ can be computed using Equations 16 and 15. Importantly, this
process does not impact the marginal distributions of $Y$ (as shown in
Equation 14), allowing for the evaluation of ground-truth MI using (17).

####  5.1.4 Ground-truth MI

It is imperative to have the ground truths or true MI values to evaluate and
compare MI estimation approaches. The synthetic classification data sets are
generated using the conditional $P_{{X}|{Y}}(\cdot)$ with the class
distribution (marginal on $Y$) defined by $P_{Y}(\cdot)$ in (13) and (14),
respectively. The maximum possible MI, i.e., the entropy of $Y$ ($H(Y)$), is
evaluated by substituting $P_{Y}(\cdot)$ into (1). If the underlying data set
is imbalanced, then $r\in(0,\tfrac{1}{M})$ indicating that the entropy of $Y$
is less than $\lg(M)$. Conversely, when the underlying data set is balanced
with $r=\tfrac{1}{M}$, the entropy of $Y$ equals $\lg(M)$. With this
information, we derive the conditional $P_{{Y}|{X}}(\cdot)$ and the marginal
on $Y$ ($P_{Y}(\cdot)$), as derived in (16) and (15), respectively. Thus, we
can estimate the ground-truth MI by substituting them in (2) as

$\displaystyle I(X;Y)$
$\displaystyle=\int\limits_{\boldsymbol{x}\in\mathcal{X}}P_{X}(\boldsymbol{x})\sum_{y\in\mathcal{Y}}P_{{Y}|{X}}({Y=y}\,|\,{X=\boldsymbol{x}})\cdot\lg\left(P_{{Y}|{X}}({Y=y}\,|\,{X=\boldsymbol{x}})\right)d\boldsymbol{x}$
$\displaystyle\quad-\sum_{y\in\mathcal{Y}}P_{Y}(Y=y)\cdot\lg\left(P_{Y}(Y=y)\right)$
$\displaystyle=\sum_{\boldsymbol{x}\in\mathcal{X}}\sum_{y\in\mathcal{Y}}P_{{Y}|{X}}({Y=y}\,|\,{X=\boldsymbol{x}})\cdot\lg\left(P_{{Y}|{X}}({Y=y}\,|\,{X=\boldsymbol{x}})\right)d\boldsymbol{x}$
$\displaystyle\quad-\sum_{y\in\mathcal{Y}}P_{Y}(Y=y)\cdot\lg\left(P_{Y}(Y=y)\right)\,.$

Using the above formula, we can get an estimate on the actual ground-truth MI
for the generated synthetic data set
$\mathcal{D}=\\{(\boldsymbol{x}_{i},y_{i})\\}_{i=1}^{N}$ as

$\displaystyle GI(\mathcal{D})$
$\displaystyle\approxeq\frac{1}{N}\sum_{i=1}^{N}P_{{Y}|{X}}({Y=y_{i}}\,|\,{X=\boldsymbol{x}_{i}})\cdot\lg\left(P_{{Y}|{X}}({Y=y_{i}}\,|\,{X=\boldsymbol{x}_{i}})\right)$
$\displaystyle\quad-\sum_{m=1}^{M}P_{Y}(Y=m)\cdot\lg\left(P_{Y}(Y=m)\right)\,.$
(17)

###  5.2 Experimental Setup

In this section, we outline the empirical evaluation process for our proposed
MI estimation approaches compared to the baselines on various synthetic data
sets, as illustrated in Figure 4.

Figure 4: Experimental setup for evaluating MI estimation approaches

Our primary goal is to assess the _generalization_ capabilities of our
approaches compared to the baselines. The evaluation criteria encompass
various factors, including the number of classes ($M$), input dimensions
($d$), class imbalance ($r$), and noise level ($\epsilon$) within synthetic
data sets with known MI. To achieve this, we simulate the systems generating
data set using perturbation and proximity techniques for different unique
configurations of $\mathsf{g}_{r},d,M,r,\epsilon$, as outlined in Table 1.

####  5.2.1 Evaluation Process

To conduct our experiments, we employ each MI estimation approach (detailed in
Section 3) on synthetic data sets outlined in Table 1. For each unique
configuration ($\mathsf{g}_{r},M,d,r,\epsilon$), we generate $10$ synthetic
data sets denoted as $\mathcal{D}_{j}$ using different seeds $j\in[10]$. We
apply each approach to these $10$ data sets and evaluate their performance
using NMAE, as shown in Figure 4.

For accurate MI estimation on each data set $\mathcal{D}_{j}$, we utilize
nested cross-validation with HPO. The objective functions employed for HPO
include BER for PC-Softmax, AutoGluon, and TabPFN, AIC for GMM, and MSE for
MINE with specific parameter ranges provided in Table 4. We split
$\mathcal{D}_{j}$ into a training data set with $70\text{\,}\mathrm{\char
37\relax}$ of instances and a test data set with the remaining
$30\text{\,}\mathrm{\char 37\relax}$. We then perform HPO with $100$ function
evaluations using Monte Carlo cross-validation (MCCV) employing $3$ splits on
the training data set, reserving $30\text{\,}\mathrm{\char 37\relax}$ of the
instances for validation. This process is denoted as “N$1$S-$3$SCV” indicating
$1$ split for test evaluation and $3$ MCCV splits for HPO using nested cross-
validation (denoted by N), as depicted in Figure 4. For AutoGluon, we let it
run for $1800$ seconds with BER as an objective using the training data set.
The best-performing model or pipeline is chosen based on the objective
function with minimum validation loss (or maximize validation accuracy).

We use the best model or pipeline to estimate MI ($\hat{I}_{j}$) using the
complete data set $\mathcal{D}_{j}$ and compare it to the corresponding
ground-truth MI ($I_{j}$) evaluated using (17). We evaluate our approach using
the NMAE metric, comparing the actual MIs ($\boldsymbol{I}$) with the
estimated MIs ($\hat{\boldsymbol{I}}$), defined in (18).

Table 1: Overview of the synthetic data sets for MI estimation experiments

MVN Perturbation & MVN Proximity Synthetic Data Sets Data set Type Generation
Method ($\mathsf{g}_{r}$) Input Dimensions ($d$) Classes ($M$) Noise
Level/Flip Percentage ($\epsilon$) Class Imbalance ($r$) Balanced NA {2, 4, …,
20} {2, 4, …, 10} {0.0, 0.1, …, 1.0} NA Binary-class Imbalanced Minority 5 2
{0.0, 0.1, …, 1.0} {0.05, 0.1, …, 0.5} Multi-class Imbalanced Minority,
Majority 5 5 {0.0, 0.1, …, 1.0} {0.02, 0.04, …, 0.2}

_Evaluation Metric_ To access the generalization performance of our approaches
fairly, especially concerning the number of classes $M$, we propose using a
normalized version of mean absolute error (MAE). Since the MI values range
between $0$ and $\lg(M)$, the MAE will likewise fall within this same range.
The normalizing factor, which is the entropy on $Y$, is computed using
$P_{Y}(\cdot)$ defined in (14) as
$H_{Y}(\mathsf{g}_{r},M,r)=-\sum_{m=1}^{M}P_{Y}(m)\cdot\lg\left(P_{Y}(m)\right)$.
For each approach and unique data set configuration
($\mathsf{g}_{r},d,M,r,\epsilon$), we obtain $10$ ground truths and estimated
MI values, as shown in Figure 4. The ground-truth MI values are represented as
a vector $\boldsymbol{I}=({I}_{1},\ldots,{I}_{10})$, and the estimated one as
$\boldsymbol{\hat{I}}=(\hat{I}_{1},\ldots,\hat{I}_{10})$. We use the NMAE
metric to evaluate the performance as

$\operatorname{m}_{\textsc{{\text{NMAE}}}}(\boldsymbol{I},\boldsymbol{\hat{I}})=\frac{1}{H_{Y}(\mathsf{g}_{r},M,r)}\left(\frac{1}{10}\sum_{j=1}^{10}\lvert\boldsymbol{I}[j]-\boldsymbol{\hat{I}}[j]\rvert\right)=\sum_{j=1}^{10}\frac{\lvert{I}_{j}-\hat{I}_{j}\rvert}{10\cdot
H_{Y}(\mathsf{g}_{r},M,r)}\,.$ (18)

_Implementation Details_ To generate synthetic data sets using the MVN
distribution, including the ground-truth MI calculation, we implemented a
script using the multivariate_normal and ortho_group functions provided by
scipy (Virtanen et al., 2020). For HPO, we used BO implemented by the
BayesSearchCV implemented by scikit-optimize (Head et al., 2021). The detailed
and documented code, including the implementation of all MI estimation
approaches and the synthetic data set generation processes, along with the
Python scripts used for experiments, can be found on our
GitHub111https://github.com/LeakDetectAI/automl_quant_il_detect.

For MI estimation approaches, we utilized two leading AutoML tools, TabPFN and
AutoGluon detailed in Section 4 to induce the Bayes predictor. Additionally,
we incorporated the implementations of five calibration techniques for the Cal
Log-Loss approaches provided by netcal (Küppers et al., 2020, 2023). We
implemented the GMM approach using the InfoSelect implementation, which we
extended to incorporate diagonal, spherical, and tied covariance
matrices222https://github.com/felipemaiapolo/infoselect/issues/5 (Polo and Da
Silva, 2020). Considering the estimation challenges faced by both the GMM and
the provided InfoSelect implementation with high-dimensional data sets, the
dimensionality reduction technique is implemented when the data set
dimensionality exceeds $100$ ($d>100$). The architectures and objective
function for MINE approach was implemented using PyTorch (Paszke et al.,
2019), while the MLP architectures including the CCE loss for the PC-Softmax
approach was implemented using Keras (Chollet et al., 2015).

###  5.3 Results

In this section, we provide a comprehensive analysis of our results on MI
estimation approaches using MVN synthetic data sets, as detailed in Table 1.
We discuss the overall performance of various approaches on system data sets
generated with perturbation and proximity techniques in Section 5.3.1. We
assess the generalization of top-performing approaches using AutoGluon and
TabPFN including the baselines across various factors, including the number of
classes ($M$), input dimensions ($d$), class imbalance ($r$), and noise levels
($\epsilon$) in Section 5.3.2. For a complete overview of generalization
capabilities across all techniques, please refer to Section C.2.

####  5.3.1 Overall Results

In Figure 5, we present a comprehensive analysis of various MI estimation
approaches applied to data sets generated using the MVN distribution to
simulate synthetic systems. We use bar charts to illustrate the mean and
standard error (SE) of the NMAE on both MVN perturbation and proximity
synthetic data sets. Our analysis covers different data set types generated by
synthetic systems, including balanced, binary-imbalanced, and multi-class
imbalanced data sets, as detailed in Table 1. For systems generating a
balanced data set, we examine two scenarios involving vulnerable synthetic
systems: one with no noise and another with a $50\text{\,}\mathrm{\char
37\relax}$ noise level. Additionally, considering the critical importance of
accurate MI estimation on non-vulnerable systems ($\epsilon=1.0$) for the
precise prediction of the existence of IL by the respective ILD approaches, we
include non-vulnerable systems in our analysis.

(a) Performance on MVN perturbation synthetic data sets

(b) Performance on MVN proximity synthetic data sets

Figure 5: Overall NMAE of MI estimation approaches on synthetic data set

5.3.1.1 MVN Perturbation Data Sets In Figure 5(a), we delve into the
performance insights of MI estimation approaches applied to synthetic data
sets generated using the MVN perturbation technique. In cases where systems
generate balanced classification data sets, we observe that the Log-Loss
estimation method, using TabPFN in combination with suitable calibration
techniques (Cal Log-Loss), stands out as a top performer. Especially, TabPFN
IR Cal Log-Loss exhibits exceptional performance, with an impressively low
mean NMAE (approximately $0.0026$ for $\epsilon=0.0$, $0.0056$ for
$\epsilon=1.0$) with low variation, with SE being approximately $0.0008$ for
$\epsilon=0.0$ and $0.0006$ for $\epsilon=1.0$. While in the case of a
vulnerable system with noise $\epsilon=5.0$, TabPFN IR Cal Log-Loss performs
worse with NMAE approximately $0.2115$ and SE being $0.0124$. In this case,
TabPFN PS Cal Log-Loss and Beta Cal Log-Loss perform very well with NMAE at
approximately $0.0666$ and SE being $0.0032$. However, it is worth noting
that, overall, the use of the HB Cal Log-Loss method with both AutoML
approaches does not yield substantial benefits, except in cases where it is
paired with AutoGluon for estimating MI in non-vulnerable systems generating
balanced data set.

In the context of systems generating binary-class classification imbalanced
data sets, the GMM baseline emerges as a standout performer, achieving the
lowest mean NMAE (approximately $0.019$) and maintaining a low SE
(approximately $0.0017$). Furthermore, the Log-Loss and IR Cal Log-Loss
approaches using TabPFN also exhibit strong performance in terms of NMAE, with
very low SEs. Conversely, when confronted with multi-class imbalanced data
sets, the IR Cal Log-Loss approach using TabPFN takes the lead by achieving
the lowest mean NMAE (approximately $0.0293$) and maintains a low SE
(approximately $0.0013$). The GMM baseline also performs well in this context,
highlighting its adaptability in estimating MI within systems generating
imbalanced data sets. Overall, TabPFN models combined with the Cal Log-Loss
approach, especially IR Cal Log-Loss, consistently demonstrate strong
performance in estimating MI in the MVN perturbation system data sets.

5.3.1.2 MVN Proximity Data Sets In Figure 5(b), we delve into the performance
insights of MI estimation approaches applied to synthetic data sets generated
using the MVN proximity technique. Notably, we find similar trends in the
performance of MI estimation approaches compared to the baselines as on the
MVN perturbation system data sets. Once again, the Log-Loss estimation method,
coupled with TabPFN and appropriate calibration techniques (Cal Log-Loss),
excels in scenarios involving systems generating balanced classification data
sets. Specifically, the TabPFN IR Cal Log-Loss method exhibits exceptional
performance for vulnerable systems with a mean NMAE of $0.0027$ for
$\epsilon=0.0$, $0.0104$ for $\epsilon=0.5$ and a remarkably low SE of
$0.0008$ for $\epsilon=0.0$ and $0.0011$ for $\epsilon=0.5$. In contrast, for
non-vulnerable systems, the TabPFN PS Cal Log-Loss and TabPFN TS Cal Log-Loss
approaches deliver the best performance, both recording mean NMAE of $0.0163$.
The TabPFN Log-Loss also exhibits strong performance in this context.

Like before, AutoGluon models display larger NMAE values, signifying the
superiority of TabPFN for this scenario. Similar to the perturbation data set,
in the case of binary-class imbalanced data sets, the GMM baseline emerges as
the top performer, achieving the lowest mean NMAE of $0.017$. TabPFN IR Cal
Log-Loss also showcases commendable performance with very low SE. In multi-
class imbalanced data sets, TabPFN IR Cal Log-Loss again maintains its
superior performance with the lowest mean NMAE of $0.01625$. The GMM baseline
also demonstrates strong performance, highlighting its versatility in
delivering accurate MI estimates across imbalanced synthetic data sets.
Overall, TabPFN with Cal Log-Loss, especially IR Cal Log-Loss, consistently
demonstrates strong performance in estimating MI for MVN perturbation system
data sets. In general, AutoGluon models display higher NMAE values, signifying
the superiority of TabPFN models for this scenario. As mentioned in Section
5.1.3, introducing noise using perturbation techniques makes the generated
data sets more complex. Thus, the overall performance of all the models in
estimating MI within the synthetic data sets generated using the proximity
technique is slightly better than that by using the perturbation technique.

5.3.1.3 Summary In summary, TabPFN, particularly IR Cal Log-Loss, consistently
demonstrates remarkable MI estimation capabilities across a range of scenarios
for both synthetic systems simulated using MVN perturbation and proximity
techniques. The GMM baseline consistently stands out as an outlier, excelling
in accurately estimating MI for imbalanced synthetic data set. Nevertheless,
it is worth noting that while AutoGluon models deliver respectable results,
they do not reach the same level of prominence as TabPFN within these data
sets. In general, we observe that using calibration techniques (Cal Log-Loss)
does not significantly impact the performance for TabPFN Log-Loss, especially
for the balanced synthetic data sets, which suggests that TabPFN produces
well-calibrated class probabilities, which is also observed in Section C.2. In
contrast, using appropriate calibration techniques (Cal Log-Loss) mostly
significantly improves the performance of AutoGluon Log-Loss for estimating MI
in vulnerable synthetic data sets, thus hinting at miscalibrated class
probabilities returned by the resulting models, which is also observed in
Section C.2. The usage of calibration techniques (Cal Log-Loss) significantly
deteriorates the performance of AutoGluon Log-Loss by overestimating the MI
due to overfitting in non-vulnerable systems, which could hamper the ILD
performance of this approach by detecting non-existent ILs and producing false
positives.

####  5.3.2 Generalization Capability Analysis

In this section, we delve into the generalization capabilities of different MI
estimation approaches, focusing on their performance relative to baseline
methods. We identify the best-performing MI estimation approaches using both
AutoGluon and TabPFN to assess their generalization capability. Our selection
process is based on a comprehensive analysis of the overall results presented
in Section 5.3.1. Specifically, we calculate mean NMAE for each data set type,
including balanced, binary-class, and multi-class imbalance. For instance, to
select the best performing TabPFN model (out of Mid-Point, Log-Loss, IR Cal
Log-Loss, PS Cal Log-Loss, Beta Cal Log-Loss, TS Cal Log-Loss, and HB Cal Log-
Loss) on balanced data sets, we evaluate the mean NMAE across all experiments
as outlined in Table 1 and then choose the best-performing approach.

We gauge the generalization capability with respect to the number of classes
($M$) and input dimensions ($d$) for a non-vulnerable system and two types of
vulnerable synthetic systems ($\epsilon=0.0$ and $\epsilon=0.5$) generating
balanced data sets. To assess generalization with respect to the number of
classes, we aggregate the NMAE across all input dimensions
($d\in\\{2,4,\ldots,20\\}$) for each possible number of classes
($M\in\\{2,4,\ldots,10\\}$) Similarly, we assess generalization with respect
to the input dimensions by aggregating NMAE across the number of classes. To
assess generalization with respect to class imbalances, we aggregate the NMAE
across different noise levels ($\epsilon\in\\{0.0,0.1,\ldots,1.0\\}$) for each
class imbalance parameter ($r$) value, such that
$r\in\\{0.05,0.1,\ldots,0.5\\}$ for binary-class imbalanced data sets and
$r\in\\{0.02,0.04,\ldots,0.2\\}$ for multi-class imbalanced data sets.
Similarly, we assess generalization with respect to the noise levels by
aggregating NMAE across the class imbalances.

For a visual representation of our filtered generalization results using a
series of line plots, we refer to Figures 6 and 7 for results on MVN
perturbation and MVN proximity data sets respectively. The best models are
selected based on the overall results shown in Figures 5(a) and 5(b),
respectively. In these figures, the Y-axis displays the mean NMAE for MI
estimates obtained by the selected approaches with respect to various factors
of system data sets. This comprehensive analysis allows us to assess the
generalization capabilities of our proposed MI estimation approaches compared
to baselines under diverse scenarios of the number of classes ($M$), input
dimensions ($d$), class imbalance ($r$), and noise level ($\epsilon$)
(referred to as flip percentage ($\epsilon$) in case of MVN perturbation data
set).

Figure 6: Generalizability of MI estimation approaches on MVN perturbation
synthetic data sets

5.3.2.1 Number of Classes ($M$) In the top row of Figures 6 and 7, we analyze
the generalization capability of MI estimation methods across data sets with
varying numbers of classes ($M$). These figures pertain to results on MVN
perturbation and MVN proximity data sets, respectively, with the X-axis
ranging from $2$ to $10$, representing the number of classes.

_Vulnerable Systems_ In vulnerable systems, AutoGluon and TabPFN consistently
exhibit high MI estimation accuracy (lower mean NMAE), with increasing the
number of classes. Overall TabPFN consistently outperforms AutoGluon with
AutoGluon showing notable improvements, especially in synthetic systems
simulated using perturbation technique with a $50\text{\,}\mathrm{\char
37\relax}$ noise level ($\epsilon=0.5$). In the case of the GMM baseline, the
MI estimation accuracy sometimes improves (mean NMAE decreases) for vulnerable
systems, especially for multi-class balanced synthetic data sets (classes more
than $2$). This observation can be attributed to the loss of information due
to aggregation. Since the number of classes ($M$) in vulnerable synthetic data
sets increases, the precision of the MI estimation by GMM baseline tends to
degrade, as illustrated in Figures 12 and 13. In contrast, the performance of
MINE and PC-Softmax baselines declines as the number of classes increases,
with MINE being the least accurate, also illustrated in Figures 12 and 13.

_Non-vulnerable Systems_ Similar to vulnerable systems, both AutoGluon and
TabPFN consistently demonstrate improved MI estimation (lower mean NMAE) with
increasing the number of classes. AutoGluon shows more significant
improvements, especially for multi-class balanced synthetic data sets (classes
more than $2$). TabPFN consistently outperforms AutoGluon for non-vulnerable
systems. In contrast, the performance of baselines deteriorates as the number
of classes increases, especially for the GMM baseline in the case of
perturbation non-vulnerable system data sets. As an exception to estimating MI
in vulnerable systems, MINE and PC-Softmax baselines perform better than GMM,
particularly for perturbation system data sets.

5.3.2.2 Input Dimensions ($d$) In the second row of Figures 6 and 7, we
analyze the generalization capability of MI estimation methods across data
sets with varying input dimensions ($d$). These figures pertain to results on
MVN perturbation and MVN proximity data sets, respectively, with the X-axis
ranging from $2$ to $20$, representing the input dimensions.

_Vulnerable Systems_ Notably, both AutoGluon and TabPFN consistently exhibit
high MI estimation accuracy (lower mean NMAE) as input dimensions increase
from $2$ to $20$. Overall TabPFN outperforms AutoGluon with AutoGluon showing
notable improvements, especially in synthetic systems simulated using
perturbation technique with a $50\text{\,}\mathrm{\char 37\relax}$ noise level
($\epsilon=0.5$). The GMM baseline shows significant deterioration in
performance with increasing input dimensions (especially from $10$ to $20$),
specifically for MVN proximity system data sets. While MINE and PC-Softmax
baselines also show deterioration in performance most of the time with
increasing input dimensions, this becomes significant, especially beyond $10$,
also illustrated in Figures 12 and 13. These findings conclude that even using
deep MLPs for MI estimation approaches sometimes does not work for high-
dimensional system data sets.

_Non-vulnerable Systems_ Similar to vulnerable systems, both AutoGluon and
TabPFN exhibit improved MI estimation accuracy (lower mean NMAE) as input
dimensions increase from $2$ to $20$, with TabPFN outperforming AutoGluon and
AutoGluon exhibiting more significant improvements. Again, the GMM baseline
shows significant deterioration in performance with increasing input
dimensions (especially from $10$ to $20$), specifically for MVN perturbation
data sets. While MINE and PC-Softmax baselines consistently perform well with
increasing input dimensions for MVN perturbation data sets, their performance
deteriorates in the case of MVN proximity data sets.

Figure 7: Generalizability of MI estimation approaches on MVN proximity
synthetic data sets

5.3.2.3 Class Imbalance ($r$) In the third row of Figures 6 and 7, we analyze
the generalization capability of MI estimation methods across data sets with
different levels of class imbalance ($r$). These figures pertain to results on
MVN perturbation and MVN proximity data sets, respectively, with the X-axis
ranging from $0.05$ to $0.5$ for binary-class imbalanced data sets and from
$0.02$ to $0.2$ for multi-class imbalanced data sets, representing the degree
of class imbalance.

_Binary-class Imbalanced_ In the systems generating imbalanced binary-class
data set, TabPFN and GMM consistently exhibit high MI estimation accuracy as
the class imbalance in the synthetic data sets decreases from $0.0$ to $0.5$.
Notably, the GMM baseline exhibits a slight advantage over TabPFN,
particularly in MVN perturbation synthetic data set. In contrast, the
performance of AutoGluon tends to deteriorate as class imbalance decreases in
most cases. Moreover, MINE and PC-Softmax baselines also exhibit performance
deterioration, with the decline becoming particularly significant beyond a
class imbalance of $0.25$, especially in MVN perturbation synthetic data set.
These intriguing observations suggest that AutoGluon, MINE, and PC-Softmax
surprisingly perform better on imbalanced binary-class synthetic data sets.
However, it’s essential to note that they still underperform compared to
TabPFN and GMM.

_Multi-class Imbalanced_ TabPFN and GMM continue to display high MI estimation
accuracy as the class imbalance in the synthetic data sets decreases ($0.0$ to
$0.2$). In this scenario, TabPFN shows a slight advantage over the GMM
baseline, especially in MVN perturbation synthetic data set. Surprisingly, the
performance of AutoGluon also improves in most cases as class imbalance
decreases. On the other hand, the performance of MINE and PC-Softmax baselines
remains relatively consistent, with only a noticeable performance drop
occurring at a class imbalance of $0.04$. These intriguing findings contrast
the observations for systems generating imbalanced binary-class synthetic data
sets, highlighting the different behavior of AutoGluon, MINE, and PC-Softmax
in handling imbalanced multi-class synthetic data sets.

5.3.2.4 Noise Level ($\epsilon$) In the last row of Figures 6 and 7, we
analyze the generalization capability of MI estimation methods across data
sets with varying noise level ($\epsilon$) in both binary-class and multi-
class imbalanced scenarios. These figures pertain to results on MVN
perturbation and MVN proximity data sets, respectively, with the X-axis
ranging from $0.0$ to $1.0$, representing the noise level.

_Binary-class Imbalanced_ In the systems generating imbalanced binary-class
data set, TabPFN and GMM baseline demonstrated high MI estimation accuracy
consistently as the noise level in the synthetic data sets increased ($0.0$ to
$1.0$). The GMM baseline shows a slight advantage over the TabPFN, especially
in MVN perturbation synthetic data set. In contrast, the performance of
AutoGluon tends to deteriorate as the noise level in the synthetic data sets
increases. On the other hand, both MINE and PC-Softmax baselines exhibited
significant improvements in performance with increasing noise levels in the
system data sets. The findings suggest that TabPFN and GMM baseline are
reliable options for estimating MI for systems generating noisy imbalanced
binary-class synthetic data sets.

_Multi-class Imbalanced_ TabPFN and GMM continue to maintain high MI
estimation accuracy as the noise level in the synthetic data sets increases
($0.0$ to $1.0$). In this scenario, TabPFN showed a slight advantage over the
GMM baseline, especially in MVN perturbation synthetic data set. The
performance of AutoGluon continued deteriorating with increasing noise levels
in the synthetic data set. MINE and PC-Softmax baselines continued to exhibit
significant improvements with increasing noise levels in the synthetic data
set. These findings highlight the robustness of TabPFN and GMM in estimating
MI for systems generating noisy imbalanced multi-class synthetic data sets.

5.3.2.5 Summary In this section, we summarize the key takeaways regarding the
generalization capabilities of the chosen approaches for balanced and
imbalanced synthetic data sets. These findings consistently apply to the MVN
perturbation and proximity synthetic data set.

_Influence of Number of Classes ( $M$) and Input Dimensions ($d$)_ Our
performance analysis using NMAE, across various numbers of classes ($M$) and
input dimensions ($d$), consistently demonstrates strong generalization
capabilities for TabPFN and AutoGluon. In contrast, the baselines display
diminished generalization capability, as the performance consistently
deteriorates when the number of classes and input dimensions increases.
Notably, GMM encounters significant challenges when dealing with high-
dimensional synthetic data sets, resulting in substantial performance
degradation beyond $10$ input dimensions, aligning with the well-known
limitations of GMM in high-dimensional spaces. Furthermore, MINE demonstrates
the lowest generalization capability concerning the number of classes and
input dimensions, which suggests that it may not be well-suited for estimating
MI between multi-dimensional inputs and one-dimensional outputs in
classification synthetic data sets. A comprehensive analysis of the
generalization capabilities of all approaches, considering both the number of
classes ($M$) and input dimensions ($d$), without aggregating results across
either dimension, is provided in Section C.2.1.

_Influence of Class Imbalance ( $r$) and Noise ($\epsilon$)_ Our performance
analysis across various class imbalances ($r$) and noise levels ($\epsilon$)
in imbalanced synthetic data sets highlights the robust generalization
capabilities of TabPFN and GMM. Notably, the strong generalization capability
of GMM can be attributed to the low dimensionality ($d=5$) of the imbalanced
synthetic data set. This characteristic suggests that GMM primarily struggles
with estimating MI in high-dimensional synthetic data sets, remaining largely
unaffected by noise and class imbalance. In contrast, AutoGluon, MINE, and PC-
Softmax baselines did not exhibit strong generalization capabilities
concerning variations in class imbalance ($r$) and noise levels ($\epsilon$)
within the imbalanced synthetic data set. A comprehensive analysis of the
generalization capabilities of all approaches, considering both the class
imbalance ($r$) and noise levels ($\epsilon$), without aggregating results
across either dimension, is provided in Section C.2.2.

_Conclusion_ TabPFN based approaches consistently demonstrates remarkable
generalization capabilities in estimating MI across various factors for both
MVN perturbation and proximity data sets.

##  6 Information Leakage Detection Approaches

In this section, we outline the ILD process in a given system, as depicted in
Figure 8(b), and present various approaches, including baseline methods, shown
in Figure 8(a). As discussed in Section 2, detecting IL in a system involves
leveraging MI and approximating the predictive performance of Bayes predictor
on the system data set $\mathcal{D}$. Consequently, we partition our ILD
approaches into two parts: one based on MI, and the other centered on
classification, explained in Sections 6.2 and 6.3 respectively. The MI-based
approach employs OTT on MI estimates, obtained through various techniques,
including Log-Loss, Mid-Point, and baseline approaches, detailed in Section 3.
For the classification-based approach, we introduce two methodologies from our
prior work (Gupta et al., 2022), employing FET and PTT to assess classifier
performance in detecting ILs. We recommend using AutoML tools, specifically
AutoGluon and TabPFN, for ILD, as outlined in Section 4. To enhance IL
detection confidence, we propose applying Holm-Bonferroni correction on
$p$-values acquired from top-$10$ performing AutoML pipelines or models
obtained through HPO, as described in Section 6.4.

###  6.1 IL Detection Process

The IL detection process begins with the acquisition of accurate MI or Bayes
predictor performance estimates obtained from the models. This is achieved
using nested $K$-fold cross-validation (KFCV) coupled with HPO. The objective
functions employed for HPO include BER for PC-Softmax, AutoGluon, and TabPFN,
AIC for GMM, and MSE for MINE with specific parameter ranges provided in Table
4. The process begins by partitioning the data set $\mathcal{D}$ generated by
the system into a training data set ($90\text{\,}\mathrm{\char 37\relax}$ of
instances) and a test data set ($10\text{\,}\mathrm{\char 37\relax}$ of
instances). We conduct HPO for $100$ function evaluations using MCCV employing
$3$ splits on the training data set. Each split reserves
$30\text{\,}\mathrm{\char 37\relax}$ of the training data for validation,
enabling an effective performance assessment of current parameter values.
Refer to Table 4 for the hyperparameters and their corresponding ranges for
the approaches. This process is annotated as “N$10$F-$3$SCV” signifying the
use of KFCV with $K=10$ ($10$F) for evaluation and $3$S MCCV splits for HPO
using nested cross-validation, as depicted in Figure 8(a). Subsequently, we
identify the top-performing models using HPO. In the case of AutoGluon, the
tool is executed for a fixed duration of $1800$ seconds on $\mathcal{D}$ to
select the best-performing pipelines. KFCV is then applied to the top-$10$
performing models or pipelines ($J=10$), generating $K=10$ estimates from
each, using the complete system data set $\mathcal{D}$, as depicted in Figure
8(b).

The next crucial step involves assessing $p$-values, derived from appropriate
statistical tests applied to the $K=10$ estimates for each of the top-$10$
performing models or pipelines to determine the presence of IL within the
system, as depicted in Figure 8(b). To compute these $p$-values, the OTT is
utilized on the $10$ MI estimates of $j$-th best-performing model or pipeline
(denoted as $\boldsymbol{\hat{I}}_{j}$), as elaborated in Section 6.2.
Additionally, the PTT is employed to compare the accuracies of the $j$-th
best-performing AutoML pipeline (denoted as $\boldsymbol{a}_{j}$) with those
of the marginal Bayes predictor accuracies (denoted by $\boldsymbol{a}_{mc}$),
as detailed in Section 6.3.1. We also apply FET on the confusion matrices of
AutoML pipelines (denoted as
$\mathcal{M}_{j}=\\{\boldsymbol{M}_{j}^{k}\\}_{k=1}^{K}$), and the results are
then aggregated using both the _mean_ and _median_ operators, as detailed in
Section 6.3.2.

To enhance the confidence in the IL detection decision, multiple $p$-values
obtained from top-$10$ performing AutoML pipelines or models through HPO using
Holm-Bonferroni correction, as described in Section 6.4. Based on our prior
work (Gupta et al., 2022), we concluded that setting the _rejection threshold_
equal to $\lfloor\tfrac{J}{2}\rfloor$, on the cut-off parameter $\tau$
provides a robust and accurate IL detection in a system. In this context, we
aggregate $J=10$ $p$-values using the Holm-Bonferroni correction by using the
condition $\tau\geq 5$ to deduce the occurrence of IL in the system. In
summary, the ILD process involves rigorous estimation, statistical assessment,
and the use of correction techniques to confidently detect the presence of IL
within a given system.

(a) Calculation of $p$-value by different ILD approaches

(b) IL detection process

Figure 8: Procedure of using ILD approaches to detect ILs in a system
generating $\mathcal{D}$

###  6.2 MI-based Approach

In this section, we will describe the process of applying OTT on the MI
estimates obtained through Log-Loss and Mid-Point, including the state-of-the-
art techniques to perform ILD in a given system. These approaches are built
upon the condition for the existence of IL in a system, as discussed in
Section 2.3. According to information theory, IL occurs in a system if the MI
between system input and output is significantly greater than $0$. To ensure
precise estimation of MI using Log-Loss and Mid-Point techniques, we suggest
utilizing AutoML tools, specifically AutoGluon and TabPFN, described in
Section 4.

_One-Sample T-Test Based Approach_ The one-sample t-test (OTT) is a
statistical test used to determine if the mean of a sample differs
significantly from a known or hypothesized population mean (Demšar, 2006). In
the context of detecting IL in the system using MI estimates, we apply this
test to evaluate whether the MI values obtained from the system significantly
deviate from the expected or null hypothesis values, in this case being $0$.
To deduce the presence of IL in a system, using the OTT, we start by
considering the null hypothesis ($H_{0}(\boldsymbol{\hat{I}}_{j}\sim 0)$) and
the alternative hypothesis ($H_{0}(\boldsymbol{\hat{I}}_{j}\gg 0)$). In this
case, the null hypothesis ($H_{0}$) would state that the MI for the system is
approximately equal to the expected mean ($\approxeq 0$), indicating no IL.
Conversely, the alternative hypothesis ($H_{1}$) suggests that there is a
significant deviation from the expected mean ($\gg 0$), implying the presence
of IL. We collect $10$ MI estimates from the system using KFCV for $K=10$,
which are treated as sample data points and denote them by the vector
$\boldsymbol{\hat{I}}_{j}=(\hat{I}_{1},\ldots,\hat{I}_{10})$. Next, we collect
the sample mean ($\mu_{\boldsymbol{\hat{I}}_{j}}$) and the standard deviation
($\sigma_{\boldsymbol{\hat{I}}_{j}}$) of these MI estimates. The sample mean
is an estimate of the population mean, while the sample standard deviation
measures the variability within our sample. The $t$-statistic quantifying the
amount of SE the sample mean is away from the expected population mean under
the null hypothesis is calculated as
$t=\frac{\mu_{\boldsymbol{\hat{I}}_{j}}-\mu_{0}}{\nicefrac{{\sigma_{\boldsymbol{\hat{I}}_{j}}}}{{\sqrt{K}}}}$,