# Efficient CDF Approximations for Normalizing Flows

Chandramouli Sastry  Andreas Lehrmann  Marcus Brubaker  Alexander Radovic

###### Abstract

Normalizing flows model a complex target distribution in terms of a bijective
transform operating on a simple base distribution. As such, they enable
tractable computation of a number of important statistical quantities,
particularly likelihoods and samples. Despite these appealing properties, the
computation of more complex inference tasks, such as the cumulative
distribution function (CDF) over a complex region (e.g., a polytope) remains
challenging. Traditional CDF approximations using Monte-Carlo techniques are
unbiased but have unbounded variance and low sample efficiency. Instead, we
build upon the diffeomorphic properties of normalizing flows and leverage the
divergence theorem to estimate the CDF over a closed region in target space in
terms of the flux across its _boundary_ , as induced by the normalizing flow.
We describe both deterministic and stochastic instances of this estimator:
while the deterministic variant iteratively improves the estimate by
strategically subdividing the boundary, the stochastic variant provides
unbiased estimates. Our experiments on popular flow architectures and UCI
benchmark datasets show a marked improvement in sample efficiency as compared
to traditional estimators.

Machine Learning, ICML

##  1 Introduction

Figure 1: Overview. We propose an efficient estimator to compute the
cumulative density $P$ induced by a normalizing flow $f$ over a complex region
$\bf V$. We show that $P$ corresponds to integrating over the boundary
$\partial{\bf W}=g(\partial{\bf V})$ in a uniform base space w.r.t. a vector
field ${\bf F}$ with constant divergence $1$, which can be expressed as an
integration over the boundary $\partial{\bf V}$ in the target space w.r.t. an
equivalent vector field ${\bf G}$.

Normalizing Flows (Kobyzev et al., 2020; Papamakarios et al., 2021) are a form
of generative model which constructs tractable probability distributions
through invertible, differentiable transformations, _i.e._ , diffeomorphisms.
They admit both efficient and exact density evaluation and sampling and, as
valid probability distributions, theoretically support a range of other
probabilistic inference tasks. However, some inference tasks remain
computationally challenging. Here we consider the task of computing cumulative
densities over arbitrary closed regions of distributions represented by a
normalizing flow. Cumulative densities are a fundamental statistical measure.
They answer the question: what is the probability of a sample in this range of
values? Example applications are particularly common in the domain of
uncertainty estimation (Liu et al., 2019; Mazaheri et al., 2020; Abdar et al.,
2021), and the evaluation of risk in financial settings (Richardson et al.,
1997; Mehta et al., 2012). Other recent examples include distributional
reinforcement learning (Sun et al., 2021), few-shot learning (Ridgeway &
Mozer, 2018), and alternative methods for training normalizing flows (Dai &
Seljak, 2021).

Despite the importance of cumulative densities, current techniques for
computing them with normalizing flows are restricted to traditional Monte-
Carlo based estimators which, while unbiased, can have unbounded variance and
poor computational efficiency. Further, such approaches fail to exploit the
inherent structure of normalizing flows, namely their construction as a
diffeomorphic transformation.

Normalizing Flows (Kobyzev et al., 2020; Papamakarios et al., 2021) are
explicit likelihood models with tractable distributions that enable efficient
and exact density evaluation and sampling. As valid probability distributions,
normalizing flows further theoretically support a range of probabilistic
inference tasks and inverse problems. Further, recent work has connected
normalizing flows to diffusion models (DBLP:journals/corr/abs-2011-13456), a
SOTA approach for generative modeling, allowing trained diffusion models to be
interpreted as a consistent continuous Normalizing Flow. However, not all
inference problems currently have computationally efficient implementations.

In this paper, we consider the task of inferring cumulative densities – over
arbitrary closed regions – as learnt by a normalizing flow. Cumulative
densities are a fundamental statistical measure, and represent the answer to a
natural question, what is the probability of a sample in this range of values?
Example applications are particularly common in the domain of uncertainty
estimation (Liu et al., 2019; Mazaheri et al., 2020; Abdar et al., 2021), and
include evaluations of risk in the stock market (Richardson et al., 1997;
Mehta et al., 2012). Other recent examples include distributional
reinforcement learning (Sun et al., 2021), few-shot learning (Ridgeway &
Mozer, 2018), and alternative methods for training normalizing flows (Dai &
Seljak, 2021). Despite this, the popular (and only) techniques for estimating
CDFs in Normalizing Flows are traditional Monte-Carlo based estimators which
while being unbiased can have unbounded variance and low sample-efficiency.

In this paper we describe a novel estimator that exploits the unique
characteristics of normalizing flows to efficiently estimate the cumulative
density in closed regions. To do this we exploit the homeomorphic property of
normalizing flow transforms to relate the cumulative density in the target
space to volume in the base space. Next we adapt the divergence theorem to
show that we only need to consider the boundary of the region in the target
space. We analyze the resulting theoretical estimator to identify best-case
performance scenarios and use this analysis to motivate an adaptive,
approximate estimator. To practically realize this estimator and evaluate its
performance, we explore its application to computation of cumulative density
over convex regions.

##### Contributions.

The contributions of this paper are as follows. First, we derive a novel
formulation of the cumulative density of a normalizing flow by relating it to
the volume in a (uniform) base space and the boundary integrals in both base
and data space based on the divergence theorem. Second, based on this we
describe an unbiased, stochastic estimator for cumulative densities over a
closed region. Next, we analyze this estimator to build an adaptive,
deterministic approximation for cumulative densities over a closed region by
strategically adding points along the boundary of the region. Finally, we
describe a comprehensive evaluation protocol to evaluate CDF estimators for
normalizing flows as a function of region sizes, architecture families and
capacities. The code to reproduce our results, including training popular
normalizing flow architectures, approximating cumulative densities with the
proposed adaptive boundary estimator and other baseline methods is publicly
available.111 https://github.com/BorealisAI/nflow-cdf-approximations

##  2 Related Work

Our proposed method is principally related to two bodies of existing work: (1)
algorithms for efficient estimation of cumulative densities for popular
distribution functions; and (2) modifications of popular normalizing flows
architectures to allow for more efficient cumulative density estimation.

Classic Algorithms. Computing cumulative densities is not tractable in higher
dimensions even for distributions such as the multivariate Gaussian which have
closed-form density functions. Consequently, approximations to volume
integrals/cumulative density functions have been an important research
problem: for example, (Genz & Bretz, 2009) discuss methods for approximating
integrals over multivariate Gaussian and t-distributions. (Botev & l’Ecuyer,
2015) and (Botev, 2017) propose a tilting method for i.i.d. sampling and
approximating cumulative densities over convex polytopes in multivariate
Gaussian and t-distributions respectively. Approximating integrals over 1D
distributions is also challenging for certain distributions; see (Lange, 1999)
for a discussion on quadrature methods and chapters 4-6 of (Press et al.,
1988) for a discussion on numerical methods for estimating 1D integrals.
(Cunningham et al., 2013) extensively study Expectation-Propagation (Minka,
2001) as an approximate integration method for computing cumulative densities
over polyhedral or hyperrectangular regions in Gaussian space; if we consider
the base space of a normalizing flow to be Gaussian, one could consider –
based on our discussion in the methods section – applying this method to the
region in Gaussian space corresponding to the region $\bf V$ in the target
space; however, it is not easy to ensure that polyhedral regions in target
space remain polyhedral in base space without restricting the expressive power
of the flow transform.

Normalizing Flows. (Cundy & Ermon, 2020) consider the general problem of
estimating integrals and propose an approximate inference procedure that
allows explicit control of the bias-variance tradeoff; the method involves
first creating partitions of the base uniform space of the normalizing flow
and then training a separate variational approximator for each partition: when
one partition is used for the entire region, the algorithm reduces to the
standard variational approximation resulting in a biased low-variance estimate
– however, for an exact estimate, the number of partitions have to be
increased and the algorithm reduces to Monte-Carlo sampling. Furthermore, it
is not easy to determine the partitioning strategy to arrive at the exact
estimate with as few samples as possible. Similar to our problem setting,
(Liang et al., 2020) consider the problem of approximating range densities in
a tabular database using auto-regressive models: in order to efficiently
marginalize over data dimensions, they train the model by randomly masking few
data dimensions, which is referred to as variable-skipping. At inference time,
they only need to sample the variables whose ranges are specified to form the
estimate; even though the training objective is indeed the maximum-likelihood
objective, the estimate obtained from the model trained with variable-skipping
technique is only an approximate value of the range densities. (Chilinski &
Silva, 2020) suggest an alternative to normalizing flows with some
similarities, where CDF estimates are available through a simple forward pass.
However, this new model family has some limitations in comparison to
normalizing flows, most notably sampling from the learned distribution is not
supported.

##  3 Background

Before we describe our approach we provide a brief review of normalizing flows
and traditional sampling-based methods of estimating their cumulative
densities.

###  3.1 Normalizing Flows

Normalizing flows transform a random variable $\bf Y$ with a simple base
distribution $p_{\bf Y}({\bf y})$ into a random variable $\bf X$ with a
complex target distribution $p_{\bf X}({\bf x})$ using a bijective mapping
$f:\mathbb{R}^{d}\longrightarrow\mathbb{R}^{d}$. Base and target space are
connected by the change-of-variables formula

$p_{\bf X}({\bf x})=p_{\bf Y}({\bf
y})\left|\det\frac{\operatorname{d}\\!{g}}{\operatorname{d}\\!{\bf
x}}\right|,$ (1)

where $g$ is the inverse of $f$. We set ${\bf Y}\sim\mathcal{U}_{[0,1]^{d}}$,
which has the important property $P({\bf y}\in{\bf W})=\mathrm{vol}({\bf W})$
for any compact set ${\bf W}\subset[0,1]^{d}$, i.e., the cumulative density
equals the enclosed volume.222We use $\mathcal{U}_{S}$ to denote the uniform
distribution with support $S$. For pre-existing flows $f^{\prime}$ with non-
uniform (e.g., Gaussian) base distribution $p_{\bf Y^{\prime}}({\bf
y}^{\prime})$, we can simply set $f:=f^{\prime}(F^{-1}_{{\bf Y}^{\prime}}({\bf
y}))$ to align them with the case above. Sampling from a normalizing flow
involves first sampling ${\bf y}\sim p_{\bf Y}$ and then applying the flow
transform $f$ to obtain ${\bf x}=f({\bf y})$.

##### Objective.

Given a compact set of interest ${\bf V}\subset\mathbb{R}^{d}$ in target
space, our goal will be the efficient computation of

$P({\bf x}\in{\bf V})=\int_{\bf V}p_{\bf X}({\bf x})\,\operatorname{d}\\!{\bf
x}.$ (2)

###  3.2 Monte-Carlo Estimation

If the volume $\bf V$ admits a tractable membership function, one could
estimate $P({\bf x}\in{\bf V})$ from $N$ random samples from the model as the
fraction of points falling into $\bf V$; we refer to this method as the
_Monte-Carlo (MC) estimate_. Similarly, if the volume $\bf V$ admits a
uniformly distributed random variable ${\bf U}\sim\mathcal{U}_{\bf V}$, one
could estimate $P({\bf x}\in{\bf V})$ as

$\scalebox{1.0}{ $\int_{\bf V}p_{\bf X}({\bf x})\,\operatorname{d}\\!{\bf
x}=\int_{\bf V}p_{\bf U}({\bf x})\cdot\frac{p_{\bf X}({\bf x})}{p_{\bf U}({\bf
x})}\,\operatorname{d}\\!{\bf x}=\mathbb{E}_{{\bf x}\sim p_{\bf
U}}\left[\frac{p_{\bf X}({\bf x})}{p_{\bf U}({\bf x})}\right]$},$

which we refer to as the _Importance Sampling (IS) estimate_. While both of
the described estimators are unbiased, they have unbounded variances for
smaller sample sizes and may require several runs to obtain accurate results
(Cunningham et al., 2013).

##  4 Efficient Estimation of Cumulative Densities

Equipped with these tools, we will now describe our technique to efficiently
compute the cumulative density of Eq.(2). We proceed in three steps: first, we
leverage the divergence theorem to relate $P({\bf x}\in{\bf V})$ to a flux
through the _transformed_ boundary $g(\partial{\bf V})$ of $\bf V$ in base
space (Section 4.1). Then, we introduce an equivalent flux through the
_untransformed_ boundary $\partial{\bf V}$ of ${\bf V}$ to enable direct
integration in target space (Section 4.2). Finally, we propose an iterative
protocol that enables fine-grained control over the quality of the
approximation (Section 4.3).

###  4.1 Cumulative Density as Boundary Flux

We will first show that computing the cumulative density over $\bf V$ in
target space is equivalent to computing the volume of ${\bf W}:=g({\bf V})$ in
base space.

##### Lemma 1.

Let ${\bf V}\subset\mathbb{R}^{d}$ be compact and ${\bf W}:=g({\bf V})$. Then
$P({\bf x}\in{\bf V})=\mathrm{vol}({\bf W})$.

_Proof._ Noting that $p({\bf x})=p({\bf
y})\left|\det\frac{\operatorname{d}\\!{g}}{\operatorname{d}\\!{\bf x}}\right|$
and $\operatorname{d}\\!{\bf
y}=\left|\det\frac{\operatorname{d}\\!{g}}{\operatorname{d}\\!{\bf
x}}\right|\operatorname{d}\\!{\bf x}$, we have

$\begin{split}P({\bf x}\in{\bf V})&=\int_{\bf V}p({\bf
x})\,\operatorname{d}\\!{\bf x}=\int_{\bf V}p({\bf
y})\left|\det\frac{\operatorname{d}\\!{g}}{\operatorname{d}\\!{\bf
x}}\right|\,\operatorname{d}\\!{\bf x}\\\ &=\int_{\bf W}p({\bf
y})\,\operatorname{d}\\!{\bf y}=\int_{\bf W}\operatorname{d}\\!{\bf
y}=\mathrm{vol}({\bf W})\ \raisebox{-2.84526pt}{\qed}\end{split}$ (3)

We can control the properties of ${\bf W}$ through $f$:

##### Definition 1 (Diffeomorphism).

A differentiable bijection $f$ on $\mathbb{R}^{d}$ is called a diffeomorphism,
if it has a differentiable inverse $g=f^{-1}$.

Importantly, diffeomorphisms map points on the boundary to points on the
boundary and points in the interior to points in the interior (Armstrong,
2013), implying $\partial{\bf W}=g(\partial{\bf V})$ for any diffeomorphic
flow $f$. This includes popular flow architectures like Glow (Kingma &
Dhariwal, 2018), MAF (Papamakarios et al., 2017), and FFJORD (Grathwohl et
al., 2019) and ensures that $\bf W$ meets the requirements for the following
theorem relating volume integrals to surface integrals.

##### Theorem 1 (Divergence Theorem).

Let ${\bf T}\subset\mathbb{R}^{d}$ be compact with piecewise smooth boundary
$\partial\bf T$. Given a vector field $\bf B$, the volume integral of the
divergence $\nabla\cdot\bf B$ over $\bf T$ and the surface integral of $\bf B$
over $\partial\bf T$ are related by

$\int_{\bf T}(\nabla\cdot{\bf B})\>\operatorname{d}\\!{\bf
T}=\int_{\partial\bf T}{\bf B}\cdot{\bf n}\>\operatorname{d}\\!{(\partial\bf
T)},$ (4)

where ${\bf n}$ is the outward-facing unit normal.

_Proof._ (Wade, 2017) ∎

Setting ${\bf A}:=\partial{\bf W}$, an application of the divergence theorem
using the vector field ${\bf F}({\bf x})=d^{-1}{\bf x}$ yields
$\mathrm{div}({\bf F})=1$ and thus

$\mathrm{vol}({\bf W})=\int_{\bf A}{\bf F}\cdot{\bf
n}\,\operatorname{d}\\!{\bf A},$ (5)

i.e., we can compute the volume of $\bf W$ in terms of a flux through the
boundary of ${\bf W}$.333We note that $\bf F$ is not unique and other choices
are possible. Now, without loss of generality, assume that $\bf V$ is given by
a simplicial polytope. The boundary of such a polytope is defined by a set of
$(d\\!-\\!1)$-simplices; for example, the boundary of a simplicial polytope in
2D is a set of $1$-simplices (i.e., line segments). We can approximate the
boundary $\bf A$ by flowing the vertices of the $(d\\!-\\!1)$-simplices
defining the boundary of $\bf V$ through $g$. Denoting the resulting set of
simplices by $\\{{\bf A}_{i}\\}_{i}$, we have

$\mathrm{vol}({\bf W})\approx\sum_{{\bf A}_{i}}{\bf F}({\bf
A}_{i}^{\\!(c)})\cdot{\bf A}_{i}^{\\!(n)},$ (6)

where ${\bf A}_{i}$ is the $i$-th simplex, ${\bf A}_{i}^{\\!(c)}$ its
centroid, and ${\bf A}_{i}^{\\!(n)}$ its outward-facing normal. The
approximation error is rooted in the assumption that the true boundary patch
in ${\bf A}$ corresponding to a transformed simplex ${\bf A}_{i}$ is linear.
The larger the set of simplices approximating the boundary $\bf A$ the better
the estimate; see Fig. 2 for a summary of this section.

Figure 2: Estimation in Base Space. A piecewise linear boundary $\bf S$
defined by points $\\{p_{i}\\}_{i}$ in target space has a complex shape $\bf
A$ in base space. The cumulative density enclosed by $\bf S$ is equivalent to
the boundary flux of $\bf F$ through $\bf A$ (Eq.(5)). As an approximation, we
can instead compute the boundary flux of $\bf F$ through the piecewise linear
boundary $\\{{\bf A}_{i}\\}_{i}$ defined by $\\{g(p_{i})\\}_{i}$ (Eq.(6)).

###  4.2 Stochastic Estimation in Target Space

One drawback of the method described in Eq.(6) is that it requires first
constructing the set of simplices $\\{{\bf A}_{i}\\}_{i}$ in base space before
we can apply the divergence theorem. In particular, we would need to recompute
the centroids ${\bf A}_{i}^{\\!(c)}$ and normals ${\bf A}_{i}^{\\!(n)}$ every
time $\\{{\bf A}_{i}\\}_{i}$ changes. Preferable is a direct integration over
the boundary ${\bf S}:=\partial{\bf V}$ in target space using an appropriate
transformation of $\bf F$.

##### Lemma 2.

Let $\operatorname{d}\\!{\bf A}$ and $\operatorname{d}\\!{\bf S}$ be the area
vectors in base space and target space, respectively. Then we have $\int_{\bf
A}{\bf F}\cdot\operatorname{d}\\!{\bf A}=\int_{\bf S}{\bf
G}\cdot\operatorname{d}\\!{\bf S}$, with ${\bf G}:=|J|J^{-1}{\bf F}$ and
$J:=\nabla g$.

_Proof._ Consider a simplex with vertices $p_{i}$ ($i=1,2,\ldots,d$) in target
space and its corresponding points $g(p_{i})$ defining a simplex $\Delta{\bf
A}\in\\{{\bf A}_{i}\\}_{i}$ in base space. $g$ can be approximated with a
Taylor expansion

$g({\bf x}+{\bf\epsilon})=g({\bf
x})+J{\bf\epsilon}+O(||{\bf\epsilon}||_{2}^{2}),$ (7)

where $J$ is the Jacobian at the centroid of the simplex. With ${\bf F}={\bf
F}(\Delta{\bf A}^{\\!(c)})$ denoting the field at the centroid, we can write
out ${\bf F}\cdot\Delta{\bf A}^{\\!(n)}$ as (see Appendix A for a discussion
on computing surface normals)

$\begin{split}{\bf F}\cdot\Delta{\bf
A}^{\\!(n)}&=\frac{1}{(d-1)!}\begin{vmatrix}{\bf F}\\\
J\left(p_{2}-p_{1}\right)+2O(||{\bf\epsilon}||_{2}^{2})\\\
J\left(p_{3}-p_{1}\right)+2O(||{\bf\epsilon}||_{2}^{2})\\\ \vdots\\\
J\left(p_{d}-p_{1}\right)+2O(||{\bf\epsilon}||_{2}^{2})\\\ \end{vmatrix}.\\\
\end{split}$ (8)

In the limit $\Delta{\bf A}\rightarrow 0$, we have $\Delta{\bf
A}^{\\!(n)}\rightarrow\operatorname{d}\\!{\bf A}$ and
$||{\bf\epsilon}||_{2}^{2}\rightarrow 0$:

$\begin{split}{\bf F}\cdot\operatorname{d}\\!{\bf
A}&=\frac{1}{(d-1)!}\left\lvert J\begin{bmatrix}J^{-1}{\bf F}\\\
p_{2}-p_{1}\\\ p_{3}-p_{1}\\\ \vdots\\\ p_{d}-p_{1}\\\
\end{bmatrix}\right\rvert={\bf G}\cdot\operatorname{d}\\!{\bf S},\end{split}$
(9)

where ${\bf G}=|J|J^{-1}{\bf F}$ ∎

For notational convenience we assume $|J|>0$ so that outward-facing surface
normals in base space remain outward-facing in target space; the general case
requires straightforward tracking of their signs. Interestingly, $\bf G$ has a
tractable divergence, paving the way for an application of the divergence
theorem in target space.

##### Lemma 3.

With ${\bf G}$ defined as in Lemma 2, we have $\mathrm{div}({\bf G})=|J|$.

_Proof._

$\begin{split}\nabla\cdot{\bf G}&=\nabla\cdot(\text{adj}(J){\bf F})\\\
&=(\nabla\cdot\text{adj}(J)){\bf F}+\text{tr}(\text{adj}(J)\nabla_{x}\\{{\bf
F}\\})\\\ &=d^{-1}\text{tr}(\text{adj}(J)J)+(\nabla\cdot\text{adj}(J)){\bf
F}\\\ &=|J|+(\nabla\cdot\text{adj}(J)){\bf F}\\\ &=|J|,\end{split}$ (10)

because $\text{tr}(\text{adj}(J)J)=d|J|$. Furthermore, (Evans, 2010) shows
$\nabla\cdot\text{adj}(J)=0$ for any $C^{2}$-map ∎

Lemma 3 allows us to close the loop and tie the field $\bf G$ back to our
original objective of computing $P({\bf x}\in{\bf V})$. Indeed, an application
of Eq.(4) and Eq.(1) shows

$\int_{\bf S}{\bf G}\cdot\operatorname{d}\\!{\bf S}=\int_{\bf
V}|J|\operatorname{d}\\!{\bf V}=\int_{\bf V}p({\bf x})\operatorname{d}\\!{\bf
V}=P({\bf x}\in{\bf V}).$ (11)

See Fig. 3 for an overview of all equivalencies derived.

####  4.2.1 Stochastic Boundary Flux Estimator

For a simplicial polytope $\bf V$ with boundary ${\bf S}$ given by non-
overlapping simplices $\\{{\bf S}_{i}\\}_{i}$, we can leverage $\bf G$ to
define a _stochastic boundary flux (BF-S)_ estimator of the cumulative density
$P({\bf x}\in{\bf V})$ using points sampled from the boundary:

$\begin{split}\int_{\bf S}{\bf G}\cdot\operatorname{d}\\!{\bf
S}&=\sum_{i}\int_{{\bf S}_{i}}{\bf G}({\bf x})\cdot{\bf
S}_{i}^{(n)}\operatorname{d}\\!{\bf x}\\\ &=\sum_{i}{\bf
S}_{i}^{(A)}\int_{{\bf S}_{i}}\frac{{\bf G}({\bf x})}{{\bf
S}_{i}^{(A)}}\cdot{\bf S}_{i}^{(n)}\operatorname{d}\\!{\bf x}\\\
&=\sum_{i}{\bf S}_{i}^{(A)}\int_{{\bf S}_{i}}\mathcal{U}_{{\bf S}_{i}}({\bf
x})\cdot{\bf G}({\bf x})\cdot{\bf S}_{i}^{(n)}\operatorname{d}\\!{\bf x}\\\
&=\sum_{i}{\bf S}_{i}^{(A)}\cdot\mathbb{E}_{{\bf x}\sim\mathcal{U}_{{\bf
S}_{i}}}[{\bf G}({\bf x})\cdot{\bf S}_{i}^{(n)}],\end{split}$ (12)

where ${\bf S}_{i}$ is the $i$-th boundary simplex, ${\bf S}_{i}^{(A)}$ its
area, and ${\bf S}_{i}^{(n)}$ its outward-facing unit normal. Alternatively,
multiplying and dividing the surface integral $\int_{\bf S}{\bf
G}\cdot\operatorname{d}\\!{\bf S}$ with the total surface area
$\left(\sum_{i}{\bf S}_{i}^{(A)}\right)$:

$\begin{split}\int_{\bf S}{\bf G}\cdot\operatorname{d}\\!{\bf
S}&=\left(\sum_{i}{\bf S}_{i}^{(A)}\right)\cdot\mathbb{E}_{{\bf
x}\sim\mathcal{U}_{\bf S}}\left[{\bf G}({\bf x})\cdot{\bf n}({\bf
x})\right],\\\ \end{split}$ (13)

where ${\bf n}({\bf x})$ is the outward-facing unit normal at $\bf x$. Note
that $\mathcal{U}_{\bf S}$ can be decomposed in terms of a categorical
distribution $c_{i}={\bf S}_{i}^{(A)}/\sum_{j}{\bf S}_{j}^{(A)}$ over the
simplices ${\bf S}_{i}$ and the uniform distributions $\mathcal{U}_{{\bf
S}_{i}}$. More specifically, we have $\mathcal{U}_{\bf S}({\bf
x})=\sum_{i}c_{i}\cdot\mathcal{U}_{{\bf S}_{i}}({\bf x})$, wherein all but one
of the terms are zero for a point ${\bf x}$ on the surface ${\bf S}$ — because
the simplices are non-overlapping and the point ${\bf x}$ belongs to exactly
one of them.

Figure 3: Equivalencies. We show equivalent expressions for the cumulative
density of interest $P({\bf x}\in{\bf V})$ and their associated spaces. The
proposed estimator is motivated by a counter-clockwise chain of arguments,
ultimately leveraging $\int_{\bf S}{\bf G}\,\operatorname{d}\\!{\bf S}$.

###  4.3 Adaptive Estimation in Target Space

While the stochastic boundary flux estimator introduced in Section 4.2.1 is
unbiased, it does not allow for a strategic placement of evaluation points in
an iterative fashion. We can turn Eq.(12) into a sequential process by first
approximating $\mathbb{E}_{{\bf x}\sim\mathcal{U}_{{\bf S}_{i}}}[{\bf G}({\bf
x})\cdot{\bf S}_{i}^{(n)}]$ with a deterministic estimate and then performing
prioritized iterative refinements.

####  4.3.1 Deterministic Boundary Flux Estimator

Given a set of boundary simplices $\\{{\bf S}_{i}\\}_{i}$, we can approximate
the expected dot-product in Eq.(12) by the average dot-product over the
_vertices_ of ${\bf S}_{i}$. Formally, this _adaptive boundary flux (BF-A)_
estimator of $P({\bf x}\in{\bf V})$ is given by

$\int_{\bf S}{\bf G}\cdot\operatorname{d}\\!{\bf S}\approx\sum_{i}{\bf
S}_{i}^{(A)}\cdot\overline{\bf G}({\bf S}_{i})\cdot{\bf S}_{i}^{(n)},$ (14)

where $\overline{\bf G}({\bf S}_{i})$ is the average ${\bf G}$-field over all
vertices of ${\bf S}_{i}$. While Eq.(14) is a biased estimate of the
cumulative density $P({\bf x}\in{\bf V})$, its deterministic nature and
inherent structure increase sample efficiency (e.g., prudent placement,
reuseable computations) and enable strategic refinements to arbitrary
precision.

####  4.3.2 Iterative Refinement

The quality of the approximation in Eq.(14) can be controlled via iterative
refinement of the initial set of boundary simplices $\\{{\bf S}_{i}\\}_{i}$.
Each refinement step consists of a splitting operation that first introduces
one new vertex at the midpoint of a boundary line segment444The
($d-1$)-simplices $\\{{\bf S}_{i}\\}_{i}$ forming the boundary $\bf S$ each
consist of ${d\choose 2}$ line segments. and then splits all simplices sharing
that edge. The key component of this splitting process is a priority queue $Q$
managing the next boundary line segment $l$ to split, where $l$ is the longest
edge of the simplex with highest priority according to a priority function
$\mathrm{Pr}$.

$\mathrm{Pr}$ is designed such that it gives high priority to simplices which
probably have high error. As a general rule, the larger the simplex the larger
the variation in the dot-products, and hence the larger the estimation bias.
Although the area ${\bf S}_{i}^{(A)}$ of a simplex squarely fits the
definition of size, it can be a misleading measure that should not be used in
isolation — see Fig. 4 for an illustrated explanation. Therefore, we also use
the total edge length as a measure of size. While this approach works well in
lower dimensions, we saw marked improvements in terms of sample efficiency in
higher dimensions after complementing these two criteria with the standard
deviation of the dot-products at the vertices. Taken together, we define the
priority of the $i$-th simplex ${\bf S}_{i}$ as

$\Pr({\bf S}_{i})={\bf S}_{i}^{(A)}\times(\sigma_{d}({\bf
S}_{i})+\epsilon)\times\sum_{e_{j}\in{\bf S}_{i}}||e_{j}||_{2}^{2},$ (15)

where $e_{j}$ is the $j$-th edge ($1$-simplex) constituting the simplex ${\bf
S}_{i}$, $\sigma_{d}({\bf S}_{i})$ is the standard deviation of the dot-
products – computed with the unit surface normal – across the simplex ${\bf
S}_{i}$, and $\epsilon$ is a small number to guard against assigning a zero
priority to large simplices whose dot-products at the vertices are all zero.
See Appendix B for a summary of the entire splitting process in pseudo code.

Figure 4: Splitting Criterion (Total Edge Length). For two simplices with
equal area, the uncertainty about the dot-product within a simplex grows with
a point’s distance to its closest vertex, motivating the use of total edge
length as a splitting criterion.

###  4.4 Analysis

####  4.4.1 Best Case Scenarios

The integral of ${\bf G}$ over a simplex ${\bf S}_{i}$ can be computed exactly
if the flow transformation $g$ is: (a) linear over the simplex ${\bf S}_{i}$;
or (b) ${\bf G}\cdot{\bf S}_{i}^{(n)}$ is constant for all ${\bf x}\in{\bf
S}_{i}$. The latter is equivalent to ${\bf F}\cdot\operatorname{d}\\!{\bf A}$
being constant, such as for hyperspheres. In both cases, we have
$\overline{\bf G}({\bf S}_{i})=\mathbb{E}_{{\bf x}\sim\mathcal{U}_{{\bf
S}_{i}}}[{\bf G}({\bf x})]$. In practice, it is rarely possible to realize
these best cases, however, if the dot-product is locally-linear over the
simplex, our biased estimate is a good approximation to the exact integral.

####  4.4.2 Runtime Complexity & Implementation Considerations

For a given simplicial polytope, the algorithm first computes the surface
normals for each of the boundary simplices, which requires computing $d$
determinants of $(d-1)\times(d-1)$ matrices for each boundary simplex. The
total complexity of this step is thus $\mathcal{O}(kd^{4})$, where $k$ is the
number of boundary simplices. Note that the surface normals do not need to be
recomputed upon splitting edges as they remain the same.

When computing ${\bf G}$, one has to first flow the input to the uniform
distribution and then compute the Jacobian-vector product of the backward-
pass; in order to avoid computing the Jacobian-vector product, we can first
compute the vector-Jacobian product ${\bf S}_{i}^{(n)\top}J^{-1}$ of the
backward-pass followed by a dot-product with $|J|{\bf F}$ to obtain ${\bf
S}_{i}^{(n)\top}|J|J^{-1}{\bf F}$. This avoids an extra backward-pass in
PyTorch, which does not support forward-mode auto-differentiation of Jacobian-
vector-products. In fact, the derivative of ${\bf
S}_{i}^{(n)\top}|J|J^{-1}{\bf F}$ with respect to ${\bf S}_{i}^{(n)\top}$
yields exactly $\bf G$.

##  5 Experiments

(a) Glow (Kingma & Dhariwal, 2018)

(b) MAF (Papamakarios et al., 2017)

(c) FFJORD (Grathwohl et al., 2019)

Figure 5: Sample Efficiency. We highlight the sample efficiency of the
proposed adaptive boundary flux (BF-A) estimator as compared to Monte-Carlo
(MC) and Importance Sampling (IS) estimators. We show absolute estimation
errors for 3 flow architectures: two discrete normalizing flows (Glow and MAF)
— composed of 7 flow-layers with 64 hidden units in their coupling nets — and
one continuous normalizing flow (FFJORD) — parameterized by a 2-layer MLP with
64 hidden units in each layer. All models were trained on a 2D slice of the
Power dataset, with volumes $\bf V$ given by convex hulls constructed from
points sampled on the surface of spheres with radius 1.0, centered at points
sampled from the flow.

We will now evaluate the proposed adaptive boundary flux estimator (BF-A)
against traditional sampling-based estimators such as the Monte-Carlo (MC) and
Importance Sampling (IS) estimators. For the purpose of evaluation, we train
normalizing flows on $k$-dimensional ($k\in[2,5]$) data derived from $4$
tabular datasets open sourced on the UCI Machine Learning Repository (Dua &
Graff, 2017) and preprocessed as in (Papamakarios et al., 2017): Power, Gas,
Hepmass, and Miniboone. We construct normalizing flows from the following
three architecture families: Glow ((Kingma & Dhariwal, 2018)), Masked
Autoregressive Flows ((Papamakarios et al., 2017)), and FFJORD ((Grathwohl et
al., 2019)). For every pair (dataset, $k$), we obtain $2$ random
$k$-dimensional slices of the dataset over which we train the normalizing
flows. In addition, we also consider how the capacities of these models affect
the performance: for the two discrete flows (Glow and MAFs), we vary the
number of hidden dimensions in the coupling networks and the number of flow
transforms stacked together; for FFJORD, we instead vary the number of hidden
dimensions in each layer of the MLP parameterizing the ODE. For constructing
discrete flows, we choose $3$, $5$ or $7$ flow layers and construct coupling
layers with $16$, $32$ or $64$ hidden units. While one Glow layer corresponds
to a sequence of (ActNorm) – (Glow Coupling) – (Invertible $1\times 1$)
transformations, one MAF layer corresponds to a sequence of (ActNorm) – (MAF
Coupling) transformations. For continuous flows, we parameterize the neural
ODE with $2$ hidden layers, each consisting of $16$, $32$ or $64$ hidden
units. In summary, we obtain a total of $168$ flow models for every
$k\in[2,5]$ by training $21$ flow models ($9$ each from the $2$ discrete flow
families and $3$ from the continuous flows) on $2$ slices derived from each of
the $4$ datasets. We include training details and the box-plots of training
log-likelihoods in Appendix C.

##### Generating Simplicial Polytopes.

The crucial component in ensuring a fair and unbiased evaluation of the
estimators is in uniformly sampling simplicial polytopes $\bf V$ in target
space such that they are not biased towards any specific regions or cumulative
densities. In order to actually sample a simplicial polytope, we construct a
convex hull around a fixed number of points in target space. To ensure that
convex hulls of varying sizes are equally explored, we sample these points
from the boundary of a sphere centered at a flow sample and control the size
of the convex hull by controlling the radius of the sphere. Convex hulls with
a cumulative density of less than $0.01$ are rejected. While the number of
hull points affects the average size and cumulative density of the enclosed
volume, it does not markedly change the average size of the simplices as
defined in Section 4.3.2 and we therefore do not consider the number of hull
points as an evaluation axis.

In our evaluation we consider $20$ points to construct the hull, trading off
lower cumulative densities against linear approximation of the spherical
boundary; also note that all $20$ points lie on the boundary of the hull and
count against the sample budget of the BF-A algorithm. Furthermore, we
consider spheres with radius $0.5$, $0.75$, or $1.0$, amounting to a coverage
radius of up to $1$ standard deviation for our normalized data. Finally, in
order to evaluate a normalizing flow model, we construct $5$ convex hulls for
each choice of the radius. Overall, we evaluate each CDF estimator over a
total of $10,\\!080$ hulls.555$10,\\!080\ \mathrm{hulls}=168\
\frac{\mathrm{models}}{\mathrm{dim}}\cdot 4\ \mathrm{\mathrm{dim}}\cdot 5\
\frac{\mathrm{hulls}}{\mathrm{models}\cdot\mathrm{radii}}\cdot 3\
\mathrm{radii}$

##### Evaluation.

For each convex hulls we compute a reference cumulative density using $2M$
importance samples and evaluate all estimators against this gold standard
based on their predicted cumulative density using a fixed sample budget $K$.
Unlike our proposed adaptive boundary flux estimator (BF-A), which is
deterministic, the Monte-Carlo (MC) and Importance Sampling (IS) estimators
are stochastic; thus, in order to accurately capture the variance of these
stochastic estimators, we collect results across $5$ runs for each hull. We
then evaluate the estimators in terms of both relative and absolute errors.
Table 1 shows the aggregate results across all $10,\\!080$ configurations with
a sample budget of $4000$ points. In the remainder of this section, we will
unroll these results and analyze our estimator along individual evaluation
axes.

Estimator | Absolute Error | Relative Error
---|---|---
IS | 0.00572±0.01585 | 0.03000±0.04632
MC | 0.00381±0.00381 | 0.04864±0.05178
BF-A (ours) | 0.00152±0.00554 | 0.01822±0.06626

Table 1: Quantitative Evaluation (Aggregate): We show the absolute and
relative errors of each estimator with a sample budget of $4000$ points,
averaged across all evaluation axes. The proposed BF-A algorithm outperforms
the IS and MC estimators in terms of both the absolute and relative error.

##### Discussion.

On average, our BF-A estimator obtains approximately $3.8\times$ and
$2.5\times$ lower absolute errors than IS and MC, respectively, for a sample
budget of $4k$ points (Table 1); likewise, we obtain $1.7\times$ and
$2.7\times$ lower relative errors than IS and MC, respectively. For the
specific case of 2D flows, we report relative errors that are better than IS
by $11.6\times$ and MC by $22.1\times$ (Table 2); also see Fig. 5 for a
qualitative comparison of the estimators for 2D flows in terms of the absolute
deviation from the true cumulative density as obtained over a random hull: we
note that in all $3$ cases the proposed algorithm converges to the true
estimate well within $500$ samples while the stochastic estimates of the MC
and IS estimators exhibit high variance. This pattern is consistent across our
experiments with 2D flows and we find that our algorithm outperforms IS and MC
by $30\times$ and $55\times$, respectively, for a sample budget of just $500$
points. A comparison across different hull radii is shown in Table 3. We note
that, while the absolute errors of all estimators increase with increasing
hull-sizes, the relative errors of MC decrease with increasing hull sizes,
while those of IS and BF-A increase with increasing hull sizes; nonetheless,
BF-A outperforms IS and MC in terms of both absolute and relative errors for
the explored hull sizes. Additional quantitative evaluations can be found in
the Appendix: in terms of different capacities (Appendix D; Tables 4-5), the
observed trends are in line with the analysis above, both for discrete and
continuous flows. We also note that varying the width of the coupling
transform (in discrete flows) does not affect the observed trends as much as
the depth of the flow transform. Finally, we observe similar trends across
varying flow architectures as well (Appendix D; Table 6); also see Fig. 8 in
Appendix D for a visual summary of all quantitative evaluations.

Our qualitative and quantitative evaluations show a marked improvement of the
proposed method over traditional sampling-based estimators. BF-A benefits from
our novel formulation of cumulative densities as boundary fluxes as well as
from the sample efficiency enabled by the priority function $\mathrm{Pr}$.
Intuitively, the estimation error of the BF-A estimator can be traced back to
how much the dot-product ${\bf G}\cdot{\bf n}$ can change for small changes in
the input $\bf x$; in fact, our motivation behind evaluating the estimators as
a function of growing dimensions, growing hull-sizes, and growing model
capacities is targeted towards testing how well our estimator can adapt to
increasing complexity. We observe that while our method continues to
outperform the sampling-based estimators with growing hull sizes (Table 3),
model capacities (Appendix D; Table 4 and Table 5), varying flow architectures
(Appendix D; Table 6), and data dimensions up to and including 4D, BF-A only
performs comparably to the sampling-based estimators for a sample budget of
$4k$ points with 5D flows (Table 2).

Dims | Estimator | Absolute Error | Relative Error
---|---|---|---
2 | IS | 0.00737±0.01981 | 0.01620±0.02459
MC | 0.00482±0.00434 | 0.03101±0.03366
BF-A (ours) | 0.00071±0.00281 | 0.00146±0.00308
3 | IS | 0.00528±0.01564 | 0.02543±0.03618
MC | 0.00363±0.00381 | 0.05247±0.05202
BF-A (ours) | 0.00134±0.00584 | 0.00710±0.03254
4 | IS | 0.00405±0.00818 | 0.04495±0.05404
MC | 0.00274±0.00251 | 0.06576±0.06162
BF-A (ours) | 0.00237±0.00780 | 0.03209±0.08405
5 | IS | 0.00376±0.00640 | 0.07233±0.08153
MC | 0.00244±0.00218 | 0.07340±0.06483
BF-A (ours) | 0.00380±0.00658 | 0.07796±0.14184

Table 2: Quantitative Evaluation (Dimension-wise): We show the absolute and
relative errors of each of the methods with a sample budget of 4000 points
averaged across all the configurations separately for each dimensionality. We
note that BF-A outperforms both IS and MC on average for 2D, 3D and 4D flows
in terms of both absolute and relative errors.

Radius | Estimator | Absolute Error | Relative Error
---|---|---|---
0.5 | IS | 0.00338±0.01066 | 0.01895±0.02430
MC | 0.00336±0.00336 | 0.05670±0.05612
BF-A (ours) | 0.00073±0.00287 | 0.01016±0.04500
0.75 | IS | 0.00480±0.01327 | 0.02771±0.04061
MC | 0.00365±0.00350 | 0.04896±0.05205
BF-A (ours) | 0.00126±0.00460 | 0.01666±0.05429
1.0 | IS | 0.00808±0.01998 | 0.03941±0.05871
MC | 0.00425±0.00428 | 0.04292±0.04761
BF-A (ours) | 0.00229±0.00728 | 0.02500±0.08450

Table 3: Quantitative Evaluation (Size-wise): We show the absolute and
relative errors of each of the methods with a sample budget of 4000 points
averaged across all the configurations separately for each radius of the hull.
On average, BF-A performs better than IS and MC in terms of both absolute and
relative error.

##### Limitations & Future Work.

Future extensions of the proposed estimator could consider using higher-order
information of the model, such as curvatures, in defining the priority
function in order to improve the scalability of the method. Incorporating
curvature information may also enable the definition of a more natural
stopping criterion than sample budget. Another interesting extension to the
adaptive estimator would be to split a batch of edges together as opposed to
splitting a single edge at a time. Additionally, it remains to be explored how
different choices of field $\bf F$ can be useful in minimizing the estimation
error for a given simplicial polytope. Finally, noting that our method works
better with FFJORD than MAF and Glow on average (Appendix D; Table 6), even
though they all obtain similar train log-likelihoods (Appendix C), one could
also explore regularization methods that encourage smoother transformations
between the target and base spaces.

## Conclusion

Cumulative density functions are as fundamentally important as probability
density functions, but cannot always be estimated efficiently and exactly. In
this work, we propose a new stochastic CDF estimator derived as an adaptation
of the divergence theorem by exploiting the diffeomorphic properties of flow
that allow us to estimate CDF as the volume of the boundary. Furthermore, we
demonstrate how this stochastic CDF estimator can be made more sample-
efficient by reusing computations and strategically acquiring new points such
that the estimate obtains a good CDF estimate with as few points as possible.
In our evaluation of this new adaptive estimator, we find that it outperforms
traditional sampling-based estimators in terms of sample-efficiencies and
accuracies. We believe that our proposed efficient CDF estimator will help
open up new application areas with Normalizing Flows.

##### Acknowledgements.

This work was supported by a Borealis AI Graduate Fellowship. The authors also
thank Vector Institute and Sageev Oore for supporting the project with
computational resources.

## References

  * Abdar et al. (2021) Abdar, M., Pourpanah, F., Hussain, S., Rezazadegan, D., Liu, L., Ghavamzadeh, M., Fieguth, P. W., Cao, X., Khosravi, A., Acharya, U. R., Makarenkov, V., and Nahavandi, S.  A review of uncertainty quantification in deep learning: Techniques, applications and challenges.  _Inf. Fusion_ , 2021.
  * Armstrong (2013) Armstrong, M. A.  _Basic topology_.  Springer Science & Business Media, 2013.
  * Botev (2017) Botev, Z. I.  The normal law under linear restrictions: simulation and estimation via minimax tilting.  _Journal of the Royal Statistical Society_ , 2017.
  * Botev & l’Ecuyer (2015) Botev, Z. I. and l’Ecuyer, P.  Efficient probability estimation and simulation of the truncated multivariate student-t distribution.  _Winter Simulation Conference (WSC)_ , 2015.
  * Chilinski & Silva (2020) Chilinski, P. M. and Silva, R.  Neural likelihoods via cumulative distribution functions.  _Uncertainty in Artificial Intelligence, UAI_ , 2020.
  * Cundy & Ermon (2020) Cundy, C. and Ermon, S.  Flexible approximate inference via stratified normalizing flows.  _Uncertainty in Artificial Intelligence, UAI_ , 2020.
  * Cunningham et al. (2013) Cunningham, J. P., Hennig, P., and Lacoste-Julien, S.  Gaussian probabilities and expectation propagation.  _Journal of Machine Learning Research_ , 2013.
  * Dai & Seljak (2021) Dai, B. and Seljak, U.  Sliced iterative normalizing flows.  _ICML Workshop on Invertible Neural Networks, Normalizing Flows, and Explicit Likelihood Models_ , 2021.
  * Dua & Graff (2017) Dua, D. and Graff, C.  UCI machine learning repository, 2017.  URL http://archive.ics.uci.edu/ml.
  * Evans (2010) Evans, L. C.  Partial differential equations.  _Graduate studies in mathematics_ , 2010.
  * Genz & Bretz (2009) Genz, A. and Bretz, F.  _Computation of multivariate normal and t probabilities_ , volume 195\.  Springer Science & Business Media, 2009.
  * Grathwohl et al. (2019) Grathwohl, W., Chen, R. T. Q., Bettencourt, J., Sutskever, I., and Duvenaud, D.  FFJORD: free-form continuous dynamics for scalable reversible generative models.  _7th International Conference on Learning Representations, ICLR_ , 2019.
  * Kingma & Dhariwal (2018) Kingma, D. P. and Dhariwal, P.  Glow: Generative flow with invertible 1x1 convolutions.  _Advances in Neural Information Processing Systems, NeurIPS_ , 2018\.
  * Kobyzev et al. (2020) Kobyzev, I., Prince, S., and Brubaker, M.  Normalizing flows: An introduction and review of current methods.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 2020.
  * Lange (1999) Lange, K.  Quadrature methods.  _Numerical analysis for statisticians_ , 1999.
  * Liang et al. (2020) Liang, E., Yang, Z., Stoica, I., Abbeel, P., Duan, Y., and Chen, X.  Variable skipping for autoregressive range density estimation.  _International Conference on Machine Learning, ICML_ , 2020.
  * Liu et al. (2019) Liu, J. Z., Paisley, J. W., Kioumourtzoglou, M., and Coull, B. A.  Accurate uncertainty estimation and decomposition in ensemble learning.  _Advances in Neural Information Processing Systems, NeurIPS_ , 2019\.
  * Mazaheri et al. (2020) Mazaheri, B., Jain, S., and Bruck, J.  Robust correction of sampling bias using cumulative distribution functions.  _Advances in Neural Information Processing Systems, NeurIPS_ , 2020\.
  * Mehta et al. (2012) Mehta, A., Neukirchen, M., Pfetsch, S., and Poppensieker, T.  Managing market risk: Today and tomorrow, 2012.
  * Minka (2001) Minka, T. P.  Expectation propagation for approximate bayesian inference.  _Uncertainty in Artificial Intelligence, UAI_ , 2001.
  * Papamakarios et al. (2017) Papamakarios, G., Murray, I., and Pavlakou, T.  Masked autoregressive flow for density estimation.  _Advances in Neural Information Processing Systems, NeurIPS_ , 2017\.
  * Papamakarios et al. (2021) Papamakarios, G., Nalisnick, E., Rezende, D. J., Mohamed, S., and Lakshminarayanan, B.  Normalizing flows for probabilistic modeling and inference.  _Journal of Machine Learning Research_ , 2021.
  * Press et al. (1988) Press, W. H., Teukolsky, S. A., Vetterling, W. T., and Flannery, B. P.  _Numerical recipes in C_.  Cambridge University Press, 1988.
  * Richardson et al. (1997) Richardson, M., Boudoukh, J., and Whitelaw, R. F.  The best of both worlds: A hybrid approach to calculating value at risk, 1997.
  * Ridgeway & Mozer (2018) Ridgeway, K. and Mozer, M. C.  Learning deep disentangled embeddings with the f-statistic loss.  _Advances in Neural Information Processing Systems, NeurIPS_ , 2018\.
  * Sun et al. (2021) Sun, W., Lee, C., and Lee, C.  DFAC framework: Factorizing the value function via quantile mixture for multi-agent distributional q-learning.  _International Conference on Machine Learning, ICML_ , 2021.
  * Wade (2017) Wade, W.  _An Introduction to Analysis_.  Pearson, 2017.

##  Appendix A Computing Surface Normals

In applying the divergence theorem, one would need to compute the surface
normal vector. In the context of this paper, we compute surface normals to the
$(d-1)$-dimensional simplices constituting the boundary: for example, we would
need to compute the surface normal to line-segments in 2D, triangles in 3D,
tetrahedrons in 4D and so on. By definition of $(d-1)$-dimensional simplices,
they lie in a plane defined by $(d-1)$ linearly independent vectors; since we
are computing normals to $(d-1)$-dimensional simplices lying in
$\mathbb{R}^{d}$, the vectors spanning the plane containing the simplex and
the normal vector to this plane form a full basis of $\mathbb{R}^{d}$. In
general, one could imagine solving a system of simultaneous equations in order
to determine a vector normal to a given set of vectors. However, in this
specific case, since we know that we are computing a normal vector to $(d-1)$
linearly independent vectors, we can just define the normal vector to a
simplex ${\bf z}$ defined by the vertex set $\\{z_{i}\\}_{i}$ in terms of the
determinant:

${\bf v}=\begin{vmatrix}{\bf e}\\\ z_{2}-z_{1}\\\ z_{3}-z_{1}\\\ \vdots\\\
z_{d}-z_{1}\end{vmatrix}$ (16)

where, $\bf e$ is the orthonormal unit basis in $\mathbb{R}^{d}$. Since the
determinant of a matrix containing two identical rows evaluates to 0, we can
easily confirm that the dot product $\forall i\in[2,d]$, ${\bf
v}\cdot(z_{i}-z_{1})=0$. Also note that at a high level, this construction
resembles the vector cross-product operation. In fact, similar to vector cross
products, the vector norm of the $\bf v$ is equal to the volume of the
parallelotope spanned by the set of $(d-1)$ linearly independent vectors
$\\{z_{i}-z_{1}\\}_{i}$. Therefore, in order to derive the surface normal
vector whose magnitude is equal to the volume of the simplex $\bf z$, we would
need to divide $\bf v$ with $(d-1)!$:

${\bf z}^{(n)}=\frac{1}{(d-1)!}\begin{vmatrix}{\bf e}\\\ z_{2}-z_{1}\\\
z_{3}-z_{1}\\\ \vdots\\\ z_{d}-z_{1}\end{vmatrix}$ (17)

Accordingly, for some vector $\bf w$, the dot-product ${\bf w}\cdot{\bf
z}^{(n)}$ can be written as:

${\bf w}\cdot{\bf z}^{(n)}=\frac{1}{(d-1)!}\begin{vmatrix}{\bf w}\\\
z_{2}-z_{1}\\\ z_{3}-z_{1}\\\ \vdots\\\ z_{d}-z_{1}\end{vmatrix}$ (18)

##  Appendix B Pseudo-codes

Algorithm 1 Split edges and update the volume

Objects and Attributes:

  * edge: A $1$ simplex object in $\mathbb{R}^{d}$ having the following attributes:

    * –

endPoints: the endPoints of the edge

    * –

simplices: set of $d-1$-dimensional simplices that contains this edge

  * simplex: A $d-1$ simplex object in $\mathbb{R}^{d}$ having the following attributes:

    * –

edges: the set of edges contained in the simplex

    * –

unitSurfaceNormal: the unit surface normal of the simplex, adjusted so that it
points outwards (or inwards) of the hull.

    * –

volume: the volume of the simplex

    * –

getVolumeElement: a method that will loop through the points and computes the
mean dot-product.

Input:

  * edgeToSplit: The selected edge to split according to the priority function described in Eq.(15).

  * edges: Hashmap of Edges

    * –

Allows $\mathcal{O}(1)$ retrieval of edge objects from the endpoints. If
object does not exist in the queue, a new one is created and returned.

  * volume: The current estimated volume

Output:

  * newVolume

1:newVolume $\leftarrow$ volume $\triangleright$ Initialize the newVolume to
be equal to the current volume

2:for simplex in edgeToSplit.simplices do

3: simplex1 $\leftarrow$ new Simplex()

4: simplex2 $\leftarrow$ new Simplex()

5: simplex1.unitSurfaceNormal $\leftarrow$ simplex.unitSurfaceNormal

6: simplex2.unitSurfaceNormal $\leftarrow$ simplex.unitSurfaceNormal

7: simplex1.volume $\leftarrow$ simplex.volume/2

8: simplex2.volume $\leftarrow$ simplex.volume/2

9: for otherEdge in simplex.edges-{edgeToSplit} do

10: newEdge1 $\leftarrow$ edges.get(otherEdge.endPoint[0],edge.midPoint)

11: newEdge2 $\leftarrow$ edges.get(otherEdge.endPoint[1],edge.midPoint)

12: if otherEdge.hasEndPoint(edgeToSplit.endPoint[0]) then

13: otherEdge.add(simplex1)

14: simplex1.edges.add(otherEdge, newEdge1, newEdge2)

15: else if otherEdge.hasEndPoint(edgeToSplit.endPoint[1]) then

16: otherEdge.add(simplex2)

17: simplex2.edges.add(otherEdge, newEdge1, newEdge2)

18: else

19: otherEdge.add(simplex1, simplex2)

20: newEdge1.add(simplex1, simplex2)

21: newEdge2.add(simplex1, simplex2)

22: simplex1.edges.add(otherEdge, newEdge1, newEdge2)

23: simplex2.edges.add(otherEdge, newEdge1, newEdge2)

24: otherEdge.remove(simplex)

25: newVolume $\leftarrow$ newVolume - simplex.getVolumeElement()

26: newVolume $\leftarrow$ newVolume + simplex1.getVolumeElement()

27: newVolume $\leftarrow$ newVolume + simplex2.getVolumeElement()

28:return newVolume

Algorithm 2 getVolumeElement of Simplex Object

Objects and Attributes:

  * edge: A $1$ simplex object in $\mathbb{R}^{d}$ having the following attributes:

    * –

endPoints: the endPoints of the edge

    * –

simplices: set of $d-1$ simplices that contains this edge

  * simplex: A $d-1$ simplex object in $\mathbb{R}^{d}$ having the following attributes:

    * –

edges: the set of edges contained in the simplex

    * –

unitSurfaceNormal: the unit surface normal of the simplex, adjusted so that it
points outwards (or inwards) of the hull.

    * –

volume: the volume of the simplex

    * –

getVolumeElement: a method that will loop through the points and computes the
mean dot-product.

Input:

  * simplex: A simplex object

  * flow: The normalizing flow

  * D: Number of dimensions

Output:

  * volume

1:volumeElement $\leftarrow$ $0$

2:for point in simplex.points() do

3: volumeElement $\leftarrow$ volumeElement +
simplex.volume*flow.dotProduct(point,simplex.unitSurfaceNormal)

4:return volumeElement

##  Appendix C Training Details and Log-likelihoods

For all normalizing flows, we train the models with a batch size of 10k and
stop when the log-likelihoods do not improve over 5 epochs. For the continuous
flows, we used the exact divergence for computing the log-determinant. We used
exp-scaling in the affine coupling layer of both MAF and Glow models – and, in
order to prevent numerical overflows, we applied a tanh nonlinearity before
the exp-scaling. Finally, we used softplus as our activation function for both
the Neural ODE and coupling networks. From Fig. 6 and Fig. 7, we observe both
the Continuous and Discrete flows obtain similar log-likelihoods and are able
to fit the training data well.

(a) 2D

(b) 3D

(c) 4D

(d) 5D

Figure 6: Log-Likelihoods (Discrete Flows).

(a) 2D

(b) 3D

(c) 4D

(d) 5D

Figure 7: Log-Likelihoods (Continuous Flows).

##  Appendix D Extended Results

In this section, we include results that evaluate our algorithm as a function
of flow-architectures and their capacities. We also include plots to visually
summarize our quantitative results.

Depth | Estimator | Absolute Error | Relative Error
---|---|---|---
3 | IS | 0.00449±0.01048 | 0.02683±0.04063
MC | 0.00359±0.00334 | 0.04844±0.05345
BF-A (ours) | 0.00115±0.00479 | 0.01719±0.07362
5 | IS | 0.00562±0.01381 | 0.03167±0.04787
MC | 0.00372±0.00336 | 0.04833±0.05216
BF-A (ours) | 0.00142±0.00536 | 0.02221±0.07773
7 | IS | 0.00585±0.01499 | 0.03113±0.04583
MC | 0.00375±0.00342 | 0.04776±0.05186
BF-A (ours) | 0.00192±0.00775 | 0.02470±0.08521

Width | Estimator | Absolute Error | Relative Error
---|---|---|---
16 | IS | 0.00552±0.01409 | 0.03095±0.04505
MC | 0.00368±0.00335 | 0.04851±0.05277
BF-A (ours) |  0.00195±0.00767 |  0.02290±0.07088
32 | IS | 0.00544±0.01357 | 0.02940±0.04310
MC | 0.00365±0.00334 | 0.04782±0.05253
BF-A (ours) |  0.00130±0.00609 |  0.01944±0.07800
64 | IS | 0.00501±0.01200 | 0.02932±0.04666
MC | 0.00374±0.00342 | 0.04820±0.05217
BF-A (ours) |  0.00125±0.00396 |  0.02184±0.08754

Table 4: Quantitative Evaluation (Capacity - Glow and MAF): We show the
absolute and relative errors of each of the methods with a sample budget of
4000 points as a function of the depth and width of discrete flows as
explained; we see that BF-A outperforms both IS and MC.

Width | Estimator | Absolute Error | Relative Error
---|---|---|---
16 | IS | 0.00631±0.02026 | 0.02975±0.04846
MC | 0.00383±0.00472 | 0.04961±0.04996
BF-A (ours) | 0.00171±0.00504 | 0.01430±0.03244
32 | IS | 0.00614±0.01808 | 0.03005±0.04966
MC | 0.00400±0.00420 | 0.04991±0.05094
BF-A (ours) | 0.00131±0.00343 | 0.01034±0.02098
64 | IS | 0.00701±0.02152 | 0.03080±0.04856
MC | 0.00432±0.00470 | 0.04914±0.05012
BF-A (ours) | 0.00170±0.00391 | 0.01132±0.02083

Table 5: Quantitative Evaluation (Capacity - FFJORD): We show the absolute and
relative errors of each of the methods with a sample budget of 4000 points as
a function of the number of hidden units of two-layer MLPs parameterizing
FFJORD as explained; we see that BF-A outperforms both IS and MC.

Flow | Method | Absolute Error | Relative Error
---|---|---|---
Glow | IS | 0.00395±0.01131 | 0.02524±0.03818
MC | 0.00357±0.00328 | 0.05072±0.05421
BF-A (ours) | 0.00155±0.00726 | 0.01963±0.07250
MAF | IS | 0.00667±0.01480 | 0.03445±0.05032
MC | 0.00379±0.00346 | 0.04568±0.05062
BF-A (ours) | 0.00145±0.00473 | 0.02311±0.08501
FFJORD | IS | 0.00649±0.02000 | 0.03021±0.04890
MC | 0.00405±0.00455 | 0.04955±0.05034
BF-A (ours) | 0.00157±0.00417 | 0.01196±0.02529

Table 6: Quantitative Evaluation (Flow-Architecture): We show the absolute and
relative errors of each of the methods with a sample budget of 4000 points as
a function of the flow-architecture; we see that BF-A outperforms both IS and
MC. Most interestingly, we find that our method is more efficient with Neural
ODEs than with Glow or MAF – our analyses show that this can be attributed to
the fact that Neural ODEs have much smoother diffeomorphic transformations
than Glow and MAF. Perhaps more importantly, the Neural ODEs also obtain
similar or better training log-likelihoods than Glow and MAF (Appendix C).

(a) Aggregate Evaluation of Estimators

(b) Evaluation of Estimators with Growing Dimensions

(c) Evaluation of Estimators with Growing Hull Sizes.

(d) Evaluation of Estimators as a function of Flow Architectures.

Figure 8: Visualizations of Quantitative Results.