# Generative Conditional Distributions by Neural (Entropic) Optimal Transport

Bao Nguyen  Binh Nguyen  Hieu Trung Nguyen  Viet Anh Nguyen

###### Abstract

Learning conditional distributions is challenging because the desired outcome
is not a single distribution but multiple distributions that correspond to
multiple instances of the covariates. We introduce a novel neural entropic
optimal transport method designed to effectively learn generative models of
conditional distributions, particularly in scenarios characterized by limited
sample sizes. Our method relies on the minimax training of two neural
networks: a generative network parametrizing the inverse cumulative
distribution functions of the conditional distributions and another network
parametrizing the conditional Kantorovich potential. To prevent overfitting,
we regularize the objective function by penalizing the Lipschitz constant of
the network output. Our experiments on real-world datasets show the
effectiveness of our algorithm compared to state-of-the-art conditional
distribution learning techniques. Our implementation can be found at
https://github.com/nguyenngocbaocmt02/GENTLE.

Conditional Distribution, Optimal Transport, Generative Models

##  1 Introduction

The conditional distribution of a response variable given the covariate
information is a key quantity for numerous tasks in machine learning and data
science. For example, knowing the conditional distribution of a patient’s
health outcome (the response) given the patient’s characteristics (the
covariate) improves risk assessment and prediction accuracy and allows
personalized interventions wherein healthcare providers can tailor treatments
based on individual patient profiles. It is also a useful tool to evaluate the
effects of a policy in economics (Díaz, 2020; Athey & Imbens, 2015). Access to
a response’s conditional probability distribution also leads to an estimate of
the treatment effect without explicit, additional experiments.

Unfortunately, learning the relationship between the response and the
covariate is challenging, especially in high-dimensional covariate space or
complex dependency settings. Many simple and interpretable methods, such as
kernel density estimation, may struggle to capture the complexity of the data-
generating conditional distribution accurately (Athey et al., 2021). On the
other hand, there are overly complex (deep neural network) models that perform
well on the training data but generalize poorly to new data; this overfitting
behavior triggers unreliable predictions and interpretations. Additionally,
performing inferences on these overly complex models can be time-consuming and
environmentally unfriendly. Alternatively, using additional expert inputs, one
can build physically informed simulation models to represent the real physical
world. However, running simulations can be time-consuming, making them ill-
suited for real-time decision-making.

Another critical obstacle underlying conditional distribution learning is that
we are not simply learning one distribution; instead, we need to learn a
collection of distributions, one conditional distribution for each possible
value of the covariate information. This challenge intensifies because, in
reality, the data is scarce. This scarcity is particularly prevalent in
healthcare, where privacy concerns or budget constraints prohibit the
collection and dissemination of samples. As an illustrative example, we will
consider the Lalonde-Dehjia-Wahba (LDW) (LaLonde, 1986; Dehejia & Wahba, 1999)
dataset, which is a widely used dataset for analyzing the average treatment
effects (Athey et al., 2021). The LDW dataset comprises 16,177 samples
aggregated from actual experiments and the Current Population Survey. Each
observed sample includes nine covariate attributes of a person and a scalar
response reporting the individual’s earnings in 1978. Since two persons may
have the exact attributes, it happens that for a specific covariate value,
there are multiple observations of the response. The histogram in Figure 1
shows the count numbers per distinct covariates for the LDW dataset. We
observe that more than 13,000 covariates have only one response. This is an
extreme case where we must learn the conditional distribution with a single
data point observed at that covariate. On the other hand, we also observe
several popular covariates with many responses: the most popular covariate has
47 responses. This again induces unbalanced sampling in which the number of
responses is not balanced across different covariate values. A naive solution
to fix this problem could be under-sampling (remove data) the more frequent
covariate and over-sampling (create synthetic data) for the rare covariate to
balance the dataset using popular methods such as SMOTE (Chawla et al., 2002),
SMOTER (Torgo et al., 2013) and SMOGN (Branco et al., 2017). Unfortunately,
these methods do not scale well for large datasets with more than thousands of
covariates.

Figure 1: Histogram count of Number of observed responses for distinct
covariate values in the LDW dataset (LaLonde, 1986). This dataset has over
13,000 covariates with only one response.

The current literature on conditional distribution learning reveals a notable
gap in methodologies that can effectively address the aforementioned
difficulties. Specifically, our experimental section shows that CWGAN (Athey
et al., 2021) suffers from mode collapse and fails to recover the data-
generating conditional distributions.

###  1.1 Problem Statement and Contributions

We consider $X$ a multi-dimensional covariate and $Y$ a one-dimensional
response. For any specific value $x$, the data-generating conditional
distribution of $Y|X=x$ is denoted by $\mu_{x}$. Our goal is to learn these
conditional distributions $\mu_{x}$ uniformly over all $x$ from a finite
dataset
$\mathcal{D}=\\{(x_{i},\\{y_{i,k}\\}_{k=1,\ldots,N_{i}})\in\mathbb{R}^{d}\times\mathbb{R}$.
Without any loss of generality, we suppose that the covariates $x_{i}$ are
distinct, and for each $x_{i}$, $\\{y_{ik}\\}_{k=1,\ldots,N_{i}}$ are $N_{i}$
independent samples drawn from $\mu_{x_{i}}$.

More importantly, we aim to learn a generative model of the conditional
distributions, which can allow us to draw new, unseen samples. This task is
more challenging than learning specific finite-dimensional statistics of the
conditional distributions, such as learning conditional expectations or
quantiles, especially under a small training sample size. We reemphasize that
conditional distribution learning is challenging: our goal is not simply
learning one distribution, but instead, a collection of distributions – one
conditional distribution for each possible value of the covariate information.
This challenge intensifies because of data scarcity. This difficulty is often
encountered in real-world deployment, for instance, in health care and
finance.

Contributions. We introduce GENTLE, a generative neural transport learning
model for conditional distributions. GENTLE is represented by a neural network
$T_{\theta}$ that takes the covariate $x$ and a uniform $(0,1)$ random
variable $U$ as input, and it outputs $T_{\theta}(x,U)$ with a distribution
that is close to $\mu_{x}$, the data-generating conditional distribution of
$Y|X=x$. Due to the low sample size, $T_{\theta}$ can overfit, and we propose
learning $\theta$ by minimizing the sum of the fitness and the overfitting
regularization trade-off. The overfitting regularizer ensures that
$T_{\theta}$ is sufficiently smooth in the covariate $x$, and we use an
entropic optimal transport distance to measure the discrepancy of the
generative distribution outputs. To this end, we create a second network
$v_{\phi}$ to approximate the conditional potential associated with the semi-
dual formulation of the entropic optimal transport distance. This leads to a
min-max optimization problem over the pair of generative-conditional potential
network parameters $(\theta,\phi)$. This min-max problem is then solved using
a state-of-the-art gradient descent-ascent algorithm.

Using real-life datasets, we empirically validate that GENTLE performs better
than state-of-the-art methods such as conditional Wasserstein GANs across
multiple performance metrics (Wasserstein metric and Kolmogorov-Smirnov
metric). In particular, even for the case of LDW-CPS data that exhibits severe
imbalances in observed responses, GENTLE still demonstrates a remarkable
ability to generate observations that follow ground truth distributions
closely, quantitatively, and qualitatively.

This paper unfolds as follows: Section 2 provides an overview of the existing
literature. Section 3 reviews the technical background on optimal transports
and the tools we will utilize. Section 4 presents our approach to generative
conditional distributions. The numerical results are provided in Section 5.

Notations. We write the capital letter $X$ for the random variable and the
case letter $x$ for a specific realization of $X$. For simplicity, we use
$Y(x)$ to denote the conditional random variable $Y|X=x$. For a measurable map
$f:\mathbb{R}\to\mathbb{R}$ and any probability measure $\mathbb{P}$ on
$\mathbb{R}$, we use $f_{\\#}\mathbb{P}$ to denote the pushforward measure of
$\mathbb{P}$ under $f$. The (standard) uniform distribution on $(0,1)$ is
denoted $\mathcal{U}(0,1)$, and we use $\mathbb{U}$ to represent the
probability measure associated with this uniform distribution. We denote
$N=\sum_{i}N_{i}$ as the total number of observations in the training dataset.

##  2 Related Works

Metamodeling is a computational approach that efficiently captures the
relationship between input variables and key output statistics like means or
quantiles. This is crucial for real-time simulation applications such as
personalized decision making (Shen et al., 2021) and online risk monitoring
(Jiang et al., 2020). Traditional metamodeling techniques, while effective,
often require pre-specification of relevant statistics, limiting their
flexibility. To address this, generative metamodeling has emerged as a
promising technique that estimates conditional distributions quickly and
dynamically. Generative metamodeling differs from traditional approaches by
simulating a wider range of possible outcomes, providing a more comprehensive
analysis of potential scenarios. Athey et al. (2021) propose using a
Wasserstein Generative Adversarial Model (WGAN, Arjovsky et al. 2017) to
simulate the distribution of outcomes conditioned on covariates, offering a
novel approach to outcome simulation. Similarly, Hong et al. (2023) propose
quantile-regression-based generative metamodeling (QRGMM), another method of
conditional distribution learning. QRGMM leverages the conditional quantile
regression technique to provide a more detailed understanding of the
distribution tails, which is crucial in risk assessment and decision-making
processes. These advancements in generative metamodeling can significantly
enhance the accuracy and applicability of real-time simulations, paving the
way for more personalized and precise decision-making tools in various fields.

Synthetic evaluation data. Previous works (van Breugel et al., 2024; Parikh et
al., 2022) leverage the deep generative model to approximate conditional
distribution for synthetic data generation. (Shrivastava & Tuzel, 2019)
parameterize a conditionally Gaussian model with a Deep Neural Network to
model error distribution for each state of the object.

Neural Optimal Transport is a recent emergent framework that proposes to learn
the optimal transport plan by deep neural networks. Various existing works in
this direction include parameterization of the OT map as the gradient of an
input convex neural network (Jacob et al., 2019; Makkuva et al., 2020; Bunne
et al., 2022, 2023); or solutions of convex-concave problems (Korotin et al.,
2022; Gushchin et al., 2022; Choi et al., 2023). However, these formulations
are primarily suited to computer vision or computational biology applications.
To the best of our knowledge, our paper is the first work attempting to
formulate the neural OT framework into a conditional distribution learning
task.

Learning and decision-making with conditional distributions also attracted
tremendous attention in artificial intelligence and operations research
communities. Nguyen et al. (2020b) studied the robust parameter estimation of
conditional distributions under the parametric exponential family of
distributions, while Nguyen et al. (2020a) studied the non-parametric
Wasserstein robust estimation of conditional statistics. Making conditional
decisions based on observed covariates is a well-known problem in finance
(Nguyen et al., 2024) and operations research (Esteban-Pérez & Morales, 2022);
interested readers are invited to the survey by Sadana et al. (2024) for
recent developments.

##  3 Background

We focus on a one-dimensional response space for $Y$. The optimal transport
distance between two distributions supported on $\mathbb{R}$ is formally
defined as follows.

######  Definition 3.1 (Optimal Transport Distance, Villani 2009).

Given two probability measures $\mu$ and $\nu$ supported on $\mathbb{R}$ with
finite second moments, the optimal transport (2-Wasserstein) distance between
them is defined as

$\mathrm{W}_{2}^{2}(\mu,\nu)\stackrel{{\scriptstyle\text{def.}}}{{=}}\inf_{\pi\in\Pi(\mu,\nu)}\int_{\mathbb{R}\times\mathbb{R}}|y-y^{\prime}|^{2}\,\mathrm{d}\pi(y,y^{\prime}),$
(1)

where $\Pi(\mu,\nu)$ denotes the set of probability distributions on
$\mathbb{R}^{2}$ with marginals $\mu$ and $\nu$, respectively.

#### Optimal transport in 1d.

In one dimension, the 2-Wasserstein distance admits a closed form expression
(Bobkov & Ledoux, 2019, Theorem 2.10): for two measures $\mu$ and $\nu$ with
cumulative density functions (CDF) $F_{\mu}$ and $F_{\nu}$, respectively, we
have

$\displaystyle\mathrm{W}_{2}^{2}(\mu,\nu)$
$\displaystyle=\int_{0}^{1}|F_{\mu}^{-1}(t)-F_{\nu}^{-1}(t)|^{2}\mathrm{d}t$
(2)
$\displaystyle=\mathbb{E}_{U\sim\mathcal{U}(0,1)}\left[|F_{\mu}^{-1}(U)-F_{\nu}^{-1}(U)|^{2}\right].$

#### Entropic-regularized optimal transport.

The recent popularity in applications of OT to machine learning is largely due
to its efficient computation by adding negative entropy of the transport plan
as a form of regularization to the linear OT program (1). This problem is
called entropic optimal transport (EOT), which is smooth and strongly convex,
therefore has a unique minimizer:

$\mathrm{W}_{2,\varepsilon}^{2}(\mu,\nu)\stackrel{{\scriptstyle\text{def.}}}{{=}}\min_{\pi\in\Pi(\mu,\nu)}\int_{\mathbb{R}^{2}}|y-y^{\prime}|^{2}\mathrm{d}\pi-\varepsilon
H(\pi).$ (3)

Cuturi (2013) showed that the minimization problem (3) can be solved by
Sinkhorn’s algorithm (a type of matrix scaling algorithm, also called
iterative proportional fitting procedure, Kullback 1968) that can be
implemented at large scale and is analytically tractable. In addition to its
computational advantages, problem (3) can also be interpreted as a specific
instance of the Schrödinger bridge problem, which has a rich history in
physics. For an in-depth theoretical introduction to EOT, readers can refer to
Nutz (2021) or Léonard (2014).

A nice property of the EOT problem (3) is that it allows an equivalent smooth
unbounded semi-dual formulation (Genevay et al., 2016, Proposition 2.1)

$\displaystyle\min_{\pi\in\Pi(\mu,\nu)}\int_{\mathbb{R}^{2}}|y-y^{\prime}|^{2}\mathrm{d}\pi-\varepsilon
H(\pi)$ (4)
$\displaystyle=\sup_{v(\cdot)\in\mathcal{C}}\int_{\mathbb{R}}v^{\varepsilon}(y)\mathrm{d}\mu+\int_{\mathbb{R}}v(y^{\prime})\mathrm{d}\nu-\varepsilon,$

where $\mathcal{C}$ is the space of all continuous functions from $\mathbb{R}$
to $\mathbb{R}$, and $v^{\varepsilon}$ is defined as

$v^{\varepsilon}(y)=-\varepsilon\log\left(\int_{\mathbb{R}}\exp\left(\frac{v(y^{\prime})-\lvert
y-y^{\prime}\rvert^{2})}{\varepsilon}\right)\mathrm{d}\nu\right).$ (5)

The function $v$ is called Kantorovich potential, along with its smoothed
$c$-transform $v^{\varepsilon}$, with $c$ being the absolute difference bi-
function. Using the semi-dual form is that we need to optimize over only one
function $v$ associating with $\nu$, and we can obtain $v^{\varepsilon}$ as
the potential associating with $\mu$.

##  4 Methodologies

We propose to learn a mapping $T:\mathbb{R}^{d}\times(0,1)\to\mathbb{R}$ such
that uniformly over all $x$, $T(x,U)=Y(x)$ in distribution when
$U\sim\mathcal{U}(0,1)$. The existence of this mapping $T(x,U)$ is guaranteed
by the following celebrated result.

######  Lemma 4.1 (Noise Outsourcing, Austin 2015).

Suppose that $\mathcal{X}$ and $\mathcal{Y}$ are standard Borel spaces and
that $(X,Y)$ is an $(\mathcal{X}\times\mathcal{Y})$-valued random variable.
Then, there are random variables $U\sim\mathcal{U}(0,1)$ coupled with
$\mathcal{X}$ and $\mathcal{Y}$ and a Borel function
$T:\mathcal{X}\times(0,1)\rightarrow\mathcal{Y}$ such that $U$ is independent
from $X$ and

$(X,Y)=(X,T(X,U))\quad\text{almost surely}.$

In general, the mapping $T$ is not unique. Following the neural OT framework
(Korotin et al., 2022; Makkuva et al., 2020; Bunne et al., 2022), we use a
neural network to parametrize a possible instance of $T(x,U)$, and we denote
this network as $T_{\theta}(x,U)$. For any specific value of the parameter
$\theta$, the network outputs a probability measure
$T_{\theta}(x,\cdot)_{\\#}\mathbb{U}$, which represents the distribution of
$T_{\theta}(x,U)$ with $U\sim\mathcal{U}(0,1)$. We propose the following
regularized loss to find the optimal map $T_{\theta}$:

$\min_{\theta}\mathrm{Fit}(\theta)+\lambda\mathrm{Reg}(\theta),$ (6)

where the first term is a measure of the fitness of the map, the second is a
regularization to prevent overfitting, and $\lambda\geq 0$ is the
regularization weight. We explain the detailed formulation of each term and
its rationale in the next sections.

###  4.1 Fitness Measure

Motivated by the OT distance in (2), a natural choice to measure the
discrepancy between the network outputs and the conditional data-generating
distributions is by integrating the discrepancy of the respective inverse CDFs
over all possible values $x$, i.e. , by evaluating

$\displaystyle\mathbb{E}_{X\sim\mathcal{D},U\sim\mathcal{U}(0,1)}[|F_{\mu_{X}}^{-1}(U)-T_{\theta}(X,U)|^{2}]$
(7) $\displaystyle=$
$\displaystyle\mathbb{E}_{X\sim\mathcal{D}}\left[\mathbb{E}_{U\sim\mathcal{U}(0,1)}[|F_{\mu_{X}}^{-1}(U)-T_{\theta}(X,U)|^{2}]\right].$

Above, $F_{\mu_{x}}$ is the cumulative distribution function of $Y(x)$, and
the outer expectation is taken over the (uniform) sampling of all distinct
values of $X$ in the training dataset $\mathcal{D}$. Unfortunately, the
inverse CDF $F_{\mu_{X}}^{-1}$ can be unbounded, and the magnitude of
$F_{\mu_{X}}^{-1}$ explodes to infinity for tail events whenever $U$ is to
close to the boundary 0 or 1. Hence, learning $T_{\theta}$ with the above
metric is not stable. To improve stability, we use the following fitness
function

$\mathrm{Fit}(\theta)\stackrel{{\scriptstyle\text{def.}}}{{=}}\mathbb{E}_{X\sim\mathcal{D}}\big{[}\mathbb{E}_{U\sim\mathcal{U}(0,1)}[|U-F_{\mu_{X}}(T_{\theta}(X,U))|^{2}]\big{]}.$
(8)

In practice, we do not have access to the real cumulative distribution of
$Y(x)$. However, we have a set of responses sampled from $Y(x)$, so we propose
to use a simple method such as Kernel Density Estimator (KDE) with the kernel
bandwidth of $h$ to have a smooth conditional density from responses of the
covariate $x$ in the training dataset. The kernel bandwidth plays a vital role
in our objective function, and tuning this parameter is important. In theory,
our KDE smoothness using the Gaussian kernel is strongly related to the
smoothed Wasserstein distance. From Nietert et al. (2021, Lemma 1 and
Corollary 1), the smooth-2-Wasserstein distance is continuous and
monotonically non-increasing in the Gaussian kernel bandwidth. This means that
the larger the bandwidth (the more smoothing applied), the smaller the fitness
term will become. In practice, we could implement an automatic kernel
bandwidth selection such as the Scott (Scott, 1979) or the Silverman
(Silverman, 2018) methods. In our implementation, the value of the bandwidth
is selected by a grid search to minimize the Wasserstein distance on the
validation set.

###  4.2 Regularization Measure

Minimizing the $\mathrm{Fit}(\theta)$ term can lead to overfitting due to the
low sample sizes for many values of $x_{i}$. To combat this problem, we add a
regularization term to promote the transfer learning across different values
of $x_{i}$. Consider the CPS dataset, we examine the relationship between the
covariate distance $\|x_{i}-x_{j}\|$ versus the empirical Wasserstein distance
between the conditional random variables $Y(x_{i})$ and $Y(x_{j})$. For better
accuracy, we select only pairs $(x_{i},x_{j})$ such that both $N_{i}$ and
$N_{j}$ are greater than 18. The scatter plot in Figure 2 indicates a positive
correlation between the covariate distance and the Wasserstein distance of the
conditional distributions. We remark that we observed the same trend in other
datasets.

Figure 2: Positive correlation between the Wasserstein distance between
$Y(x_{i})$ and $Y(x_{j})$ and the covariate distance $\|x_{i}-x_{j}\|$. Only
covariates $x$ with more than 18 observations are selected.

From this empirical evidence, we can assume that if two covariates $x_{i}$ and
$x_{j}$ are similar, their respective conditional distributions of $Y(x_{i})$
and $Y(x_{j})$ should also be similar. We emphasize that relying on this
assumption to design the regularization term does not restrict the
generalization of our framework. This assumption, in general, holds for the
case where the data-generating process follows the form $Y=f(X)+\text{noise}$,
where $f$ is a continuous function. This is a classic setting in distribution
learning, a more general form of linear assumption. Moreover, classic
algorithms for supervised and unsupervised learning, such as kNN or
k-clustering, also require this assumption to work. We can, therefore, assume
a Lipschitz condition that there exists a finite positive number $L$ such that

$W_{2}^{2}(\mu_{x_{i}},\mu_{x_{j}})\leq L\|x_{i}-x_{j}\|\qquad\forall i,j,$

where $\mu_{x}$ is the data-generating conditional probability measure for
$Y(x)$. If the generative network $T_{\theta}$ matches sufficiently well the
conditional distributions, then it is reasonable to translate the above
Lipschitz assumption into a regularization term: let $\mathcal{E}$ be a set of
pairs $(x_{i},x_{j})$ for which we would like to impose Lipschitz conditions,
we can regularize using the sum of the pairwise entropic OT distance:

$\sum_{(x_{i},x_{j})\in\mathcal{E}}\mathrm{W}^{2}_{2,\varepsilon}\left(T_{\theta}(x_{i},\cdot)_{\\#}\mathbb{U},T_{\theta}(x_{j},\cdot)_{\\#}\mathbb{U}\right).$
(9)

The semi-dual form of the EOT (4) leads to

$\displaystyle\sum_{(x_{i},x_{j})\in\mathcal{E}}\leavevmode\nobreak\
\Big{(}\max_{v_{ij}(\cdot)\in\mathcal{C}}\int_{\mathbb{R}}v_{ij}^{\varepsilon}(y)\mathrm{d}\mathbb{P}_{\theta,j}+\int_{\mathbb{R}}v_{ij}(y)\mathrm{d}\mathbb{P}_{\theta,i}-\varepsilon\Big{)}.$

Above, we have used the shorthand
$\mathbb{P}_{\theta,i}=T_{\theta}(x_{i},\cdot)_{\\#}\mathbb{U}$, the
pushforward measure of the uniform measure under the map
$T_{\theta}(x_{i},\cdot)$. Each maximization problem inside the sum is an
infinite-dimensional optimization problem, in which we need to search over all
potentials $v_{ij}$ continuous. To alleviate this computational burden, a
naive approach is to use a neural network $v_{\phi}(x_{i},x_{j},\cdot)\approx
v_{ij}(\cdot)$, where $\phi$ is a learnable parameter. By omitting the
constant $\varepsilon$ and switching the order of the summation and the
maximization operator, we can regularize using

$\displaystyle\max_{\phi}\\!\\!\\!\sum_{(x_{i},x_{j})\in\mathcal{E}}\\!\\!\\!\left(\int_{\mathbb{R}}v_{\phi}^{\varepsilon}(x_{i},x_{j},y)\mathrm{d}\mathbb{P}_{\theta,j}\\!+\\!\int_{\mathbb{R}}v_{\phi}(x_{i},x_{j},y)\mathrm{d}\mathbb{P}_{\theta,i}\right).$

A downside of this parametrization is that it requires both $x_{i}$ and
$x_{j}$ as input to $v_{\phi}$ to output the potential. This leads to
redundancy in the parametrization. To combat this problem, we propose a single
input network $v_{\phi}(x_{i},y)$ to approximate $v_{ij}(y)$, the potential
associating with $\mathbb{P}_{\theta,i}$. Further, we compute the smoothed
transform $v_{\phi}^{\varepsilon}(x_{i},y)$

$\displaystyle
v_{\phi}^{\varepsilon}(x_{i},y)\stackrel{{\scriptstyle\text{def.}}}{{=}}$
$\displaystyle-\varepsilon\log\left(\int_{\mathbb{R}}\exp\left(\frac{v_{\phi}(x_{i},y^{\prime})-\lvert
y-y^{\prime}\rvert^{2})}{\varepsilon}\right)\mathrm{d}\mathbb{P}_{\theta,i}\right),$

and use $v_{\phi}^{\varepsilon}(x_{i},y)$ to approximate
$v_{ij}^{\varepsilon}(y)$, the potential associating with
$\mathbb{P}_{\theta,j}$. The above single input parametrization no longer
takes $x_{j}$ as input; thus, to make the network well-defined, we need
disambiguation by imposing the following condition: for each distinct value
$x_{i}$, there is at most one pair $(x_{i},x_{j})$ that belongs to
$\mathcal{E}$. If this condition on $\mathcal{E}$ holds, then we can uniquely
identify the $x_{j}$ such that $v_{\phi}^{\varepsilon}(x_{i},y)$ is serving as
the potential.

Thus, for a valid set $\mathcal{E}$, we propose the following regularization
term:

$\mathrm{Reg}(\theta)\stackrel{{\scriptstyle\text{def.}}}{{=}}\max_{\phi}\leavevmode\nobreak\
\mathcal{R}(\theta,\phi),$ (10a) where $\mathcal{R}(\theta,\phi)$ is
$\displaystyle\mathcal{R}(\theta,\phi)\stackrel{{\scriptstyle\text{def.}}}{{=}}$
(10b)
$\displaystyle\quad\sum_{(x_{i},x_{j})\in\mathcal{E}}\leavevmode\nobreak\
\Big{(}\int_{\mathbb{R}}v_{\phi}^{\varepsilon}(x_{i},y)\mathrm{d}\mathbb{P}_{\theta,j}+\int_{\mathbb{R}}v_{\phi}(x_{i},y)\mathrm{d}\mathbb{P}_{\theta,i}\Big{)}.$

In this way, $\mathrm{Reg}(\theta)$ is the sum of the semi-dual neural
entropic OT distances of all pairwise distributions specified by
$\mathcal{E}$.

To conclude this section, we elaborate on a reasonable approach to building
$\mathcal{E}$. Ideally, $\mathcal{E}$ should capture neighborhood information
about which pair $(x_{i},x_{j})$ is close to each other; further, it should
satisfy our disambiguation criterion that for each distinct value $x_{i}$,
there is at most one pair $(x_{i},x_{j})$ that belongs to $\mathcal{E}$. A
reasonable choice of $\mathcal{E}$ that is also computationally efficient is
to build $\mathcal{E}$ as the directed edge set of a minimum spanning tree of
all $x$. Here, the spanning tree minimizes the sum of distance
$\|x_{i}-x_{j}\|$ between two connected nodes, and we can employ Kruskal’s
algorithm (Kruskal, 1956). For the direction of the edge: at each iteration,
when we add $x_{i}$ to the incumbent minimum spanning tree, then $x_{i}$ is
chosen as the tail of the edge, and the other endpoint of the adding edge
(which already belongs to the incumbent tree) is chosen as the head of the
edge. This construction will satisfy the disambiguation criterion for the set
$\mathcal{E}$. We emphasize that the graph-based approach to identify which
$x_{j}$ is associated with $x_{i}$ in the regularization term is not
restricted only to a discrete set of covariates. This is because the graph-
based approach is only employed during the training when the number of
covariates is finite. In inference phases, the model can give the
distributional prediction to any covariates.

######  Remark 4.2.

In (8), we use the unregularized OT distance as a discrepancy function for the
fit term $\mathrm{Fit}(\theta)$ because therein, we are measuring the
discrepancy between a fixed target (a uniform distribution $U$) and a quantity
that depends on the generative network output $T_{\theta}(X,U)$ and the data-
generating CDFs $F_{\mu_{X}}$. This contrasts with the regularization term
$\mathrm{Reg}(\theta)$, where we need to compare two distributions that are
both outputs of the generative network $T_{\theta}$. Using the unregularized
OT distance for the $\mathrm{Reg}(\theta)$ will require evaluating the CDFs of
$T_{\theta}$. To alleviate this difficulty, we use the semi-dual form of EOT,
where we can parametrize a single conditional potential network to evaluate
all the pairwise EOT distances between covariate pairs in $\mathcal{E}$.

###  4.3 Min-Max Optimization with Smooth Gradient Descent-Ascent

By substituting the form of the regularization term (10) into (6), we obtain
the min-max optimization problem:

$\min_{\theta}\max_{\phi}\leavevmode\nobreak\
\mathrm{Fit}(\theta)+\lambda\mathcal{R}(\theta,\phi).$ (11)

Naively optimizing problem (11) will result in a common pitfall in
optimization of training GANs: gradient descent-ascent (GDA) on nonconvex-
nonconcave objective is, in general, hard to converge (Heusel et al., 2017;
Lin et al., 2020; Zheng et al., 2023). To alleviate this problem, we employ
the state-of-the-art first-order algorithm for GDA (Zheng et al., 2023) with a
convergence guarantee. First, for parameters $r_{1}\neq r_{2}$,
$r_{1},r_{2}>0$, we add two extra regularization (smoothing) terms with
auxiliary variables $p$ and $q$ into the loss function. This leads to the
augmented objective

$\displaystyle\mathcal{L}(\theta,\phi,p,q)$ (12)
$\displaystyle\stackrel{{\scriptstyle\text{def.}}}{{=}}\mathrm{Fit}(\theta)+\lambda\mathcal{R}(\theta,\phi)+\frac{r_{1}}{2}\|\theta-p\|_{2}^{2}-\frac{r_{2}}{2}\|\phi-q\|_{2}^{2}.$

Including the quadratic terms enhances the smoothness of the primal and dual
update, facilitating a more balanced trade-off between the gradient updates in
our nonconvex-nonconcave setting. In practice, we must evaluate this loss with
training samples, detailed in Algorithm 1. Our training procedure is
summarized in Algorithm 2. This doubly-smoothed GDA update is guaranteed to
converge to a stationary point (Zheng et al., 2023, Theorem 2) as long as the
max term $\mathcal{R}(\theta,\phi)$ satisfies Kurdyka-Łojasiewicz condition (a
nonsmooth extension of the more popular Polyak-Łojasiewicz condition).

In general, we resort to a validation procedure to choose the appropriate
parameters $h,\lambda,\varepsilon,r_{1},r_{2}$, and the optimization
algorithm’s learning rates $(\alpha,\beta,\gamma,\delta)$. Details on the
parameters grid are given in the Appendix A.

Algorithm 1 Empirical loss evaluation

1: Input: data $(x_{i},\hat{F}_{x_{i}})$, network $T_{\theta}$ and $v_{\phi}$,
auxiliary variables $p,q$, batch size $B$, set of pair covariates
$\mathcal{E}$, regularization parameter $\varepsilon,\lambda,r_{1},r_{2}$

2: Output: loss estimate $\widehat{\mathcal{L}}(\theta,\phi,p,q)$

3: Initialize: take a batch of $B$ training nodes
$\\{x_{i}^{(b)}\\}_{b\in[B]}$

4: for each $(x_{i}^{(b)},x_{j}^{(b)})\in\mathcal{E}$ do

5: Sample $U^{(b)}\sim\mathcal{U}[0,1]$

6: for $m=1,\cdots,M$ do

7: Sample $U^{(b,m)}\sim\mathcal{U}[0,1]$

8: end for

9: end for

10:
$\widehat{\mathrm{Fit}}(\theta)\leftarrow\frac{1}{B}\sum_{b}\left[U^{(b)}-\hat{F}_{x_{i}^{(b)}}(T_{\theta}(x_{i}^{(b)},U^{(b)}))\right]^{2}_{2}$
(following (8))

11: Compute following (10):
$\hat{\mathcal{R}}(\theta,\phi)\leftarrow\frac{1}{M\times
B}\sum_{(x_{i}^{(b)},x_{j}^{(b)},U^{(b,m)})}\\\
\left[v_{\phi}^{\varepsilon}(x_{i}^{(b)},T_{\theta}(x_{j}^{(b)},U^{(b,m)}))+v_{\phi}(x_{i}^{(b)},T_{\theta}(x_{i}^{(b)},U^{(b,m)}))\right]$

12: Compute following (12):  $\widehat{\mathcal{L}}(\theta,\phi,p,q)\\\
\leftarrow\widehat{\mathrm{Fit}}(\theta)+\lambda\hat{\mathcal{R}}(\theta,\phi)+\frac{r_{1}}{2}\lVert{\theta-p}\rVert^{2}_{2}-\frac{r_{2}}{2}\lVert{\phi-q}\rVert_{2}^{2}$

Algorithm 2 Training of the Minimax loss (12)

1: Input: data $\mathcal{D}=(x_{i},\\{y_{ik}\\}_{k=1,\ldots,N_{i}})$, network
$T_{\theta}$ and $v_{\phi}$, learning rate $\alpha,\beta$, extrapolation
parameters $\gamma,\delta$

2: Initialize: $\theta^{0},\phi^{0},p^{0},q^{0}$

3: for each number of training iteration $t$ do

4:
$\theta^{t+1}\leftarrow\theta^{t}-\alpha\nabla_{\theta}\widehat{\mathcal{L}}(\theta^{t},\phi^{t},p^{t},q^{t})$

5:
$\phi^{t+1}\leftarrow\phi^{t}+\beta\nabla_{\phi}\widehat{\mathcal{L}}(\theta^{t+1},\phi^{t},p^{t},q^{t})$

6: $p^{t+1}\leftarrow p^{t}+\gamma(\theta^{t+1}-p^{t})$

7: $q^{t+1}\leftarrow q^{t}+\delta(\phi^{t+1}-q^{t})$

8: end for

##  5 Experiments

Datasets. To benchmark our proposed method GENTLE, we use the following two
datasets and the same preprocessing steps as in Athey et al. (2021).

The LDW-CPS dataset, constructed by LaLonde (1986); Dehejia & Wahba (1999), is
widely used in studies of Average Treatment Effects. This dataset contains
16177 samples from the real-world experiment and the Current Population
Survey. Each covariate $x_{i}$ includes the individual’s eight attributes: two
earnings measure (earning 74, earning 75), two indicators for ethnicity
(black, hispanic), marital status (married), age (age), two education measures
(education, nodegree), and a binary variable for whether this individual
receives treatment (is_treated). The response $y_{i}$ is earning 78, this
individual’s yearly income in 1978. This dataset has limited samples for each
covariate, shown in Figure 1. To achieve better testing evaluation and model
selection, we choose covariates based on their frequency for the test and
validation sets. Specifically, covariates with a frequency higher than 30 are
included in the test set, while those with a frequency higher than 20 but less
than or equal to 30 are included in the validation set. According to this
split strategy, the training set contains 15,255 samples corresponding to
13,745 covariates, the validation set contains 488 samples corresponding to 19
covariates, and the test set contains 434 samples corresponding to 20
covariates.

We also apply frameworks to a practical simulator called the Esophageal Cancer
Markov Chain Simulation Model (ECM).111Code is publicly available at
https://simopt.github.io/ECSim This simulator enables us to simulate patients’
quality-adjusted life years (QALYs) from their current treatment until death.
The approach to generate data for training and testing follows the authors in
Hong et al. (2023). Each covariate $x_{i}$ consists of five attributes which
are Barrett, aspirinEffect, statinEffect, drugIndex, and initialAge. The
response feature for this dataset is QALY (quality-adjusted life years).
Because we have access to the simulator, we can sample infinite responses for
each covariate. We follow the dataset generation in Hong et al. (2023) in
which they use the simulator to uniformly generate $2\times 10^{4}$ different
covariates for the train set, and each covariate has one response. The
validation set has $200$ covariates that are different from those in the train
set, while the test set also has $200$ covariates that are different from both
the train and validation sets. Each covariate has $10^{4}$ responses.

Baselines. We compare our method against state-of-the-art baselines, including
CWGAN (Athey et al., 2021), WGAN-GP (Gulrajani et al., 2017), MGAN (Baptista
et al., 2020), and CDSB (Shi et al., 2022). We keep the architecture of
$T_{\theta}(x,U)$ and $v_{\phi}$ simple with seven fully connected layers with
ReLU activation.

We denote $D_{\text{train}}$, $D_{\text{val}}$, and $D_{\text{test}}$ as the
set of covariates in the training, validation, and test sets, respectively.
While the covariates in the training set are used for training the model
parameters, the ones in the validation test are utilized to select
hyperparameters. For each distinct $x\in D_{\text{test}}$, the set of actual
respective responses is denoted as $Y_{GT}(x)$, while we use GENTLE, CWGAN,
WGAN-GP, MGAN, and CDSB to generate $K=10^{4}$ approximating points for
$Y_{\text{GT}}(x)$, denoted as $Y_{\text{GENTLE}}(x)$, $Y_{\text{CWGAN}}(x)$,
$Y_{\text{WGAN-GP}}(x)$, $Y_{\text{MGAN}}(x)$, $Y_{\text{CDSB}}(x)$
respectively. We then compute the average Wasserstein distance and Kolmogorov-
Smirnov statistic (Kolmogorov, 1933) over $D_{\text{test}}$ between
$Y_{\text{GT}}(x)$ and the outputs $Y_{\text{GENTLE}}(x)$,
$Y_{\text{CWGAN}}(x)$, $Y_{\text{WGAN-GP}}(x)$, $Y_{\text{MGAN}}(x)$, and
$Y_{\text{CDSB}}(x)$, as shown in Table 1.

Table 1: Wasserstein distance and K-S statistic on generated $\hat{Y}(x)$ using different methods vs. corresponding ground truth. The symbol $\pm$ represents the standard deviation. Dataset | Method |  WD($\downarrow$) |  KS($\downarrow$)
---|---|---|---
LDW-CPS |  | GENTLE
---
| 3767.62
---
($\pm$ 1337.56)
0.48 $\pm$ 0.05
CWGAN |  | 13371.99
---
($\pm$ 3637.66)
0.88 $\pm$ 0.09
MGAN |  | 4096.75
---
($\pm$ 1004.5)
0.91 $\pm$ 0.19
WGAN-GP |  | 7828.66
---
($\pm$ 2306.00)
0.99 $\pm$ 0.02
CDSB |  | 4970.66
---
($\pm$ 1094.74)
0.56 $\pm$ 0.14
ECM | GENTLE | 0.97 $\pm$ 0.14 | 0.21 $\pm$ 0.02
| CWGAN |  4.42 $\pm$ 2.01 |  0.94 $\pm$ 0.05
| MGAN |  1.09 $\pm$ 0.59 |  0.55 $\pm$ 0.13
| WGAN-GP |  1.59 $\pm$ 0.77 |  0.66 $\pm$ 0.27
| CDSB |  2.04 $\pm$ 0.43 |  0.47 $\pm$ 0.12

Results. We plot the qualitative results of estimated densities for different
covariates $X$ in Figure 3 and 4. Overall, synthetic observations generated by
GENTLE exhibit a strong resemblance to the ground truth distributions in both
datasets, while the ones generated by baselines deviate significantly from the
data-generating distributions. For each covariate, the distributions generated
by CWGAN seem to collapse into the mean of the training data distribution,
while those of MGAN and WGAN-GP tend to shrink towards the mean of the
conditional data distribution of the covariate. Although CDSB does not suffer
from this collapse tendency, it still fails to match the data-generating
distributions accurately. This is consistent with the quantitative results in
Table 1, which show that the data generated by GENTLE have not only smaller
Wasserstein distance (WD) and Kolmogorov-Smirnov (KS) values but also smaller
standard deviations. This suggests GENTLE’s superior performance compared to
the other baselines.

###  5.1 Ablation Study on Training Loss Terms

We explore the impact of removing specific components from the loss function
(12) by considering two ablation versions. The first version omits the
smoothing term in the loss function $\mathcal{L}$ by setting $r_{1}=r_{2}=0$,
while the second version omits the regularization term by setting $\lambda=0$.
When we exclude the regularization term, the model generates uniformly
distributed responses in the LDW-CPS dataset, as illustrated in Figure 6. We
hypothesize that the absence of $\mathrm{Reg}(\theta)$ compromises the model’s
ability to enforce the Lipschitz condition between the proximity in the
covariates and responses space. This absence significantly hampers the
method’s performance in the dataset with an unbalanced observation of
covariates like LDW-CPS. This problem does not exist in the ECM dataset which
covariates are uniformly drawn from the simulation, so the performance of this
version in ECM just slightly decreases. On the other hand, removing the
smoothing term leads the model to get trapped in bad local optima, learning
only simple two-peak patterns in both two datasets, as depicted in Figure 7
and Figure 8. This becomes more apparent when examining Table 2, where the KS
statistic of this version surpasses that of the main method. This improvement
may be attributed to the model’s inclination to match high-density regions
exclusively, lacking overall generalization. In summary, the results in Table
2 highlight that both ablation versions exhibit significantly worse
performance than the official method, particularly regarding Wasserstein
distance. This underscores the importance of the regularization and smoothing
terms in maintaining the model’s robustness and enhancing its overall
performance.

Table 2: Wasserstein distance and K-S statistic on generated $\hat{Y}(x)$ using different ablation versions of GENTLE. The version GENTLE w/o SMO does not use the smoothing term while GENTLE w/o REG does not use the regularization term. The symbol $\pm$ represents the standard deviation. Dataset | Method |  WD($\downarrow$) |  KS($\downarrow$)
---|---|---|---
LDW-CPS | GENTLE |  | 3767.62
---
($\pm$ 1337.56)
0.48 $\pm$ 0.05
| GENTLE w/o
---
reg
| 13625.83
---
($\pm$ 562.10)
0.57 $\pm$ 0.04
| GENTLE w/o
---
smo
| 6075.02
---
($\pm$ 1340.75)
0.41 $\pm$ 0.04
ECM | GENTLE |  | 0.97 $\pm$ 0.14
---
0.21 $\pm$ 0.02
| GENTLE w/o
---
reg
|  1.13 $\pm$ 0.47
---
0.29 $\pm$ 0.04
| GENTLE w/o
---
smo
|  1.80 $\pm$ 1.28
---
0.30 $\pm$ 0.05

##  6 Conclusions

This paper introduces a neural generative model for conditional distributions,
whose parameters can be learned from limited data. Our method demonstrates
remarkable efficacy in generating distributions over response spaces from
given inputs. This is exemplified by successfully applying conditional dataset
synthesis tasks on LDW-CPS and ECM datasets, where our framework exhibits
superior performances.

An exciting extension of the current work could involve exploring the
application of our proposed method to dynamic or time-series data.
Incorporating temporal dependencies could enhance the model’s ability to
capture evolving conditional distributions over time, enabling more accurate
predictions in dynamic environments. Additionally, extending the setting of
our approach beyond one-dimensional responses could further broaden its scope
and impact. Finally, exploring interpretability techniques to elucidate the
underlying mechanisms driving the learned conditional distributions could
enhance understanding of the model’s predictions in real-world decision-making
scenarios.

Figure 3: The qualitative results of methods on the LDW-CPS dataset. The
density graph ‘GT’, ‘GENTLE’ (ours), ‘CWGAN’, ‘WGAN-GP’, ‘MGAN’, ‘CDSB’ is
constructed by applying kernel density estimate (KDE) on $Y_{\text{GT}}(x)$,
$Y_{\text{GENTLE}}(x)$, $Y_{\text{CWGAN}}(x)$, $Y_{\text{WGAN-GP}}(x)$,
$Y_{\text{MGAN}}(x)$, $Y_{\text{CDSB}}(x)$, respectively. Each subfigure is
for a different value $x$ on the test set. We observe that GENTLE can produce
the true distribution much more effectively than baselines. Figure 4: The
qualitative results of methods on the ECM dataset. The density graph ‘GT’,
‘GENTLE’ (ours), ‘CWGAN’, ‘WGAN-GP’, ‘MGAN’, ‘CDSB’ is constructed by applying
kernel density estimate (KDE) on $Y_{\text{GT}}(x)$, $Y_{\text{GENTLE}}(x)$,
$Y_{\text{CWGAN}}(x)$, $Y_{\text{WGAN-GP}}(x)$, $Y_{\text{MGAN}}(x)$,
$Y_{\text{CDSB}}(x)$, respectively. Each subfigure is for a different value
$x$ on the test set. We observe that GENTLE can produce the true distribution
much more effectively than baselines.

Acknowledgments. The work of Binh Nguyen is supported by Singapore’s Ministry
of Education grant A-0004595-00-00. Viet Anh Nguyen gratefully acknowledges
the generous support from the CUHK’s Improvement on Competitiveness in Hiring
New Faculties Funding Scheme and the CUHK’s Direct Grant Project Number
4055191.

## Impact Statement

Our paper presents an advancement in the field of machine learning and also
has potential societal implications. By improving the understanding and
capabilities of conditional distribution modeling, our framework can
contribute to various real-world applications. For instance, it can aid in
personalized treatment planning in healthcare by modeling conditional
distributions of patient responses to different therapies. In finance, it can
enhance risk assessment by predicting conditional distributions of asset
returns under varying market conditions. By enabling more accurate and nuanced
predictions in diverse domains, our research can positively impact decision-
making processes, resource allocation, and, ultimately, the well-being of
individuals and society.

## References

  * Arjovsky et al. (2017) Arjovsky, M., Chintala, S., and Bottou, L.  Wasserstein generative adversarial networks.  In _International Conference on Machine Learning_ , pp. 214–223. PMLR, 2017.
  * Athey & Imbens (2015) Athey, S. and Imbens, G. W.  Machine learning for estimating heterogeneous causal effects.  Technical report, 2015.
  * Athey et al. (2021) Athey, S., Imbens, G. W., Metzger, J., and Munro, E.  Using Wasserstein generative adversarial networks for the design of Monte Carlo simulations.  _Journal of Econometrics_ , 2021.
  * Austin (2015) Austin, T.  Exchangeable random measures.  In _Annales de l’IHP Probabilités et Statistiques_ , volume 51, pp. 842–861, 2015.
  * Baptista et al. (2020) Baptista, R., Hosseini, B., Kovachki, N. B., and Marzouk, Y.  Conditional sampling with monotone GANs: From generative models to likelihood-free inference.  _arXiv preprint arXiv:2006.06755_ , 2020.
  * Bobkov & Ledoux (2019) Bobkov, S. and Ledoux, M.  _One-Dimensional Empirical Measures, Order Statistics, and Kantorovich Transport distances_ , volume 261.  American Mathematical Society, 2019.
  * Branco et al. (2017) Branco, P., Torgo, L., and Ribeiro, R. P.  SMOGN: A pre-processing approach for imbalanced regression.  In _First International Workshop on Learning with Imbalanced Domains: Theory and Applications_ , pp. 36–50. PMLR, 2017.
  * Bunne et al. (2022) Bunne, C., Krause, A., and Cuturi, M.  Supervised training of conditional Monge maps.  _Advances in Neural Information Processing Systems_ , 35:6859–6872, 2022.
  * Bunne et al. (2023) Bunne, C., Stark, S. G., Gut, G., Del Castillo, J. S., Levesque, M., Lehmann, K.-V., Pelkmans, L., Krause, A., and Rätsch, G.  Learning single-cell perturbation responses using neural optimal transport.  _Nature Methods_ , 20(11):1759–1768, 2023.
  * Chawla et al. (2002) Chawla, N. V., Bowyer, K. W., Hall, L. O., and Kegelmeyer, W. P.  Smote: synthetic minority over-sampling technique.  _Journal of artificial intelligence research_ , 16:321–357, 2002.
  * Choi et al. (2023) Choi, J., Choi, J., and Kang, M.  Generative modeling through the semi-dual formulation of unbalanced optimal transport.  _arXiv preprint arXiv:2305.14777_ , 2023.
  * Cuturi (2013) Cuturi, M.  Sinkhorn distances: Lightspeed computation of optimal transport.  _Advances in Neural Information Processing Systems_ , 26, 2013.
  * Dehejia & Wahba (1999) Dehejia, R. H. and Wahba, S.  Causal effects in nonexperimental studies: Reevaluating the evaluation of training programs.  _Journal of the American statistical Association_ , 94(448):1053–1062, 1999.
  * Díaz (2020) Díaz, I.  Machine learning in the estimation of causal effects: Targeted minimum loss-based estimation and double/debiased machine learning.  _Biostatistics_ , 21(2):353–358, 2020.
  * Esteban-Pérez & Morales (2022) Esteban-Pérez, A. and Morales, J. M.  Distributionally robust stochastic programs with side information based on trimmings.  _Mathematical Programming_ , 195(1):1069–1105, 2022.
  * Genevay et al. (2016) Genevay, A., Cuturi, M., Peyré, G., and Bach, F.  Stochastic optimization for large-scale optimal transport.  In _Proceedings of the 30th International Conference on Neural Information Processing Systems_ , pp. 3440–3448, 2016.
  * Gulrajani et al. (2017) Gulrajani, I., Ahmed, F., Arjovsky, M., Dumoulin, V., and Courville, A. C.  Improved training of Wasserstein GANs.  _Advances in Neural Information Processing Systems_ , 30, 2017.
  * Gushchin et al. (2022) Gushchin, N., Kolesov, A., Korotin, A., Vetrov, D., and Burnaev, E.  Entropic neural optimal transport via diffusion processes.  _arXiv preprint arXiv:2211.01156_ , 2022.
  * Heusel et al. (2017) Heusel, M., Ramsauer, H., Unterthiner, T., Nessler, B., and Hochreiter, S.  GANs trained by a two time-scale update rule converge to a local Nash equilibrium.  _Advances in Neural Information Processing Systems_ , 30, 2017.
  * Hong et al. (2023) Hong, L. J., Hou, Y., Zhang, Q., and Zhang, X.  Learning to simulate: Generative metamodeling via quantile regression.  _arXiv preprint arXiv:2311.17797_ , 2023.
  * Jacob et al. (2019) Jacob, L., She, J., Almahairi, A., Rajeswar, S., and Courville, A.  W2GAN: Recovering an optimal transport map with a GAN, 2019.  URL https://openreview.net/forum?id=BJx9f305t7.
  * Jiang et al. (2020) Jiang, G., Hong, L. J., and Nelson, B. L.  Online risk monitoring using offline simulation.  _INFORMS Journal on Computing_ , 32(2):356–375, 2020.
  * Kolmogorov (1933) Kolmogorov, A.  Sulla determinazione empirica di una legge di distribuzione.  _Giorn Dell’inst Ital Degli Att_ , 4:89–91, 1933.
  * Korotin et al. (2022) Korotin, A., Selikhanovych, D., and Burnaev, E.  Neural optimal transport.  In _The Eleventh International Conference on Learning Representations_ , 2022.
  * Kruskal (1956) Kruskal, J. B.  On the shortest spanning subtree of a graph and the traveling salesman problem.  _Proceedings of the American Mathematical Society_ , 7(1):48–50, 1956.
  * Kullback (1968) Kullback, S.  Probability densities with given marginals.  _The Annals of Mathematical Statistics_ , 39(4):1236–1243, 1968.
  * LaLonde (1986) LaLonde, R. J.  Evaluating the econometric evaluations of training programs with experimental data.  _The American Economic Review_ , pp. 604–620, 1986.
  * Léonard (2014) Léonard, C.  A survey of the Schrödinger problem and some of its connections with optimal transport.  _Discrete and Continuous Dynamical Systems-Series A_ , 34(4):1533–1574, 2014.
  * Lin et al. (2020) Lin, T., Jin, C., and Jordan, M.  On gradient descent ascent for nonconvex-concave minimax problems.  In _International Conference on Machine Learning_ , pp. 6083–6093. PMLR, 2020.
  * Makkuva et al. (2020) Makkuva, A., Taghvaei, A., Oh, S., and Lee, J.  Optimal transport mapping via input convex neural networks.  In _International Conference on Machine Learning_ , pp. 6672–6681. PMLR, 2020.
  * Nguyen et al. (2020a) Nguyen, V. A., Zhang, F., Blanchet, J., Delage, E., and Ye, Y.  Distributionally robust local non-parametric conditional estimation.  _Advances in Neural Information Processing Systems_ , 33:15232–15242, 2020a.
  * Nguyen et al. (2020b) Nguyen, V. A., Zhang, X., Blanchet, J., and Georghiou, A.  Distributionally robust parametric maximum likelihood estimation.  In _Advances in Neural Information Processing Systems_ , volume 33, pp. 7922–7932, 2020b.
  * Nguyen et al. (2024) Nguyen, V. A., Zhang, F., Wang, S., Blanchet, J., Delage, E., and Ye, Y.  Robustifying conditional portfolio decisions via optimal transport.  _Operations Research_ , 2024.
  * Nietert et al. (2021) Nietert, S., Goldfeld, Z., and Kato, K.  Smooth $p$-Wasserstein distance: Structure, empirical approximation, and statistical applications.  In _International Conference on Machine Learning_ , pp. 8172–8183. PMLR, 2021.
  * Nutz (2021) Nutz, M.  Introduction to Entropic Optimal Transport.  _Lecture notes, Columbia University_ , 2021.
  * Parikh et al. (2022) Parikh, H., Varjao, C., Xu, L., and Tchetgen, E. T.  Validating causal inference methods.  In _International Conference on Machine Learning_ , pp. 17346–17358. PMLR, 2022.
  * Sadana et al. (2024) Sadana, U., Chenreddy, A., Delage, E., Forel, A., Frejinger, E., and Vidal, T.  A survey of contextual optimization methods for decision-making under uncertainty.  _European Journal of Operational Research_ , 2024.
  * Scott (1979) Scott, D. W.  On optimal and data-based histograms.  _Biometrika_ , 66(3):605–610, 1979.
  * Shen et al. (2021) Shen, H., Hong, L. J., and Zhang, X.  Ranking and selection with covariates for personalized decision making.  _INFORMS Journal on Computing_ , 33(4):1500–1519, 2021.
  * Shi et al. (2022) Shi, Y., De Bortoli, V., Deligiannidis, G., and Doucet, A.  Conditional simulation using diffusion Schrödinger bridges.  In _Uncertainty in Artificial Intelligence_ , pp. 1792–1802. PMLR, 2022.
  * Shrivastava & Tuzel (2019) Shrivastava, A. and Tuzel, O.  Learning conditional error model for simulated time-series data.  In _CVPR Workshops_ , pp. 91–94, 2019.
  * Silverman (2018) Silverman, B. W.  _Density Estimation for Statistics and Data Analysis_.  Routledge, 2018.
  * Torgo et al. (2013) Torgo, L., Ribeiro, R. P., Pfahringer, B., and Branco, P.  SMOTE for regression.  In _Portuguese Conference on Artificial Intelligence_ , pp. 378–389. Springer, 2013.
  * van Breugel et al. (2024) van Breugel, B., Seedat, N., Imrie, F., and van der Schaar, M.  Can you rely on your model evaluation? Improving model evaluation with synthetic test data.  _Advances in Neural Information Processing Systems_ , 36, 2024.
  * Villani (2009) Villani, C.  _Optimal Transport: Old and New_ , volume 338.  Springer, 2009.
  * Zheng et al. (2023) Zheng, T., Zhu, L., So, A. M.-C., Blanchet, J., and Li, J.  Universal gradient descent ascent method for nonconvex-nonconcave minimax optimization.  In _Advances in Neural Information Processing Systems_ , 2023.

##  Appendix A Additional Experimental Details

###  A.1 Parameter Tuning

We heuristically fix the learning rates $\alpha=\beta=0.001$ as a reasonable
value for almost every simple multilayer perception model. For other
parameters, we conducted a grid search on minimizing the WD metric of the
validation set. We found that the best combination for our parameters in the
LDW-CPS dataset is KDE bandwidth of $h=0.3$, $\varepsilon=1.0$, $\lambda=0.4$,
$r_{1}=3.0$, $r_{2}=2.0$, $\gamma=0.5$, $\delta=0.7$. Regarding the ECM
dataset, the best combination is $h=0.2$, $\varepsilon=1.0$, $\lambda=0.4$,
$r_{1}=3.0$, $r_{2}=2.0$, $\gamma=0.5$, $\delta=0.7$. We fix these parameter
values for all our experiments in the main paper.

###  A.2 Comparison between GENTLE and Baselines on Other Metrics

In this section, we report two additional metrics, including Mean Squared
Error (MSE) and $R^{2}$ score, to evaluate the performance of GENTLE and
baselines. Tables 3 and 4 represent the average results of algorithms across
ten different seed numbers in terms of these two metrics. We observe that
GENTLE outperforms baselines on both MSE and $R^{2}$ criteria. Although MGAN
has comparable results to GENTLE in terms of MSE and $R^{2}$, the qualitative
results shown in Figures 3 and 4 and the statistical metrics of WD and KS in
Tables 1 reveals its trend of collapsing to the mean of the ground truth
distribution.

Table 3: MSE and $R^{2}$ on generated $\hat{Y}(x)$ using different methods vs. corresponding ground truth on the LDW-CPS dataset. The symbol $\pm$ represents the standard deviation. Method | MSE ($\times 10^{6}$)($\downarrow$) |  $R^{2}$($\uparrow$)
---|---|---
GENTLE | 3.70 $\pm$ 0.73 | 0.92 $\pm$ 0.05
CWGAN | 108.43 $\pm$ 42.93 | -1.83 $\pm$ 0.13
MGAN | 3.82 $\pm$ 0.71 | 0.92 $\pm$ 0.09
WGAN-GP | 59.64 $\pm$ 9.65 | -0.59 $\pm$ 0.07
CDSB | 8.54 $\pm$ 1.21 | 0.77 $\pm$ 0.08
Table 4: MSE and $R^{2}$ on generated $\hat{Y}(x)$ using different methods vs. corresponding ground truth on the ECM dataset. The symbol $\pm$ represents the standard deviation. Method | MSE ($\downarrow$) |  $R^{2}$($\uparrow$)
---|---|---
GENTLE | 0.62 $\pm$ 0.32 | 0.96 $\pm$ 0.05
CWGAN | 24.30 $\pm$ 4.36 | 0.00 $\pm$ 0.13
MGAN | 1.24 $\pm$ 0.22 | 0.91 $\pm$ 0.10
WGAN-GP | 2.11 $\pm$ 0.33 | 0.86 $\pm$ 0.08
CDSB | 2.01 $\pm$ 0.39 | 0.89 $\pm$ 0.04

###  A.3 Additional Experimental Results

In Figure 5, we plot the value of $T_{\theta}(x,U)$ for an increasing grid of
$U\sim\mathcal{U}(0,1)$. We observe that the output of $T_{\theta}$ also
increased with $U$, empirically confirming that $T_{\theta}$ is an increasing
function, which conforms with the fact that $T_{\theta}(x,U)$ is an
approximation of an increasing function of $U$.

(a) GENTLE trained on LDW-CPS

(b) GENTLE trained on ECM

Figure 5: Empirical evidence for the monotonically increasing property of the
learned network $T_{\theta}(x,U)$ in the variable $U$. Figure 6: Qualitative
result of the GENTLE version without regularization term on the LDW-CPS
dataset. Each subfigure represents the results of a different covariate value
in the test set. The learned distribution of this version tends to span over
the range of the response term uniformly. Figure 7: Qualitative result of the
GENTLE version without smoothing term on the LDW-CPS dataset. Each subfigure
represents the results of a different covariate value in the test set. Figure
8: Qualitative result of the GENTLE version without smoothing term on the ECM
dataset. Each subfigure represents the results of a different covariate value
in the test set.


# Phish-Defence: Phishing Detection Using Deep Recurrent Neural Networks

Aman Rangapur, Tarun Kanakam and Dhanvanthini P Aman Rangapur is with the
School of Computer Science and Engineering, VIT-AP University, Andhra Pradesh,
India, 522237. Email<EMAIL_ADDRESS>Kanakam is with the School
of Computer Science and Engineering, VIT-AP University, Andhra Pradesh, India,
522237, Email<EMAIL_ADDRESS>Jubilson is with the School
of Computer Science and Engineering, VIT-AP University, Andhra Pradesh, India,
522237, Email<EMAIL_ADDRESS>
###### Abstract

In the growing world of the internet, the number of ways to obtain crucial
data such as passwords and login credentials, as well as sensitive personal
information has expanded. Page impersonation, often known as phishing, is one
method of obtaining such valuable information. Phishing is one of the most
straightforward forms of cyber-attack for hackers, as well as one of the
simplest for victims to fall for. It can also provide hackers with everything
they need to access their targets’ personal and corporate accounts. Such
websites do not provide a service but instead gather personal information from
users. In this paper, we achieved state-of-the-art accuracy in detecting
malicious URLs using recurrent neural networks. Unlike previous studies, which
looked at online content, URLs, and traffic numbers, we merely look at the
text in the URL, which makes it quicker and catches zero-day assaults. The
network has been optimized so that it may be utilized on tiny devices like
Mobiles, Raspberry Pi without sacrificing the inference time.

##  I Introduction

The rapid growth and adoption of new technologies, such as smart devices and
5G connectivity, has resulted in the rapid growth of internet-based services.
Some of the services are essential for day-to-day activities. The increasing
number of people who utilize web-based services demonstrates their popularity
over the years. The goal is to increase the availability and accessibility of
web-based services used frequently so that users do not face any
inconvenience[1, 2]. Nonetheless, the unlimited access and usage of web-based
services or technology on the internet offers potential for cyber-attacks
because of non-uniform cyberspace regulations or processes [3].

There has been incredible growth in the number of illegal activities with
increased malicious and destructive content in the networks over the past
several decades. According to statistics, financial transactions, online
gaming services, and social media are among the most popular web-based
services with a large user base. These services are highly prone to cyber
assaults. The way to protect the user from cyberattacks is called
Cybersecurity [4]. Cybersecurity is responsible for safeguarding, prohibiting,
and recovering all the assets that use the internet from cyber assaults.

Phishing is a social-engineering-based cyber-attack that aims to steal
personal data such as credit card information, login credentials, and
passwords. Most of the time, the hacker registers a counterfeit domain
address, meticulously creates a website that resembles the organization’s
actual website, then sends a mass email to thousands, instructing recipients
to click on a link on the fake website. These assaults primarily result in the
loss of crucial data, financial loss, intellectual property loss, and
deterioration of trust and national security [5, 6]. Several cybersecurity
professionals and researchers have proposed and built a variety of anti-
phishing techniques for identifying phishing websites [7]. Elmer Lastdanger
proposed in his study that “Phishing is a scalable act of deception whereby
impersonation is used to obtain information from a target” [8].

Machine learning is a field within computer science that differs from
traditional computing. In traditional computing, algorithms follow the rules
and get executed, whereas Machine Learning algorithms allow computers to train
on data and use statistical analysis to output the value. In addition, neural
networks and genetic algorithms make the code very efficient. Neural networks
are a means of machine learning, in which a computer learns to perform some
task by analyzing training data; it’s a process that mimics how the human
brain operates. A neural network generates the best possible result without
redesigning the output criteria.

The motivation of this research is to obtain state-of-the-art accuracy using
Recurrent Neural Networks, which is efficient in detecting malicious websites.
Furthermore, our model has been optimized and can be utilized on tiny devices
like Raspberry Pi without sacrificing the inference time.

##  II Background study

The constant development of web services and e-commerce platforms has
encouraged many phishers and criminals to develop new ways to exploit and
deceive novice users into sending their financial information [9, 10].

Phishing is by far the most common cybercrime as of 2022, with the FBI’s
Internet Crime Complaint Centre reporting more than twice as many phishing
cases than any other kind of cybercrime. In a web phishing attack using email,
the phisher deceives web users by developing a fake website to steal the
financial and personal information of the users. Web phishing attacks can be
achieved in many ways [10, 11]. Initially, the phisher makes a phishing
website, which looks like the original website in its appearance. Eventually,
many emails containing hyperlinks to the phishing website are continually sent
to users by the phisher, which requests validating or updating credentials and
financial information to deceive the victims [12]. Finally, the target victims
are redirected to the phishing website when the hyperlinks are clicked. When
the victim enters the information through the phishing website, the phisher
can fully control the victim’s financial information. As a result, many
financial and identity thefts can be executed after successful web phishing
attacks [10, 13, 11]. Online social networks, emails, messengers, blogs,
forums, mobile applications, and voice over IP (VoIP) are different
communication platforms. Phishing attacks may target these communication
systems. Figure 1 shows the life cycle of a web phishing attack.

Figure 1: Life cycle of a web phishing attack.

Phishing attempts against software-as-a-service (SaaS), webmail users were the
most common in 2021. During the COVID-19 pandemic, the number of phishing
attacks increased significantly. Approximately 70% of healthcare facilities
with less than 500 employees are targeted, indicating that these smaller
facilities have probably weaker security systems due to low-security budgets.
It’s also been observed that 75% of phishing websites use SSL (Secure Sockets
Layer) security, which means that a website that uses the SSL protocol is not
always legitimate. The countermeasures in tackling phishing attacks are
categorized as follows.

###  II-A Traditional methods

####  II-A1 Blacklisting and whitelisting

Blacklisting and whitelisting are phishing URL lists in which the specified
URL is compared to a preset phishing URL. The disadvantage of this methodology
is that the blacklist cannot cover all phishing sites because a newly created
phishing site takes a long time to get updated in the list. This delay in
adding the suspect site to the list might be enough for the attackers to
achieve their goals [1, 9]. Whitelisting methods are contrary to blacklisting
methods. It determines if a URL is legitimate by comparing it to the
“whitelist” Database. The “whitelist” Database mainly contains a list of well-
known legitimate URLs [10].

####  II-A2 Visual Similarity

By analyzing the appearance of sites, Visual Similarity techniques may
distinguish phishing sites from genuine ones. Attackers embed images, Flash,
ActiveX, and Java Applet instead of HTML text to avoid phishing detection. As
enrollment items appear on phishing site webpages, visual similarity-based
detection algorithms may immediately distinguish them. It classifies URLs
based on visual similarity strategies using HTML Document Object Model tree,
CSS similarity, visual perception, visual characteristics, pixel-based, and
hybrid approaches [14].

###  II-B Non-Traditional methods

####  II-B1 Machine Learning and Deep Learning

To create phishing detection models, traditional machine learning methods
(also known as shallow learning) such as Support Vector Machines (SVM),
Decision Trees (DT), and Random Forests (RF) can be used.

The Hybrid Ensemble Feature Selection (HEFS) presented by Chiew K. L. et al.
[15] is a feature selection structure for a machine learning-based phishing
detection system. In the first step of HEFS, an original Cumulative
Distribution Function gradient (CDF-g) calculation is used to generate
critical feature sub-sets, which are then used as an input to a data
perturbation collection to generate optional feature sub-sets. The subsequent
stage infers pattern features from the optional feature subsets by utilizing a
function disturbance set. Compared to SVM, Naive Bayes, C4.5, JRip, and PART
classifiers, HEFS performs better when coordinated with the Random Forest
classifier.

Deep Learning, which comes under machine learning, focuses on building deep
networks with various layer types such as convolutional, pooling, dropout, and
fully connected layers. Deep learning-based phishing detection models are
built using algorithms like Convolutional Neural Networks (CNN), Recurrent
Neural Networks (RNN), Autoencoders.

Wei W. et al. [16] suggested a deep neural network with convolutional layers
for identifying phishing sites based solely on the URL address content.
Despite previous research that looked at URLs, traffic statistics, and online
content, they only looked at the URL text. The network they demonstrated had
been appropriately updated, so it can now be used even on mobile phones
without affecting its performance. As a result, the approach is faster and can
detect zero-day threats.

####  II-B2 Heuristic based methods

In Heuristic-based methods, minor characteristics of a website are extracted
to determine if the URL is malicious or legitimate. In contrast to
blacklisting approaches, heuristic-based solutions may detect new phishing
websites on a regular basis[17]. Intelligent machine learning classifiers are
trained using certain phishing and authentic websites as training data and
then utilized to detect newly developed phishing websites accurately [11].

Rao R. S. et al. [18] developed a heuristic technique with TWSVM (twin support
vector machine) classifier to detect malicious enrolled phishing websites.
Moreover, websites are facilitated on arrangement servers. Their research has
focused on detecting phishing sites hosted on arrangement domains by
inspecting the sign-in page and the website’s main page. The hyperlink and URL
– based characteristics are used to identify maliciously enlisted phishing
websites.

####  II-B3 Hybrid Techniques

A combination of different approaches is used to create a better model in
terms of accuracy and precision. A DTOF-ANN (”Decision Tree and Optimal
Features based Artificial Neural Network”), which is a neural-network phishing
identification model based on a decision tree and ideal feature selection, was
proposed by Zhu E. et al. [19] to tackle over-fitting. The initial step was to
improve the traditional K-medoids clustering technique by selecting a
consistent number of introduction centers to eliminate copy dots in popular
datasets. In Addition to DTOF-ANN, optimal feature selection based on
innovative features, such as decision trees and neighborhood search
techniques, is meant to eliminate the ineffective and undesirable parts.
Finally, the optimal architecture of the neural network classifier is
constructed by modifying boundaries and training the neural network classifier
with the specified ideal features.

##  III Related Work

Researchers have studied phishing attacks repeatedly, but most of them were
not entirely perfect. Despite getting solutions and accuracy, it required
complicated calculations, which made it difficult to use and took a good
amount of inference time to detect.

Using a variety of characteristics, da Silva et al. [20] developed a phishing
prediction model. The suggested model assesses static characteristics like
Keywords and patterns in the phishing URL. A qualitative examination of these
characteristics that do not define elements such as relationships and
similarities between qualities was undertaken in addition to the quantitative
data in the research.

A neural network model based on decision trees and optimal feature selection
was suggested by Zhu et al. [19]. Using the classic K-medoids clustering
process, they created an incremental selection strategy to eliminate duplicate
points from public datasets. An optimum feature selection method was created
to exclude the negative and useless traits. Tan et al. [21] developed graph
theory-based anti-phishing techniques. The suggested technique’s initial stage
is deleting hyperlinks from the web page in question and bringing in relevant
local web pages.

For phishing detection applications, Zabihimayvan and Doran used Fuzzy Rough
Set (FRS) theory to pick significant characteristics from the dataset [12].
Rough Set (RS) theory is extended into the Fuzzy Rough Set (FRS) theory.
Websites A and B are both phishing websites. Their features A and B have the
same value, so RS is a method to find a decision boundary by calculating the
equality of each data point based on certain features and the same classes,
such as websites A and B being both phishing websites and their features A and
B having the same value. RS is appropriate for the original dataset used in
this paper, in which the features are used as a discrete value, i.e., they are
a collection of 1, 0, -1 elements. However, after the dataset executes the
nominalization process, the value of the feature is transferred to a
continuous number from 0 to 1, and the FRS strategy is applied.

In 2021, El-Rashidy proposed a new method for selecting characteristics for an
online phishing detection model [22]. There are two parts to the feature
selection process. The first phase determines the impact of each feature’s
absence by training the random forest model on a new dataset with one feature
removed. A feature queue graded from high to low accuracy is created after
removing each element in the loop. The model is then trained and tested by
starting with one feature, adding a new feature from the ranking feature list
each time to create a dataset, calculating the accuracy of each time, and
ultimately finding the feature subset with the most excellent accuracy. This
strategy is used to find the most useful feature subset. However, because each
new dataset must go through the algorithm training and testing procedure,
there is a lot of computational complexity and time involved.

Yang et al. [23] proposed a new technique for phishing detection based on an
inverted matrix online sequential over-learning machine that characterizes a
website using three characteristics. To reduce matrix inversion, they employed
the Sherman Morrison Woodbury equation. They introduced the online queue
extreme learning machine to modify the training model.

Subasi and Kremic [24] also developed a system for detecting phishing
websites. They classified websites as legitimate or phishing using a variety
of machine learning algorithms. They presented two types of learners (Adaptive
Boosting (AdaBoost) and Multiboost) to help with anti-phishing and attack
detection.

For phishing and Botnet attack detection, De La Torre Parra et al. [25]
suggested a cloud-based deep learning architecture. The model consists of two
basic security mechanisms Distributed Convolutional Neural Network (DCNN)
model that works together to detect phishing and application layer Distributed
Denial of Service (DDoS) attacks and a cloud-based temporary long short-term
memory (LSTM) network model hosted on the backend to detect botnet attacks and
feed CNN to detect phishing attacks.

The primary focus of Yi et al. [26] was developing a deep learning framework
for detecting phishing websites. The researchers created two types of web
phishing features original features and interactive features. These
characteristics are employed in a Deep Belief Networks-based detection model
(DBN). The DBN-based detection model produced promising results during
experiments utilizing natural IP streams.

Wei et al. [27] suggested a lightweight deep learning technique for detecting
fraudulent URLs, allowing for a real-time and energy-efficient phishing
detection system. They demonstrated that the suggested technology could
identify phishing in real time using an energy-efficient integrated single
board computer using website URLs. Table I depicts other existing machine
learning models for phishing detection.

In this paper, we have developed and assessed web phishing detection models
using Recurrent Neural Networks such as LSTM and GRU to achieve maximum
accuracy and precision without compromising inference time for detecting
malicious websites on small devices.

TABLE I: Existing Machine Learning Models for Phishing Detection

Model or Algorithm | Type | Dataset | Challenges | Limitations | Accuracy
---|---|---|---|---|---
Random Forest (RF_1) [28] | Single | ISCXURL-2016 |  | Achieved high
---
accuracy and low
response time
without relying on
third-party services
and using limited
features extracted
from a URL.
| Neither used multiple
---
different datasets to
train the model,
nor compared the results
or evaluated the
robustness of
the model.
99.57%
Adaboost [29] | Single |  | Websites (PhishTank,
---
MillerSmiles, Google
Search), size of the
dataset not
mentioned, each
instance has
30 features.
| The proposed model
---
used Weka 3.6,
Python, and
MATLAB.
| Didn’t compare the results
---
and evaluate the
robustness of
the model.
98.30%
Random Forest (RF_2) [30] | Single |  | Websites (PhishTank,
---
OpenPhish, Alexa,
online payment gateway)
5223 instances, 2500
phishing URLs,
2723 legitimate URLs.
Has 20 features.
| 20 features are manually
---
extracted, some features
need to be obtained
by calling a third-party
service, and some
features need to
parse the website’s
HTML source code.
| Did not use multiple
---
different datasets to train
the model, compare the
results, or to evaluate
the robustness of the
model. The experimental
dataset is small.
99.50%
PSL1 \+ PART [31] | Hybrid |  | Websites (PhishTank,
---
Relbank)
30,500 original
instances,
20,500 phishing
URLs,
10,000 legitimate
URLs. Has 18 features.
| Extracted 3000
---
comprehensive
features and applied
different set of
parameters to ML
models to compare
the experimental
results.
| The legitimate URLs
---
in the dataset are all
related to banks, and
some features are
limited to e-banking
websites.
99.30%
| Recurrent Neural
---
Network (RNN) +
Convolutional
Neural Network
(CNN) [32]
Deep learning |  | Websites (PhishTank,
---
Alexa, etc),
490,408 instances,
245,385 phishing URLs,
245,023 legitimate URLs.
| Dataset is large-scale.
---
The first one to use
deep learning
model to detect
malicious URLs.
| The maximum length
---
of the URL is
255 characters.
The phishing website
URLs did not have
relevant semantics.
95.79%
| Random forest +
---
Neural network +
bagging [33].
Hybrid | UCI |  | No previous research
---
focuses on using a
feedforward NN
ensemble learning.
| Insufficient data
---
sources and lack of a
feature extraction
process.
97.40%
| Auto encoder +
---
NIOSELM [34]
Hybrid |  | Websites (PhishTank,
---
Alexa, DMOZ),
60,000 legitimate URLs,
5000 phishing URLs.
Has 56 features.
| The dataset is
---
imbalanced.
| The detection
---
accuracy may not be
the best compared
with the existing
methods.
94.60%

##  IV Methodology

Long Short-Term Memory [35, 36] was introduced to address the long-term
dependency problem. LSTM has a hidden state in which an LSTM block would
output and an internal cell that maintains the information across a temporal
context. The cell stores the long-term information, and the LSTM can erase,
write and read information from that cell based on whatever a context defines.
In step t, there is a hidden state $h_{t}$ and a cell state $C_{t}$. The
selection of the information erased or written, or read is controlled by three
corresponding gates each for erasing, writing, and reading respectively. These
gate values are dynamic and are learned and computed based on the input at a
particular time step and the hidden state from the previous time step.

Figure 2: Vanilla RNN architecture.

Figure. 2 is a single neural network with inputs $X_{t-1}$, $X_{t}$, $X_{t+1}$
at $t-1$, $t$ and $t+1$ that has a $tanh$ activation function. The input at a
particular time step $h_{t}$, as well as $h_{t-1}$ (output of the previous RNN
block), which is provided as input, and $tanh$ is applied that gives us
$h_{t}$. It is given as output and input for the next RNN block at the next
time step. LSTM has a similar structure, and it is a chain of repeating
modules where the same block applies at every step. The structure of each
block contains four different layers that interact with each other and are not
recurrent layers.

TABLE II: Variables of LSTM Unit

Variables of LSTM Unit
---
$x_{t}\in\mathbb{R}^{d}$ | input vector to LSTM Unit
$f_{t}\in(0,1)^{h}$ | forget gate’s activation vector
$i_{t}\in(0,1)^{h}$ | input/update forget gate’s activation vector
$o_{t}\in(0,1)^{h}$ | output forget gate’s activation vector
$h_{t}\in(-1,1)^{h}$ | hidden state vector also known as output vector
$\tilde{c}_{t}\in(-1,1)^{h}$ | cell input activation vector
$c_{t}\mathbb{R}^{h}$ | cell state vector
$W\in\mathbb{R}^{h\times d},U\in\mathbb{R}^{h\times d}~{}and~{}b\in\mathbb{R}^{h}$ | weight matrices and bias vector parameters

Figure 3: LSTM architecture.

The sigmoid function is used in these three blocks because the output should
lie between 0 & 1\. The cell state ($C_{t}$) can remove and add information to
the cell state, which is regulated by gates. The following points show the
working of LSTM.

  * •

Forget gate layer uses a sigmoid function to decide the information to be
thrown away from the cell state.

$f_{t}=\sigma(U_{f}h_{t-1}+x_{t}W_{f}+b_{f})$ (1)

  * •

Input gate layer(sigmoid function) decides the values to be updated, and the
$tanh$ function creates the new vector value, which is added to the state.

$i_{t}=\sigma(x_{t}W_{i}+h_{t-1}U_{i}+b_{i})$ (2)
$\tilde{C}_{t}=tanh(x_{t}W_{C}+h_{t-1}U_{c}+b_{c})$ (3)

  * •

Multiply the old cell with the cell which contains information that needs to
be dropped, and sum it with the dot product of the result from the previous
step.

$C_{t}=\sigma(f_{t}\times C_{t-1}+i_{t}\times\tilde{C}_{t})$ (4)

  * •

A filtered version of output is formed in the cell state; run the sigmoid
layer to decide the part of the cell states for output. Apply $tanh$ activator
to the cell state and multiply it with the output of the sigmoid gate, where
the value lies between 0-1.

$o_{t}=\sigma(x_{t}W_{o}+h_{t-1}U_{o}+b_{o})$ (5)

$h_{t}=o_{t}\times tanh(C_{t})$ (6)

Gated Recurrent Unit [37] was proposed as an alternative to LSTM, which is
more straightforward; it combines forget and input gates into a single update
gate as well as merges the cell state and hidden state. Figure 4 shows GRU
architecture.

Figure 4: GRU architecture.

The update gate controls what parts of the hidden state are updated and
preserved.

$z_{t}=\sigma(x_{t}W^{z}+h_{t-1}U^{z}+b_{z})$ (7)

Reset gate controls which parts of the previous hidden state are used to
compute the fresh subject. New hidden state(reset gate) content selects useful
parts of the previous hidden state. It uses the current input to compute new
hidden content.

$\tilde{h}_{t}=tan(r_{t}\times h_{t-1}U+x_{t}W+b)$ (8)

The hidden state (update gate) controls what is kept from the previous hidden
state and what is updated to new hidden state content.

$h_{t}=(1-z_{t})\times\tilde{h}_{t}+z_{t}\times h_{t}-1$ (9)

If reset gate is set to all 1’s and update gate to all 0’s then in reset gate
$h_{t}$ will be $tanh(W*[h_{t-1}],X_{t})$ and in update gate, $h_{t}$ becomes,
$h_{t-1}$ respectively. Input and forget gates of LSTM are coupled by an
update gate in GRU; the reset gate in GRU is applied directly to the previous
hidden state.

TABLE III: Variables of GRU Unit

Variables of GRU Unit
---
$x_{t}$ | input vector
$h_{t}$ | output vector
$\tilde{h}_{t}$ | candidate activation vector
$z_{t}$ | update gate vector
$r_{t}$ | reset gate vector
$W,U~{}and~{}b$ | parameter matrices and vector

##  V Proposed Model

###  V-A Preparation of Dataset

For effective detection of malicious URLs, the dataset should contain recent
URLs which are malicious, for recognizing fresh features to train the model.
Attackers will change the production of phishing links by anti-phishing
regulations and procedures that have been released. Anti-phishing models and
algorithms must also be improved based on new phishing data. Furthermore, the
training dataset’s quantity and validity significantly impact the performance
of machine learning-based solutions. The performance of deep learning models
increases with the variety of content in the training dataset. Hence, it is
advised that phishing URLs and legitimate URLs should be extracted from data
repositories. The dataset used in this paper consists of numerous legitimate
and malicious URLs which are taken from PhishTank, OpenPhish, and Common
Crawl. It consists of 46839 instances, and we merely looked into the text in
the URL and extracted features to train the model. The dataset was split into
75% and 25% for training and testing, respectively.

###  V-B Network Architecture and Training Parameters

Figure 5: Comparison of Accuracy at different epochs. Figure 6: Comparison of
Accuracy at different epochs.

All the features captured from the URL are fed to the LSTM layer with an
orthogonal recurrent initializer. The embedded matrix is used as weights,
further combined with dense layers with a sigmoid activation function. A
dropout of 0.5 is added between the LSTM layer and dense layer. Similarly, the
features are captured from the URL and fed to the GRU layer with an orthogonal
recurrent initializer and sigmoid recurrent activation, combined with dense
layers with a softmax activation function. A dropout of 0.2 is added between
the dense layers. Figure 5 and 6 depict the network architecture of Phish-
Defence of LSTM and GRU, respectively. The trained model is converted to a
tflite model for producing faster inference time on small edge devices like
Raspberry Pi.

The model is trained for 40 epochs with a batch size of 500 using the $tanh$
activation function for the LSTM model and the sigmoid activation function for
the GRU model. Binary cross entropy loss with Adam optimizer is attached to
the dense layer for classifying legitimate and malicious URLs in both LSTM and
GRU models. Binary cross-entropy loss is defined as

$loss=-\frac{1}{N}\sum_{i=1}^{N}y_{i}\bullet
log\left(p\left(y_{i}\right)\right)+\left(1-y_{i}\right)\bullet
log\left(1-p\left(y_{i}\right)\right)$ (10)

##  VI Results

As described in Section V-B the network architectures were trained with a
batch size of 500 with 70 training steps for 40 epochs, and the training and
testing loss started oscillating. The Adam optimization method was utilized
with an initial learning rate of 1e-3. This learning rate was decreased by a
learning rate scheduler, which checks for lack of reduction in training loss
for five epochs and decreases the learning rate by a factor of 0.1. The
reduction was stopped when the learning rate was reduced to an absolute
minimum of 1e-5. Early stopping was set at 6 epochs, which ends the training
process when the training loss oscillates for more than 6 epochs, yielding a
maximum validation accuracy of 90.50% for the LSTM model.

TABLE IV: Results of Phish-Defense

Model | Precision | Inference time | F-Score | Recall
---|---|---|---|---
Random Forest (RF_1) [28] | 98.30 | 2.30 | 99.01 | 99.65
Adaboost [29] | 99.00 | 2.80 | 98.80 | 98.60
Random Forest (RF_2) [30] | 99.40 | 1.95 | 99.40 | -
PSL 1 \+ PART [31] | 98.70 | 3.10 | 99.00 | 99.30
| Recurrent Neural
---
Network (RNN) +
Convolutional
Neural Network
(CNN) [32]
97.33 | 2.03 | 95.52 | 93.78
| Random forest +
---
Neural network +
bagging [33]
96.00 | 1.90 | 97.00 | 98.10
| Auto encoder +
---
NIOSELM [34]
97.80 | 2.20 | 95.52 | 98.30
PD-LSTM | 91.02 | 0.60 | 96.40 | 97.54
PD-GRU | 99.08 | 0.53 | 99.24 | 98.97

Figure 7: Accuracy vs Epochs for GRU model. Figure 8: Loss vs Epochs for GRU
model. Figure 9: Comparison of accuracy with existing models. Figure 10:
Comparison of inference time with existing models.

Later, we built and trained a fresh RNN model using GRU and achieved a state-
of-the-art accuracy of 98.51%. Figure 7 and 8 shows the accuracy vs epochs and
loss vs epochs graph for Gated Recurrent Unit model. This accuracy was
achieved using the softmax activation function with a 0.25 dropout. After
performing a significance test, the resulting trained model was found to
perform with an inference time of 0.8 seconds on Raspberry Pi-4. Figure 9 and
10 illustrate the accuracy and inference time in the form of a line graph.
Table IV shows the results of Phish-Defense compared to other machine learning
algorithms.

##  VII Conclusion

This study recommended a web phishing detection approach using deep recurrent
neural networks to predict phishing websites. In phishing site prediction, the
most influential features and optimal weights of website features were
tokenized, and these features were used to train the RNN to make better
predictions of phishing websites. The experimental results demonstrated that
the proposed phishing websites prediction approaches based on RNN increased
the classification performance using fewer features. Additionally, the
proposed method approaches produced the competitive classification accuracy of
phishing websites among all other classifiers with all the feature selection
methods used in this study with the lowest inference time. Thus, the proposed
approach based on RNNs can also be used on smaller devices like Raspberry-Pi
and Arduino to successfully predict phishing websites to contribute and
provide more confidence for online commerce and business customers.

## Acknowledgments

This research is carried in the Artificial Intelligence and Robotics (AIR)
Research Centre, VIT-AP University. We also thank the management for
motivating and supporting AIR Research Centre, VIT-AP University in building
this project.

## References

  * [1] R. Mohammad, F. Thabtah, and T. Mccluskey, “Predicting phishing websites based on self-structuring neural network,” Neural Computing and Applications, vol. 25, pp. 443–458, 08 2013.
  * [2] G. Vrbančič, I. Fister, and V. Podgorelec, “Swarm intelligence approaches for parameter setting of deep learning neural network: Case study on phishing websites classification,” in Proceedings of the 8th International Conference on Web Intelligence, Mining and Semantics, WIMS ’18, (New York, NY, USA), Association for Computing Machinery, 2018.
  * [3] R. Verma and A. Das, “What’s in a url: Fast feature extraction and malicious url detection,” in Proceedings of the 3rd ACM on International Workshop on Security And Privacy Analytics, IWSPA ’17, (New York, NY, USA), p. 55–63, Association for Computing Machinery, 2017.
  * [4] A. Basit, M. Zafar, X. Liu, A. R. Javed, Z. Jalil, and K. Kifayat, “A comprehensive survey of ai-enabled phishing attacks detection techniques,” Telecommunication Systems, 01 2021.
  * [5] A. Aleroud and L. Zhou, “Phishing environments, techniques, and countermeasures: A survey,” Computers & Security, vol. 68, 04 2017.
  * [6] Z. Ramzan and C. Wueest, “Phishing attacks: Analyzing trends in 2006,” in CEAS 2007 - The Fourth Conference on Email and Anti-Spam, 2-3 August 2007, Mountain View, California, USA, 2007.
  * [7] A. Ubing, S. Kamilia, A. Abdullah, N. Zaman, and M. Supramaniam, “Phishing website detection: An improved accuracy through feature selection and ensemble learning,” International Journal of Advanced Computer Science and Applications, vol. 10, 01 2019.
  * [8] E. Lastdrager, “Achieving a consensual definition of phishing based on a systematic review of the literature,” Crime Science, vol. 3, p. 9, 09 2014\.
  * [9] M. He, S.-J. Horng, P. Fan, K. Khan, R.-S. Run, J.-L. Lai, R.-J. Chen, and A. Sutanto, “An efficient phishing webpage detector,” Expert Syst. Appl., vol. 38, pp. 12018–12027, 09 2011.
  * [10] R. Mohammad, F. Thabtah, and T. Mccluskey, “Tutorial and critical analysis of phishing websites methods,” Computer Science Review, vol. 17, 05 2015.
  * [11] W. Ali, “Phishing website detection based on supervised machine learning with wrapper features selection,” International Journal of Advanced Computer Science and Applications, vol. 8, 01 2017.
  * [12] M. Zabihimayvan and D. Doran, “Fuzzy rough set feature selection to enhance phishing attack detection,” in 2019 IEEE International Conference on Fuzzy Systems (FUZZ-IEEE), pp. 1–6, 2019.
  * [13] N. Abdelhamid, A. Ayesh, and F. Thabtah, “Phishing detection based associative classification data mining,” Expert Systems with Applications, vol. 41, p. 5948–5959, 10 2014.
  * [14] A. Jain and B. B. Gupta, “Phishing detection: Analysis of visual similarity based approaches,” Security and Communication Networks, vol. 2017, pp. 1–20, 01 2017.
  * [15] K. L. Chiew, C. L. Tan, K. Wong, K. Yong, and W. Tiong, “A new hybrid ensemble feature selection framework for machine learning-based phishing detection system,” Information Sciences, vol. 484, pp. 153–166, 05 2019.
  * [16] W. Wei, Q. Ke, J. Nowak, M. Korytkowski, R. Scherer, and M. Woźniak, “Accurate and fast url phishing detector: A convolutional neural network approach,” Computer Networks, vol. 178, p. 107275, 04 2020.
  * [17] R. Mohammad, F. Thabtah, and T. Mccluskey, “Predicting phishing websites based on self-structuring neural network,” Neural Computing and Applications, vol. 25, pp. 443–458, 08 2013.
  * [18] R. Rao, A. Pais, and P. Anand, “A heuristic technique to detect phishing websites using twsvm classifier,” Neural Computing and Applications, vol. 33, pp. 1–20, 06 2021.
  * [19] E. Zhu, Y. Ju, Z. Chen, F. Liu, and X. Fang, “Dtof-ann: An artificial neural network phishing detection model based on decision tree and optimal features,” Applied Soft Computing, vol. 95, p. 106505, 06 2020.
  * [20] E. Feitosa, C. Silva, and V. Garcia, “Heuristic-based strategy for phishing prediction: A survey of url-based approach,” Computers & Security, 09 2019\.
  * [21] C. L. Tan, K. L. Chiew, K. Yong, S. Sze, J. Abdullah, and Y. Sebastian, “A graph-theoretic approach for the detection of phishing webpages,” Computers & Security, vol. 95, p. 101793, 05 2020.
  * [22] M. El-Rashidy, “A smart model for web phishing detection based on new proposed feature selection technique,” Menoufia Journal of Electronic Engineering Research, vol. 30, pp. 97–104, 01 2021.
  * [23] L. Yang, J. Zhang, X. Wang, Z. Li, Z. Li, and Y. He, “An improved elm-based and data preprocessing integrated approach for phishing detection considering comprehensive features,” Expert Systems with Applications, vol. 165, p. 113863, 08 2020.
  * [24] A. Subasi and E. Kremic, “Comparison of adaboost with multiboosting for phishing website detection,” Procedia Computer Science, vol. 168, pp. 272–278, 01 2020.
  * [25] G. Parra, P. Rad, K.-K. R. Choo, and N. Beebe, “Detecting internet of things attacks using distributed deep learning,” Journal of Network and Computer Applications, vol. 163, p. 102662, 04 2020.
  * [26] P. yi, Y. Guan, F. Zou, Y. Yao, W. Wang, and T. Zhu, “Web phishing detection using a deep learning framework,” Wireless Communications and Mobile Computing, vol. 2018, pp. 1–9, 09 2018.
  * [27] B. Wei, R. Hamad, L. Yang, X. He, H. Wang, B. Gao, and W. L. Woo, “A deep-learning-driven light-weight phishing detection sensor,” Sensors, vol. 19, p. 4258, 09 2019.
  * [28] B. B. Gupta, K. Yadav, I. Razzak, K. Psannis, A. Castiglione, and X. Chang, “A novel approach for phishing urls detection using lexical based machine learning in a real-time environment,” Computer Communications, vol. 175, pp. 47–57, 2021.
  * [29] A. Odeh, I. Keshta, and E. Abdelfattah, “Phiboost- a novel phishing detection model using adaptive boosting approach,” Jordanian Journal of Computers and Information Technology (JJCIT), vol. 07, pp. 65–74, 03 2021.
  * [30] E. Gandotra and D. Gupta, “Improving spoofed website detection using machine learning,” Cybernetics and Systems, vol. 52, no. 2, pp. 169–190, 2021\.
  * [31] P. Barraclough, G. Fehringer, and J. Woodward, “Intelligent cyber-phishing detection for online,” Computers & Security, vol. 104, p. 102123, 2021\.
  * [32] W. Wang, F. Zhang, X. Luo, and S. Zhang, “Pdrcnn: Precise phishing detection with recurrent convolutional neural networks,” Security and Communication Networks, vol. 2019, pp. 1–15, 10 2019.
  * [33] A. Zamir, H. Khan, T. Iqbal, N. Yousaf, F. Aslam, A. Anjum, and M. Hamdani, “Phishing web site detection using diverse machine learning algorithms,” The Electronic Library, vol. ahead-of-print, 01 2020.
  * [34] L. Yang, J. Zhang, X. Wang, Z. Li, Z. Li, and Y. He, “An improved elm-based and data preprocessing integrated approach for phishing detection considering comprehensive features,” Expert Syst. Appl., vol. 165, p. 113863, 2021\.
  * [35] H. Sak, A. Senior, and F. Beaufays, “Long short-term memory recurrent neural network architectures for large scale acoustic modeling,” Proceedings of the Annual Conference of the International Speech Communication Association, INTERSPEECH, pp. 338–342, 01 2014.
  * [36] S. Hochreiter and J. Schmidhuber, “Long short-term memory,” Neural computation, vol. 9, pp. 1735–80, 12 1997.
  * [37] R. Dey and F. M. Salem, “Gate-variants of gated recurrent unit (gru) neural networks,” in 2017 IEEE 60th International Midwest Symposium on Circuits and Systems (MWSCAS), pp. 1597–1600, 2017.