# Generative Modeling with Flow-Guided Density Ratio Learning

Alvin Heng1,*, Abdul Fatir Ansari1,2,$\dagger$, Harold Soh1,3
1Dept. of Computer Science, National University of Singapore
2AWS AI Labs
3Smart Systems Institute, National University of Singapore

###### Abstract

We present Flow-Guided Density Ratio Learning (FDRL), a simple and scalable
approach to generative modeling which builds on the stale (time-independent)
approximation of the gradient flow of entropy-regularized $f$-divergences
introduced in DG$f$low. In DG$f$low, the intractable time-dependent density
ratio is approximated by a stale estimator given by a GAN discriminator. This
is sufficient in the case of sample refinement, where the source and target
distributions of the flow are close to each other. However, this assumption is
invalid for generation and a naive application of the stale estimator fails
due to the large chasm between the two distributions. FDRL proposes to train a
density ratio estimator such that it learns from progressively improving
samples during the training process. We show that this simple method
alleviates the density chasm problem, allowing FDRL to generate images of
dimensions as high as $128\times 128$, as well as outperform existing gradient
flow baselines on quantitative benchmarks. We also show the flexibility of
FDRL with two use cases. First, unconditional FDRL can be easily composed with
external classifiers to perform class-conditional generation. Second, FDRL can
be directly applied to unpaired image-to-image translation with no
modifications needed to the framework. Code is publicly available at
https://github.com/ajrheng/FDRL.

11footnotetext:  Correspondence at<EMAIL_ADDRESS>Work done while at the National University of Singapore, prior to joining
Amazon.

##  1 Introduction

Figure 1: Left: Illustration of FDRL’s training setup for the $\tau$ training
iteration. For clarity, we emphasize the choices of the $f$-divergence and
$g$, the Bregman divergence function, as part of training. Right: the various
applications of FDRL, ranging from unconditional image generation to class-
conditional generation by composition with external classifiers and unpaired
image-to-image translation.

Deep generative modeling has made significant strides in recent years.
Advances in techniques such as generative adversarial networks [11] and
diffusion models [13, 32] have enabled the generation of high-quality,
photorealistic images and videos. Recently, gradient flows have emerged as an
interesting alternative approach to generative modeling. Unlike the
aforementioned techniques, gradient flows have a unique interpretation where
the objective is to find the path of steepest descent between two
distributions. Further, where diffusion models typically require a tractable
terminal distribution such as Gaussian, the source and target distributions in
gradient flows can be varied for different tasks in generative modeling. For
instance, the gradient flow performs data synthesis when the source
distribution is a simple prior and the target is the data distribution. When
the two distributions are from different domains of the same modality, the
gradient flow performs domain translation.

Among gradient flow methods, Wasserstein gradient flows (WGF) have recently
become a popular specialization for a variety of applications [20, 8, 1, 9].
WGFs model the gradient dynamics on the space of probability measures with
respect to the Wasserstein metric; these methods aim to construct the optimal
path between two probability measures — a source distribution
$q({\mathbf{x}})$ and a target distribution $p({\mathbf{x}})$, where the
notion of optimality refers to the path of steepest descent in Wasserstein
space.

Despite the immense potential of gradient flow methods in generative modeling
[1, 9, 8], applications to high-dimensional image synthesis remain limited,
and its applications to other domains, such as class-conditional generation
and image-to-image translation have been underexplored. A key difficulty is
that existing methods to solve for the gradient flow often involve complex
approximation schemes. For instance, the JKO scheme [14] necessitates
estimation of the 2-Wasserstein distance and divergence functionals, while
particle simulation methods, which involve simulating an equivalent stochastic
differential equation (SDE), requires solving for the time-dependent marginal
distribution over the flow.

A recent work on adopting gradient flows for sample refinement, DG$f$low [1],
proposes a stale (time-independent) estimate given by the discriminator of a
pretrained GAN, as the density ratio estimate necessary in the corresponding
SDE. While this dramatically simplifies the particle simulation approach, in
practice this only works when $q({\mathbf{x}})$ and $p({\mathbf{x}})$ are
close together, such as between the GAN generated images and the true data
distribution. When the two distributions are far apart, such as the case of
image synthesis, where $q({\mathbf{x}})$ is a simple prior distribution,
estimating the density ratio $q({\mathbf{x}})/p({\mathbf{x}})$ becomes a
trivial problem due to the large density chasm [26], and use of this naive
stale estimator in simulating the SDE fails to generate realistic images.

We therefore pose the question: _is it possible to extend the refinement
scheme of DG $f$flow to the more challenging task of sample generation?_ We
answer this in the affirmative by proposing Flow-Guided Density Ratio Learning
(FDRL), a new training approach to address the problem of estimating the
density ratio between the prior and the data distribution, by progressively
training a density ratio estimator with evolving samples. FDRL operates
exclusively in the data space, and does not require training additional
generators, unlike related particle methods [9]. Compared to gradient flow
baselines, we achieve the best quantitative scores in image synthesis and are
the first to scale to image dimensions as high as $128\times 128$. In
addition, we show that our framework can be seamlessly applied to two other
critical tasks in generative modeling: class-conditional generation and
unpaired image-to-image translation, which have been under-explored in the
context of gradient flow methods.

##  2 Background

We provide a brief overview of gradient flows and density ratio estimation.
For a more comprehensive introduction to gradient flows, please refer to [28].
A thorough overview of density ratio estimation can be found in [33].

#### Wasserstein Gradient Flows.

Consider Euclidean space equipped with the familiar $L_{2}$ distance metric
$(\mathcal{X},\|\cdot\|_{2})$. Given a function
$F:\mathcal{X}\rightarrow\mathbb{R}$, the curve
$\\{{\mathbf{x}}(t)\\}_{t\in\mathbb{R}^{+}}$ that follows the direction of
steepest descent is called the gradient flow of $F$:

${\mathbf{x}}^{\prime}(t)=-\nabla F({\mathbf{x}}(t)).$ (1)

In generative modeling, we are interested in sampling from the underlying
probability distribution of a given dataset. Hence, instead of Euclidean
space, we consider the space of probability measures with finite second
moments equipped with the 2-Wasserstein metric
$(\mathcal{P}_{2}(\Omega),\mathcal{W}_{2})$. Given a functional
$\mathcal{F}:\mathcal{P}_{2}(\Omega)\rightarrow\mathbb{R}$ in the
2-Wasserstein space, the gradient flow of $\mathcal{F}$ is the steepest
descent curve of $\mathcal{F}$. Such curves are called Wasserstein gradient
flows.

#### Density Ratio Estimation via Bregman Divergence.

Let $q({\mathbf{x}})$ and $p({\mathbf{x}})$ be two distributions over
$\mathcal{X}\in\mathbb{R}^{d}$ where we have access to i.i.d samples
${\mathbf{x}}_{q}\sim q({\mathbf{x}})$ and ${\mathbf{x}}_{p}\sim
p({\mathbf{x}})$. The goal of density ratio estimation (DRE) is to estimate
the true density ratio
$r^{*}({\mathbf{x}})=\frac{q({\mathbf{x}})}{p({\mathbf{x}})}$ based on samples
${\mathbf{x}}_{q}$ and ${\mathbf{x}}_{p}$.

We will focus on density ratio fitting under the Bregman divergence (BD),
which is a framework that unifies many existing DRE techniques [33, 34]. Let
$g:\mathbb{R}^{+}\rightarrow\mathbb{R}$ be a twice continuously differentiable
convex function with a bounded derivative. The BD seeks to quantify the
discrepancy between the estimated density ratio $r_{\theta}$ and the true
density ratio $r^{*}$:

$BD_{g}(r^{*}||r_{\theta})=\mathbb{E}_{p}[g(r^{*}({\mathbf{x}}))-g(r_{\theta}({\mathbf{x}}))+\partial
g(r_{\theta}({\mathbf{x}}))r_{\theta}({\mathbf{x}})]\\\
-\mathbb{E}_{q}[\partial g(r_{\theta}({\mathbf{x}}))].$ (2)

In practice, we estimate the expectations in Eq. 2 using Monte Carlo samples.
We can also drop the term $\mathbb{E}_{p}[g(r^{*}({\mathbf{x}}))]$ during
optimization as it does not depend on the model $r_{\theta}$. The minimizer of
Eq. 2, which we denote $\theta^{*}$, satifies
$r_{\theta^{*}}(x)=r^{*}({\mathbf{x}})=q({\mathbf{x}})/p({\mathbf{x}})$.
Moving forward, we will use the hatted symbol $\mathbb{\widehat{E}}_{p}$ to
refer to Monte Carlo estimates for ease of notation. We consider two common
instances of $g$ in this paper, namely $g(y)=\frac{1}{2}(y-1)^{2}$ and
$g(y)=y\log y-(1+y)\log(1+y)$, which correspond to the Least-Squares
Importance Fitting (LSIF) and Logistic Regression (LR) objectives. We provide
the full form of the BD for these choices of $g$ in Appendix C.

##  3 Discriminator Gradient Flow

Ansari et al. propose DG$f$low [1] for sample refinement by simulating the
SDE:

$d{\mathbf{x}}_{t}=-\nabla_{\mathbf{x}}f^{\prime}(q_{t}({\mathbf{x}}_{t})/p({\mathbf{x}}_{t}))dt+\sqrt{2\gamma}d{\mathbf{w}}_{t},$
(3)

where $d{\mathbf{w}}_{t}$ denotes the standard Wiener process.
$p({\mathbf{x}})$ represents the target distribution of refinement, while
$q_{t}({\mathbf{x}})$ represents the time-evolved source distribution
$q_{0}({\mathbf{x}})$. In theory, as $t\rightarrow\infty$, the stationary
distribution of the particles ${\mathbf{x}}_{t}$ from Eq. 3 converges to the
target distribution $p({\mathbf{x}})$. The SDE is the equivalent particle
system of the corresponding Fokker-Planck equation, which is the gradient flow
of the functional

$\mathcal{F}^{f}_{p}(q)=\int
p({\mathbf{x}})f(q({\mathbf{x}})/p({\mathbf{x}}))d{\mathbf{x}}+\gamma\int
q({\mathbf{x}})\log q({\mathbf{x}})d{\mathbf{x}},$ (4)

where $f:\mathbb{R}^{+}\rightarrow\mathbb{R}$ is a twice-differentiable convex
function with $f(1)=0$. Eq. 3 can therefore be understood as the flow which
minimizes Eq. 4. The first term in Eq. 4 is the $f$-divergence, which measures
the discrepancy between $q({\mathbf{x}})$ and $p({\mathbf{x}})$. Popular
$f$-divergences include the Kullback-Leibler (KL), Pearson-$\chi^{2}$
divergence and Jensen-Shannon (JS) divergence. Meanwhile, the second term is a
differential entropy term that improves expressiveness of the flow. We list
the explicit forms of $f$ that will be considered in this work in Table 2 in
the Appendix.

In DG$f$low, the source distribution is the distribution of samples from the
GAN generator, while the target distribution is the true data distribution.
The time-varying density ratio $q_{t}({\mathbf{x}})/p({\mathbf{x}})$ is
approximated by a stale estimate $r_{\theta}({\mathbf{x}})$

$d{\mathbf{x}}_{t}=-\nabla_{\mathbf{x}}f^{\prime}(r_{\theta}({\mathbf{x}}_{t}))dt+\sqrt{2\gamma}d{\mathbf{w}}_{t},$
(5)

where $r_{\theta}$ is represented by the GAN discriminator. The assumption
underlying this approach is that the distribution of generated samples from
the GAN is sufficiently close to the data distribution, so that the density
ratio $q_{t}({\mathbf{x}})/p({\mathbf{x}})$ is approximately constant
throughout the flow.

As we are unable to compute the marginal distributions $q_{t}({\mathbf{x}})$
or estimate the time-dependent density ratio
$q_{t}({\mathbf{x}})/p({\mathbf{x}})$ in Eq. 3 a priori, we are motivated to
extend the stale formulation of DG$f$low from refinement to the more difficult
problem of synthesis, i.e., generating samples from noise. We show that in
such a scenario, the stationary distribution of Eq. 5 shares the same MLE as
the data distribution.

######  Lemma 1.

Let $\gamma=1$ and assume that $q({\mathbf{x}})\sim U[a,b]$, where $a,b$ are
chosen appropriately (e.g., [-1,1] for image pixels). Then the stationary
distribution of Eq. 5, $\rho_{\infty}({\mathbf{x}})$, has the same maximum
likelihood estimate as $p({\mathbf{x}})$,
$\operatorname*{arg\,max}_{\mathbf{x}}\rho_{\infty}({\mathbf{x}})=\operatorname*{arg\,max}_{\mathbf{x}}p({\mathbf{x}})$.

Based on Lemma 1, we hypothesize that samples drawn from Eq. 5 have high
likelihoods under the data distribution. This motivates us to train a stale
estimator $r_{\theta}({\mathbf{x}})$ where $q({\mathbf{x}})$ is a simple prior
and $p({\mathbf{x}})$ is the true data distribution.

###  3.1 Where The Stale Estimator Breaks Down

We thus wish to modify the source distribution such that $q({\mathbf{x}})$ is
a simple prior, such as a Gaussian, while keeping $p({\mathbf{x}})$ to be the
true data distribution. Given samples from $q({\mathbf{x}})$ and
$p({\mathbf{x}})$, we can leverage the Bregman divergence to train a stale
estimator and simulate Eq. 5. However, we observe empirically that this fails
catastrophically as the density ratio estimation problem becomes trivial and
the Bregman loss diverges. We convey the key intuitions with a toy example of
simple 2D Gaussians.

Consider an SDE of the form

$d{\mathbf{x}}_{t}=-\nabla_{\mathbf{x}}f^{\prime}(q({\mathbf{x}}_{t})/p({\mathbf{x}}_{t}))dt+\sqrt{2\gamma}d{\mathbf{w}}_{t},$
(6)

where $q({\mathbf{x}})$ and $p({\mathbf{x}})$ are both simple 2D Gaussian
densities. This differs from Eq. 3 as the distribution $q({\mathbf{x}}_{t})$
is time-independent, and fixed as the source Gaussian distribution. We would
like to investigate how such an SDE, with $q/p$ represented by a stale
estimator $r_{\theta}$, behaves in transporting particles from
$q({\mathbf{x}})$ to $p({\mathbf{x}})$.

Although we have access to the analytical form of the density ratio
$q({\mathbf{x}})/p({\mathbf{x}})$ in this example, we choose to train a
multilayer perceptron (MLP) to estimate the density ratio using the Bregman
divergence Eq. 2. This is to evaluate the general case where we may only have
samples and not the analytic densities. We investigate two scenarios by
varying the distance between $q$ and $p$, as measured by their KL divergence
(or more generally, the $f$-divergences). $q({\mathbf{x}})$ is fixed at
$\mathcal{N}(\mathbf{0},0.1\mathbb{I})$, where bold numbers represent 2D
vectors. We set $p({\mathbf{x}})\sim\mathcal{N}(\mathbf{1},0.1\mathbb{I})$ in
the first case and $p({\mathbf{x}})\sim\mathcal{N}(\mathbf{6},0.1\mathbb{I})$
in the second.

We first train separate MLPs to estimate the density ratio, using
${\mathbf{x}}_{q}$ and ${\mathbf{x}}_{p}$ drawn from the respective
distributions. After convergence, we simulate a finite number of steps of Eq.
5 on samples drawn from $q({\mathbf{x}})$, with the density ratio given by the
MLP. We characterize the transition of particles through Eq. 5 as
parameterized by a generalized kernel

$\tilde{q}({\mathbf{x}})=\int
q({\mathbf{x}}^{\prime})M_{\theta}({\mathbf{x}}|{\mathbf{x}}^{\prime})d{\mathbf{x}}^{\prime},$
(7)

where $M_{\theta}$ is the transition kernel parameterized by the MLP. The
results are shown in Fig. 2. In (a), where the two distributions are close
together, the particles from $q({\mathbf{x}})$ are successfully transported to
$p({\mathbf{x}})$. Meanwhile in (b), even after flowing for many more steps,
the particles have not converged to the target distribution. Flowing for more
steps in (b) results in the particles drifting even further from the target
(and out of the image frame). This suggests that the stale estimator is only
appropriate in the case where the two distributions are not too far apart,
which indeed corresponds to the assumption of DG$f$low that the GAN generated
images are already close to the data distribution.

This observation can be attributed to two complementary explanations. Firstly,
the Bregman divergence is optimized using samples in a Monte Carlo fashion;
when the two distributions are separated by large regions of vanishing
likelihood, there are insufficient samples in these regions for the model to
learn an accurate estimate of the density ratio. Consequently, the inaccurate
estimate causes the particles to drift away from the correct direction of the
target as it crosses these regions, as seen in Fig. 2(b).

The second explanation may be attributed to the density chasm problem [26].
When two distributions are far apart, a binary classifier (which is equivalent
to a density ratio estimator, see appendix Sec. E) can obtain near perfect
accuracy while learning a relatively poor estimate of the density ratio. For
example, in the second case above, the neural network can trivially separate
the two modes with a multitude of straight lines, thus failing to estimate the
density ratio accurately.

(a)

(b)

Figure 2: Toy experiments by simulating the flow Eq. 5 with an MLP density
ratio estimator. The source distribution $q({\mathbf{x}})$ is set as
$\mathcal{N}(\mathbf{0},0.1\mathbb{I})$. Blue particles represent the source
particles, and orange particles represent the same particles after flowing for
$K$ steps. (a) We set
$p({\mathbf{x}})\sim\mathcal{N}(\mathbf{1},0.1\mathbb{I})$ and flow for $K=15$
steps. (b) We set $p({\mathbf{x}})\sim\mathcal{N}(\mathbf{6},0.1\mathbb{I})$
and flow for $K=400$ steps. We can clearly see that particles in (a) have
converged to the target distribution, while particles in (b) have not,
demonstrating the density chasm problem.

The eventual goal of this work is the generative modeling of images, where we
set our source distribution to be a simple prior distribution, and the target
distribution to be the distribution of natural images. In such cases involving
much higher-dimensional distributions, the density chasm problem will be
dramatically more pronounced, causing the stale estimator to fail to generate
realistic images. As mentioned previously, in early experiments where we
attempt to learn a density ratio estimator with $q({\mathbf{x}})$ being a
diagonal Gaussian or uniform prior and $p({\mathbf{x}})$ represented by
samples of natural images (e.g. CIFAR10), the Bregman loss diverges early
during training, leading to failed image samples.

###  3.2 Data-Dependent Priors

One approach to address the density chasm problem for images is to devise a
prior $q({\mathbf{x}})$ that is as close to the target distribution as
possible, while remaining easy to sample from. Inspired by generation from
seed distributions with robust classifiers [29], we first adopted a data-
dependent prior (DDP) by fitting a multivariate Gaussian to the training
dataset:

$q({\mathbf{x}})=\mathcal{N}(\bm{\mu}_{D},\bm{\Sigma}_{D}),\textrm{where}\\\
\bm{\mu}_{D}=\mathbb{E}_{D}[{\mathbf{x}}],\ \ \
\bm{\Sigma}_{D}=\mathbb{E}_{D}[({\mathbf{x}}-\bm{\mu}_{D}]^{T}({\mathbf{x}}-\bm{\mu}_{D})]$
(8)

where subscript $D$ represents the dataset. Samples from this distribution
contain low frequency features of the dataset, such as common colors and
outlines, which we visualize in Fig. 11 of the appendix. We proceeded to train
a CNN with the Bregman divergence in the same manner as the toy experiments,
with ${\mathbf{x}}_{q}$ drawn from the DDP and ${\mathbf{x}}_{p}$ from the
dataset. However, we found that the Bregman loss still diverged early during
training — the DDP alone is insufficient to bridge the density chasm. In the
next section, we describe a novel approach to cross the chasm using _flow-
guided_ training, which when used with DDP achieves significantly better
generative performance.

##  4 Flow-Guided Density Ratio Learning

We propose FDRL, a training scheme that is able to scale to the domain of
high-dimensional images, while retaining the simplicity of the stale
formulation of DG$f$low. Instead of merely training $r_{\theta}$ to estimate
$q({\mathbf{x}})/p({\mathbf{x}})$, our method proposes to estimate
$\tilde{q}({\mathbf{x}})/p({\mathbf{x}})$ by flowing samples from
$q({\mathbf{x}})$ at each training step with Eq. 5. This approch is based on
observations that $\tilde{q}({\mathbf{x}})$ progressively approaches the
target distribution as $r_{\theta}$ is being trained. In this formulation,
$\tilde{q}({\mathbf{x}})$ is no longer static but changes with $r_{\theta}$ as
$\theta$ is updated. Therefore, unlike the toy examples of Sec. 3.1, the
density ratio estimator is being trained on samples that improve with
training.

We formulate our training setup as follows. Consider a given training step
$\tau$, where our model has parameters $\theta_{\tau}$ from the gradient
update of the previous iteration. We switch notations slightly and denote
$q^{\prime}({\mathbf{x}})$ as our source distribution, which is a simple prior
(such as the DDP). We propose to first draw samples ${\mathbf{x}}_{0}\sim
q^{\prime}({\mathbf{x}})$, and simulate Eq. 5 for $K$ steps with the stale
estimator $r_{\theta_{\tau}}$ using the Euler-Maruyama discretization:

${\mathbf{x}}_{k+1}=-\eta\nabla_{\mathbf{x}}f^{\prime}(r_{\theta_{\tau}})({\mathbf{x}}_{k})+\nu\bm{\xi}_{k}$
(9)

where $\bm{\xi}_{k}\sim\mathcal{N}(0,I)$, $\eta$ is the step size and
$k\in[0,K-1]$. We label the resultant particles as ${\mathbf{x}}_{K}$. The
${\mathbf{x}}_{K}$ are drawn from the distribution
$\tilde{q}_{\tau}({\mathbf{x}}_{K})$ given by

$\tilde{q}_{\tau}({\mathbf{x}}_{K})=\int
q^{\prime}({\mathbf{x}}^{\prime})M_{\theta_{\tau}}({\mathbf{x}}_{K}|{\mathbf{x}}^{\prime})d{\mathbf{x}}$
(10)

where $M_{\theta_{\tau}}({\mathbf{x}}_{K}|{\mathbf{x}})$ is the transition
kernel of simulating Eq. 9 with parameters $\theta_{\tau}$. We optimize
$r_{\theta_{\tau}}({\mathbf{x}})$ for the $\tau$ training iteration using the
Bregman divergence

$\mathcal{L}(\theta_{\tau})=\mathbb{\widehat{E}}_{p}[\partial
g(r_{\theta_{\tau}}({\mathbf{x}}))r({\mathbf{x}})-g(r_{\theta_{\tau}}({\mathbf{x}}))]-\mathbb{\widehat{E}}_{\tilde{q}_{\tau}}[\partial
g(r_{\theta_{\tau}}({\mathbf{x}}))]$ (11)

where expectation over $p({\mathbf{x}})$ can be estimated using samples from
the dataset. We update the parameters for the next training step $\tau+1$ as
$\theta_{\tau+1}\leftarrow\theta_{\tau}-\alpha\nabla_{\theta}\mathcal{L}(\theta_{\tau})$
with learning rate $\alpha$. The process repeats until convergence. We provide
pseudocode in Algorithm 1, and an illustration of the training setup in Fig.
1. We illustrate our training process in Fig. 3, where we utilize the toy
example from Fig. 2(b) which the stale estimator $r_{\theta}=q^{\prime}/p$
failed previously. We can see from the top row that as $\tau$ increases, the
orange samples from $\tilde{q}_{\tau}$, which are obtained by flowing the blue
samples from $q^{\prime}$ for fixed $K$ steps, progressively approach the
target distribution. In the bottom row, we plot the progression of the mean of
$\tilde{q}_{\tau}$ as training progresses, which shows the iterative approach
towards the target distribution. Together with the experiments for high-
dimensional image synthesis in Sec. 6, these results empirically justify that
FDRL indeed overcomes the density chasm problem, allowing the stale estimator
$r_{\theta}({\mathbf{x}})$ to be used for sample generation.

Figure 3: Plot of $\tilde{q}_{\tau}$ as training progresses in the toy example
of Fig. 2(b). The total number of training steps is $T=1000$, such that
$\tau\in[0,T-1]$. In the top row, blue particles are samples from $q^{\prime}$
while orange particles are samples from $\tilde{q}_{\tau}$ at specific
training steps. Left: $\tau=10$, center: $\tau=40$, right: $\tau=1000$. In the
bottom row, we plot the trajectory of the mean of $\tilde{q}_{\tau}$ as
training progresses.

###  4.1 Convergence Distribution of FDRL

In this section, we analyze the convergence distribution after training has
converged. With $T$ training steps and $\tau\in[0,T-1]$, let us consider what
happens during the final training iteration $T-1$. The model has parameters
$\theta_{T-1}$ from the gradient update of the previous iteration. Following
Algorithm 1, we first sample ${\mathbf{x}}_{K}$ from the distribution
$\tilde{q}_{T-1}({\mathbf{x}}_{K})=\int
q^{\prime}({\mathbf{x}}^{\prime})M_{\theta_{T-1}}({\mathbf{x}}_{K}|{\mathbf{x}}^{\prime})d{\mathbf{x}}^{\prime}$
by running $K$ steps of Eq. 9. With samples
${\mathbf{x}}_{K}\sim\tilde{q}_{T-1}({\mathbf{x}})$ and ${\mathbf{x}}_{p}\sim
p({\mathbf{x}})$, we optimize the BD and take the gradient step
$\theta_{T}\leftarrow\theta_{T-1}-\alpha\nabla_{\theta}\mathcal{L}$. Assuming
convergence, this means by definition that the loss gradient of the final step
$\nabla_{\theta}\mathcal{L}=0$. Thus, with no further gradient updates
possible, our model is parameterized by $\theta_{T}$ after convergence. We can
now analyze what happens when sampling from our trained model with parameters
$\theta_{T}$.

In our framework, we sample from our trained model by flowing for a total of
$K+\kappa$ steps. This sampling can be viewed as two stages ($K$ “bridging”
steps, followed by $\kappa$ refinement steps). To elaborate, in the first
stage, we sample as before by running Eq. 9 for $K$ steps using our trained
model $r_{\theta_{T}}$. The samples from stage 1 will be drawn from
$\tilde{q}_{T}({\mathbf{x}}_{K})=\int
q^{\prime}({\mathbf{x}}^{\prime})M_{\theta_{T}}({\mathbf{x}}_{K}|{\mathbf{x}}^{\prime})d{\mathbf{x}}^{\prime}$.
Based on empirical evidence from our toy example in Fig. 3 and image
experiments in Sec. 6.2, $\tilde{q}_{T}({\mathbf{x}}_{K})$ is closer to the
data distribution $p({\mathbf{x}})$. Interestingly, we can introduce a second
refinement stage. Recall that convergence in the Bregman divergence,
$\nabla_{\theta}\mathcal{L}=0$, implies that our density ratio estimator has
converged to the value $r_{\theta_{T}}({\mathbf{x}})=\tilde{q}_{T-1}/p$. Since
after the first stage, we have samples from $\tilde{q}_{T}({\mathbf{x}})$, and
an estimator
$r_{\theta_{T}}({\mathbf{x}})=\tilde{q}_{T-1}/p\approx\tilde{q}_{T}/p$
(assuming the change in $\theta$ in the final training step is small), we can
perform sample refinement in exactly the formalism of DG$f$low by running the
same Eq. 9 for a further $\kappa$ steps. We provide pseudocode for sampling in
Algorithm 2 and justify the two-stage sampling scheme with ablations in Sec.
6.3.

Algorithm 1 Training

repeat

Sample ${\mathbf{x}}_{p}\sim p({\mathbf{x}}),{\mathbf{x}}_{0}\sim
q^{\prime}({\mathbf{x}})$

for $j\leftarrow 1,K$ do

Obtain ${\mathbf{x}}_{K}$ from ${\mathbf{x}}_{0}$ by simulating Eq. 9.

end for

Update $\theta$ according to
$\begin{multlined}\nabla_{\theta}[g^{\prime}(r_{\theta}({\mathbf{x}}_{p}))r({\mathbf{x}}_{p})-g(r_{\theta}({\mathbf{x}}_{p}))-g^{\prime}(r_{\theta}({\mathbf{x}}_{K}))]\end{multlined}\nabla_{\theta}[g^{\prime}(r_{\theta}({\mathbf{x}}_{p}))r({\mathbf{x}}_{p})-g(r_{\theta}({\mathbf{x}}_{p}))-g^{\prime}(r_{\theta}({\mathbf{x}}_{K}))]$

until converged

Algorithm 2 Sampling

Sample ${\mathbf{x}}_{0}\sim q^{\prime}({\mathbf{x}})$

for $j\leftarrow 1,K+\kappa$ do

Obtain ${\mathbf{x}}_{K+\kappa}$ from ${\mathbf{x}}_{0}$ by simulating Eq. 9.

end for

return ${\mathbf{x}}_{K+\kappa}$

Figure 4: Samples from FDRL-DDP on CIFAR10 $32^{2}$, CelebA $64^{2}$ and LSUN
Church $128^{2}$ using LSIF-$\chi^{2}$. More results using various BD
objectives and $f$-divergences can be found in Appendix.

##  5 Related Works

Gradient flows are a general framework for constructing the steepest descent
curve of a given functional, and consequently have been used in optimizing a
variety of distance metrics, ranging from the $f$-divergence [10, 9, 1, 8],
maximum mean discrepancy [2, 21], Sobolev distance [22] and related forms of
the Wasserstein distance [18].

Two well-known methodologies to simulate Wasserstein gradient flows are
population density approaches such as the JKO scheme [14] and particle-based
approaches. The former requires computation of the 2-Wasserstein distance and
a free energy functional, which are generally intractable. As such, several
works have proposed approximations leveraging Brenier’s theorem [20] or the
Fenchel conjugate [8].

On the other hand, particle approaches involve simulating the SDE Eq. 3.
DG$f$low leverages GAN discriminators as a stale density ratio estimator for
sample refinement. Our method is inspired by the stale formulation of
DG$f$low, but we extend it to the general case of sample generation. In fact,
DG$f$low is a specific instance of FDRL where we fix the prior $q^{\prime}$ to
be the implicit distribution defined by the GAN generator. Euler Particle
Transport (EPT) [9] similarly trains a deep density ratio estimator using the
Bregman divergence. However, EPT was only shown to be effective in the latent
space and thus necessitated the training of an additional generator network,
whereas our method is able to scale directly in the data space.

Energy-based models (EBMs) adopt Langevin dynamics to sample from an
unnormalized distribution [15, 37]. The training scheme of FDRL bears
similarity with that of Short-run EBM [24], where a sample is drawn by
restarting the Langevin chain and running $K$ gradient steps. However, FDRL is
fundamentally distinct from EBMs; the former is trained with density ratio
estimation, while the latter is trained with maximum likelihood. Further, it
can be shown that Langevin dynamics is in fact a special case of the KL
gradient flow [14, 17].

##  6 Experiments

The goal of our experiments is to compare FDRL against gradient flow
baselines. However, for completeness, we also include results from other
modeling approaches, including EBMs, GANs, and Diffusion models.

###  6.1 Setup

We test FDRL with unconditional generation on 32$\times$32 CIFAR10,
64$\times$64 CelebA and 128$\times$128 LSUN Church datasets, class-conditional
generation on CIFAR10 and image-to-image translation on 64$\times$64 Cat2dog
dataset. We label experiments with the DDP as FDRL-DDP, and ablation
experiments with a uniform prior as FDRL-UP. We use modified ResNet
architectures for all experiments in this study. See Appendix F for more
details.

###  6.2 Unconditional Image Generation

Table 1: FID scores for 32$\times$32 CIFAR10 and 64$\times$64 CelebA. Lower is
better.

Method | CIFAR10 | CelebA
---|---|---
EBM-based |  |
IGEBM [6] | 40.58 | -
SR-EBM [24] | 44.50 | 23.02
LEBM [25] | 70.15 | 37.87
CoopNets [36] | 33.61 | -
VAEBM [35] | 12.19 | 5.31
Diffusion-based |  |
NCSN [30] | 25.32 | 26.89
DDPM [13] | 3.17 | -
GAN-based |  |
WGAN-GP [12] | 55.80 | 30.00
SNGAN [19] | 21.70 | -
Gradient flow-based |  |
EPT [10] | 24.60 | -
JKO-Flow [8] | 23.70 | -
$\ell$-SWF [5] | 59.70 | 38.30
Ours |  |
FDRL-DDP (LSIF-$\chi^{2}$) | 22.28 | 17.77
FDRL-DDP (LR-KL) | 22.61 | 18.09
FDRL-DDP (LR-JS) | 22.80 | -
FDRL-DDP (LR-logD) | 22.82 | -
FDRL-UP (LSIF-$\chi^{2}$ | 30.23 | -
FDRL-UP (LR-KL) | 31.99 | -

In Fig. 4, we show uncurated samples of FDRL-DDP on different combinations of
$g$ and $f$-divergences, where the combinations are chosen due to numerical
compatibility (see Appendix C). Visually, our model is able to produce high-
quality samples on a variety of datasets up to resolutions of $128\times 128$,
surpassing existing gradient flow techniques [9, 8]. Samples with other
$f$-divergences can be found in Appendix H. In Table 1, we show the FID scores
of FDRL-DDP in comparison with relevant baselines utilizing different
generative approaches.

We emphasize that our experiments are geared towards comparisons with gradient
flow baselines, in which FDRL-DDP outperforms both EPT and JKO-Flow, as well
as the non-parametric method $\ell$-SWF. Our method outperforms all baseline
EBMs except VAEBM, which combines a VAE and EBM for improved performance. We
also outperform WGAN and the score-based NCSN, while matching SNGAN’s
performance. Comparing method approaches, modern diffusion models (DDPM)
outperform EBM and gradient flow methods, including FDRL. Interesting future
work could investigate reasons for this gap, which may lead to improvements
for EBMs and gradient flows.

To intuit the flow process, we show how samples evolve with the LSUN Church
dataset in Fig. 8 of the appendix, where the prior contains low frequency
features such as the rough silhouette of the building, and the flow process
generates the higher frequency details to create a realistic image.

We also visualize samples obtained by interpolating in the prior space for
CelebA in Fig. 9 of the appendix. The model smoothly interpolates in the
latent space despite using the DDP, indicating a semantically relevant latent
representation that is characteristic of a valid generative model. Finally, to
verify that our model has not merely memorized the dataset, we show nearest
neighbor samples of the generated images in the training set in Appendix G.
The samples produced by FDRL-DDP are distinct from their closest training
samples, telling us that FDRL-DDP is capable of generating new and diverse
samples beyond the training data.

###  6.3 Ablations

Figure 5: FID as a function of the total flow length on CIFAR10 for
LSIF-$\chi^{2}$ DDP when sampling from a model trained with $K=100$.

#### Data-Dependent Prior vs Uniform Prior

We motivate the use of the data-dependent prior in conjunction with FDRL by
conducting ablations with $q^{\prime}$ being a uniform prior (UP),
${\mathbf{x}}_{0}\sim U[-1,1]$. We train UP runs for 20$\%$ more training
steps than DDP runs, but keep all other hyperparameters unchanged for control.
We include qualitative samples of FDRL-UP in Fig. 21 in the appendix. Table 1
reports the quantitative scores. Visually, FDRL-UP produces diverse and
appealing samples even with a uniform prior. However, the quantitative scores
are significantly poorer than using the DDP, which validates our hypothesis in
Sec. 3.2.

#### Length of Flow vs Image Quality

Fig. 5 illustrates the relationship between total flow length and FID.
Although the model is trained with $K=100$, sampling with $K$ steps at testing
does not produce the best image quality, as sampling with an additional
$\kappa=20$ steps yields the best FID score. This supports the two-stage
sampling approach discussed in Sec. 4.1. First, $K=100$ steps of the flow are
run to sample images from $\tilde{q}_{T}({\mathbf{x}}_{K})$, followed by
$\kappa=20$ additional steps of the same Eq. 9 using the estimator
$r{\theta_{T}}\approx\tilde{q}_{T}/p$ for sample refinement, as done in
DG$f$low.

###  6.4 Conditional Image Generation

Figure 6: Class-conditional samples obtained by composing an unconditional
FDRL with a pretrained robust classifier.

We can sample from target densities different from what FDRL was trained on by
simply composing density ratios. Consider a multiclass classifier which
classifies a given image into one of $N$ classes. We show in Appendix E that
we can express such classifiers as a density ratio
$p(y=n|{\mathbf{x}})=N^{-1}p({\mathbf{x}}|y=n)/p({\mathbf{x}})$. Thus, we can
obtain a conditional density ratio estimator
$r_{\theta}({\mathbf{x}}|y=n)=q_{T}({\mathbf{x}})/p({\mathbf{x}}|y=n)$ by
composing our unconditional estimator $r_{\theta}({\mathbf{x}})$ with the
classifier output (see Appendix E):

$r_{\theta}({\mathbf{x}}|y=n)=\frac{1}{N}r_{\theta}({\mathbf{x}})p(y=n|{\mathbf{x}})^{-1}.$
(12)

When $r_{\theta}({\mathbf{x}}|y=n)$ is used in the flow in Eq. 9, we perform
class-conditional generation. This is conceptually similar to the idea
proposed in [32, 4], where an unconditional score model
$\nabla_{\mathbf{x}}\log p_{t}({\mathbf{x}}(t))$ is composed with a time-
dependent classifier $\nabla_{\mathbf{x}}\log p_{t}(y|{\mathbf{x}}(t))$ to
form a class-conditional model. However, whereas [32] requires separately
training a time-dependent classifier, our formulation allows the use
pretrained classifiers off-the-shelf, with _no further retraining_. Inspired
by earlier work on image synthesis with robust classifiers [29], we found that
using a pretrained adversarially-robust classifier was necessary in obtaining
useful gradients for generation. We show our results in Fig. 6, where each row
represents conditional samples of each class in the CIFAR10 dataset.

###  6.5 Unpaired Image-to-image Translation

Figure 7: Image-to-image translation process from cat to dog images using
FDRL.

FDRL can also be seamlessly applied to unpaired image-to-image-translation
(I2I). We simply fix the prior distribution $q^{\prime}({\mathbf{x}})$ to a
source image domain and $p({\mathbf{x}})$ to a target domain. We then train
the model in exactly the same manner as unconditional generation using
Algorithm 1.

The I2I model is tested on the Cat2dog dataset [16]. From Fig. 7, FDRL is able
to smoothly translate images of cats to dogs while preserving relevant
features, such as the background colors, pose and facial tones. For instance,
a cat with light fur is translated to a dog with light fur. Quantitatively,
the I2I baseline CycleGAN [39] achieves better FID scores than FDRL (Appendix
Table 4). However, like many I2I methods [16, 3, 38, 23], CycleGAN relies on
specific inductive biases, such as dual generators and discriminators, and
cycle-consistency loss. Future work could explore incorporating these biases
into FDRL to improve translation performance.

##  7 Conclusion

We propose FDRL, a method that enables the stale approximation of the gradient
flow to be applied to generative modeling of high-dimensional images. Beyond
unconditional image generation, the FDRL framework can be seamlessly adapted
to other key applications, including class-conditional generation and unpaired
image-to-image translation. Future work could focus on theoretically
characterizing the exact distribution induced by the stale SDE of Eq. 5 and
further improving our understanding of its convergence properties.

##  8 Acknowledgements

This research is supported by the National Research Foundation Singapore and
DSO National Laboratories under the AI Singapore Programme (AISG Award No:
AISG2-RP-2020-016).

## References

  * [1] Abdul Fatir Ansari, Ming Liang Ang, and Harold Soh.  Refining deep generative models via discriminator gradient flow.  In ICLR, 2021.
  * [2] Michael Arbel, Anna Korba, Adil Salim, and Arthur Gretton.  Maximum mean discrepancy gradient flow.  Advances in Neural Information Processing Systems, 32, 2019.
  * [3] Yunjey Choi, Youngjung Uh, Jaejun Yoo, and Jung-Woo Ha.  Stargan v2: Diverse image synthesis for multiple domains.  In Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, pages 8188–8197, 2020.
  * [4] Prafulla Dhariwal and Alexander Nichol.  Diffusion models beat gans on image synthesis.  Advances in Neural Information Processing Systems, 34:8780–8794, 2021.
  * [5] Chao Du, Tianbo Li, Tianyu Pang, Shuicheng Yan, and Min Lin.  Nonparametric generative modeling with conditional and locally-connected sliced-wasserstein flows.  arXiv preprint arXiv:2305.02164, 2023.
  * [6] Yilun Du and Igor Mordatch.  Implicit generation and generalization in energy-based models.  arXiv preprint arXiv:1903.08689, 2019.
  * [7] Logan Engstrom, Andrew Ilyas, Shibani Santurkar, and Dimitris Tsipras.  Robustness (python library), 2019.
  * [8] Jiaojiao Fan, Amirhossein Taghvaei, and Yongxin Chen.  Variational wasserstein gradient flow.  arXiv preprint arXiv:2112.02424, 2021.
  * [9] Yuan Gao, Jian Huang, Yuling Jiao, Jin Liu, Xiliang Lu, and Zhijian Yang.  Deep generative learning via euler particle transport.  In Mathematical and Scientific Machine Learning, pages 336–368. PMLR, 2022.
  * [10] Yuan Gao, Yuling Jiao, Yang Wang, Yao Wang, Can Yang, and Shunkang Zhang.  Deep generative learning via variational gradient flow.  In International Conference on Machine Learning, pages 2093–2101. PMLR, 2019.
  * [11] Ian Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair, Aaron Courville, and Yoshua Bengio.  Generative adversarial nets.  Advances in neural information processing systems, 27, 2014.
  * [12] Ishaan Gulrajani, Faruk Ahmed, Martin Arjovsky, Vincent Dumoulin, and Aaron C Courville.  Improved training of wasserstein gans.  Advances in neural information processing systems, 30, 2017.
  * [13] Jonathan Ho, Ajay Jain, and Pieter Abbeel.  Denoising diffusion probabilistic models.  Advances in Neural Information Processing Systems, 33:6840–6851, 2020.
  * [14] Richard Jordan, David Kinderlehrer, and Felix Otto.  The variational formulation of the fokker–planck equation.  SIAM journal on mathematical analysis, 29(1):1–17, 1998.
  * [15] Yann LeCun, Sumit Chopra, Raia Hadsell, M Ranzato, and F Huang.  A tutorial on energy-based learning.  Predicting structured data, 1(0), 2006.
  * [16] Hsin-Ying Lee, Hung-Yu Tseng, Jia-Bin Huang, Maneesh Singh, and Ming-Hsuan Yang.  Diverse image-to-image translation via disentangled representations.  In Proceedings of the European conference on computer vision (ECCV), pages 35–51, 2018.
  * [17] Chang Liu, Jingwei Zhuo, and Jun Zhu.  Understanding mcmc dynamics as flows on the wasserstein space.  In International Conference on Machine Learning, pages 4093–4103. PMLR, 2019.
  * [18] Antoine Liutkus, Umut Simsekli, Szymon Majewski, Alain Durmus, and Fabian-Robert Stöter.  Sliced-wasserstein flows: Nonparametric generative modeling via optimal transport and diffusions.  In International Conference on Machine Learning, pages 4104–4113. PMLR, 2019.
  * [19] Takeru Miyato, Toshiki Kataoka, Masanori Koyama, and Yuichi Yoshida.  Spectral normalization for generative adversarial networks.  arXiv preprint arXiv:1802.05957, 2018.
  * [20] Petr Mokrov, Alexander Korotin, Lingxiao Li, Aude Genevay, Justin M Solomon, and Evgeny Burnaev.  Large-scale wasserstein gradient flows.  Advances in Neural Information Processing Systems, 34:15243–15256, 2021.
  * [21] Youssef Mroueh and Truyen Nguyen.  On the convergence of gradient descent in gans: Mmd gan as a gradient flow.  In International Conference on Artificial Intelligence and Statistics, pages 1720–1728. PMLR, 2021.
  * [22] Youssef Mroueh, Tom Sercu, and Anant Raj.  Sobolev descent.  In The 22nd International Conference on Artificial Intelligence and Statistics, pages 2976–2985. PMLR, 2019.
  * [23] Weili Nie, Arash Vahdat, and Anima Anandkumar.  Controllable and compositional generation with latent-space energy-based models.  Advances in Neural Information Processing Systems, 34:13497–13510, 2021.
  * [24] Erik Nijkamp, Mitch Hill, Song-Chun Zhu, and Ying Nian Wu.  Learning non-convergent non-persistent short-run mcmc toward energy-based model.  Advances in Neural Information Processing Systems, 32, 2019.
  * [25] Bo Pang, Tian Han, Erik Nijkamp, Song-Chun Zhu, and Ying Nian Wu.  Learning latent space energy-based prior model.  Advances in Neural Information Processing Systems, 33:21994–22008, 2020.
  * [26] Benjamin Rhodes, Kai Xu, and Michael U Gutmann.  Telescoping density-ratio estimation.  Advances in neural information processing systems, 33:4905–4916, 2020.
  * [27] Hannes Risken and JH Eberly.  The fokker-planck equation, methods of solution and applications.  Journal of the Optical Society of America B Optical Physics, 2(3):508, 1985.
  * [28] Filippo Santambrogio.  $\\{$Euclidean, metric, and Wasserstein$\\}$ gradient flows: an overview.  Bulletin of Mathematical Sciences, 7(1):87–154, 2017.
  * [29] Shibani Santurkar, Andrew Ilyas, Dimitris Tsipras, Logan Engstrom, Brandon Tran, and Aleksander Madry.  Image synthesis with a single (robust) classifier.  Advances in Neural Information Processing Systems, 32, 2019.
  * [30] Yang Song and Stefano Ermon.  Generative modeling by estimating gradients of the data distribution.  Advances in Neural Information Processing Systems, 32, 2019.
  * [31] Yang Song and Stefano Ermon.  Improved techniques for training score-based generative models.  Advances in neural information processing systems, 33:12438–12448, 2020.
  * [32] Yang Song, Jascha Sohl-Dickstein, Diederik P Kingma, Abhishek Kumar, Stefano Ermon, and Ben Poole.  Score-based generative modeling through stochastic differential equations.  arXiv preprint arXiv:2011.13456, 2020.
  * [33] Masashi Sugiyama, Taiji Suzuki, and Takafumi Kanamori.  Density ratio estimation in machine learning.  Cambridge University Press, 2012.
  * [34] Masashi Sugiyama, Taiji Suzuki, and Takafumi Kanamori.  Density-ratio matching under the bregman divergence: a unified framework of density-ratio estimation.  Annals of the Institute of Statistical Mathematics, 64(5):1009–1044, 2012.
  * [35] Zhisheng Xiao, Karsten Kreis, Jan Kautz, and Arash Vahdat.  Vaebm: A symbiosis between variational autoencoders and energy-based models.  arXiv preprint arXiv:2010.00654, 2020.
  * [36] Jianwen Xie, Yang Lu, Ruiqi Gao, and Ying Nian Wu.  Cooperative learning of energy-based model and latent variable model via mcmc teaching.  In Proceedings of the AAAI Conference on Artificial Intelligence, volume 32, 2018.
  * [37] Jianwen Xie, Yang Lu, Song-Chun Zhu, and Yingnian Wu.  A theory of generative convnet.  In International Conference on Machine Learning, pages 2635–2644. PMLR, 2016.
  * [38] Yang Zhao and Changyou Chen.  Unpaired image-to-image translation via latent energy transport.  In Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, pages 16418–16427, 2021.
  * [39] Jun-Yan Zhu, Taesung Park, Phillip Isola, and Alexei A Efros.  Unpaired image-to-image translation using cycle-consistent adversarial networks.  In Proceedings of the IEEE international conference on computer vision, pages 2223–2232, 2017.

Figure 8: Illustration of the flow process for LSUN Church, starting from a
sample from the data-dependent prior on the leftmost column. Figure 9:
Interpolation results between leftmost and rightmost samples with CelebA.

##  Appendix A Proofs

See 1

###### Proof.

The proof draws inspiration from the Langevin equation. Therefore, let us
first consider the standard overdamped Langevin SDE:

$d{\mathbf{x}}_{t}=-\nabla f({\mathbf{x}}_{t})dt+\sqrt{2}d{\mathbf{w}}_{t}$
(13)

where $f$ is Lipschitz. The associated Fokker-Planck equation is given by [27]

$\partial_{t}\rho_{t}=\nabla\cdot(\rho\nabla f)+\Delta\rho.$ (14)

where $\Delta$ is the Laplace operator. It is straightforward to show that the
stationary distribution of Eq. 13 is given by

$\rho_{\infty}({\mathbf{x}})=\frac{e^{-f({\mathbf{x}})}}{\int
e^{-f({\mathbf{x}})}d{\mathbf{x}}}=\frac{e^{-f({\mathbf{x}})}}{Z}.$ (15)

We simply verify by substituting Eq. 15 into Eq. 14:

$\displaystyle\partial_{t}\rho_{\infty}$
$\displaystyle=\nabla\cdot\left(\frac{e^{-f({\mathbf{x}})}}{Z}\nabla
f\right)+\nabla\cdot\nabla\left(\frac{e^{-f({\mathbf{x}})}}{Z}\right)$
$\displaystyle=\nabla\cdot\left(\frac{e^{-f({\mathbf{x}})}}{Z}\nabla
f\right)-\nabla\cdot\left(\frac{e^{-f({\mathbf{x}})}}{Z}\nabla f\right)$
$\displaystyle=0.$

Now, let us return to the stale estimate SDE of Eq. 5 with $\gamma=1$ and
$q({\mathbf{x}})\sim U[a,b]=\frac{1}{b-a}=C$:

$d{\mathbf{x}}_{t}=-\nabla
f^{\prime}(C/p({\mathbf{x}}_{t}))dt+\sqrt{2}d{\mathbf{w}}_{t}.$ (16)

Since $f$ is convex and twice-differentiable, $f^{\prime}$ is a monotonically
increasing function. Using what we showed above, the stationary distribution
of Eq. 16 is thus

$\rho_{\infty}({\mathbf{x}})=\frac{e^{-f^{\prime}(C/p({\mathbf{x}}))}}{\int
e^{-f(C/p({\mathbf{x}}))}d{\mathbf{x}}}=\frac{e^{-f^{\prime}(C/p({\mathbf{x}}))}}{Z}.$
(17)

Then we can show that

$\displaystyle\operatorname*{arg\,max}_{\mathbf{x}}\rho_{\infty}({\mathbf{x}})$
$\displaystyle=\operatorname*{arg\,max}_{\mathbf{x}}\frac{e^{-f^{\prime}(C/p({\mathbf{x}}))}}{Z}$
$\displaystyle=\operatorname*{arg\,max}_{\mathbf{x}}e^{-f^{\prime}(C/p({\mathbf{x}}))}$
$\displaystyle=\operatorname*{arg\,max}_{\mathbf{x}}-f^{\prime}(C/p({\mathbf{x}}))$
($e$ is monotonically increasing)
$\displaystyle=\operatorname*{arg\,min}_{\mathbf{x}}f^{\prime}(C/p({\mathbf{x}}))$
$\displaystyle=\operatorname*{arg\,min}_{\mathbf{x}}C/p({\mathbf{x}})$
($f^{\prime}$ is monotonically increasing)
$\displaystyle=\operatorname*{arg\,max}_{\mathbf{x}}p({\mathbf{x}})$
$\displaystyle\text{($p({\mathbf{x}})>0$)}.$

∎

##  Appendix B Toy Datasets

Figure 10: Comparison of different FDRL pairings of Bregman and
$f$-divergences on the 2DSwissroll dataset.

To affirm that samples generated by FDRL indeed converge to the target
distribution, we train FDRL on the synthetic 2DSwissroll dataset. The density
ratio estimator is parameterized by a simple feedforward multilayer
perceptron. We train the model to flow samples from the prior
$q^{\prime}({\mathbf{x}})=\mathcal{N}(\mathbf{0},\mathbb{I})$ to the target
distribution, which we sample from the make_swiss_roll function in scikit-
learn. We plot the results in Fig. 10, from which we can see that the model
indeed recovers $p({\mathbf{x}})$ successfully for all combinations of $f$ and
$g$.

##  Appendix C Bregman Divergence and $f$-divergence Pairing

We elucidate the full forms of the Bregman divergences here. The LSIF
objective $g(y)=\frac{1}{2}(y-1)^{2}$ is given as

$\mathcal{L}_{LSIF}(\theta)=\frac{1}{2}\mathbb{\widehat{E}}_{p}[r_{\theta}({\mathbf{x}})]^{2}-\mathbb{\widehat{E}}_{q_{\tau}}[r_{\theta}({\mathbf{x}})],$
(18)

while the LR objective $g(t)=y\log y-(1+y)\log(1+y)$ is given as

$\mathcal{L}_{LR}(\theta)=-\mathbb{\widehat{E}}_{p}\left[\log\frac{1}{1+r_{\theta}({\mathbf{x}})}\right]-\mathbb{\widehat{E}}_{q_{\tau}}\left[\log\frac{r_{\theta}({\mathbf{x}})}{1+r_{\theta}({\mathbf{x}})}\right].$
(19)

When computing the LR objective, we find that we run into numerical stability
issues when letting $r_{\theta}({\mathbf{x}})$ be the output of an
unconstrained neural network and subsequently taking the required logarithms
in Eq. 19. To circumvent this issue, we let $r_{\theta}({\mathbf{x}})$ be
expressed as the exponential of the neural network’s output, i.e., the output
of the neural network is $\log r_{\theta}({\mathbf{x}})$. This formulation
naturally lends itself to the flow of the KL, JS and logD divergences, whose
first derivatives $f^{\prime}$ that is required in Eq. 9 are also logarithmic
functions of $r_{\theta}({\mathbf{x}})$, as seen from Table. 2. We can thus
utilize numerically stable routines in existing deep learning frameworks,
avoiding the need for potentially unstable operations like exponentiations
(see Appendix D for details). As such, we pair LR with the aforementioned
divergences and abbreviate the combinations as LR-KL, LR-JS, LR-logD. We did
not run into such stability issues for the LSIF objective (Eq. 18) as the
model learns to automatically output a non-negative scalar over the course of
training, hence for LSIF we allow the neural network to estimate
$r_{\theta}({\mathbf{x}})$ directly and pair it with the Pearson-$\chi^{2}$
divergence. We abbreviate this pairing as LSIF-$\chi^{2}$.

##  Appendix D Stable Computation of LR and $f$-divergences

As mentioned in Sec. C, computing the logarithm of unconstrained neural
networks leads to instabilities in the training process. This is a problem
when computing the LR objective in Eq. 19 and the various first derivatives of
$f$-divergences. We can circumvent this by letting $r_{\theta}({\mathbf{x}})$
be expressed as the exponential of the neural network and use existing stable
numerical routines to avoid intermediate computations that lead to the
instabilities (for example, computing logarithms and exponentials directly).
Let us express the neural network output as
$NN_{\theta}({\mathbf{x}})\triangleq\log r_{\theta}({\mathbf{x}})$. The LR
objective can then be rewritten as

$\displaystyle\mathcal{L}_{LR}(\theta)$
$\displaystyle=-\mathbb{\widehat{E}}_{p}\left[\log\frac{1}{1+r_{\theta}({\mathbf{x}})}\right]-\mathbb{\widehat{E}}_{q_{\tau}}\left[\log\frac{r_{\theta}({\mathbf{x}})}{1+r_{\theta}(x)}\right]$
(20)
$\displaystyle=-\mathbb{\widehat{E}}_{p}\left[\texttt{LS}(-NN_{\theta}(x))\right]-\mathbb{\widehat{E}}_{q_{\tau}}\left[\texttt{LS}(NN_{\theta}(x))\right]$
(21)

where $\texttt{LS}(x)=\log\frac{1}{1+\exp(-x)}$, the log-sigmoid function,
which has stable implementations in modern deep learning libraries.

Similarly for the $f$-divergences whose first derivatives involve logarithms,
we can calculate them stably as

$\displaystyle f^{\prime}_{KL}(r_{\theta}({\mathbf{x}}))$ $\displaystyle=\log
r_{\theta}({\mathbf{x}})+1=NN_{\theta}({\mathbf{x}})+1$ (22) $\displaystyle
f^{\prime}_{JS}(r_{\theta}({\mathbf{x}}))$
$\displaystyle=\log\frac{2r_{\theta}({\mathbf{x}})}{1+r_{\theta}({\mathbf{x}})}=\log
2+\texttt{LS}(NN_{\theta}({\mathbf{x}}))$ (23) $\displaystyle
f^{\prime}_{logD}(r_{\theta}({\mathbf{x}}))$
$\displaystyle=\log(r_{\theta}({\mathbf{x}})+1)+1=-\texttt{LS}(-NN_{\theta}({\mathbf{x}}))+1.$
(24)

Table 2: $f$-divergences and their first derivatives $f^{\prime}$. $f$-divergence | $f$ | $f^{\prime}$
---|---|---
Pearson-$\chi^{2}$ | $(r-1)^{2}$ | $2(r-1)$
KL | $r\log r$ | $\log r+1$
JS | $r\log r-(r+1)\log\frac{r+1}{2}$ | $\log\frac{2r}{r+1}$
log D | $(r+1)\log(r+1)-2\log 2$ | $\log(r+1)+1$

##  Appendix E Classifiers are Density Ratio Estimators

To perform class-conditional generation in the FDRL framework, we would like
to estimate the density ratio of a certain class over the data distribution:
$p({\mathbf{x}}|y=n)/p({\mathbf{x}})$. With Bayes rule, we can write this as

$\displaystyle\frac{p({\mathbf{x}}|y=n)}{p({\mathbf{x}})}$
$\displaystyle=\frac{p(y=n|{\mathbf{x}})p({\mathbf{x}})/p(y=n)}{p({\mathbf{x}})}$
(25) $\displaystyle=\frac{p(y=n|{\mathbf{x}})}{p(y=n)}.$ (26)

The denominator term $p(y=n)$ can be viewed as a constant, e.g. assume the $N$
classes are equally distributed, then $p(y=n)=1/N$. Therefore, we have that
the class probability given by the softmax output of a classifier is actually
a density ratio:

$Np(y=n|{\mathbf{x}})=\frac{p({\mathbf{x}}|y=n)}{p({\mathbf{x}})}.$ (27)

We can use this equation to convert an unconditional FDRL to a class-
conditional generator. Recall the stale approximation of the gradient flow:

$d{\mathbf{x}}_{t}=-\nabla_{\mathbf{x}}f^{\prime}(r_{\theta}({\mathbf{x}}_{t}))dt+\sqrt{2\gamma}d{\mathbf{w}}_{t}$
(28)

Consider the case where we have a trained unconditional FDRL estimator
$r_{\theta}$. We can multiply the inverse of the classifier output with
$r_{\theta}({\mathbf{x}}_{t})=\tilde{q}_{\tau}({\mathbf{x}}_{t})/p({\mathbf{x}}_{t})$
to get a density ratio between $\tilde{q}_{\tau}({\mathbf{x}}_{t})$ and the
conditional data distribution $p({\mathbf{x}}_{t}|y=n)$:

$\displaystyle r_{\theta}({\mathbf{x}}_{t})p(y=n|{\mathbf{x}}_{t})^{-1}$
$\displaystyle=\frac{\tilde{q}_{\tau}({\mathbf{x}}_{t})}{p({\mathbf{x}}_{t})}\frac{Np({\mathbf{x}}_{t})}{p({\mathbf{x}}_{t}|y=n)}$
$\displaystyle=N\frac{\tilde{q}_{\tau}({\mathbf{x}}_{t})}{p({\mathbf{x}}_{t}|y=n)}.$

That is, we took our unconditional model and converted it to a conditional
generative model by composing it with a pretrained classifier. To get the
correct class-conditional density ratio that can be used for generation, we
should therefore compute

$r_{\theta}({\mathbf{x}}_{t}|y=n)=\frac{1}{N}r_{\theta}({\mathbf{x}}_{t})p(y=n|{\mathbf{x}}_{t})^{-1}$
(29)

and use this conditional density ratio estimator in our sampling method of
Algorithm Eq. 2.

##  Appendix F Experimental Details

Table 3: Network structures for the density ratio estimator
$r_{\theta}({\mathbf{x}})$.

CIFAR10
---
3$\times$3 Conv2d, 128
3 $\times$ ResBlock 128
ResBlock Down 256
2 $\times$ ResBlock 256
ResBlock Down 256
2 $\times$ ResBlock 256
ResBlock Down 256
2 $\times$ ResBlock 256
Global Mean Pooling
Dense $\rightarrow$ 1

CelebA 64
---
3$\times$3 Conv2d, 64
ResBlock Down 64
ResBlock Down 128
ResBlock 128
ResBlock Down 256
ResBlock 256
ResBlock Down 256
ResBlock 256
Global Mean Pooling
Dense $\rightarrow$ 1

LSUN Church 128
---
3$\times$3 Conv2d, 64
ResBlock Down 64
ResBlock Down 128
ResBlock Down 128
ResBlock 128
ResBlock Down 256
ResBlock 256
ResBlock Down 256
ResBlock 256
Global Mean Pooling
Dense $\rightarrow$ 1

Cat2dog 64
---
3$\times$3 Conv2d, 64
ResBlock Down 64
ResBlock Down 128
Self Attention 128
ResBlock Down 128
ResBlock Down 256
Global Mean Pooling
Dense $\rightarrow$ 1

#### Unconditional Image Generation.

For all datasets, we perform random horizontal flip as a form of data
augmentation and normalize pixel values to the range [-1, 1]. For CelebA, we
center crop the image to 140$\times$140 before resizing to 64$\times$64\. For
LSUN Church, we resize the image to 128$\times$128 directly. We use the same
training hyperparameters across all three datasets, which is as follows. We
train all models with 100000 training steps with the Adam optimizer with a
batch size of 64, excepet for CIFAR10 with uniform prior where we trained for
120000 steps. We use a learning rate of $1\times 10^{-4}$ and decay twice by a
factor of 0.1 when there are 20000 and 10000 remaining steps. We set the
number of flow steps to $K=100$ at training time. When sampling, we set
$\kappa=20$ for all DDP experiments and $\kappa=10$ for UP experiments. Other
flow hyperparameters include a step size of $\eta=3$ and noise factor
$\nu=10^{-2}$. The specific ResNet architectures are given in Table 3. We
update model weights using an exponential moving average [31] given by
$\bm{\theta}^{\prime}\leftarrow m\bm{\theta}^{\prime}+(1-m)\bm{\theta}_{i}$,
where $\bm{\theta}_{i}$ is the parameters of the model at the $i$-th training
step, and $\bm{\theta}^{\prime}$ is an independent copy of the parameters that
we save and use for evaluation. We set $m=0.998$. We use the LeakyReLU
activation for our networks with a negative slope of 0.2. We experimented with
spectral normalization and self attention layers for unconditional image
generation, but found that training was stable enough such that they were not
worth the added computational cost. The FID results in Table 1 are obtained by
generating 50000 images, and testing the results against the training set for
both CIFAR10 and CelebA.

#### Conditional Generation with Robust Classifier.

The unconditional model used for conditional generation is the same model
obtained from the section above. The pretrained robust classifier checkpoint
is obtained from the robustness111https://github.com/MadryLab/robustness
Python library [7]. It is based on a ResNet50 architecture and is trained with
$L_{2}$-norm perturbations of $\varepsilon=1$.

We choose the LR-KL variant for our results in Fig. 6. This means that our
conditional flow is given by

$\displaystyle{\mathbf{x}}_{k+1}$
$\displaystyle={\mathbf{x}}_{k}-2\eta\nabla_{{\mathbf{x}}}\log\left(r_{\theta}({\mathbf{x}}_{k})*\frac{1}{N}p(y=n|{\mathbf{x}}_{k})^{-1}\right)+\nu\bm{\xi}_{k}$
(30) $\displaystyle={\mathbf{x}}_{k}-2\eta\nabla_{{\mathbf{x}}}\left(\log
r_{\theta}({\mathbf{x}}_{k})-\phi\log
p(y=n|{\mathbf{x}}_{k})\right)+\nu\bm{\xi}_{k}$ (31)

where in the second line we introduce $\phi$ as a parameter that scales the
magnitude of the classifier’s gradients so they are comparable to the
magnitude of FDRL’s gradients. We use $\phi=0.1$.

#### Unpaired Image-to-image Translation.

The Cat2dog dataset contains 871 Birman cat images and 1364 Samoyed and Husky
dog images. 100 of each are set aside as test images. We first resize the
images to 84$\times$84 before center cropping to 64$\times$64\. Due to the
relatively small size of the dataset, we use a shallower ResNet architecture
as compared to CelebA $64^{2}$ despite the same resolution (Table 3) to
prevent overfitting. We also utilize spectral normalization and a self
attention layer at the 128-channel level to further boost stability. We set
$K=100$ during training and $K=110$ at test time. As we observed the model
tends to diverge late in training, we limit the number of training steps to
40000, with a decay factor of 0.1 applied to the learning rate at steps 20000
and 30000. All other hyperparameters are kept identical to the experiments on
unconditional image generation. We report results for LSIF-$\chi^{2}$,
although we have experimented with LR-KL and found performance to be similar.
The FID result in Table. 4 is obtained by translating the 100 test cat images,
and testing the results against the 100 test dog images.

Table 4: FID scores for image-to-image translation with the Cat2dog dataset. Model | FID $\downarrow$
---|---
FDRL | 108.10
CycleGAN | 51.79

(a)

(b)

(c)

Figure 11: Samples drawn from the data-dependent priors of (a) CIFAR10
$32^{2}$, (b) CelebA $64^{2}$ and (c) LSUN Church $128^{2}$.

##  Appendix G Nearest Neighbors

(a)

(b)

Figure 12: Nearest neighbor images for CIFAR10 as measured by $L_{2}$ distance
in (a) the feature space of an Inception V3 network pretrained on ImageNet and
(b) data space. The column to the left of the red line are samples from FDRL
LSIF-$\chi^{2}$. The images to the right of the line are the 10 nearest
neighbors in the training dataset.

(a)

(b)

Figure 13: Nearest neighbor images for CelebA as measured by $L_{2}$ distance
in (a) the feature space of an Inception V3 network pretrained on ImageNet and
(b) data space. The column to the left of the red line are samples from FDRL
LSIF-$\chi^{2}$. The images to the right of the line are the 10 nearest
neighbors in the training dataset.

##  Appendix H Uncurated Samples FDRL-DDP

Figure 14: Uncurated samples of CIFAR10 LSIF-$\chi^{2}$. Figure 15: Uncurated
samples of CIFAR10 LR-KL. Figure 16: Uncurated samples of CIFAR10 LR-JS.
Figure 17: Uncurated samples of CIFAR10 LR-logD. Figure 18: Uncurated samples
of CelebA LSIF-$\chi^{2}$. Figure 19: Uncurated samples of CelebA LR-KL.
Figure 20: Uncurated samples of LSUN Church LSIF-$\chi^{2}$.

##  Appendix I Uncurated Samples FDRL-UP

(a) LSIF-$\chi^{2}$ uniform prior.

(b) LR-KL uniform prior.

Figure 21: Uncurated CIFAR10 samples with FDRL-UP.


Further correspondence, please send an email to<EMAIL_ADDRESS>
* These authors contributed equally.

# Unsupervised Anomaly Detection of Paranasal Anomalies in the Maxillary Sinus

Debayan Bhattacharya  Finn Behrendt  Benjamin Tobias Becker  Dirk Beyersdorff
Elina Petersen  Marvin Petersen  Bastian Cheng  Dennis Eggert  Christian Betz
Anna Sophie Hoffmann  Alexander Schlaefer

###### Abstract

Deep learning (DL) algorithms can be used to automate paranasal anomaly
detection from Magnetic Resonance Imaging (MRI). However, previous works
relied on supervised learning techniques to distinguish between normal and
abnormal samples. This method limits the type of anomalies that can be
classified as the anomalies need to be present in the training data. Further,
many data points from normal and anomaly class are needed for the model to
achieve satisfactory classification performance. However, experienced
clinicians can segregate between normal samples (healthy maxillary sinus) and
anomalous samples (anomalous maxillary sinus) after looking at a few normal
samples. We mimic the clinicians ability by learning the distribution of
healthy maxillary sinuses using a 3D convolutional auto-encoder (cAE) and its
variant, a 3D variational autoencoder (VAE) architecture and evaluate cAE and
VAE for this task. Concretely, we pose the paranasal anomaly detection as an
unsupervised anomaly detection problem. Thereby, we are able to reduce the
labelling effort of the clinicians as we only use healthy samples during
training. Additionally, we can classify any type of anomaly that differs from
the training distribution. We train our 3D cAE and VAE to learn a latent
representation of healthy maxillary sinus volumes using L1 reconstruction
loss. During inference, we use the reconstruction error to classify between
normal and anomalous maxillary sinuses. We extract sub-volumes from larger
head and neck MRIs and analyse the effect of different fields of view on the
detection performance. Finally, we report which anomalies are easiest and
hardest to classify using our approach. Our results demonstrate the
feasibility of unsupervised detection of paranasal anomalies from MRIs with an
AUPRC of 85% and 80% for cAE and VAE, respectively.

###### keywords:

paranasal anomaly, unsupervised anomaly detection, autoencoder, VAE

##  1 INTRODUCTION

Anomalies occurring in the paranasal sinuses are commonly reported in patients
who undergo neuroradiological assessment of the head using diagnostic imaging
[1]. These incidental findings pose clinical challenges[2] and we have limited
knowledge on the importance of these reported findings on the general
population. To this end, numerous studies have been done to analyse the
significance of these findings [3, 4, 5, 6, 7]. Most of these works are
population studies involving large sample sizes which are manually annotated
by clinicians. In a three year retrospective study, it was observed that
malignant tumours and inverted papillomas were classified as nasal polyps with
a misdiagnosis rate of 5.63% and 8.45% respectively. Therefore, computer aided
diagnosis systems (CADx) [8, 9, 10] have been proposed with the idea of
working in conjunction with the clinician to reduce the misdiagnosis rate of
paranasal anomalies. However, all these works rely on supervised learning and
consider at most one anomaly. Apart from requiring large labelled datasets,
supervised learning models also need labelled data that are representative of
the classes for accurate prediction [11]. Further, the type of anomaly to be
classified has to be decided beforehand. In our case, we consider three
anomalies, namely: (i) mucosal thickening (ii) polyps (iii) cysts. This is
particularly challenging in our case where the considered anomalies are known
to have high intra-class morphological variations [12, 13, 14] and have
unequal occurrences in our dataset.

In light of the aforementioned points, we are motivated to perform
Unsupervsised Anomaly Detection (UAD) using autoencoders. Autoencoders are a
common choice for data compresison and outlier detection [15]. In medical
imaging, brain anomaly detection and segmentation [16] has gained popularity
over the years. The underlying concept of reconstruction-based UAD is that the
autoencoder learns to compress and reconstruct only normal images during
training. The assumption is that during testing, reconstruction errors will be
low for normal images whereas anomalous images will have a large
reconstruction error as the autoencoder will fail to properly reconstruct
anomalous regions. In our case, we use cAE and VAE learn to compress and
reconstruct healthy maxillary sinus (MS) volumes. Through our approach, we
derive multiple benefits, namely: (i) our autoencoder becomes indifferent to
the anomaly distribution thereby allowing detection of more than one anomaly,
(ii) we reduce the labelling effort of the clinicians as the training dataset
does not require anomalous samples, (iii) we are able to generate a heat map
based on the reconstruction error between the original volume and the
reconstructed volume. The heat maps visualize the region of the potential
anomaly and highlight it. This may prove to be beneficial to the clinician
when making a diagnosis.

In summary, our contributions are three-fold. First, we pose the paranasal
anomaly detection problem as a UAD problem and thereby become indifferent to
the anomaly distribution. Second, we systematically evaluate our cAE and VAE
approach on MS volumes with different field of view. Third, we report which
anomalies are easiest and the hardest to classify using our UAD method.

##  2 METHODS

###  2.1 Dataset and Implementation Details

Our labelled dataset consists of head and neck MRIs of 199 patients. Each MRI
is a fluid attenuated inversion recovery (FLAIR) MRI. Our labelled dataset is
part of the Hamburg City Health Study [17]. Out of the 199 patients, 93
patients exhibit one or multiple anomalies in at least left or right MS. Two
Ear, Nose and Throat (ENT) surgeons and one ENT specialised radiologist
confirmed the diagnosis of the observed pathology. We group the 3 anomalies
into a single class called ”anomaly” and the normal MS are categorized into
”normal” class. Altogether, we have 269 normal MS volumes and 130 anomalous MS
volumes . Each MRI has a resolution of 173 x 319 x 319 voxels along the
saggital, coronal and axial directions respectively with each voxel of size
0.53 $mm$ x 0.75 $mm$ x 0.75 $mm$.

Preprocessing : We randomly selected a FLAIR MRI as the fixed MRI and
performed rigid registration on the remaining MRIs. This was followed by
resampling to a dimension of 128 $\times$ 128 $\times$ 128\. We extracted two
sub-volumes from the resampled head and neck MRIs, one for each MS. The
extracted sub-volumes were of sizes: 33 $\times$ 47 $\times$ 45 (small), 46
$\times$ 57 $\times$ 55 (medium) and 53 $\times$ 67 $\times$65 (large). Owing
to the symmetry of left and right MS, we horizontally flipped the coronal
planes of right MS to give it the appearance of left MS for each patient.
Finally, these sub-volumes were reshaped to a standard size of 64 $\times$ 64
$\times$ 64 voxels for the 3D cAE and VAE. All the MS volumes were normalised
to a range of 0 to 1. Our preprocessing pipeline is illustrated in figure 1
(a,b)

Data split: Our training set contains 172 normal MS volumes , validation set
contains 43 normal MS volumes and 52 anomalous MS volumes and test set
contains 54 normal MS volumes and 78 anomalous MS volumes. We perform a three-
fold cross validation split for all our experiments.

Implementation Details: We use PyTorch [18] and PyTorch Lightning [19] for all
our experiments. We use batch size $N$ = 16 and latent dimension $n_{z}$ = 512
for the cAE. We use a learning rate of $1e^{-4}$ and Adam optimizer [20] with
default parameters to train our cAE and VAE. We run all our experiments for
100 epochs.

###  2.2 Deep Learning Methods

Figure 1: (a) Extraction of left and right MS from head and neck MRI (b)
Flipping of the coronal plane of right MS (c) Cyst in the right MS (d) Polyp
in the left MS (e) Cyst in the left MS (f) MS (small) showing mucosal
thickening (g) MS (medium) showing polyp (h) MS (large) showing cyst (i) TOP:
Our cAE architecture with latent vector of size 512 used as bottleneck.
BOTTOM: Our VAE architecture. In both the cAE and VAE decoders, we perform 3D
convolution followed by trilinear upsampling. ReLU is used as non-linear
activation function in both the networks. (j)-(m) Images extracted from
original , reconstructed , intensity difference and heat map volumes. (j) is
an image extracted from normal MS volume. (k),(l),(m) are mucosal thickening,
polyps and cyst anomalies respectively. The red markings denote anomalies and
green circles denote normal MS.

Similar to Bengs et al [21], we extend our architecture to 3D as it has shown
to improve the detection performance in 3D MRI scans. Our cAE and VAE
architectures are shown in figure 1 (i) . There are two stages to our
approach. First, we train our cAE and VAE to learn the distribution of $X_{h}$
where $X_{h}$ represents healthy MS volumes. For training our cAE, we use the
L1 reconstruction loss as described below:

$L_{1}=\sum_{k=1}^{N}|x^{k}-\hat{x}^{k}|$ (1)

Here, $x^{k}$ and $\hat{x}^{k}$ denote the $k$-th MS volume and reconstructed
MS volume respectively. N denotes the mini-batch size. For training our VAE,
we use L1 reconstruction loss and KL Divergence. While training, the VAE
learns a mean $\mu_{z}$ and variance $\sigma_{z}$ from which a sample is drawn
and reconstructed. The loss used to train our VAE is shown below:

$L_{VAE}=L_{1}+\lambda_{KL}D_{KL}(q(z|x)||p(z))$ (2)

Here, $D_{KL}(.||.)$ represents the Kullback–Leibler divergence between the
parameterized latent distribution $q(z|x)$ $\sim$ $N(\mu_{z},\sigma_{z})$ and
the prior p(z) which follows a multivariate normal distribution.
$z\in\mathbb{R}^{n_{z}}$ represents the latent vector. $\lambda_{KL}$ is a
Lagrangian multiplier and we have set it to 1 for our experiments. VAE
projects the the input MS volume to $q(z|x)$ and KL-Divergence loss attempts
to bring it close to a prior $p(z)$ [22].

Second, we use the trained cAE and VAE to reconstruct the MS volumes in the
validation set. For each MS volume, we calculate the L1 and L2 reconstruction
loss denoted as $t_{L1}$ and $t_{L2}$ respectively. We choose the optimal
thresholds by plotting the precision recall curve and select the threshold
with the highest F1 score. During inference, $\hat{x}^{k}$ with $L_{1}>$
$t_{L1}$ and $L_{2}>$ $t_{L2}$ is classified as anomalous MS volume. Here,
$L_{2}$ is the L2 reconstruction loss defined as follows:

$L_{2}=\sum_{k=1}^{N}(x^{k}-\hat{x}^{k})^{2}$ (3)

Additionally, for the MS volumes classified as anomalous, a further analysis
is done by calculating voxel-wise intensity difference
$D_{k}=|x^{k}-\hat{x}^{k}|$ after which a median filter of kernel size 5 is
applied on it to remove sporadic reconstruction errors. Finally, a heat map is
rendered for the individual slices along the coronal, axial and sagittal
planes as shown in figure 1 (j, k, l, m). The regions in red denote the
regions which the cAE failed to reconstruct. Our qualitative results indicate
an overlap between the anomalous regions and the poorly reconstructed regions
of MS volumes.

##  3 Results

Table 1: Anomaly Detection Performance on two thresholds $t_{L1}$ and $t_{L2}$
where positive labels are assigned to the anomalous class.

Method | MS Size | Precision | Recall | F1 | AUPRC
---|---|---|---|---|---
|  | $t_{L1}$ | $t_{L2}$ | $t_{L1}$ | $t_{L2}$ | $t_{L1}$ | $t_{L2}$ | $t_{L1}$ | $t_{L2}$
VAE | small | 0.69 | 0.76 | 0.63 | 0.62 | 0.64 | 0.68 | 0.76 | 0.80
VAE | medium | 0.63 | 0.63 | 0.84 | 0.91 | 0.72 | 0.75 | 0.70 | 0.75
VAE | large | 0.61 | 0.64 | 0.81 | 0.86 | 0.70 | 0.73 | 0.65 | 0.69
cAE | small | 0.75 | 0.81 | 0.74 | 0.66 | 0.74 | 0.73 | 0.81 | 0.83
cAE | medium | 0.73 | 0.77 | 0.62 | 0.74 | 0.67 | 0.75 | 0.80 | 0.85
cAE | large | 0.68 | 0.74 | 0.82 | 0.73 | 0.74 | 0.73 | 0.73 | 0.78

Table 2: Accuracy per anomaly on the test set reported in percentage (%).
Here, # refers to the number of correctly classified samples divided by the
total samples from a particular category.

Method | MS Size | Normal (%/#) | Mucosal Thickening (%/#) | Polyps (%/#) | Cysts (%/#)
---|---|---|---|---|---
VAE | small | 0.48 (26/54) | 0.75 (22/29) | 0.82 (28/34) | 0.73 (11/15)
cAE | medium | 0.61 (33/54) | 0.62 (18/29) | 0.91 (31/34) | 0.8 (12/15)

The anomaly detection performance is shown in table 1. We consider the Area
Under Precision Recall Curve (AUPRC) to evaluate our classifiers as we have an
unbalanced test set. We rank the classifiers based on AUPRC. We report the
mean values of the mentioned metrics. In terms of AUPRC, the VAE that uses
small MS volume and cAE that uses medium MS volume are the best performing
classifiers.

From table 1, we notice that the all the VAEs have relatively lower AUPRC in
comparison to cAEs. Furthermore, our results show that considering L2 loss
when computing the anomaly score and using $t_{L2}$ as threshold leads to
better performance for all cAEs and VAEs. In table 2, we report the accuracy
per anomaly in terms of percentage and number of occurances for the best
performing VAE and cAE from table 1. We observe that polyps are the easiest to
classify, followed by cysts. Mucosal thickening anomaly is the hardest to
classify for both VAE and cAE.

##  4 Discussion and Conclusion

From the results in table 1 we observe that using $t_{L2}$ as threshold leads
to better performance. This can be attributed to the fact that L2 loss
penalizes voxel-wise intesity outliers more heavily than L1 loss. Therefore,
even small regions of poor reconstruction in the MS volume can amount to high
overall reconstruction error. Additionally, our accuracy percentage per
anomaly reported in table 2 shows that mucosal thickening anomalies are the
most difficult to classify. We believe this to be the case because unlike
polyps and cysts which occur as visible masses in the MRI (See figure 1
(l),(m)), mucosal thickening anomalies have more subtle appearances as these
anomalies mostly cause inflammation of the mucosal walls. Therefore, unless
the inflammations are too noticeable, they almost have the appearance of a
healthy MS. We also observe that classifiers using small and medium MS volumes
have the best AUPRCs. This can be attributed to the fact that cropping large
volumes lead to inclusion of unnecessary surrounding anatomical structures
outside of the MS. These additional anatomical structures effect the
reconstruction error and thereby, we end up selecting sub-optimal thresholds.
Our work has some limitations, one being the accuracy of healthy MS volume
detection needs to be higher. We think this can be achieved by better
localisation and cropping strategies of the MS in the head and neck MRI and by
labelling more healthy MS volumes. Second, we have not experimented with
autoencoders with skip connections or with more parameters and studied its
effect on reconstruction error.

In conclusion, we evaluate UAD for paranasal anomaly detection. Previous
methods [8, 9, 10] have used supervised learning methods and as a result are
constrained to classify the anomalies that are included in the training
distribution. Through our UAD approach, our models learn the healthy MS volume
distribution $X_{h}$ thereby reducing the labelling effort of the clinicians.
Also, we are able to detect multiple anomalies. Further, we render heat maps
of poor reconstruction. This can provide valuable insights to clinicians while
making a diagnosis.

## References

  * [1] Wilson, R., Kuan Kok, H., Fortescue-Webb, D., Doody, O., Buckley, O., and Torreggiani, W. C., “Prevalence and seasonal variation of incidental mri paranasal inflammatory changes in an asymptomatic irish population,” Irish medical journal 110(9), 641 (2017).
  * [2] Hansen, A. G., Helvik, A.-S., Nordgård, S., Bugten, V., Stovner, L. J., Håberg, A. K., Gårseth, M., and Eggesbø, H. B., “Incidental findings in mri of the paranasal sinuses in adults: a population-based study (hunt mri),” BMC ear, nose, and throat disorders 14(1), 13 (2014).
  * [3] Tarp, B., Fiirgaard, B., Christensen, T., Jensen, J. J., and Black, F. T., “The prevalence and significance of incidental paranasal sinus abnormalities on mri,” Rhinology 38(1), 33–38 (2000).
  * [4] Rak, K. M., Newell, J. D., Yakes, W. F., Damiano, M. A., and Luethke, J. M., “Paranasal sinuses on mr images of the brain: significance of mucosal thickening,” AJR. American journal of roentgenology 156(2), 381–384 (1991).
  * [5] Stenner, M. and Rudack, C., “Diseases of the nose and paranasal sinuses in child,” GMS current topics in otorhinolaryngology, head and neck surgery 13, Doc10 (2014).
  * [6] Rege, I. C. C., Sousa, T. O., Leles, C. R., and Mendonça, E. F., “Occurrence of maxillary sinus abnormalities detected by cone beam ct in asymptomatic patients,” BMC oral health 12, 30 (2012).
  * [7] Cooke, L. D. and Hadley, D. M., “Mri of the paranasal sinuses: incidental abnormalities and their relationship to symptoms,” The Journal of laryngology and otology 105(4), 278–281 (1991).
  * [8] Kim, Y., Lee, K. J., Sunwoo, L., Choi, D., Nam, C.-M., Cho, J., Kim, J., Bae, Y. J., Yoo, R.-E., Choi, B. S., Jung, C., and Kim, J. H., “Deep learning in diagnosis of maxillary sinusitis using conventional radiography,” Investigative radiology 54(1), 7–15 (2019).
  * [9] Jeon, Y., Lee, K., Sunwoo, L., Choi, D., Oh, D. Y., Lee, K. J., Kim, Y., Kim, J.-W., Cho, S. J., Baik, S. H., Yoo, R.-E., Bae, Y. J., Choi, B. S., Jung, C., and Kim, J. H., “Deep learning for diagnosis of paranasal sinusitis using multi-view radiographs,” Diagnostics (Basel, Switzerland) 11(2) (2021).
  * [10] Liu, G. S., Yang, A., Kim, D., Hojel, A., Voevodsky, D., Wang, J., Tong, C. C. L., Ungerer, H., Palmer, J. N., Kohanski, M. A., Nayak, J. V., Hwang, P. H., Adappa, N. D., and Patel, Z. M., “Deep learning classification of inverted papilloma malignant transformation using 3d convolutional neural networks and magnetic resonance imaging,” International forum of allergy & rhinology (2022).
  * [11] Althnian, A., AlSaeed, D., Al-Baity, H., Samha, A., Dris, A. B., Alzakari, N., Abou Elwafa, A., and Kurdi, H., “Impact of dataset size on classification performance: An empirical evaluation in the medical domain,” Applied Sciences 11(2) (2021).
  * [12] Tos, M., Larsen, P. L., Larsen, K., and Cayé-Thomasen, P., [Nasal Polyps ], 103–125, Springer Berlin Heidelberg, Berlin, Heidelberg (2000).
  * [13] Janner, S. F. M., Caversaccio, M. D., Dubach, P., Sendi, P., Buser, D., and Bornstein, M. M., “Characteristics and dimensions of the schneiderian membrane: a radiographic analysis using cone beam computed tomography in patients referred for dental implant surgery in the posterior maxilla,” Clin Oral Implants Res 22, 1446–1453 (Mar. 2011).
  * [14] Hung, K., Hui, L., Yeung, A. W. K., Wu, Y., Hsung, R. T.-C., and Bornstein, M. M., “Volumetric analysis of mucous retention cysts in the maxillary sinus: A retrospective study using cone-beam computed tomography,” Imaging Sci Dent 51, 117–127 (Jan. 2021).
  * [15] Pang, G., Shen, C., Cao, L., and Hengel, A. V. D., “Deep learning for anomaly detection: A review,” ACM Comput. Surv. 54 (mar 2021).
  * [16] Baur, C., Denner, S., Wiestler, B., Navab, N., and Albarqouni, S., “Autoencoders for unsupervised anomaly segmentation in brain mr images: A comparative study,” Medical Image Analysis , 101952 (2021).
  * [17] Jagodzinski, A., Johansen, C., Koch-Gromus, U., Aarabi, G., Adam, G., Anders, S., Augustin, M., der Kellen, R. B., Beikler, T., Behrendt, C.-A., Betz, C. S., Bokemeyer, C., Borof, K., Briken, P., Busch, C.-J., Büchel, C., Brassen, S., Debus, E. S., Eggers, L., Fiehler, J., Gallinat, J., Gellißen, S., Gerloff, C., Girdauskas, E., Gosau, M., Graefen, M., Härter, M., Harth, V., Heidemann, C., Heydecke, G., Huber, T. B., Hussein, Y., Kampf, M. O., von dem Knesebeck, O., Konnopka, A., König, H.-H., Kromer, R., Kubisch, C., Kühn, S., Loges, S., Löwe, B., Lund, G., Meyer, C., Nagel, L., Nienhaus, A., Pantel, K., Petersen, E., Püschel, K., Reichenspurner, H., Sauter, G., Scherer, M., Scherschel, K., Schiffner, U., Schnabel, R. B., Schulz, H., Smeets, R., Sokalskis, V., Spitzer, M. S., Terschüren, C., Thederan, I., Thoma, T., Thomalla, G., Waschki, B., Wegscheider, K., Wenzel, J.-P., Wiese, S., Zyriax, B.-C., Zeller, T., and Blankenberg, S., “Rationale and design of the hamburg city health study,” European Journal of Epidemiology 35, 169–181 (Feb 2020).
  * [18] Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J., Chanan, G., Killeen, T., Lin, Z., Gimelshein, N., Antiga, L., Desmaison, A., Köpf, A., Yang, E., DeVito, Z., Raison, M., Tejani, A., Chilamkurthy, S., Steiner, B., Fang, L., Bai, J., and Chintala, S., “Pytorch: An imperative style, high-performance deep learning library.”
  * [19] Falcon et al., W., “Pytorch lightning,” GitHub. Note: https://github.com/PyTorchLightning/pytorch-lightning 3 (2019).
  * [20] Kingma, D. P. and Ba, J., “Adam: A method for stochastic optimization,” (2014).
  * [21] Bengs, M., Behrendt, F., Krüger, J., Opfer, R., and Schlaefer, A., “Three-dimensional deep learning with spatial erasing for unsupervised anomaly segmentation in brain mri,” CARS , 1–11 (2021).
  * [22] Kingma, D. P. and Welling, M., “Auto-encoding variational bayes,” (2013).