# Reliable Estimation of KL Divergence using a Discriminator in Reproducing
Kernel Hilbert Space

Sandesh Ghimire, Aria Masoomi, Jennifer Dy
Department of Electrical and Computer Engineering
Northeastern University
Boston, MA, USA
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>

###### Abstract

Estimating Kullback–Leibler (KL) divergence from samples of two distributions
is essential in many machine learning problems. Variational methods using
neural network discriminator have been proposed to achieve this task in a
scalable manner. However, we noted that most of these methods using neural
network discriminators suffer from high fluctuations (variance) in estimates
and instability in training. In this paper, we look at this issue from
statistical learning theory and function space complexity perspective to
understand why this happens and how to solve it. We argue that the cause of
these pathologies is lack of control over the complexity of the neural network
discriminator function and could be mitigated by controlling it. To achieve
this objective, we 1) present a novel construction of the discriminator in the
Reproducing Kernel Hilbert Space (RKHS), 2) theoretically relate the error
probability bound of the KL estimates to the complexity of the discriminator
in the RKHS space, 3) present a scalable way to control the complexity (RKHS
norm) of the discriminator for a reliable estimation of KL divergence, and 4)
prove the consistency of the proposed estimator. In three different
applications of KL divergence – estimation of KL, estimation of mutual
information and Variational Bayes – we show that by controlling the complexity
as developed in the theory, we are able to reduce the variance of KL estimates
and stabilize the training.

##  1 Introduction

Estimating Kullback–Leibler (KL) divergence from data samples is an essential
component in many machine learning problems including Bayesian inference,
calculation of mutual information or methods using information theoretic
objectives. Variational formulation of Bayesian Inference requires KL
divergence computation, which could be challenging when we only have finite
samples from two distributions. Similarly, computation of information
theoretic objectives like mutual information requires computation of KL
divergence between the joint and the product of marginals.

KL divergence estimation from samples was studied thoroughly by Nguyen et al.
[1] using a variational technique, convex optimization and RKHS norm
regularization, while also providing theoretical guarantees and insights.
However, their technique requires handling the whole dataset at once and is
not scalable. Many modern models need to use KL divergence with large scale
data, and often with neural networks, for example total correlation
variational autoencoder (TC-VAE) [2], adversarial variational Bayes (AVB) [3],
information maximizing GAN (InfoGAN) [4], and amortized MAP [5] all need to
compute KL divergence in a deep learning setup. These large scale models have
imposed new requirements on KL divergence estimation like scalability (able to
handle large amount of data samples) and minibatch compatibility (compatible
with minibatch-based optimization).

Methods like Nguyen et al. [1] are not suitable in the large scale setup.
These modern needs were later met by modern neural network based methods such
as variational divergence minimization (VDM) [6], mutual information neural
estimation (MINE) [7], and discriminator based KL estimation with GAN-type
objective [8, 5]. A key attribute of these methods is that they are based on
updating a neural-net based discriminator to estimate KL divergence from a
subset of samples making them scalable and minibatch compatible. We, however,
noticed that even in simple examples, these methods exhibited pathologies like
unreliability (high fluctuation of estimates) or instability during training
(KL estimates blowing up). Similar observations of instability of VDM and MINE
have also been reported in the literature [8, 9].

Why are these techniques unreliable? In this paper, we attempt to understand
the core problem in the KL estimation using discriminator network. We look at
it from the perspective of statistical learning theory and discriminator
function space complexity and draw insights. Based on these insights, we
propose that these fluctuations are a consequence of not controlling the
smoothness and the complexity of the discriminator function space. Measuring
and controlling the complexity of function space itself becomes a difficult
problem when the discriminator is a deep neural network. Note that naive
approaches to bound complexity by the number of parameters would neither be
guaranteed to yield meaningful bound [10], nor be easy to implement.

Therefore, we present the following contributions to resolve these challenges.
First, we propose a novel construction of the discriminator function using
deep network such that it lies in a smooth function space, the Reproducing
Kernel Hilbert Space(RKHS). By utilizing the learning theory and the
complexity analysis of the RKHS space, we bound the probability of the error
of KL-divergence estimates in terms of the radius of RKHS ball and kernel
complexity. Using this bound, we propose a scalable way to control the
complexity by penalizing the RKHS norm. This additional regularization of the
complexity is still linear, (${O(m)}$) in time complexity with the number of
data samples. Then, we prove consistency of the proposed KL estimator using
ideas from empirical process theory. Experimentally, we demonstrate that the
proposed way of controlling complexity significantly improves KL divergence
estimation and significantly reduce the variance. In mutual information
estimation, our method is competitive with the state-of-the-art method and in
Variational Bayesian application, our method stabilizes training of MNIST
dataset leading to sharp reconstruction.

##  2 Related Work

Nguyen et al. [1] used variational method to estimate KL divergence from
samples of two distribution using convex risk minimization (CRM). They used
the RKHS norm as a way to both measure and penalize the complexity of the
variational function. However, their work required handling all data at once
and solving a convex optimization problem which has time complexity in the
order of $O(m^{3})$ and space complexity in the order of $O(m^{2})$ . Ahuja
[11] used similar convex formulation in RKHS space and found it difficult to
scale. VDM reformulated the f-Divergence objective using Fenchel duality and
used a neural network to represent the variational function [6]. Although
close in concept to [1], it is scalable since it uses a separate discriminator
network and adversarial optimization. It, however, did not control the
complexity of the neural-net function, and faced issues with stability.

One area of modern application of KL-divergence estimation is in computing
mutual information, which is useful in applications such as stabilizing GANs
[7]. MINE [7] also optimized a lower bound to KL divergence (Donsker-Varadhan
representation). Similar to VDM, MINE used a neural network as the dual
variational function: it is thus scalable, but without complexity control and
is unstable. Another use of KL divergence is scalable variational inference
(VI) as shown in AVB [8]. VI requires KL divergence estimation between the
posterior and the prior, which becomes nontrivial when a sample based scalable
estimation is required. AVB solved it using GAN-type adversarial formulation
and a neural network discriminator. Similarly, [5] used GAN-type adversarial
formulation to obtain KL divergence in amortized inference.

Chen et al. [2] proposed TC-VAE to improve disentanglement by penalizing the
KL divergence between the marginal latent distribution and the product of
marginals in each dimension. The KL divergence was computed by a minibatch-
based sampling strategy that gives a biased estimate. Our work is close to
Song et al. [9] who investigated the high variance in existing mutual
information estimators and found that clipping the discriminator output is
helpful in reducing variance. In our work, we take a principled way to connect
variance to the complexity of discriminator function space and constrain it by
penalizing its RKHS norm instead. None of the existing works considered
looking at the discriminator function space, connecting its complexity to the
unreliable KL-divergence estimation, or mitigating the problem by controlling
the complexity.

##  3 Reproducing Kernel Hilbert Space

Let $\mathcal{H}$ be a Hilbert space of functions $f:\mathcal{X}\to{\rm
I\\!R}$ defined on non-empty space $\mathcal{X}$. It is a Reproducing Kernel
Hilbert Space (RKHS) if the evaluation functional,
$\delta_{x}:\mathcal{H}\to{\rm I\\!R}$, $\delta_{x}:f\mapsto f(x)$, is linear
continuous $\forall x\in\mathcal{X}$. Every RKHS, $\mathcal{H}_{K}$, is
associated with a unique positive definite kernel,
$K:\mathcal{X}\times\mathcal{X}\to{\rm I\\!R}$, called the reproducing kernel
[12], such that it satisfies:
1\. $\forall x\in\mathcal{X},K(.,x)\in\mathcal{H}_{K}$ 2\. $\forall
x\in\mathcal{X},\forall f\in\mathcal{H}_{K},\hskip 2.84544pt\langle
f,K(.,x)\rangle_{\mathcal{H}_{K}}=f(x)$

RKHS is often studied using a specific integral operator. Let
$\mathcal{L}_{2}(d\rho)$ be a space of functions $f:\mathcal{X}\to{\rm I\\!R}$
that are square integrable with respect to a Borel probability measure $d\rho$
on $\mathcal{X}$, we define an integral operator
$\mathscr{L}_{K}:\mathcal{L}_{2}(d\rho)\to\mathcal{L}_{2}(d\rho)$ [13, 14]:
$(\mathscr{L}_{K}f)(x)=\int_{\mathcal{X}}f(y)K(x,y)d\rho(y)$ This operator
will be important in constructing a function in RKHS and in computing sample
complexity.

##  4 Problem Formulation and Contribution

GAN-type Objective for KL Estimation: Let $p(x)$ and $q(x)$ be two probability
density functions in space $\mathcal{X}$ and we want to estimate their KL
divergence using finite samples from each distribution in a scalable and
minibatch compatible manner. As shown in [8, 5], this can be achieved by using
a discriminator function. First, a discriminator $f:\mathcal{X}\to{\rm I\\!R}$
is trained with the objective:

$\displaystyle
f^{*}=\underset{f}{\operatorname{argmax}}[{E_{p(x)}\log\sigma(f(x))+E_{q(x)}\log(1-\sigma(f(x)))}]$
(1)

where $\sigma$ is the Sigmoid function given by
$\sigma(x)=\frac{e^{x}}{1+e^{x}}$. Then it can be shown [8, 5] that the KL
divergence $KL(p(x)||q(x))$ is given by: $KL(p(x)||q(x))=E_{p(x)}[f^{*}(x)]$

Sources of Error: Eq. (1) is ambiguous in the sense that it is silent about
the discriminator function space over which the optimization is carried out.
Typically, a neural network is used as the discriminator. This implies that we
are considering the space of functions represented by the neural network of
given architecture as the hypothesis space, over which the maximization occurs
in eq. (1). Hence, we must rewrite eq. (1) as

$\displaystyle f^{*}_{h}=\underset{f\in
h}{\operatorname{argmax}}[{E_{p(x)}\log\sigma(f(x))+E_{q(x)}\log(1-\sigma(f(x)))}]$
(2)

where $h$ is the discriminator function space. Furthermore, we also
approximate integrals in eq. (2) with the Monte Carlo estimate using finite
number of samples, say $m$, from the distribution $p$ and $q$.

$\displaystyle f^{m}_{h}=\underset{f\in
h}{\operatorname{argmax}}\Big{[}{\frac{1}{m}\sum_{x_{i}\sim
p(x_{i})}\log\sigma(f(x_{i}))+\frac{1}{m}\sum_{x_{j}\sim
q(x_{j})}\log(1-\sigma(f(x_{j})))}\Big{]}$ (3)

Similarly, we write KL estimate obtained from, respectively, infinite and
finite samples as:

$\displaystyle KL(f)=E_{p(x)}[f(x)],\hskip
11.38092ptKL_{m}(f)=\frac{1}{m}\sum_{x_{i}\sim p(x_{i})}[f(x)]$ (4)

Each of these steps introduce some error in our estimate. We can now start our
analysis by first decomposing the total estimation error as:

$\displaystyle
KL_{m}(f^{m}_{h})-KL(f^{*})=\underbrace{KL_{m}(f^{m}_{h})-KL(f^{m}_{h})}_{\text{Deviation-
from-mean
error}}+\underbrace{KL(f^{m}_{h})-KL(f^{*}_{h})}_{\text{Discriminator induced
error}}+\underbrace{KL(f^{*}_{h})-KL(f^{*})}_{Bias}$ (5)

This equation decomposes total estimation error into three terms: 1) deviation
from the mean error, 2) error in KL estimate by the discriminator due to using
finite samples in optimization eq. (3), and 3) bias when the considered
function space does not contain the optimal function. Here, we concentrate on
quantifying the probability of deviation-from-mean error which is directly
related to observed variance of the KL estimate.

Summary of Technical Contributions: Since the deviation is the difference
between a random variable and its mean, we can bound the probability of this
error using concentration inequality and the complexity of the function space
of $f^{m}_{h}$. To use smooth function space, we propose to construct a
function out of neural networks such that it lies on RKHS (Section 5). Then,
we bound the probability of deviation-from-mean error through the covering
number of the RKHS space (Section 6.1), then control complexity (Section 6.2)
and prove consistency of the proposed estimator (Section 7).

##  5 Constructing $f$ in RKHS

The following theorem due to [15] paves a way for us to construct a neural
function in RKHS.

######  Theorem 1.

[[15] Appendix A] A function $f\in\mathcal{L}_{2}(d\rho)$ is in Reproducing
Kernel Hilbert Space, $\mathcal{H}_{K}$, if and only if it can be expressed as

$\displaystyle\forall
x\in\mathcal{X},f(x)=\int_{\mathcal{W}}g(w)\psi(x,w)d\tau(w),$ (6)

for a certain function $g:\mathcal{W}\to\mathbb{R}$ such that
$||g||^{2}_{\mathcal{L}_{2}(d\tau)}<\infty$. The RKHS norm of $f$ satisfies
$||f||^{2}_{\mathcal{H}_{K}}\leq||g||^{2}_{\mathcal{L}_{2}(d\tau)}$ and the
kernel $K$ is given by

$\displaystyle K(x,t)=\int_{\mathcal{W}}\psi(x,w)\psi(t,w)d\tau(w)$ (7)

Theorem 1 not only gives us a condition when a square integrable function is
guaranteed to lie in RKHS, it also provides us with a recipe to construct a
function in RKHS. We use this theorem with the neural networks as $\psi$ and
$g$. We sample $w\sim\mathcal{N}(0,\gamma\text{I})$ and pass it through two
neural networks, $\psi$ and $g$, where $\psi$ takes $x$ and $w$ as two
arguments and $g$ takes only $w$ as an argument. More precisely, we consider
$\psi(x,w)=\phi_{\theta}(x)^{T}w$. The kernel $K$, as defined in eq. (7), can
be obtained as:

$\displaystyle K_{\theta}(x^{*},t^{*})$
$\displaystyle=\int_{\mathcal{W}}\phi_{\theta}(x^{*})^{T}ww^{T}\phi_{\theta}(t^{*})d\tau(w)=\gamma\phi_{\theta}(x^{*})^{T}\phi_{\theta}(t^{*})$
(8)

where $E_{w\sim\mathcal{N}(0,\gamma\text{I})}[ww^{T}]=\gamma\text{I}$. We
sometimes denote the kernel $K$ by $K_{\theta}$ to emphasize that it is a
function of neural network parameters, $\theta$.

Traditionally, kernel $K$ remains fixed and the norm of the function $f$
determines the complexity of the function space. In our formulation, both the
RKHS kernel and its norm with respect to the kernel change during training
since the kernel depends on neural network parameters, $\theta$. Therefore,
the challenge is to tease out how neural parameters, $\theta$, affect the
deviation-from-mean error in eq. (5).

##  6 Error Analysis and Control

Assumptions: Before starting our analysis, we list assumptions upon which our
theory is based.
A1. The input domains $\mathcal{X}$ and $\mathcal{W}$ are compact.
A2. The functions $\phi_{\theta}$ and $g$ are Lipschitz continuous with
Lipschitz constants $L_{\phi}$ and $L_{g}$ respectively.
A3. Higher order derivatives $D_{x}^{\alpha}K(x,t)$ up to some high order
$\tau=h/2$ of kernel $K$ exist.

Assumptions A1 is satisfied in our experiments since we consider a bounded set
in $\mathbb{R}^{n}$ and $\mathbb{R}^{D}$ as our domains. Similarly, A2 is
satisfied since we enforce Lipschitz continuity of $\phi$ and $g$ by using
spectral normalization [16]. Assumption A3 is a bit subtle. By the definition
of $K$ in eq.(8), higher order derivative of $K$ exists iff higher order
derivative of $\phi_{\theta}$ exists. This is readily satisfied by deep
networks with smooth activation functions, and is true everywhere except at
origin for ReLU activation. Using the boundedness of the input domain and
Lipschitz continuity, we show the following:

######  Proposition 1.

Under the assumptions A1, A2, we have $\underset{K_{\theta}}{sup}$
$K_{\theta}(x,t)<\infty$ and $||g||^{2}_{\mathcal{L}_{2}(d\tau)}<\infty$.

###  6.1 Bounding the Error Probability of KL Estimates

Bounding the probability of deviation-from-mean error (eq. (5)) is tricky
since, in our case, the kernel is not fixed and we are also optimizing over
them. We bound it in two steps: 1) we derive a bound for a fixed kernel, 2) we
take supremum of this bound over all the kernels parameterized by $\theta$.

For a fixed kernel, we first bound the probability of deviation-from-mean
error in terms of the covering number in Lemma 1. We then use an estimate of
the covering number of RKHS due to [14] to relate the bound to kernel
$K_{\theta}$ in Theorem 2, identifying the role of neural networks in this
error bound.

######  Lemma 1.

Let $f^{m}_{\mathcal{H}_{K}}$ be the optimal discriminator function in an RKHS
$\mathcal{H}_{K}$ which is M-bounded. Let
${KL}_{m}(f^{m}_{\mathcal{H}_{K}})=\frac{1}{m}\sum_{i}f^{m}_{\mathcal{H}_{K}}(x_{i})$
and $KL(f^{m}_{\mathcal{H}_{K}})=E_{p(x)}[f^{m}_{\mathcal{H}_{K}}(x)]$ be the
estimate of KL divergence from m samples and that by using the true
distribution $p(x)$ respectively. Then the probability of error at some
accuracy level, $\epsilon$, is lower-bounded as:

$\displaystyle\text{Prob.}($
$\displaystyle|{KL}_{m}(f^{m}_{\mathcal{H}_{K}})-{KL}(f^{m}_{\mathcal{H}_{K}})|\leq\epsilon)\geq
1-2\mathcal{N}(\mathcal{H}_{K},\frac{\epsilon}{4\sqrt{S_{K}}})\exp(-\frac{m\epsilon^{2}}{4M^{2}})$

where $\mathcal{N}(\mathcal{H}_{K},\eta)$ denotes the covering number of an
RKHS space $\mathcal{H}_{K}$ with disks of radius $\eta$, and
$S_{K}=\underset{x,t}{sup}$ ${K(x,t)}$ which we refer to as kernel complexity.

###### Proof Sketch.

We cover RKHS with discs of radius $\eta=\frac{\epsilon}{4\sqrt{S_{K}}}$.
Within this radius, the deviation does not change too much. So, we can bound
deviation probability at the center of disc and apply union bound over all the
discs. To bound deviation probability at the center, we apply Hoeffding’s
inequality and applying union bound simply leads to counting number of discs
which is exactly the covering number. See supplementary materials for the full
proof. ∎

Lemma 1 bounds the probability of error in terms of the covering number of the
RKHS space. Note that the radius of the disc is inversely related to $S_{K}$
which indicates how complex the RKHS space defined by the kernel $K_{\theta}$
is. Here $K_{\theta}$ depends on the neural network parameters $\theta$.
Therefore, we denote $S_{K}$ as a function of $\theta$ as $S_{K}(\theta)$ and
term it kernel complexity. Next, we use Lemma 2 due to [14] to obtain an error
bound in estimating KL divergence with finite samples in Theorem 2.

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

Let $K:\mathcal{X}\times\mathcal{X}\to{\mathbb{R}}$ be a
$\mathcal{C}^{\infty}$ Mercer kernel and the inclusion
$I_{K}:\mathcal{H}_{K}\xhookrightarrow{}\mathcal{C}(\mathcal{X})$ be the
compact embedding defined by $K$ to the Banach space
$\mathcal{C}(\mathcal{X})$. Let $B_{R}$ be the ball of radius $R$ in RKHS
$\mathcal{H}_{K}$. Then $\forall\eta>0,R>0,h>n$, we have

$\displaystyle\ln\mathcal{N}(I_{K}(B_{R}),\eta)\leq\left(\frac{RC_{h}}{\eta}\right)^{\frac{2n}{h}}$
(9)

where $\mathcal{N}$ gives the covering number of the space $I_{K}(B_{R})$ with
discs of radius $\eta$, and $n$ represents the dimension of the input space
$\mathcal{X}$. $C_{h}$ is given by $C_{h}=C_{s}\sqrt{||\mathscr{L}_{s}||}$
where $\mathscr{L}_{s}$ is a linear embedding from square integrable space
$\mathcal{L}_{2}(d\rho)$ to the Sobolev space $H^{h/2}$ and $C_{s}$ is a
constant.

To prove Lemma 2 [14], the RKHS space is embedded in the Sobolev Space
$H^{h/2}$ using $\mathscr{L}_{s}$ and then the covering number of the Sobolev
space is used. Thus the norm of $\mathscr{L}_{s}$ and the degree of Sobolev
space, $h/2$, appears in the covering number of a ball in $\mathcal{H}_{K}$.
In Theorem 2, we use Lemma 1 and 2 to bound the estimation error of KL
divergence.

######  Theorem 2.

Let ${KL}(f^{m}_{\mathcal{H}})$ and ${KL}_{m}(f^{m}_{\mathcal{H}})$ be the
estimates of KL divergence obtained by using true distribution $p(x)$ and $m$
samples respectively as described in Lemma 1, then the probability of error in
the estimation at the error level $\epsilon$ is given by:

$\displaystyle\text{Prob.}($
$\displaystyle|{KL}_{m}(f^{m}_{\mathcal{H}})-{KL}(f^{m}_{\mathcal{H}})|\leq\epsilon)\geq
1-2\exp\Bigg{[}\left(\frac{4RC_{p}\sqrt{S_{p}||\mathscr{L}_{p}||}}{\epsilon}\right)^{\frac{2n}{h}}-\frac{m\epsilon^{2}}{4M^{2}}\Bigg{]}$

where
$C_{p}\sqrt{S_{p}||\mathscr{L}_{p}||}=\underset{K_{\theta}}{sup}C_{s}\sqrt{S_{K}(\theta)||\mathscr{L}_{s}||}$,
i.e. $C_{p},S_{p},\mathscr{L}_{p}$ correspond to a kernel for which the bound
is maximum.

###### Proof.

We prove this in two steps: First we obtain an error bound for a fixed kernel
space and apply supremum over all $\theta$. For any RKHS
$\mathcal{H}_{K_{\theta}}$, with fixed kernel $K_{\theta}$, we have

$\displaystyle\text{Prob.}($
$\displaystyle|{KL}_{m}(f^{m}_{\mathcal{H}_{K_{\theta}}})-{KL}(f^{m}_{\mathcal{H}_{K_{\theta}}})|\geq\epsilon)\leq
2\exp\Bigg{[}\left(\frac{4RC_{s}\sqrt{S_{K}(\theta)||\mathscr{L}_{s}||}}{\epsilon}\right)^{\frac{2n}{h}}-\frac{m\epsilon^{2}}{4M^{2}}\Bigg{]}$
(10)

We prove this error bound as follows. Lemma 2 gives the covering number of an
RKHS ball of radius $R$, which we apply to Lemma 1. We fix the radius of discs
to $\eta=\frac{\epsilon}{4\sqrt{S_{K}}}$ in Lemma 1 and substitute
$C_{h}=C_{s}\sqrt{||\mathscr{L}_{s}(\theta)||}$ to obtain eq.(10).

Since we are continuously changing $\theta$ during training, the kernel also
changes. Hence, to find the upper bound over all possible kernels, we take the
supremum over all kernels.

$\displaystyle\text{Prob.}(|{KL}_{m}(f^{m}_{\mathcal{H}})-{KL}(f^{m}_{\mathcal{H}})|\geq\epsilon)$
$\displaystyle\leq\underset{K_{\theta}}{sup}\hskip
5.69046pt\text{Prob.}(|{KL}_{m}(f^{m}_{\mathcal{H}_{K_{\theta}}})-{KL}(f^{m}_{\mathcal{H}_{K_{\theta}}})|\geq\epsilon)$
(11) $\displaystyle\leq
2\exp\Bigg{[}\left(\frac{4RC_{p}\sqrt{S_{p}||\mathscr{L}_{p}||}}{\epsilon}\right)^{\frac{2n}{h}}-\frac{m\epsilon^{2}}{4M^{2}}\Bigg{]}$
(12)

where $S_{p}=S_{K}(\theta_{p})$ and
$\mathscr{L}_{p}=\mathscr{L}_{K}(\theta_{p})$, i.e., $S_{p}$ and
$\mathscr{L}_{p}$ correspond to kernel complexity and Sobolev operator norm
corresponding to optimal kernel $K_{\theta_{p}}$ that extremizes eq. (11).
Theorem statement readily follows from eq. (12) ∎

Theorem 2 shows that the error increases exponentially with the radius of the
RKHS space, $R$, complexity of the kernel $S_{K}(\theta_{p})$, and the norm of
the Sobolev space embedding operator $||\mathscr{L}_{p}||$. The Sobolev
embedding operator, $\mathscr{L}_{p}$, is a mapping from
$\mathcal{L}_{2}(d\rho)$ to the Sobolev space $H^{h/2}$. It can be shown [14]
that the operator norm can be bounded as
$||\mathscr{L}_{p}||\leq\rho(\mathcal{X})\sum_{|\alpha|\leq
h/2}\underset{x,t\in\mathcal{X}}{sup}(D^{\alpha}_{x}K_{\theta_{p}}(x,t))^{2},$
where $\rho$ is the measure of the input space $\mathcal{X}$. Therefore, the
norm $||\mathscr{L}_{p}||$ directly measures smoothness of $K_{\theta_{p}}$ in
terms of norm of its derivative in addition to the supremum value of $K$,
while $S_{K}(\theta_{p})$ only depends on the supremum value of
$K_{\theta_{p}}$.

###  6.2 Complexity Control

From Theorem 2, we see that the error probability could be decreased by
decreasing $R,||\mathscr{L}_{p}||$ and $S_{K}(\theta_{p})$. Using argument
similar to the proof of Proposition 1, we can show that the Lipschitz
constraint on $\phi_{\theta}$ also affects $S_{K}$ and may affect
$||\mathscr{L}_{p}||$. In our experiments, however, we fix the Lipschitz
constraints during optimization and do not change $S_{K}$ and
$||\mathscr{L}_{p}||$ dynamically. Here, we focus on the norm, $R$ from
Theorem 2. To obtain the optimal discriminator $f^{m}_{h}$, we optimize the
following objective with an extra penalization of the upper bound, i.e.
$||g||$ on the RKHS norm of $f$:

$\displaystyle f^{m}_{h}=\underset{f\in
h}{\operatorname{argmax}}{\frac{1}{m}\sum_{x_{i}\sim
p(x_{i})}\log\sigma(f(x_{i}))+\frac{1}{m}\sum_{x_{j}\sim
q(x_{j})}\log(1-\sigma(f(x_{j})))}-\frac{\lambda_{0}}{m}||g||_{\mathcal{L}_{2}(d\tau)}^{2}$
(13)

The regularization term prevents the radius of RKHS ball from growing,
maintaining a low error probability. Optimization of eq. (13) w.r.t. neural
network parameters $\theta$ allows dynamic control of the complexity of the
discriminator function on the fly in a scalable and efficient way. Note that,
computation of $||g||_{\mathcal{L}_{2}(d\tau)}$ requires randomly sampling
$w\sim\mathcal{N}(0,\gamma\textbf{I})$ and passing through neural network $g$
independent of the data $x_{i},x_{j}$. Therefore, if the computational
complexity of optimization is O(m), it will remain the same after
incorporating this additional term, i.e. regularization does not increase
asymptotic time complexity which is linear with the number of samples, m.

Figure 1: a) Top scatter plot compares KL divergence estimates between a
method using Neural network discriminator without complexity control (red) and
that using RKHS discriminator with compelxity control (blue); b) In the
bottom, we show the effect of varying the regularization parameter $\lambda$
on bias and variance while using the RKHS discriminator with complexity
control as in eq.(13).

##  7 Variance and Consistency of the Estimate

###  7.1 Variance Analysis

Theorem 2 gives an upper bound on the probability of error. Intuitively, the
variance and probability of error behave similarly for many distributions,
i.e. higher variance might indicate higher probability of error. Below we
quantify this intuition for a Gaussian distributed estimate:

######  Theorem 3.

Let $X=KL_{m}(f^{m}_{\mathcal{H}})$ be the estimated KL divergence using m
samples as described in Theorem 2. Assuming that $X$ follows a Gaussian
distribution $X\sim\mathcal{N}(\mu,\sigma)$, we can obtain an upper bound on
this variance of the estimate as follows:

$\displaystyle\sigma$
$\displaystyle\leq\frac{\epsilon}{\text{erf}^{-1}\Big{[}-4\exp\Big{[}\left(\frac{4RC_{p}\sqrt{S_{p}||L_{p}||}}{\epsilon}\right)^{\frac{2n}{h}}-\frac{m\epsilon^{2}}{4M^{2}}\Big{]}+1\Big{]}}$
(14)

where erf is the Gauss error function and is a monotonic function.

Obviously, this relation applies only to Gaussian distributed estimate, a
strong assumption. However, Theorem 3 is presented for illustrative purpose.
It suggests that by decreasing $R$, the radius of the RKHS ball, the variance
of the estimate could be decreased. Experimentally, we observe that the
variance decreases as we penalize the RKHS norm more, consistent with the
spirit of Theorem 3.

###  7.2 Consistency of Estimates

Here we show that the regularized objective leads to a consistent estimation.

######  Theorem 4.

Let $f^{*}$ and $f^{m}$ be optimal discriminators as described in eq. (1) and
eq. (13) respectively, and the KL estimate is given by
$KL(f)=E_{p(x)}[f(x)],\hskip 11.38092ptKL_{m}(f)=\frac{1}{m}\sum_{x_{i}\sim
p(x_{i})}[f(x)]$. Then, in the limiting case as $m\to\infty$,
$|KL_{m}(f^{m}_{h})-KL(f^{*})|\to 0$.

###### Proof Sketch.

The difference between the true KL divergence and the estimated KL divergence
can be divided into three terms as shown in eq. (5). We assume that our
function space is rich enough to contain the true solution, driving bias to
zero. From Theorem 2, we see that in the limiting case of $m\to 0$, the
deviation-from-mean error goes to $0$. Therefore, the key step that remains to
be shown is that the discriminator induced error (second term in eq.(5)) also
goes to 0 as $m\to\infty$.

It can be shown if we can prove that the optimal discriminator in eq. (13)
approaches the optimal discriminator in eq. (2). To prove this, we show that
the argument being maximized by $f^{m}_{h}$ approaches the argument being
maximized by $f_{h}^{*}$ in the limiting case. To show this, we need to show
that the function space, $\log\sigma f$, is Glivenko Cantelli [17], which we
prove in following steps:
1\. We show that $f$ is Lipschitz continuous by definition and due to
Lipschitz continuity of $\phi_{\theta}$. Then we show that $\log\sigma f$ is
Lipschitz continuous if $f$ is Lipschitz continuous.
2\. Then we show that for a class of functions with Lipschitz constant $L$,
the metric entropy, $\log N$, can be obtained in terms of $L$ and entropy
number of the bounded input space, $\mathcal{X}$.
3\. Since the metric entropy does not grow with the number of samples $m$, we
show that $\frac{1}{m}\log N\to 0$ which lets us show that $\log\sigma f$
belongs to Glivenko Cantelli class of functions by using Theorem 2.4.3 from
[17]. See supplementary material for the complete proof. ∎

##  8 Experimental Results

We present results on three applications of KL divergence estimation: 1. KL
estimation between simple Guaussian distributions, 2. Mutual information
estimation, 3. Variational Bayes. In our experiments, the RKHS discriminator
is constructed with $\psi$ and $g$ networks as described in Section 5, where
the network $\psi$ is very close to a regular neural network. In two
experiments, we compare our results with the models using regular neural net
discriminator to ensure that the difference in performance between RKHS and
regular neural network is not due to architectural difference.

Figure 2: Comparing our method with CPC [18], convex risk minimization(NWJ)
[1] and SMILE [9] regarding mutual information estimation between two
variables.

##### KL Estimation between Two Gaussians

We assume that we have finite sets of samples from two distributions. We
further assume that we are required to apply minibatch based optimization. We
consider estimating KL divergence between two Gaussian distributions in 2D,
where we know the analytical KL divergence between the two distributions as
the ground truth. We consider three different pairs of distributions
corresponding to true KL divergence values of $1.3,13.8$ and $38.29$,
respectively and use $m=5000$ samples from each distribution to estimate KL in
the finite case. We repeat the estimation experiments with random
initialization 30 times and report the mean, standard deviation, scatter and
box plots.

Fig. 1 top row compares the estimation of KL divergence with regular neural
net and RKHS discriminator with complexity control based on eq. (13). With our
proposed RKHS discriminator, the KL estimates are significantly more reliable
and accurate: error reduced from 0.5 to 0.04, 5.8 to 1.07 and 60.6 to 9.7 and
variance reduced from 0.2 to 0.002, 223 to 4.4 and 3521 to 33 for true KL 1.3,
13.8 and 38.29 respectively. In Fig. 1 bottom row, we investigate our
complexity control method on the effect of varying the regularization
parameter $\lambda=\lambda_{0}/m$. As expected, increasing regularization
parameter penalizes more on the RKHS norm and therefore reduces variance. This
is consistent with our theory. Regarding bias, however, as we increase the
$\lambda$, the bias decreases and then starts to increase. Hence, one needs to
strike a balance between bias and variance while choosing $\lambda$.

##### Mutual Information Estimation

Computation of mutual information is a direct use case of KL divergence
computation. We replicate the experimental setup of [19, 9] to estimate mutual
information between $(x,y)$ drawn from 20-d Gaussian distributions, where the
mutual information is increased by step size of 2 from 2 to 10. We compare the
performance of our method with traditional KL divergence computation methods
like contrastive predictive coding (CPC) [18], convex risk minimization (NWJ)
[1] and SMILE [9]. In Fig.2, our method with RKHS discriminator (with
$\lambda=1e^{-5}$) performs better than CPC [18] and NWJ [1], and is
competitive with the state-of-the-art, SMILE [9]. In the bottom row, we also
show the effect of regularization parameter $\lambda$ in our method. Similar
to the previous experiment, increasing the regularization parameter decreases
the variance and increases the bias. It is consistent with our theoretical
insights about the effect of reducing RKHS norm on variance.

Figure 3: (a) Comparison of MNIST digit reconstruction using AVB autoencoder
model [8]. Trace of KL divergence and reconstruction loss in AVB model with
Neural network discriminator (b) and RKHS discriminator in (c).

##### Adversarial Variational Bayes

Variational Bayes requires KL divergence estimation. When we do not have
access to analytical form of the posterior/prior distributions, but only have
access to the samples, we need to estimate KL divergence from samples.
Adversarial Variational Bayes (AVB) [8] presents a way to achieve this using a
discriminator network. We adopt this setup and demonstrate that the training
becomes unstable if we do not constrain the complexity of the discriminator.
First, we train AVB on MNIST dataset with a simple neural network
discriminator architecture. As the training progresses, the KL divergence
blows up after about 500 epochs (Fig. 3(b)) and the reconstruction starts to
get worse (Fig. 3(a)). We modify the same architecture according to our
construction such that it lies in RKHS and then penalize the RKHS norm as in
eq. (13). It stabilizes the training for a large number of epochs and the
reconstruction does not deteriorate as the training progresses, resulting into
sharp reconstruction (Fig. 3(a)). We want to clarify that this instability in
training neural net discriminator is present if we use a basic discriminator
architecture. It does not mean that there exists no other method to design a
stable neural net discriminator. In fact, AVB [8] presents a discriminator
that adds additional inner product structure to stabilize the discriminator
training. Our point here is that we can stabilize the training by ensuring
that the discriminator lies in a well behaved function space (the RKHS) and
controlling its complexity, consistent with our theory.

##  9 Limitations, Discussion and Conclusion

##### Limitations:

The proposed construction of neural function in RKHS exhibits good properties
of both the deep learning and kernel methods. However, it requires
constructing two separate deep networks, $\psi$ and $g$. It makes our model a
bit bulky and also requires more parameter due to additional $g$. Moreover,
currently our RKHS discriminator’s output is scalar; generalizing this
function to a multivariable output could make our model bulkier and increase
parameters even more. Second limitation is the requirement of higher order
derivative of kernel $K$ in assumption A3. While this requirement is satisfied
if smooth activation function is used in $\phi_{\theta}$, for activations like
ReLU or LeakyReLU, the derivatives exist everywhere except at the origin. In
these cases, we need to carefully investigate if we can use subgradients to
define operator norm $||\mathscr{L}_{p}||$.

##### Discussion and Conclusion:

We have shown that using a regular neural network as a discriminator in
estimating KL divergence results in unreliable estimation if the complexity of
the function space is not controlled. We then showed a solution by
constructing a discriminator function in RKHS space using neural networks and
penalizing its complexity in a scalable way. Although the idea to use RKHS
norm to penalize complexity is not new (see for example [1]), it is not clear
how to use this idea directly on the function $f$. In traditional kernel
methods, algorithms often do not work with RKHS function $f$ directly, but
rather work with kernel matrix, $K$ by using, for example, the Representer
Theorem [20]. In the case of big data, working with the big kernel matrix is
computationally expensive although some methods have been proposed to speed up
the computation, like Random Fourier Feature [21]. We propose a different view
by directly constructing a function in RKHS space, which led us to scalable
algorithm while incorporating the advantages of neural networks. Moreover, our
representation could also be seen as an improvement over RFF by using neural
basis, $\psi$, instead of Fourier basis. The idea of constructing a neural-net
function in RKHS and complexity control could also be useful in stabilizing
GANs in general. Currently, the most successful way to stabilize GANs is to
enforce smoothness by gradient penalization [22, 23, 24]. On the light of the
present analysis, gradient penalty could also be thought as a way to control
the complexity of the discriminator.

## References

  * [1] X. Nguyen, M. J. Wainwright, and M. I. Jordan, “Estimating divergence functionals and the likelihood ratio by convex risk minimization,” IEEE Transactions on Information Theory, vol. 56, no. 11, pp. 5847–5861, 2010.
  * [2] T. Q. Chen, X. Li, R. B. Grosse, and D. K. Duvenaud, “Isolating sources of disentanglement in variational autoencoders,” in Advances in Neural Information Processing Systems, pp. 2610–2620, 2018.
  * [3] L. Mescheder, A. Geiger, and S. Nowozin, “Which training methods for gans do actually converge?,” in International Conference on Machine Learning, pp. 3481–3490, 2018.
  * [4] X. Chen, Y. Duan, R. Houthooft, J. Schulman, I. Sutskever, and P. Abbeel, “Infogan: Interpretable representation learning by information maximizing generative adversarial nets,” in Advances in neural information processing systems, pp. 2172–2180, 2016.
  * [5] C. K. Sønderby, J. Caballero, L. Theis, W. Shi, and F. Huszár, “Amortised map inference for image super-resolution,” ICLR, 2017.
  * [6] S. Nowozin, B. Cseke, and R. Tomioka, “f-gan: Training generative neural samplers using variational divergence minimization,” in Advances in neural information processing systems, pp. 271–279, 2016.
  * [7] M. I. Belghazi, A. Baratin, S. Rajeshwar, S. Ozair, Y. Bengio, A. Courville, and D. Hjelm, “Mutual information neural estimation,” in International Conference on Machine Learning, pp. 531–540, 2018.
  * [8] L. Mescheder, S. Nowozin, and A. Geiger, “Adversarial variational bayes: Unifying variational autoencoders and generative adversarial networks,” in International Conference on Machine Learning (ICML), 2017.
  * [9] J. Song and S. Ermon, “Understanding the limitations of variational mutual information estimators,” in International Conference on Learning Representations, 2020.
  * [10] C. Zhang, S. Bengio, M. Hardt, B. Recht, and O. Vinyals, “Understanding deep learning requires rethinking generalization,” arXiv preprint arXiv:1611.03530, 2016.
  * [11] K. Ahuja, “Estimating kullback-leibler divergence using kernel machines,” in 2019 53rd Asilomar Conference on Signals, Systems, and Computers, pp. 690–696, IEEE, 2019.
  * [12] A. Berlinet and C. Thomas-Agnan, Reproducing Kernel Hilbert Spaces in Probability and Statistics.  Springer US, 2011.
  * [13] F. Bach, “On the equivalence between kernel quadrature rules and random feature expansions,” The Journal of Machine Learning Research, vol. 18, no. 1, pp. 714–751, 2017.
  * [14] F. Cucker and S. Smale, “On the mathematical foundations of learning,” Bulletin of the American mathematical society, vol. 39, no. 1, pp. 1–49, 2002\.
  * [15] F. Bach, “Breaking the curse of dimensionality with convex neural networks,” The Journal of Machine Learning Research, vol. 18, no. 1, pp. 629–681, 2017\.
  * [16] T. Miyato, T. Kataoka, M. Koyama, and Y. Yoshida, “Spectral normalization for generative adversarial networks,” in International Conference on Learning Representations, 2018.
  * [17] A. W. Van Der Vaart and J. A. Wellner, “Weak convergence,” in Weak convergence and empirical processes, Springer, 1996.
  * [18] A. v. d. Oord, Y. Li, and O. Vinyals, “Representation learning with contrastive predictive coding,” arXiv preprint arXiv:1807.03748, 2018.
  * [19] B. Poole, S. Ozair, A. Van Den Oord, A. Alemi, and G. Tucker, “On variational bounds of mutual information,” in International Conference on Machine Learning, pp. 5171–5180, PMLR, 2019.
  * [20] B. Schölkopf, A. J. Smola, F. Bach, et al., Learning with kernels: support vector machines, regularization, optimization, and beyond.  MIT press, 2002.
  * [21] A. Rahimi, B. Recht, et al., “Random features for large-scale kernel machines.,” in Neural Information Processing Systems, vol. 3, p. 5, Citeseer, 2007.
  * [22] M. Arjovsky, S. Chintala, and L. Bottou, “Wasserstein generative adversarial networks,” in International Conference on Machine Learning, pp. 214–223, 2017.
  * [23] I. Gulrajani, F. Ahmed, M. Arjovsky, V. Dumoulin, and A. C. Courville, “Improved training of wasserstein gans,” in Advances in Neural Information Processing Systems 30 (I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett, eds.), pp. 5767–5777, Curran Associates, Inc., 2017.
  * [24] M. Bińkowski, D. J. Sutherland, M. Arbel, and A. Gretton, “Demystifying MMD GANs,” in International Conference on Learning Representations, 2018\.
  * [25] I. Goodfellow, J. Pouget-Abadie, M. Mirza, B. Xu, D. Warde-Farley, S. Ozair, A. Courville, and Y. Bengio, “Generative adversarial nets,” in Advances in neural information processing systems, pp. 2672–2680, 2014.
  * [26] A. Kontorovich and R. Weiss, “Maximum margin multiclass nearest neighbors,” International Conference on Machine Learning, vol. abs/1401.7898, 2014.
  * [27] L.-A. Gottlieb, A. Kontorovich, and R. Krauthgamer, “Efficient regression in metric spaces via approximate lipschitz extension,” in Similarity-Based Pattern Recognition (E. Hancock and M. Pelillo, eds.), (Berlin, Heidelberg), pp. 43–58, Springer Berlin Heidelberg, 2013.
  * [28] E. J. McShane, “Extension of range of functions,” Bulletin of the American Mathematical Society, vol. 40, no. 12, pp. 837 – 842, 1934.
  * [29] H. Whitney, “Analytic extensions of differentiable functions defined in closed sets,” Transactions of the American Mathematical Society, vol. 36, no. 1, pp. 63–89, 1934.
  * [30] R. Krauthgamer and J. R. Lee, “Navigating nets: Simple algorithms for proximity search,” in Proceedings of the fifteenth annual ACM-SIAM symposium on Discrete algorithms, pp. 798–807, 2004.

##  Appendix A Problem Formulation and Contribution

##### GAN-type Objective for KL Estimation

Let $f$ be a discriminator, $f:\mathcal{X}\to{\rm I\\!R}$. Let $p(x)$ and
$q(x)$ be two probability density functions defined over the space
$\mathcal{X}$. First, we train a discriminator as:

$\displaystyle
f^{*}=\underset{f}{\operatorname{argmax}}[{E_{p(x)}\log\sigma(f(x))+E_{q(x)}\log(1-\sigma(f(x)))}]$
(15)

where $\sigma$ is the Sigmoid function given by
$\sigma(x)=\frac{e^{x}}{1+e^{x}}$. Then the KL divergence $KL(p(x)||q(x))$ is
given by:

$\displaystyle KL(p(x)||q(x))=E_{p(x)}[f^{*}(x)]$ (16)

###### Proof.

The proof is based on similar proofs in [8, 5] and presented here for the sake
of completeness.

We rewrite the objective as :

$\displaystyle\int p(x)\log\sigma(f(x))+q(x)\log(1-\sigma(f(x)))dx$ (17)

This integral is maximum with respect to $f$ if and only if the integrand is
maximal for every $x$. As argued in the Proposition 1 of [25], the function

$\displaystyle t\mapsto a\log(t)+b\log(1-t)$ (18)

attains its maximum at $t=\frac{a}{a+b}$ showing that,

$\displaystyle\sigma(f^{*}(x))=\frac{p(x)}{p(x)+q(x)}$ (19)

Plugging the expression for Sigmoid function, we obtain,

$\displaystyle f^{*}(x)=\frac{p(x)}{q(x)}$ (20)

Therefore, by the definition of KL divergence, we have:

$\displaystyle KL(p(x)||q(x))=E_{p(x)}[\frac{p(x)}{q(x)}]=E_{p(x)}[f^{*}(x)]$
(21)

∎

##  Appendix B Error Analysis and Control

We start with the set of assumptions based on which our theory is developed.

  1. A1.

The input domains $\mathcal{X}$ and $\mathcal{W}$ are compact.

  2. A2.

The functions $\phi_{\theta}$ and $g$ are Lipschitz continuous with Lipschitz
constant $L_{\phi}$ and $L_{g}$ respectively.

  3. A3.

Higher order derivatives $D_{x}^{\alpha}K(x,t)$ of kernel $K$ exist up to some
high order $\tau=h/2$ .

######  Proposition 2.

Under the assumptions A1, A2, we have
i) $\underset{x,t}{sup}$ $K_{\theta}(x,t)<\infty$, and
ii) $||g||^{2}_{\mathcal{L}_{2}(d\tau)}<\infty$.

###### Proof.

i) By the definition
$K_{\theta}(x,t)=\gamma\langle\phi_{\theta}(x),\phi_{\theta}(t)\rangle$. Using
Cauchy Schwartz,

$\displaystyle K_{\theta}(x,t)$
$\displaystyle\leq\gamma||\phi_{\theta}(x)||||\phi_{\theta}(t)||$ (22)
$\displaystyle\leq\gamma L_{\phi}||x||L_{\phi}||t||$ (23)
$\displaystyle<\infty$ (24)

where we used the fact that $\mathcal{X}$ is bounded, and therefore, $||x||$
and $||t||$ are finite.
ii) By definition,

$\displaystyle||g||^{2}_{\mathcal{L}_{2}(d\tau)}$ $\displaystyle=\int
g(w)^{2}d\tau(w)$ (25) $\displaystyle\leq\int L_{g}^{2}||w||^{2}d\tau(w)$ (26)
$\displaystyle=L_{g}^{2}tr(C_{w})$ (27)

where $C_{w}$ is the uncentered covariance matrix of the Gaussian distributed
w. Therefore, we immediately obtain
$||g||^{2}_{\mathcal{L}_{2}(d\tau)}<\infty$. ∎

These results are useful in constructing a function $f$ in RKHS in Theorem 1
(Section 5) of the main paper.

###  B.1 Bounding the Error Probability of KL Estimates

We bound the deviation-from-mean error in two steps: 1) we derive a bound for
a fixed kernel, 2) we take supremum of this bound over all the kernels
parameterized by $\theta$.

For a fixed kernel, we first bound the probability of deviation-from-mean
error in terms of the covering number in Lemma 1. Then, we use an estimate of
the covering number of RKHS due to [14] to obtain a bound of error probability
in terms of the kernel $K_{\theta}$ in Lemma 3. Note that, Lemma 3 is proved
for a fixed kernel $K_{\theta}$, where $\theta$ is fixed. Then finally in
Theorem 2, we take supremum over all kernels $K_{\theta}$s to obtain a bound
on error probability on a space of functions with all possible kernels.

######  Lemma 1.

Let $f^{m}_{\mathcal{H}_{K}}$ be the optimal discriminator function in a RKHS
$\mathcal{H}_{K}$ which is M-bounded. Let
${KL}_{m}(f^{m}_{\mathcal{H}_{K}})=\frac{1}{m}\sum_{i}f^{m}_{\mathcal{H}_{K}}(x_{i})$
and $KL(f^{m}_{\mathcal{H}_{K}})=E_{p(x)}[f^{m}_{\mathcal{H}_{K}}(x)]$ be the
estimate of KL divergence from m samples and that by using true distribution
$p(x)$ respectively. Then the probability of error at some accuracy level,
$\epsilon$ is lower-bounded as:

$\displaystyle\text{Prob.}($
$\displaystyle|{KL}_{m}(f^{m}_{\mathcal{H}_{K}})-{KL}(f^{m}_{\mathcal{H}_{K}})|\leq\epsilon)\geq
1-2\mathcal{N}(\mathcal{H}_{K},\frac{\epsilon}{4\sqrt{S_{K}}})\exp(-\frac{m\epsilon^{2}}{4M^{2}})$

where $\mathcal{N}(\mathcal{H}_{K},\eta)$ denotes the covering number of a
RKHS space $\mathcal{H}_{K}$ with disks of radius $\eta$, and
$S_{K}=\underset{x,t}{sup}$ ${K(x,t)}$ which we refer as kernel complexity

###### Proof.

Let $\ell_{z}(f)=E_{p(x)}[f(x)]-\frac{1}{m}\sum_{i}f(x_{i})$ denotes the error
in the estimate such that we want to bound $|\ell_{z}(f)|$. We have,

$\displaystyle\ell_{z}(f_{1})-\ell_{z}(f_{2})=E_{p(x)}[f_{1}(x)-f_{2}(x)]-\frac{1}{m}\sum_{i}f_{1}(x_{i})-f_{2}(x_{i})$

We know $E_{p(x)}[f_{1}(x)-f_{2}(x)]\leq||f_{1}-f_{2}||_{\infty}$ and
$\frac{1}{m}\sum_{i}f_{1}(x_{i})-f_{2}(x_{i})\leq||f_{1}-f_{2}||_{\infty}$.
Using the triangle inequality, we obtain
$|\ell_{z}(f_{1})-\ell_{z}(f_{2})|\leq 2||f_{1}-f_{2}||_{\infty}$. Now,
consider $f\in\mathcal{H}_{K}$, then,

$\displaystyle|f(x)|=|\langle
K_{x},f\rangle|\leq||f||||K_{x}||=||f||\sqrt{K(x,x)}$ (28)

This implies the RKHS space norm and $\ell_{\infty}$ norm of a function are
related by

$\displaystyle||f||_{\infty}\leq\sqrt{S_{K}}||f||_{\mathcal{H}_{K}}$ (29)

Hence, we have:

$\displaystyle|\ell_{z}(f_{1})-\ell_{z}(f_{2})|\leq
2\sqrt{S_{K}}||f_{1}-f_{2}||_{\mathcal{H}_{K}}$ (30)

The idea of the covering number is to cover the whole RKHS space
$\mathcal{H}_{K}$ with disks of some fixed radius $\eta$, which helps us bound
the error probability in terms of the number of such disks. Let
$\mathcal{N}(\mathcal{H}_{K},\eta)$ be such disks covering the whole RKHS
space. Then, for any function $f$ in $\mathcal{H}_{K}$, we can find some disk,
$D_{j}$ with centre $f_{j}$, such that
$||f-f_{j}||_{\mathcal{H}_{K}}\leq\eta$. If we choose
$\eta=\frac{\epsilon}{2\sqrt{S_{K}}}$, then from eq.(30), we obtain,

$\displaystyle\underset{f\in D_{j}}{sup}{|\ell_{z}(f)|\geq
2\epsilon}\implies|\ell_{z}(f_{j})|\geq\epsilon$ (31)

Using the Hoeffding’s inequality,
$\text{Prob.}(|\ell_{z}(f_{j})|\geq\epsilon)\leq
2e^{-\frac{m\epsilon^{2}}{2M^{2}}}$ and eq.(31),

$\displaystyle\text{Prob.}(\underset{f\in D_{j}}{sup}{|\ell_{z}(f)|\geq
2\epsilon})\leq 2e^{-\frac{m\epsilon^{2}}{2M^{2}}}$ (32)

Applying union bound over all the disks, we obtian,

$\displaystyle\text{Prob.}(\underset{f\in\mathcal{H}}{sup}{|\ell_{z}(f)|\geq
2\epsilon})\leq
2\mathcal{N}(\mathcal{H},\frac{\epsilon}{2\sqrt{S_{K}}})e^{-\frac{m\epsilon^{2}}{2M^{2}}}$
(33)
$\displaystyle\text{Prob.}(\underset{f\in\mathcal{H}}{sup}{|\ell_{z}(f)|\leq\epsilon})\geq
1-2\mathcal{N}(\mathcal{H},\frac{\epsilon}{4\sqrt{S_{K}}})e^{-\frac{m\epsilon^{2}}{4M^{2}}}$

which proves the lemma.

On M-boundedness of $f^{m}_{\mathcal{H}_{K}}$
To prove the lemma, we assumed that $f^{m}_{\mathcal{H}_{K}}$ is M bounded. To
see why this is reasonable, from eq.29, we have
$||f^{m}_{\mathcal{H}_{K}}||_{\infty}\leq\sqrt{S_{K}}||f^{m}_{\mathcal{H}_{K}}||_{\mathcal{H}_{K}}\leq\sqrt{S_{K}}||g||_{\mathcal{L}_{2}(d\rho)}$.
Therefore, $f^{m}_{\mathcal{H}_{K}}$ is bounded if $S_{K}$ and
$||g||_{\mathcal{L}_{2}(d\rho)}$ are bounded, which is true by Proposition 1.
∎

######  Remark 1.

We derived the error bound based on the Hoeffding’s inequality by assuming
that our only knowledge about $f$ is that it is bounded. If we have other
knowledge, for example, if we know the variance of $f$, we could use
Bernstein’s inequality instead of Hoeffding’s inequality with minimal change
to the proof. To the extent we are interested in the contribution of neural
network in error bound, however, there is not much gain by using one
inequality or the other. Hence, we stick with Hoeffding’s inequality and note
other possibilities.

######  Remark 2.

Note that in Lemma 1, the radius of disks are inversely related to the the
quantity, $S_{K}$, meaning that if $S_{K}$ is high, we would need large number
of disks to fill the RKHS space. Hence, it denotes a quantity that reflects
the complexity of the RKHS space. We, therefore, term it kernel complexity.
Also in eq. 29 and the discussion about the M-boundedness, we see that the
maximum value $|f(x)|$ depends on $S_{K}$, again providing insight into how
$S_{K}$ may control both maximum fluctuation and the boundedness.

Lemma 1 bounds the probability of error in terms of the covering number of the
RKHS space. Next, we use Lemma 2 due to [14] to obtain an error bound in
estimating KL divergence with finite samples in Theorem 2.

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

Let $K:\mathcal{X}\times\mathcal{X}\to{\mathbb{R}}$ is a
$\mathcal{C}^{\infty}$ Mercer kernel and the inclusion
$I_{K}:\mathcal{H}_{K}\xhookrightarrow{}\mathcal{C}(\mathcal{X})$ is the
compact embedding defined by $K$ to the Banach space
$\mathcal{C}(\mathcal{X})$ . Let $B_{R}$ be the ball of radius $R$ in RKHS
$\mathcal{H}_{K}$. Then $\forall\eta>0,R>0,h>n$, we have

$\displaystyle\ln\mathcal{N}(I_{K}(B_{R}),\eta)\leq\left(\frac{RC_{h}}{\eta}\right)^{\frac{2n}{h}}$
(34)

where $\mathcal{N}$ gives the covering number of the space $I_{K}(B_{R})$ with
disks of radius $\eta$, and $n$ represents the dimension of inputs space
$\mathcal{X}$. $C_{h}$ is given by

$\displaystyle C_{h}=C\sqrt{||\mathscr{L}_{K}||}$ (35)

where $L_{K}$ is a linear embedding from square integrable space
$\mathcal{L}_{2}(d\rho)$ to the Sobolev space $H^{h/2}$ and $C$ is a constant.

To prove Lemma 2, the RKHS space is embedded in the Sobolev Space $H^{h/2}$
using $\mathscr{L}_{K}$ and then covering number of Sobolev space is used.
Thus the norm of $\mathscr{L}_{K}$ and the degree of Sobolev space, $h/2$,
appears in the covering number of a ball in $\mathcal{H}_{K}$. In Lemma 3, we
use this Lemma to bound the estimation error of KL divergence.

######  Lemma 3.

Let ${KL}(f^{m}_{\mathcal{H}_{K_{\theta}}})$ and
${KL}_{m}(f^{m}_{\mathcal{H}_{K_{\theta}}})$ be the estimates of KL divergence
obtained by using true distribution $p(x)$ and $m$ samples respectively and
using a fixed kernel, ${K_{\theta}}$ as described in Lemma 1, then the
probability of error in the estimation at the error level $\epsilon$ is given
by:

$\displaystyle\text{Prob.}($
$\displaystyle|{KL}_{m}(f^{m}_{\mathcal{H}_{K_{\theta}}})-{KL}(f^{m}_{\mathcal{H}_{K_{\theta}}})|\geq\epsilon)\leq
2\exp\Bigg{[}\left(\frac{4RC_{s}\sqrt{S_{K}(\theta)||\mathscr{L}_{s}(\theta)||}}{\epsilon}\right)^{\frac{2n}{h}}-\frac{m\epsilon^{2}}{4M^{2}}\Bigg{]}$

###### Proof.

Lemma 2 gives the covering number of a ball of radius $R$ in an RKHS space. In
Lemma 1, if we consider the hypothesis space to be a ball of radius $R$, we
can apply Lemma 2 in it. Additionally, since we fix the radius of disks to be
$\eta=\frac{\epsilon}{4\sqrt{S_{K}}}$ in Lemma 1, we obtain,

$\displaystyle\text{Prob.}(|{KL}_{m}(f^{m}_{\mathcal{H}_{K_{\theta}}})-{KL}(f^{m}_{\mathcal{H}_{K_{\theta}}})|\geq\epsilon)\leq
2\exp\Big{[}\left(\frac{4\sqrt{S_{K_{\theta}}}RC_{h}}{\epsilon}\right)^{\frac{2n}{h}}-\frac{m\epsilon^{2}}{4M^{2}}\Big{]}$

Substituting $C_{h}=C\sqrt{||\mathscr{L}_{K_{\theta}}||}$, we obtain,

$\displaystyle\text{Prob.}(|{KL}_{m}(f^{m}_{\mathcal{H}_{K}})-{KL}(f^{m}_{\mathcal{H}_{K_{\theta}}})|\geq\epsilon)\leq
1-2\exp\Bigg{[}\left(\frac{4RC\sqrt{S_{K_{\theta}}||\mathscr{L}_{K_{\theta}}||}}{\epsilon}\right)^{\frac{2n}{h}}-\frac{m\epsilon^{2}}{4M^{2}}\Bigg{]}$
(37)

∎

######  Theorem 2.

Let ${KL}(f^{m}_{\mathcal{H}})$ and ${KL}_{m}(f^{m}_{\mathcal{H}})$ be the
estimates of KL divergence obtained by using true distribution $p(x)$ and $m$
samples respectively as described in Lemma 1, then the probability of error in
the estimation at the error level $\epsilon$ is given by:

$\displaystyle\text{Prob.}($
$\displaystyle|{KL}_{m}(f^{m}_{\mathcal{H}})-{KL}(f^{m}_{\mathcal{H}})|\leq\epsilon)\geq
1-2\exp\Bigg{[}\left(\frac{4RC_{p}\sqrt{S_{p}||\mathscr{L}_{p}||}}{\epsilon}\right)^{\frac{2n}{h}}-\frac{m\epsilon^{2}}{4M^{2}}\Bigg{]}$

where $C_{p}\sqrt{S_{p}||\mathscr{L}_{p}||}=\underset{K_{\theta}}{sup}\hskip
5.69046ptC_{s}\sqrt{S_{K}(\theta)||\mathscr{L}_{s}||}$, i.e.
$C_{p},S_{p},\mathscr{L}_{p}$ correspond to a kernel for which the bound is
maximum.

###### Proof.

Lemma 3 gives an error bound for a fixed kernel, $K_{\theta}$. To find an
upper bound over all possible kernels, we take the supremum over all kernels.

$\displaystyle\text{Prob.}(|{KL}_{m}(f^{m}_{\mathcal{H}})-{KL}(f^{m}_{\mathcal{H}})|\geq\epsilon)$
$\displaystyle\leq\underset{K_{\theta}}{sup}\hskip
5.69046pt\text{Prob.}(|{KL}_{m}(f^{m}_{\mathcal{H}_{K_{\theta}}})-{KL}(f^{m}_{\mathcal{H}_{K_{\theta}}})|\geq\epsilon)$
(38) $\displaystyle\leq
2\exp\Bigg{[}\left(\frac{4RC_{p}\sqrt{S_{p}||\mathscr{L}_{p}||}}{\epsilon}\right)^{\frac{2n}{h}}-\frac{m\epsilon^{2}}{4M^{2}}\Bigg{]}$
(39)

where $S_{p}=S_{K}(\theta_{p})$ and
$\mathscr{L}_{p}=\mathscr{L}_{K}(\theta_{p})$, i.e., $S_{p}$ and
$\mathscr{L}_{p}$ correspond to kernel complexity and Sobolev operator norm
corresponding to optimal kernel $K_{\theta_{p}}$ that extremizes eq. (38).
Theorem statement readily follows from eq. (39) ∎

##  Appendix C Variance and Consistency of the Estimate

###  C.1 Variance Analysis

######  Theorem 3.

Let $X=KL_{m}(f^{m}_{\mathcal{H}})$ be the estimated KL divergence using m
samples as described in Theorem 2. Assuming that $X$ follows a Gaussian
distribution $X\sim\mathcal{N}(\mu,\sigma)$, we can obtain an upper bound on
this variance of the estimate as follows:

$\displaystyle\sigma$
$\displaystyle\leq\frac{\epsilon}{\sqrt{2}\text{erf}^{-1}\Big{[}-4\exp\Big{[}\left(\frac{4RC_{p}\sqrt{S_{p}||\mathscr{L}_{p}||}}{\epsilon}\right)^{\frac{2n}{h}}-\frac{m\epsilon^{2}}{4M^{2}}\Big{]}+1\Big{]}}$
(40)

where erf is the Gauss error function

$\displaystyle\text{erf}(x)=\frac{2}{\sqrt{\pi}}\int_{0}^{x}e^{-t^{2}}dt$ (41)

and it is a monotonic function.

###### Proof.

$X$ follows a Gaussian distribution with mean $\mu$ and variance $\sigma$. Let
its cumulative distribution function be $\Phi_{\mu,\sigma}$. By definition,

$\displaystyle P(X\leq\hat{x})=\Phi_{\mu,\sigma}(\hat{x})$ (42) $\displaystyle
P(X\geq\hat{x})=1-\Phi_{\mu,\sigma}(\hat{x})$ (43) $\displaystyle
P(X-\mu\geq\epsilon)=1-\Phi_{\mu,\sigma}(\mu+\epsilon)$ (44)

Since two sided probability is higher than one sided, we have,

$\displaystyle P(X-\mu\geq\epsilon)$ $\displaystyle\leq
P(|X-\mu|\geq\epsilon)$ (45) $\displaystyle\leq
2\exp\Bigg{[}\left(\frac{4RC_{s}\sqrt{S_{K}||\mathscr{L}_{s}||}}{\epsilon}\right)^{\frac{2n}{h}}-\frac{m\epsilon^{2}}{4M^{2}}\Bigg{]}$
(46)

where we used Theorem 2. Using eq.44, we have,

$\displaystyle 1-\Phi_{\mu,\sigma}(\mu+\epsilon)\leq
2\exp\Bigg{[}\left(\frac{4RC_{s}\sqrt{S_{K}||\mathscr{L}_{s}||}}{\epsilon}\right)^{\frac{2n}{h}}-\frac{m\epsilon^{2}}{4M^{2}}\Bigg{]}$
(47)

For a Gaussian distribution, we can use the following expression for the
cumulative distribution function,

$\displaystyle\Phi_{\mu,\sigma}(\hat{x})=\frac{1}{2}\Big{[}1+\text{erf}\big{(}\frac{\hat{x}-\mu}{\sigma\sqrt{2}}\big{)}\Big{]}$
(48)

where erf is the Gauss error function. Using this in the eq.47,

$\displaystyle 1-\text{erf}\big{(}\frac{\epsilon}{\sigma\sqrt{2}}\big{)}\leq
4\exp\Bigg{[}\left(\frac{4RC_{s}\sqrt{S_{K}||\mathscr{L}_{s}||}}{\epsilon}\right)^{\frac{2n}{h}}-\frac{m\epsilon^{2}}{4M^{2}}\Bigg{]}$
(49)
$\displaystyle\text{erf}\big{(}\frac{\epsilon}{\sigma\sqrt{2}}\big{)}\geq-4\exp\Bigg{[}\left(\frac{4RC_{s}\sqrt{S_{K}||\mathscr{L}_{s}||}}{\epsilon}\right)^{\frac{2n}{h}}-\frac{m\epsilon^{2}}{4M^{2}}\Bigg{]}+1$
(50)

Since the function erf is invertible within domain (-1,1), we have,

$\displaystyle\frac{\epsilon}{\sigma\sqrt{2}}$
$\displaystyle\geq\text{erf}^{-1}\Bigg{[}-4\exp\Big{[}\left(\frac{4RC_{s}\sqrt{S_{K}||\mathscr{L}_{s}||}}{\epsilon}\right)^{\frac{2n}{h}}-\frac{m\epsilon^{2}}{4M^{2}}\Big{]}+1\Bigg{]}$
(51) $\displaystyle\sigma$
$\displaystyle\leq\frac{\epsilon}{\sqrt{2}\text{erf}^{-1}\Big{[}-4\exp\Big{[}\left(\frac{4RC_{s}\sqrt{S_{K}||\mathscr{L}_{s}||}}{\epsilon}\right)^{\frac{2n}{h}}-\frac{m\epsilon^{2}}{4M^{2}}\Big{]}+1\Big{]}}$
(52)

∎

###  C.2 Consistency of Estimates

######  Theorem 4.

Let $f^{*}$ and $f^{m}_{h}$ and $f^{*}_{h}$ be optimal discriminators defined
as

$\displaystyle f^{*}$
$\displaystyle=\underset{f}{\operatorname{argmax}}[{E_{p(x)}\log\sigma(f(x))+E_{q(x)}\log(1-\sigma(f(x)))}]$
(53) $\displaystyle f^{*}_{h}$ $\displaystyle=\underset{f\in
h}{\operatorname{argmax}}[{E_{p(x)}\log\sigma(f(x))+E_{q(x)}\log(1-\sigma(f(x)))}]$
(54) $\displaystyle f^{m}_{h}$ $\displaystyle=\underset{f\in
h}{\operatorname{argmax}}\Big{[}{\frac{1}{m}\sum_{x_{i}\sim
p(x_{i})}\log\sigma(f(x_{i}))+\frac{1}{m}\sum_{x_{j}\sim
q(x_{j})}\log(1-\sigma(f(x_{j})))}\Big{]}-\frac{\lambda_{0}}{m}||g||_{\mathcal{L}_{2}(d\tau)}^{2}$
(55)

and the KL estimate is given by $KL(f)=E_{p(x)}[f(x)],\hskip
11.38092ptKL_{m}(f)=\frac{1}{m}\sum_{x_{i}\sim p(x_{i})}[f(x)]$. Then, in the
limiting case as $m\to\infty$, $|KL_{m}(f^{m}_{h})-KL(f^{*})|\to 0$.

###### Proof.

Estimation error can be divided into three terms as

$\displaystyle
KL_{m}(f^{m}_{h})-KL(f^{*})=\underbrace{KL_{m}(f^{m}_{h})-KL(f^{m}_{h})}_{\text{Deviation-
from-mean
error}}+\underbrace{KL(f^{m}_{h})-KL(f^{*}_{h})}_{\text{Discriminator induced
error}}+\underbrace{KL(f^{*}_{h})-KL(f^{*})}_{Bias}$ (56)

Therefore,

$\displaystyle\small|KL_{m}(f^{m}_{h})-KL(f^{*})|$
$\displaystyle\leq|{KL_{m}(f^{m}_{h})-KL(f^{m}_{h})}|+|{KL(f^{m}_{h})-KL(f^{*}_{h})}|$
$\displaystyle+|{KL(f^{*}_{h})-KL(f^{*})}|$ (57)

To show that the total error goes to zero, we show that each term on the right
goes to zero. The last term is the bias and we assume that the RKHS space
$h=\mathcal{H}$ we consider consists the true solution, $f^{*}$. Hence the
bias goes to zero.

Using Theorem 2, it is immediately clear that the first term,
$|{KL_{m}(f^{m}_{h})-KL(f^{m}_{h})}|$ approaches zero in the limiting case as
$m\to\infty$.

The only remaining is the second term, $|{KL(f^{m}_{h})-KL(f^{*}_{h})}|$. In
Theorem 5 we show that this term also goes to zero as $m\to 0$. ∎

######  Theorem 5.

Let $f^{*}_{h}$ and $f^{*}_{h}$ be the optimal discriminators as defined in
eq. (54) and eq. (55), and the KL divergence estimate using discriminators
learned using finite and infinite samples be
$KL(f^{m}_{h})=\int[f^{m}_{h}(x)]p(x)dx$ and
$KL(f^{*}_{h})=\int[f^{*}_{h}(x)]p(x)dx$, where, Then, in the limiting case,
we have

$\underset{m\to\infty}{lim}|KL(f^{m}_{h})-KL(f^{*}_{h})|=0$

###### Proof.

$\displaystyle|KL(f^{m}_{h})-KL(f^{*}_{h})|$
$\displaystyle=|\int[f^{m}_{h}(x)-f^{*}_{h}(x)]p(x)dx|$
$\displaystyle\leq\underset{x}{sup}{|f^{m}_{h}(x)-f^{*}_{h}(x)|}=||f^{m}_{h}(x)-f^{*}_{h}(x)||_{\infty}$

Therefore, we can show
$\underset{m\to\infty}{lim}KL(f^{m}_{h})-KL(f^{*}_{h})=0$ if
$\underset{m\to\infty}{lim}||f^{m}_{h}(x)-f^{*}_{h}(x)||_{\infty}=0$, that is,
if the function $f^{m}_{h}(x)$ converges uniformly to function $f^{*}_{h}(x)$
in the limiting case.

The two maximizer functions are given by

$\displaystyle f^{*}_{h}$ $\displaystyle=\underset{f\in
h}{\text{argmax}}[{E_{p(x)}\log\sigma(f(x))+E_{q(x)}\log(1-\sigma(f(x)))}]$
(58) $\displaystyle f^{m}_{h}$ $\displaystyle=\underset{f\in
h}{\text{argmax}}\Big{[}{\frac{1}{m}\sum_{x_{i}\sim
p(x_{i})}\log\sigma(f(x_{i}))+\frac{1}{m}\sum_{x_{j}\sim
q(x_{j})}\log(1-\sigma(f(x_{j})))}\Big{]}-\frac{\lambda_{0}}{m}||g||^{2}$ (59)

As a first step in showing that $f^{m}_{h}$ uniformly approaches $f^{*}_{h}$,
we first show that
$\underset{m\to\infty}{lim}\frac{\lambda_{0}}{m}||g||^{2}=0$ in Lemma 4.

Then, to prove the rest, let us denote,

$\displaystyle G_{m}(f)$ $\displaystyle={\frac{1}{m}\sum_{x_{i}\sim
p(x_{i})}\log\sigma(f(x_{i}))+\frac{1}{m}\sum_{x_{j}\sim
q(x_{j})}\log(1-\sigma(f(x_{j})))}$ $\displaystyle G(f)$
$\displaystyle={E_{p(x)}\log\sigma(f(x))+E_{q(x)}\log(1-\sigma(f(x)))}$

In Lemma 5, we prove that functionals $G(f)$ and $G_{m}(f)$ are concave with
respect to function $f$. In the light of these two lemmas, we argue

$\displaystyle\underset{m\to\infty}{lim}||f^{m}_{h}(x)-f^{*}_{h}(x)||_{\infty}=0\hskip
5.69046pt\text{if}\hskip
5.69046pt\underset{m\to\infty}{lim}\underset{f}{\sup}|G_{m}(f)-G(f)|=0$ (60)

Next, we show $\underset{m\to\infty}{lim}\underset{f}{\sup}|G_{m}(f)-G(f)|=0$
as follows. We have,

$\displaystyle|G_{m}-G|=$ $\displaystyle|{\frac{1}{m}\sum_{x_{i}\sim
p(x_{i})}\log\sigma(f(x_{i}))+\frac{1}{m}\sum_{x_{j}\sim
q(x_{j})}\log(1-\sigma(f(x_{j})))}$
$\displaystyle-{E_{p(x)}\log\sigma(f(x))+E_{q(x)}\log(1-\sigma(f(x)))}|$ (61)
$\displaystyle\leq$ $\displaystyle|\frac{1}{m}\sum_{x_{i}\sim
p(x_{i})}\log\sigma(f(x_{i}))-E_{p(x)}\log\sigma(f(x))|$
$\displaystyle+|\frac{1}{m}\sum_{x_{j}\sim
q(x_{j})}\log(1-\sigma(f(x_{j})))-E_{q(x)}\log(1-\sigma(f(x)))|$ (62)
$\displaystyle\therefore\underset{m\to\infty}{lim}\underset{f}{\sup}$
$\displaystyle|G_{m}(f)-G(f)|\leq\underset{m\to\infty}{lim}\underset{f}{\sup}|\frac{1}{m}\sum_{x_{i}\sim
p(x_{i})}\log\sigma(f(x_{i}))-E_{p(x)}\log\sigma(f(x))|$ $\displaystyle\hskip
28.45274pt+\underset{m\to\infty}{lim}\underset{f}{\sup}|\frac{1}{m}\sum_{x_{j}\sim
q(x_{j})}\log(1-\sigma(f(x_{j})))-E_{q(x)}\log(1-\sigma(f(x)))|$ (63)

Both the terms on right hand side go to zero if $\log\circ\sigma\circ f$ is in
a Glivenko Cantelli class of functions using Empirical Process Theory [17],
which we prove in Lemma 6. That completes the proof. ∎

######  Lemma 4.

$\underset{m\to\infty}{lim}\frac{\lambda_{0}}{m}||g||^{2}=0$

###### Proof.

$||g||_{\mathcal{L}_{2}(d\rho)}$ is bounded because $g$ is Lipschitz
continuous and its domain is bounded. Since, $||g||_{\mathcal{L}_{2}(d\rho)}$
is bounded, we immediately obtain the required statement. ∎

######  Lemma 5.

The functional $G(f)$ is concave with respect to function $f$ in the following
sense: $\theta_{1}G(f_{1})+\theta_{2}G(f_{2})\leq
G(\theta_{1}f_{1}+\theta_{2}f_{2})$ for any $\theta_{1},\theta_{2}\in(0,1)$
such that $\theta_{1}+\theta_{2}=1$. The same is true for $G_{m}(f).$

###### Proof.

$\displaystyle\theta_{1}G(f_{1})+\theta_{2}G(f_{2})$
$\displaystyle=\theta_{1}\Big{[}\int p(x)\log\sigma(f_{1}(x))dx+\int
q(x)\log(1-\sigma(f_{1}(x)))dx\Big{]}$ $\displaystyle+\theta_{2}\Big{[}\int
p(x)\log\sigma(f_{2}(x))dx+\int q(x)\log(1-\sigma(f_{2}(x)))dx\Big{]}$ (64)
$\displaystyle=\int
p(x)\Big{[}\theta_{1}\log\sigma(f_{1}(x))dx+\theta_{2}\log\sigma(f_{2}(x))dx\Big{]}$
$\displaystyle+\int
q(x)\Big{[}\theta_{1}\log\sigma(-f_{1}(x))dx+\theta_{2}\log\sigma(-f_{2}(x))dx\Big{]}$
(65) $\displaystyle\leq\int
p(x)\log\sigma[\theta_{1}f_{1}(x)+\theta_{2}f_{2}(x)]dx$ $\displaystyle+\int
q(x)\log\sigma[-(\theta_{1}f_{1}(x)+\theta_{2}f_{2}(x))]dx$ (66)
$\displaystyle=G(\theta_{1}f_{1}+\theta_{2}f_{2})$ (67)

where we used the fact that $\log(1-\sigma(f(x)))=\log\sigma(-f(x))$ (this is
straightforward using definition of Sigmoid function, $\sigma$) in line 65. In
line C.2, we used the fact that $\log\sigma$ is a concave function (see Lemma
8). ∎

######  Lemma 6.

$\log\circ\sigma\circ f$ is a Glivenko Cantelli class of function.

###### Proof.

In Lemma 7, we show that, by definition, $f$ is Lipschitz continuous with some
Lipschitz constant $L_{f}$. In Lemma 8 we show that if $f$ is a Lipschitz
continuous function from $\mathcal{X}$ to $(-\infty,\infty)$ with Lipschitz
constant, $L_{f}$, then $\log\sigma f$ is a function from $\mathcal{X}$ to
$(-\infty,0)$ with same Lipschitz constant $L_{f}$. Hence, $v=\log\sigma f$ is
a a function from $\mathcal{X}$ to $(-r,0)$. Note that since $\mathcal{X}$ is
bounded and $f$ is Lipschitz continuous from $\mathcal{X}$ to $\mathbb{R}$, we
can always find some $r$ such that $v$ maps from $\mathcal{X}$ to $(-r,0)$.

Now, we show that $v=\log\sigma f$ is Glivenko Cantelli by entropy number. Let
$\mathcal{V}=\\{v:v=\log(\sigma(f)),f\in\mathcal{F}\\}$. In Lemma 10, we use
theorem from [17] to show that $\mathcal{V}$ is Glivenko Cantelli if and only
if

$\displaystyle\frac{1}{m}\log
N(\epsilon,\mathcal{V}_{M},\ell_{1}(\mathbb{P}_{m}))\overset{\mathbb{P}}{\to}0,$
(68)

for any $M>0,\epsilon$, where $\mathcal{V}_{M}$ is the class of functions
$v\textbf{1}\\{E\leq M\\}$ where $v$ ranges over $\mathcal{V}$ and $E$ is an
envelope function to $\mathcal{V}$. Since we proved that
$\log(\sigma(f)(x))<0$ for any $x$, we can choose $E=v_{0}(x)=\textbf{0}$ as a
constant function that is an envelope to $\mathcal{V}$. For any $M>0$,
therefore, $1\\{E\leq M\\}=1$ trivially and $\mathcal{V}_{M}=\mathcal{V}$.
Hence, we just need to show

$\displaystyle\frac{1}{m}\log
N(\epsilon,\mathcal{V},\ell_{1}(\mathbb{P}_{m}))\overset{\mathbb{P}}{\to}0$
(69)

In Lemma 9, we show that the entropy number of such a function is given by

$\displaystyle\log\mathcal{N}(\epsilon,\mathcal{V},\ell_{1}(\mathbb{P}_{m}))\leq\left(\frac{16L.diam(\mathcal{X})}{\epsilon}\right)^{ddim(\mathcal{X})}\log\left(\frac{4r}{\epsilon}\right)$
(70)

and therefore is bounded and independent of the sample size $m$. Hence,
$\frac{1}{m}\log N(\epsilon,\mathcal{V},\ell_{1}(\mathbb{P}_{m}))$ goes to
$0$. ∎

######  Lemma 7.

The function $f$ defined in Theorem 1 on the main paper as:

$\displaystyle f(x)=\int_{\mathcal{W}}g(w)\psi(x,w)d\tau(w),$ (71)

where $\psi(x,w)=\phi_{\theta}(x)^{T}w$ and the function $\phi_{\theta}$ is
Lipschitz continuous with Lipschitz constant $L_{\phi}$. Then, the function
$f$ is Lipschitz continuous with some Lipschitz constant, $L_{f}$.

###### Proof.

By the definition,

$\displaystyle f(x)=\langle g(w),\psi(x,w)\rangle_{\mathcal{L}_{2}(d\tau)}$
(72)

For any two points $x_{1}$ and $x_{2}$,

$\displaystyle|f(x_{1})-f(x_{2})|$ $\displaystyle=\langle
g(w),\psi(x_{1},w)-\psi(x_{2},w)\rangle_{\mathcal{L}_{2}(d\tau)}$ (73)
$\displaystyle\leq||g(w)||_{\mathcal{L}_{2}(d\tau)}||\psi(x_{1},w)-\psi(x_{2},w)||_{\mathcal{L}_{2}(d\tau)}$
(74)

where we used Cauchy Schwartz. Now, taking the difference in $\psi$, it can be
written as

$\displaystyle||\psi(x_{1},w)-\psi(x_{2},w)||_{\mathcal{L}_{2}(d\tau)}$
$\displaystyle=\sqrt{\int[\psi(x_{1},w)-\psi(x_{2},w)]^{2}d\tau(w)}$ (75)
$\displaystyle=\sqrt{\int[(\phi_{\theta}(x_{1})-\phi_{\theta}(x_{2}))^{T}w]^{2}d\tau(w)}$
(76)
$\displaystyle\leq\sqrt{\int||\phi_{\theta}(x_{1})-\phi_{\theta}(x_{2})||^{2}||w||^{2}d\tau(w)}$
(77)

where we again used Cauchy Schwartz in the last line since
$[(\phi_{\theta}(x_{1})-\phi_{\theta}(x_{2}))^{T}w]$ is an inner product in
$\mathbb{R}^{D}$ where D is the dimension of $w$. Since $\phi_{\theta}$ is
Lipschitz continuous with Lipschitz constant $L_{\phi}$, we have

$||\phi_{\theta}(x_{1})-\phi_{\theta}(x_{2})||\leq L_{\phi}||x_{1}-x_{2}||$

Using this inequality in eq.77, we obtain

$\displaystyle||\psi(x_{1},w)-\psi(x_{2},w)||_{\mathcal{L}_{2}(d\tau)}$
$\displaystyle\leq L_{\phi}||x_{1}-x_{2}||\sqrt{\int||w||^{2}d\tau(w)}$ (78)
$\displaystyle=L_{\phi}||x_{1}-x_{2}||\sqrt{tr(C_{w})}$ (79)

where, $C_{w}$ is the uncentered covariance matrix of Gaussian distributed
$w$. Plugging eq.(79) in eq.(74), we obtain

$\displaystyle|f(x_{1})-f(x_{2})|\leq||g(w)||_{\mathcal{L}_{2}(d\tau)}L_{\phi}\sqrt{tr(C_{w})}||x_{1}-x_{2}||$
(80)

Since, we have that $||g(w)||_{\mathcal{L}_{2}(d\tau)}<\infty$ (see Lemma 4),
we have proved that $f$ is Lipschitz continuous with Lipschitz constant given
by $L_{f}\leq||g(w)||_{\mathcal{L}_{2}(d\tau)}L_{\phi}\sqrt{tr(C_{w})}$. ∎

######  Lemma 8.

The function $\log\circ\sigma$ exhibits following properties:
i) It is a concave function with its derivative always between $0$ and $1$
ii) If the Lipschitz constant of $f$ is $L_{f}$, so is the Lipschitz constant
of $\log\circ\sigma\circ f$

###### Proof.

i) Let us denote $u(x)=\log(\sigma(x))$. Then, we have,

$\displaystyle u(x)=\log\frac{e^{x}}{1+e^{x}}=x-\log(1+e^{x})$ (81)
$\displaystyle\therefore$ $\displaystyle
u^{\prime}(x)=1-\frac{e^{x}}{1+e^{x}}=\frac{1}{1+e^{x}}$ (82)
$\displaystyle\therefore$ $\displaystyle 0<u^{\prime}(x)<1,\hskip
14.22636pt\forall x\in(-\infty,\infty)$ (83)

which proves that the derivative is between $0$ and $1$. To show that $u(x)$
is concave, it is sufficient to note that its second derivative is always
negative.

ii) Let us use notation $u=\log(\sigma)$, and let $f_{2}=f(x_{2})$,
$f_{1}=f(x_{1})$, $u_{2}=u(f(x_{2}))$, $u_{1}=u(f(x_{1}))$. Since the maximum
derivative of $u$ is upper bounded by 1, $u$ as a function of $f$ has
Lipschitz constant $1$ and therefore, we can write

$\displaystyle u_{2}-u_{1}=u(f_{2})-u(f_{1})\leq f_{2}-f_{1}$
$\displaystyle=f(x_{2})-f(x_{1})$ (84) $\displaystyle\leq
L_{f}||x_{2}-x_{1}||$ (85)

where the last inequality is because $f$ is Lipschitz continuous with
Lipschitz constant $L_{f}$. This proves that the Lipschitz constant of
$\log\circ\sigma\circ f$ is also $L_{f}$.

∎

######  Lemma 9.

Let $\mathcal{F}_{L}$ be the space of L-Lipschitz functions mapping the metric
space ($\mathcal{X},\rho$) to [0,r]. Let $ddim(\mathcal{X})$ and
$diam(\mathcal{X})$ denote the doubling dimension and diameter of
$\mathcal{X}$ respectively. Then,
i) the covering numbers of $\mathcal{F}_{L}$ can be estimated in terms of the
covering numbers of $\mathcal{X}$:

$\displaystyle\mathcal{N}(\epsilon,\mathcal{F}_{L},||.||_{\infty})\leq\left(\frac{4r}{\epsilon}\right)^{\mathcal{N}(\epsilon/{8L},\mathcal{X},||.||_{\infty})}$
(86)

ii) the entropy number of $\mathcal{F}_{L}$ can be estimated as:

$\displaystyle\log\mathcal{N}(\epsilon,\mathcal{F}_{L},||.||_{\infty})\leq\left(\frac{16L.diam(\mathcal{X})}{\epsilon}\right)^{ddim(\mathcal{X})}\log\left(\frac{4r}{\epsilon}\right)$
(87)

iii) the entropy number with respect to
$\ell_{1}(\mathbb{P}_{m})=\int|f|d\mathbb{P}_{m}=\frac{1}{m}\sum_{k}|f(x_{k})|$
defined with respect to the $m$ input points, is the same as (ii), i.e.

$\displaystyle\log\mathcal{N}(\epsilon,\mathcal{F}_{L},\ell_{1}(\mathbb{P}_{m}))\leq\left(\frac{16L.diam(\mathcal{X})}{\epsilon}\right)^{ddim(\mathcal{X})}\log\left(\frac{4r}{\epsilon}\right)$
(88)

where $\mathbb{P}_{m}$ is an empirical probability measure with respect to $m$
inputs points in $\mathcal{X}$.

###### Proof.

The proof is adapted from [26] Lemma 2 and [27] Lemma 6, and modified to
handle range $[0,r]$.

i) We first cover the domain $\mathcal{X}$ by $N$ balls
${U_{1},U_{2},...,U_{|N|}}$, where
$N=\mathcal{N}(\epsilon/{8L},\mathcal{X},||.||_{\infty})$ is the covering
number of $\mathcal{X}$, $N=\\{x_{i}\in U_{i}\\}^{|N|}_{i=1}$ is a set of
center points of $|N|$ balls and $\epsilon^{\prime}=\epsilon/{8L}$ is the
radius of the covering balls.

Now, our strategy is to construct an $\epsilon$ cover
$\hat{F}=\\{\hat{f_{1}},...,\hat{f}_{|\hat{F}|}\\}$ for $\mathcal{F}_{L}$ with
respect to $||.||_{\infty}$. To do so, at every point $x_{i}\in N$, we choose
the value of $\hat{f}(x_{i})$ to be some multiple of
$2L\epsilon^{\prime}=\frac{\epsilon}{4}$, while maintaining
$||\hat{f}||_{Lip}\leq 2L$. We then construct a 2L-Lipschitz extension for
$\hat{f}$ from $N$ to all over $\mathcal{X}$(note that such an extension
always exists, see [28, 29]).

With this construction, we can show that every $f\in\mathcal{F}_{L}$ is close
to some $\hat{f}\in\hat{F}$ in the sense that
$||f-\hat{f}||_{\infty}\leq\epsilon$. To show this, note the following:

$\displaystyle|f(x)-\hat{f}(x)|$
$\displaystyle\leq|f(x)-f(x_{N})|+|f(x_{N})-\hat{f}(x_{N})|+|\hat{f}(x_{N})-\hat{f}(x)|$
(89) $\displaystyle\leq L.\rho(x,x_{N})+\epsilon/4+2L.\rho(x,x_{N})$ (90)
$\displaystyle\leq\epsilon$ (91)

where the inequality in eq.90 is due to the fact that $f$ is $L$-Lipschitz and
$\hat{f}$ is $2L$-Lipschitz and since we have covered the input space
$\mathcal{X}$, each $x$ is within $\epsilon^{\prime}$ of some $x_{N}$. Also
note that for every $f(x_{N})$ we can find $\hat{f}(x_{N})$ within some radius
$\epsilon/4$; this is because we choose $f(x_{N})$ to be some multiple of
$2L\epsilon^{\prime}$. Finally, we need to compute the cardinality of
$\hat{F}$, i.e. $|\hat{F}|$. For any $x_{i}\in|N|$, $\hat{f}$ can take one of
the multiple of $2L\epsilon^{\prime}$ values. Hence, there are
$r/2L\epsilon^{\prime}$ such possibilities as the range is $[0,r]$. Since
there are $|N|$ such possibilities for $x_{i}$, the upper bound on all
possible function values $\hat{f}$ is
$(\frac{r}{2L\epsilon^{\prime}})^{|N|}=(\frac{4r}{\epsilon})^{|N|}$, which
proves the first statement after plugging in the value of $|N|$.

ii) Taking logarithm of the result in i)

$\displaystyle\log\mathcal{N}(\epsilon,\mathcal{F}_{L},||.||_{\infty})\leq{\mathcal{N}(\epsilon/{8L},\mathcal{X},||.||_{\infty})}\log\left(\frac{4r}{\epsilon}\right)$
(92)

The covering number of the input space, $\mathcal{X}$ in terms of doubling
dimension, $ddim(\mathcal{X})$ and diameter, $diam(\mathcal{X})$ can be
written as [30]:

$\displaystyle{\mathcal{N}(\epsilon,\mathcal{X},||.||_{\infty})}\leq\left(\frac{2diam(\mathcal{X})}{\epsilon}\right)^{ddim(\mathcal{X})}$
(93)

Plugging this expression in eq.(92), we obtain the required expression.

iii) The result in i) is with respect to $||.||_{\infty}$. In eq.(90), we
showed that for any $f\in\mathcal{F}_{L}$ there is some $\hat{f}\in\hat{F}$
within a radius of $\epsilon$ such that $||f-\hat{f}||_{\infty}\leq\epsilon$.
Here, we show that this also implies that
$||f-\hat{f}||_{\ell_{1}(\mathbb{P}_{m})}\leq\epsilon$. We show this as
follows:

$\displaystyle||f-\hat{f}||_{\ell_{1}(\mathbb{P}_{m})}$
$\displaystyle=\frac{1}{m}\sum_{k=1}^{m}|f(x_{k})-\hat{f}(x_{k})|$ (94)
$\displaystyle\leq\frac{1}{m}\sum_{k=1}^{m}\epsilon=\epsilon$ (95)

Therefore, the entropy number with respect to ${\ell_{1}(\mathbb{P}_{m})}$
metric is same as the entropy number with respect to the $||.||_{\infty}$,
which proves our third claim. ∎

######  Lemma 10 ([17] Theorem 3.5. ).

Let $\mathcal{V}$ be a class of measurable functions with envelope $E$ such
that $P(E)<\infty$. Let $\mathcal{V}_{M}$ be the class of functions
$v.\textbf{1}\\{E\leq M\\}$ where $v$ ranges over $\mathcal{V}$. Then,
$\mathcal{V}$ is a Glivenco Cantelli class of functions, i.e. it satisfies

$\displaystyle\underset{v\in\mathcal{V}}{sup}\hskip
2.84544pt|\mathbb{P}_{m}v-Pv|$ (96)

, if and only if

$\displaystyle\frac{1}{m}\log
N(\epsilon,\mathcal{V}_{M},L_{1}(\mathbb{P}_{m}))\overset{\mathbb{P}}{\to}0,$
(97)

for every $\epsilon>0$ and $M>0$, where $Pv=\int vdP$ and
$\mathbb{P}_{m}v=\frac{1}{m}\sum_{k}v(x_{k})$.

##  Appendix D Experimental Results

##### Code:

The code will be publicly released.

###  D.1 Two Gaussian

####  D.1.1 Architecture and Implementation

RKHS Discriminator Architecture (Pytorch Code)

⬇

class RKHS_Net(nn.Module):

def __init__(self, dim =10, mid_dim1=20, mid_dim2=20, mid_dim3=20, D=50, gamma
=1, metric = ’rbf’, lip=5, g_lip =5):

super(RKHS_Net, self).__init__()

self.gamma = torch.FloatTensor([gamma])

self.metric = metric

self.D = D

self.act = nn.ReLU()

self.lin1 = spectral_norm( nn.Linear(dim, mid_dim1), k =g_lip)

self.lin2 = spectral_norm( nn.Linear(mid_dim1, mid_dim2), k =g_lip)

self.lin3 = spectral_norm( nn.Linear(mid_dim2 , mid_dim3), k =g_lip)

self.lin4 = spectral_norm( nn.Linear(mid_dim3, 1), k =g_lip)

self.g = nn.Sequential(self.lin1,

self.act,

self.lin2,

self.act,

self.lin3,

self.act,

self.lin4

)

self.lin_phi1 = spectral_norm(nn.Linear(2, mid_dim1), k=lip)

self.lin_phi2 = spectral_norm(nn.Linear(mid_dim1, mid_dim2), k=lip)

self.lin_phi3 = spectral_norm(nn.Linear(mid_dim2, mid_dim3), k=lip)

self.lin_phi4 = spectral_norm(nn.Linear(mid_dim3, dim), k=lip)

self.phi = nn.Sequential(self.lin_phi1,

self.act,

self.lin_phi2,

self.act,

self.lin_phi3,

self.act,

self.lin_phi4

)

def forward(self, y):

x=self.phi(y)

d = x.shape[1]

if self.metric ==’rbf’:

w= torch.sqrt(2*self.gamma)*torch.randn(size=(self.D,d))

w=w.to(x.device)

psi = ((torch.matmul(x,w.permute(1,0))
))*(torch.sqrt(2/torch.FloatTensor([self.D])).to(x.device))

w_a = w

g= self.g(w_a)

f = (psi*g.permute(1,0)).mean(1)

g_norm =(g**2).mean()

return f, g_norm

Simple Neural Network Discriminator Architecture ( Pytorch Code)

⬇

class DNet_basic(nn.Module):

def __init__(self, input_dim, mid_dim1, mid_dim2, output_dim, lip_constraint =
False, lip = 5):

super(DNet_basic, self).__init__()

self.act = nn.ReLU()

self.lin1 = nn.Linear(input_dim, mid_dim1)

self.lin2 = nn.Linear(mid_dim1, mid_dim2)

self.lin3 = nn.Linear(mid_dim2, mid_dim2)

self.lin4 = nn.Linear(mid_dim2, output_dim)

# self.sigmoid=nn.Sigmoid()

self.phi = nn.Sequential(self.lin1,

self.act,

self.lin2,

self.act,

self.lin3,

self.act,

self.lin4

)

def forward(self, x):

t = self.phi(x)

return t

##### Discrete approximation:

Both the discriminators have stacked Fully connected layers and activation
function. In the proposed RKHS discriminator, we have an additional network
self.g which we use to approximate the continuous integral
$f(x)=\int_{\mathcal{W}}g(w)\psi(x,w)d\tau(w)$ with the following discrete
approximation:

$\displaystyle
f(x)=\frac{1}{D}\sum_{k=1}^{D}g(w_{k})(\phi^{T}w_{k}\sqrt{\frac{2}{D}})$ (98)

where $w$ is sampled from a Normal distribution with variance $\gamma$. In our
experiments $D=500$ was sufficient. Note that the Neural network discriminator
is similar to $\phi^{T}w$, except that $w$ is not randomly sampled and there
is no $g$ network.

##### Lipschitz constraints:

: To enforce Lipschitz constraints on network $g$ and $\phi$ consistent with
our assumptions and theoretical results, we use spectral normalization in the
RKHS discriminator while it is absent in the basic Neural network
discriminator.

####  D.1.2 Data and Hyperparameters

Data: Since this is a toy experiment, data were generated locally using
pytorch command randn to sample from Gaussian distribution.

Learning rate: $5\times 10^{-3}$ (both models)
No. of samples from each distribution: $2500$ (both models)
Minibatch size: 50 (both models)
$\boldsymbol{\lambda:}0.005$ (RKHS disc.)
Hyperparameter selection: (RKHS disc.) The hyperparameters like learning rate
and $\lambda$ were selected by first estimating KL divergence at a mid value
like $13$. Then, same value was used in all experiments.

####  D.1.3 Computational Resources and Time

Running one experiment of KL divergence calculation takes 74 s for the basic
algorithm while it takes 245 s for the proposed method in a single GeForce GTX
1080 Ti GPU with 11GB memory.

###  D.2 Mutual Information Estimation

####  D.2.1 Models, Architecture and Implementation

RKHS Discriminator Architecture (Pytorch Code)

⬇

class ConcatLipFeatures(nn.Module):

def __init__(self, dim, hidden_dim, layers, activation,lip, gamma =1, metric =
’rbf’, D=500, mid_dim=5, g_lip =2, **extra_kwargs):

super(ConcatLipFeatures, self).__init__()

self.gamma = torch.FloatTensor([gamma])

self.metric = metric

self.D = D

self.act = nn.ReLU()

self.lin1 = spectral_norm( nn.Linear(hidden_dim, mid_dim), k = g_lip)

self.lin2 = spectral_norm( nn.Linear(mid_dim, mid_dim), k = g_lip)

self.lin3 = spectral_norm( nn.Linear(mid_dim, mid_dim), k = g_lip)

self.lin4 = spectral_norm( nn.Linear(mid_dim, 1), k = g_lip)

self.g = nn.Sequential(self.lin1,

self.act,

self.lin2,

self.act,

self.lin3,

self.act,

self.lin4

)

# output of this layer is d dim features

self.rkhs_layer = feature_perceptron(dim * 2, hidden_dim, 1, layers,
activation, lip)

def forward(self, x, y):

batch_size = x.size(0)

# Tile all possible combinations of x and y

x_tiled = torch.stack([x] * batch_size, dim=0)

y_tiled = torch.stack([y] * batch_size, dim=1)

# xy is [batch_size * batch_size, x_dim + y_dim]

xy_pairs = torch.reshape(torch.cat((x_tiled, y_tiled), dim=2), [

batch_size * batch_size, -1])

# Compute features for each x_i, y_j pair.

phi = self.rkhs_layer(xy_pairs)

d = phi.shape[1]

if self.metric == ’rbf’:

w = torch.sqrt(2 * self.gamma) * torch.randn(size=(self.D, d))

w = w.to(x.device)

psi = ((torch.matmul(phi, w.permute(1, 0)))) * (torch.sqrt(2 /
torch.FloatTensor([self.D])).to(x.device))

w_a = w # torch.cat((w,u.permute(1,0)),1)

g = self.g(w_a)

f = (psi * g.permute(1, 0)).mean(1)

g_norm = (g ** 2).mean()

return f, g_norm

Simple Neural Network Discriminator Architecture (Pytorch Code)

⬇

class ConcatCritic(nn.Module):

def __init__(self, dim, hidden_dim, layers, activation, **extra_kwargs):

super(ConcatCritic, self).__init__()

# output is scalar score

self._f = mlp(dim * 2, hidden_dim, 1, layers, activation)

def forward(self, x, y):

batch_size = x.size(0)

# Tile all possible combinations of x and y

x_tiled = torch.stack([x] * batch_size, dim=0)

y_tiled = torch.stack([y] * batch_size, dim=1)

# xy is [batch_size * batch_size, x_dim + y_dim]

xy_pairs = torch.reshape(torch.cat((x_tiled, y_tiled), dim=2), [

batch_size * batch_size, -1])

# Compute scores for each x_i, y_j pair.

scores = self._f(xy_pairs)

return torch.reshape(scores, [batch_size, batch_size]).t()

Similar to the previous experiment, the RKHS discriminator and the Neural
network discriminator are similar in core design. The main difference lies in
that the RKHS discriminator has this inner product construction same as eq.(1)
in previous subsection. To achieve this construction, the RKHS discriminator
an additional network, self.g and enforces Lipschitz constraint through
spectral normalization, which are absent in simple Neural network
discriminator.

####  D.2.2 Data and Hyperparameters

Data: The experimental setup and data generation follow
https://github.com/ermongroup/smile-mi-estimator.
Common for all methods
batch size: 64
no. of layers: 2
hidden dim: 256
no. of iterations: 40000
learning rate: $5\times 10^{-4}$

Specific to the proposed method
$\gamma:$5
Lipschitz constant enforced, $L_{\phi}$ (layer wise): 5
Lipschitz constant enforced, $L_{g}$ (layer wise): 5

####  D.2.3 Computational Resources and Time

GPU: GeForce RTX 2080 Ti 11 GB

Below, we report time taken by each method to complete an experiment to obtain
mutual information between two 20-d Gaussian distributed random variables
using $40,000$ samples from each distribution and mutual information
increasing stepwise.

Table 1: Time taken to complete one experiment CPC | NWJ | SMILE | Ours (RKHS disc.)
---|---|---|---
52 s | 48 s | 52 s | 63 s

####  D.2.4 Existing Assets

We used the code from the repo https://github.com/ermongroup/smile-mi-
estimator to generate data as well as run baseline mutual information methods.
This code corresponds to the Song et al. [9].

###  D.3 Adversarial Variational Bayes

####  D.3.1 Models, Architecture and Implementation

RKHS Discriminator Architecture (Pytorch Code)

⬇

class Discriminator_RKHS(nn.Module):

def __init__(self, x_dim, h_dim, z_dim, lip = 5, g_lip = 5, dim = 10, mid_dim1
= 20, mid_dim2 = 20, mid_dim3 = 20, D=100, gamma =1, metric = ’rbf’):

super(Discriminator_RKHS, self).__init__()

self.metric = metric

self.gamma = torch.FloatTensor([gamma])

self.D = D

self.phi = nn.Sequential(

spectral_norm(nn.Linear(x_dim + z_dim, h_dim), k = lip),

nn.LeakyReLU(),

spectral_norm(nn.Linear(h_dim, h_dim), k = lip),

nn.LeakyReLU(),

spectral_norm(nn.Linear(h_dim, h_dim), k = lip),

nn.LeakyReLU(),

spectral_norm(nn.Linear(h_dim, h_dim), k = lip),

nn.LeakyReLU(),

spectral_norm(nn.Linear(h_dim, int(h_dim/4)), k = lip)

)

self.act = nn.ReLU()

self.lin1 = spectral_norm(nn.Linear(int(h_dim/4), mid_dim1), k=g_lip)

self.lin2 = spectral_norm(nn.Linear(mid_dim1, mid_dim2), k=g_lip)

self.lin3 = spectral_norm(nn.Linear(mid_dim2, mid_dim3), k=g_lip)

self.lin4 = spectral_norm(nn.Linear(mid_dim3, 1), k=g_lip)

self.g = nn.Sequential(self.lin1,

self.act,

self.lin2,

self.act,

self.lin3,

self.act,

self.lin4

)

def weight_init(self, mean, std):

for m in self._modules:

normal_init(self._modules[m], mean, std)

def forward(self, y, z):

y = y.view(y.shape[0], -1)

y = torch.cat([y, z], 1)

x =self.phi(y)

d = x.shape[1]

if self.metric == ’rbf’:

w = torch.sqrt(2 * self.gamma) * torch.randn(size=(self.D, d))

w = w.to(x.device)

psi = ((torch.matmul(x, w.permute(1, 0)))) * (torch.sqrt(2 /
torch.FloatTensor([self.D])).to(x.device))

w_a = w

g = self.g(w_a)

f = (psi * g.permute(1, 0)).mean(1)

g_norm = (g ** 2).mean()

return f, g_norm

Simple Neural Network Discriminator Architecture (Pytorch Code)

⬇

class Discriminator_simple(nn.Module):

def __init__(self, x_dim, h_dim, z_dim):

super(Discriminator_simple, self).__init__()

self.net = nn.Sequential(

nn.Linear(x_dim + z_dim, h_dim),

nn.LeakyReLU(),

nn.Linear(h_dim, h_dim),

nn.LeakyReLU(),

nn.Linear(h_dim, h_dim),

nn.LeakyReLU(),

nn.Linear(h_dim, h_dim),

nn.LeakyReLU(),

nn.Linear(h_dim, int(h_dim/4)),

nn.LeakyReLU(),

nn.Linear(int(h_dim/4), 1)

)

def weight_init(self, mean, std):

for m in self._modules:

normal_init(self._modules[m], mean, std)

def forward(self, x, z):

x = x.view(x.shape[0], -1)

x = torch.cat([x, z], 1)

out =self.net(x)

# x = x + torch.sum(z ** 2, 1)

return out

####  D.3.2 Data and Hyperparameters

Data: Standard MNIST dataset is used.

Learning rate: $10^{-3}$ (both models)
Minibatch size: 1024 (both models)
Hidden dim of encoder/decoder: 800 (both)
Hidden dim discriminator: 1024 (both)
$\boldsymbol{\lambda:}1$ (RKHS disc.)

####  D.3.3 Computational Resources and Time

GPU: GeForce GTX 1080 Ti 11GB
Time taken to train MNIST for 1000 epochs using AVB with simple Neural net
discriminator: 11.3 hrs
Time taken to train MNIST for 1000 epochs using AVB with RKHS discriminator:
14.7 hrs

####  D.3.4 Existing Assets

We followed the official implementation of Adversarial Variational Bayes [8]
at https://github.com/LMescheder/AdversarialVariationalBayes

##  Appendix E Societal Impacts

We discuss possible negative impacts in two categories: 1) Impact of
theoretical contribution, 2) Impact of applications

##### Societal Impact of theoretical contribution:

The main theoretical contribution of the paper is its connection between
reliable/stable estimation and complexity analysis of the discriminator
function space. In its general form, this contribution does not, by itself,
pose any negative societal impact. Rather, it is about stabilizing algorithms.
So, it contributes towards more robust and stable algorithms, and may help in
developing more secure applications. We do not foresee any negative societal
impacts in safety and security of human beings and automatic systems, human
rights, human livelihood or economic security, environment. We do not see it
causing theft, harassment, fraud, bias or discrimination.

##### Societal Impact of possible applications:

As demonstrated in the experiment section, this work can be applied to
information theoretic applications that require mutual information or KL
divergence estimation. For example, it has been used in generative modeling
like variational autoencoder, variational Bayes or in stabilizing generative
adversarial networks (GANs). These generative modeling techniques are, by
themselves, quite general and can have numerous applications, including the
ones with negative impacts. By helping in accurate estimation of KL divergence
and by providing theoretical analysis, this work is contributing to develop
stronger generative models and by extension could be indirectly helping in
their negative uses. In that aspect, we appeal everyone using the algorithms
and ideas in this paper to be thoughtful and responsible in their use.