# Reliable Estimation of Kullback–Leibler Divergence by Controlling
Discriminator Complexity in the Reproducing Kernel Hilbert Space

Sandesh Ghimire, Prashnna K Gyawali, Linwei Wang
College of Computing and Information Sciences
Rochester Institute of Technology
Rochester, NY 14623, USA
{sg9872, pkg2182<EMAIL_ADDRESS>

###### Abstract

Several scalable sample-based methods to compute the Kullback–Leibler (KL)
divergence between two distributions have been proposed and applied in large-
scale machine learning models. While they have been found to be unstable, the
theoretical root cause of the problem is not clear. In this paper, we study a
generative adversarial network based approach that uses a neural network
discriminator to estimate KL divergence. We argue that, in such case, high
fluctuations in the estimates are a consequence of not controlling the
complexity of the discriminator function space. We provide a theoretical
underpinning and remedy for this problem by first constructing a discriminator
in the Reproducing Kernel Hilbert Space (RKHS). This enables us to leverage
sample complexity and mean embedding to theoretically relate the error
probability bound of the KL estimates to the complexity of the discriminator
in RKHS. Based on this theory, we then present a scalable way to control the
complexity of the discriminator for a reliable estimation of KL divergence. We
support both our proposed theory and method to control the complexity of the
RKHS discriminator through controlled experiments.

##  1 Introduction

Calculating Kullback–Leibler (KL) divergence from data samples is an essential
component in many machine learning problems that involve Bayesian inference or
the calculation of mutual information. In small data regime, this problem has
been studied using variational technique and convex optimization [16]. In the
presence of ever-increasing data, several neural network models have been
proposed which require estimation of KL divergence such as total correlation
variational autoencoder (TC-VAE) [8], adversarial variational Bayes (AVB)
[14], information maximizing GAN (InfoGAN) [9], and amortized MAP [18]. These
large scale models have imposed the following new requirements on estimating
KL divergence: 1\. Scalability: The estimation algorithm should be able to
compute KL divergence from a large amount of data samples. 2\. Minibatch
compatibility: The algorithm should be compatible with minibatch-based
optimization and allow backpropagation (or other ways of optimizing the rest
of the network) based on the estimated value of KL divergence.

These needs make classic methods such as [16] impractical, but were met by
modern neural network based methods such as variational divergence
minimization (VDM) [17], mutual information neural estimation (MINE) [4], and
GAN-based KL estimation [15, 18]. A key attribute of these methods is that
they are based on updating a neural-net discriminator function to estimate KL
divergence from a subset of samples, which makes them scalable and minibatch
compatible. We, however, noted that even in simple toy examples, these methods
tended to be either unreliable (high fluctuation of estimates, as in GAN based
approach by [15]), or unstable (discriminator yields infinity, as in MINE and
VDM) (see Table.1). This behavior exacerbated when increasing the size of the
discriminator. Similar observations of instability of VDM and MINE have been
reported in the literature [15, 19].

In this paper, we attempt to provide a theoretical underpinning for the core
problem of the large fluctuation in the GAN based estimation of KL divergence.
We approach this problem from the perspective of sample complexity, and
propose that these fluctuations are a consequence of not controlling the
complexity of the discriminator function. This direction has not been explored
in existing works, and it faces the open question of how to properly measure
the complexity of the large function space represented by neural networks.
Note that naive approaches to bound complexity by the number of parameters
would neither be guaranteed to yield tight bound, nor be easy to implement
because it requires dynamically changing the size of the network during
optimization.

We introduce the following contributions to resolve this challenge. First, to
be able to compute the complexity of the discriminator function space, we
propose a novel construction of the discriminator such that it lies in a
smooth function space, the Reproducing Kernel Hilbert Space(RKHS). Leveraging
sample complexity analysis and mean embedding of RKHS, we then bound the
probability of the error of KL-divergence estimates in terms of the complexity
of RKHS space. This further allows us to theoretically substantiate our main
proposition that not controlling the complexity of the discriminator may lead
to high fluctuation in estimation. Finally, we propose a scalable way to
control the complexity of the discriminator based on the obtained error
probability bound. In controlled experiments, we demonstrate that failing to
control the complexity of the discriminator function leads to fluctuation in
KL divergence estimates, and that the proposed method decreases such
fluctuations.

##  2 Related Works

Nguyen et al [16] used variational function to estimate KL divergence from
samples of two distribution using convex risk minimization (CRM). They used
the RKHS norm of the variational function 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
that could not be scaled. VDM reformulates the f-Divergence objective using
Fenchel duality and uses a neural network to represent the variational
function [17]. It is in concept close to [16], while the use of neural network
and adversarial optimization made the estimation scalable. It however did not
control the complexity of the neural-net function, resulting in unstable
estimations.

One area of modern application of KL-divergence estimation is in computing
mutual information which, as shown in MINE [4], is useful in applications such
as stabilizing GANs or realizing the information bottleneck principle. MINE
also optimizes a lower bound, but tighter, to KL divergence (Donsker-Varadhan
representation). Similar to VDM, MINE uses a neural network as the dual
variational function: it is thus scalable, but without complexity control and
unstable.

Another use of KL divergence is scalable variational inference (VI) as shown
in AVB [15]. VI requires KL divergence estimation between the posterior and
the prior, which becomes nontrivial when an expressive posterior distribution
is used and requires sample based scalable estimation. AVB solved it using GAN
based adversarial formulation and a neural network discriminator. Similarly,
[18] used GAN based adversarial formulation to obtain KL divergence in
amortized inference.

To disentangle latent representations in VAE, [8] proposed TC-VAE which
penalized the KL divergence between marginal latent distribution and the
product of marginals in each dimension. This KL divergence was computed by
minibatch based sampling strategy that gives a biased estimate. None of the
existing works considered the theoretical underpinning of unreliable KL-
divergence estimates, or mitigating the problem by controlling the complexity
of the discriminator function.

##  3 Preliminaries

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 reproducing
kernel [5], such that it satisfies:

  1. 1.

$\forall x\in\mathcal{X},K(.,x)\in\mathcal{H}_{K}$ (Membership property)

  2. 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)$ (Reproducing property)

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
$L_{K}:\mathcal{L}_{2}(d\rho)\to\mathcal{L}_{2}(d\rho)$ [3, 10]:
$(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.

Mean Embedding in RKHS: Let $f:\mathcal{X}\to{\rm I\\!R}$ be a function in
RKHS, $\mathcal{H}_{K}$, and $p$ be a Borel probability measures on
$\mathcal{X}$. If $E_{x\sim p}\sqrt{K(x,x)}<\infty$, then we have
$\mu_{p}\in\mathcal{H}_{K}$ called the mean embedding of the distribution $p$
and defined as [20, 5, 11]: $E_{x\sim p}f=\langle
f,\mu_{p}\rangle_{\mathcal{H}_{K}}$ The condition for the existence of mean
embedding is readily satisfied since we assume $\underset{x,t}{sup}\hskip
2.84544pt{K(x,t)<\infty}$.

##  4 Problem Formulation and Contribution

GAN-based Estimation of KL Divergence: 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 [15, 18], 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^{*}=\max_{f}[{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 [15, 18] that the KL
divergence $KL(p(x)||q(x))$ is given by:

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

Sources of Error: 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). We thus must rewrite eq.(1) as

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

where $h$ is the discriminator function space. Furthermore, we have only a
finite number of samples, say $m$, from the distribution $p$ and $q$. Then,
under finite sample, the optimum discriminator is

$\displaystyle f^{m}_{h}=\max_{f\in h}\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{]}$ (4)

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)]$ (5)

With these definitions, we can now write the error in estimation 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}$ (6)

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.(4), and 3) bias when the considered
function space does not contain optimal function. We leave quantification of
second and third term as future work. Here, we concentrate on quantifying the
probability of deviation-from-mean error which is directly related to observed
variance of the KL estimate.

Overview of Contributions: Note that the deviation is the difference between a
random variable and its mean. Based on this observation, we can bound the
probability of this error using concentration inequality and the complexity of
function space of $f^{m}_{h}$. This requires overcoming the open challenge of
measuring the complexity of neural networks function space. To this end, we
propose to construct a function out of neural network such that it lies on
RKHS. This is our first contribution (Section 5). Then, we proceed to bound
the probability of deviation-from-mean error through the covering number of
the RKHS space. Lemma 1 and Theorem 2 are our contribution (Section 6). Then
we provide insight into how the optimization of eq.(4) might affect
discriminator function space $\mathcal{H}_{K}$. Using ideas from mean
embedding, we prove Lemma 3 and Theorem 3 and provide a geometric insight
(Section 7). This allows us to present a complete story of how the
optimization setup might encourage increase in the complexity of
$\mathcal{H}_{K}$ and how to control it (Section 8).

##  5 Constructing $f$ in RKHS

To construct a function in RKHS, we use an operator $T$ related to integral
operator $L_{K}$ by $L_{K}=TT^{*}$ [3]. The following theorem due to [2] paves
a way for us to construct a neural function in RKHS.

######  Theorem 1.

[[2] Appendix A] A function $f\in\mathcal{L}_{2}(d\rho)$ is in Reproducing
Kernel Hilbert Space, say $\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),$ (7)

for a certain function $g:\mathcal{W}\to 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)$ (8)

Theorem 1 gives us a condition when a square integrable function is guaranteed
to lie in RKHS. We simply choose $g(w)=\textbf{1}$, a constant unit function
over the domain $\mathcal{W}$. This means that we can convert a square
integrable neural network function $f:\mathcal{X}\to{\rm I\\!R}$ into a
function in RKHS, if we make some weights in the neural network stochastic and
average over them. Here, we make the last layer of the neural network to be
drawn from Gaussian distribution, whose parameters are learnt during training.
More precisely, we consider $\psi(x,w)=\phi_{\theta}(x)^{T}w$, where
$\phi_{\theta}(x)$ denotes neural network transformation until the last layer,
and $w$ is the last linear layer sampled from Gaussian distribution. While in
principle any layer could be made stochastic, we chose this architecture to
reduce the computational cost of sampling. The kernel $K$, as defined in
eq.(8), 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)=\phi_{\theta}(x^{*})^{T}(\bar{w}\bar{w}^{T}+\Sigma)\phi_{\theta}(t^{*})$
(9)

where $\bar{w}$ and $\Sigma$ denote the mean and covariance of $w$. We
sometimes denote the kernel $K$ by $K_{\theta}$ to emphasize that it is a
function of neural network parameters, $\theta$.

With this construction, our discriminator function $f$ lies in RKHS denoted by
$\mathcal{H}_{K}$. With $g(w)=\textbf{1}$, it is easy to verify that
$||g||^{2}_{L_{2}(d\tau)}=1$ since $w$ is sampled from a normal distribution.
The inequality in Theorem 1 gives us $||f||^{2}_{\mathcal{H}_{K}}\leq 1$. It
is interesting that the RKHS norm of function $f$ is upper-bounded by 1.
Traditionally, kernel $K$ remains fixed and the norm of the function $f$
determines the complexity of the function space. For example, [16] penalized
the $||f||_{\mathcal{H}_{K}}$ as a way to control the function space while
estimating KL divergence. In our RKHS formulation of neural networks, the
nature of the problem has changed: $||f||_{\mathcal{H}_{K}}$ cannot increase
beyond 1, but the RKHS itself changes during training since it is determined
by the kernel that depends on neural parameters $\theta$. Therefore, the
challenge becomes teasing out how neural parameters $\theta$ affects the
complexity of the discriminator function space and how that affects the
deviation-from-mean error in eq.(6).

##  6 Bounding the Error Probability of KL Estimates

In this section, 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 [10] 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 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 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. 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 [10] to obtain an error
bound in estimating KL divergence with finite samples in Theorem 2.

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

Let $K:\mathcal{X}\times\mathcal{X}\to{\rm I\\!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}}$
(10)

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 inputs space
$\mathcal{X}$. $C_{h}$ is given by $C_{h}=C_{s}\sqrt{||L_{s}||}$ where $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, the RKHS space is embedded in the Sobolev Space $H^{h/2}$
using $L_{s}$ and then the covering number of the Sobolev space is used. Thus
the norm of $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 this
Lemma to bound the estimation error of KL divergence.

######  Theorem 2.

Let ${KL}(f^{m}_{\mathcal{H}_{K}})$ and ${KL}_{m}(f^{m}_{\mathcal{H}_{K}})$ 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}_{K}})-{KL}(f^{m}_{\mathcal{H}_{K}})|\leq\epsilon)\geq
1-2\exp\Bigg{[}\left(\frac{4RC_{s}\sqrt{S_{K}(\theta)||L_{s}||}}{\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 a RKHS space. If
we consider the hypothesis space to be a ball of radius $R$ in Lemma 1 , we
can apply Lemma 2 in it. We fix the radius of discs to
$\eta=\frac{\epsilon}{4\sqrt{S_{K}}}$ in Lemma 1 and substitute
$C_{h}=C_{s}\sqrt{||L_{s}||}$ to obtain the required result. ∎

Theorem 2 shows that the error increases exponentially with the radius of the
RKHS space, complexity of the kernel $S_{K}(\theta)$, and the norm of Sobolev
space embedding $||L_{s}||$. Since we have $||f||_{\mathcal{H}_{K}}\leq 1$, we
can consider our hypothesis space to be a ball of radius 1. To bound
$||L_{s}||$, we need to compute higher order derivatives of $K(x,t)$, which we
leave as future work. This allows us to focus on kernel complexity
$S_{K}(\theta)$, which is exponentially related to the probability of
deviation-from-mean error.

Note that to bound the deviation-from-mean error, we used union bound and
therefore, the bound does not explicitly depend on the function
$f^{m}_{\mathcal{H}_{K}}$, but only depends on the complexity $S_{K}$ of the
function space $\mathcal{H}_{K}$. However, the optimization of discriminator
(eq.(4)) also impacts the complexity $S_{K}(\theta)$. To understand this
effect, in the next section, we present an upper bound on the objective in
eq.(4), and give some geometric intuition connecting the optimization
objective with the kernel complexity $S_{K}(\theta)$. Using this intuition, we
further argue that the optimization of eq.(4) may encourage increment in the
complexity, $S_{K}(\theta)$, thereby increasing the probability of deviation
from the mean.

##  7 Mean Embedding Upper Bound

In addition to deriving complexity bound, another advantage of using RKHS is
that it allows us to use mean embedding representation. This helps us derive
some geometrical insights into the maximization objective in eq.(4), on which
we give an upper bound in Theorem 3.

######  Theorem 3.

Let $f\in\mathcal{H}_{K_{\theta}}$ be a function in RKHS
$\mathcal{H}_{K_{\theta}}$. Then we have the following upper bound on the
objective of KL divergence estimation:

$\displaystyle{}{\frac{1}{m}\sum_{x_{i}\sim
p(x)}\log\sigma(f(x_{i}))+\frac{1}{m}\sum_{x_{j}\sim
q(x)}\log(1-\sigma(f(x_{j})))}\leq\log\sigma[\langle\mu^{m}_{p}-\mu^{m}_{q},f\rangle_{\mathcal{H}_{K}}]$
(11)

and the KL divergence is given by $\text{KL}=\langle\mu^{m}_{p},f\rangle$
where $\mu^{m}_{p}$ and $\mu^{m}_{q}$ represent mean embedding of $m$ samples
from distributions $x_{i}\sim p(x)$ and $x_{j}\sim q(x)$ with respect to
$\mathcal{H}_{K}$.

The following Lemma is useful to prove this theorem.

######  Lemma 3.

${E_{p(x)}\log\sigma(f(x))+E_{q(y)}\log(1-\sigma(f(y)))}\leq\log\sigma[E_{p(x)}(f(x))-E_{q(y)}(f(y))]$

###### Proof Sketch.

It is proved by using $1-\sigma(f)=\sigma(-f)$ and using Jensen’s inequality
since $\log\sigma$ is a concave function. See supplementary material for full
proof. ∎

###### Proof of Theorem 3.

If $f$ lies in the RKHS, then there exists some $\mu^{m}_{p}$ and
$\mu^{m}_{q}$ such that

$\displaystyle\frac{1}{m}\sum_{x_{i}\sim
p(x)}f(x_{i})=\langle\mu^{m}_{p},f\rangle_{\mathcal{H}_{K}},\hskip
5.69046pt\frac{1}{m}\sum_{x_{j}\sim
q(x)}f(x)=\langle\mu^{m}_{q},f\rangle_{\mathcal{H}_{K}}$ (12)

Applying eq.(12) to the Lemma 3 for finite samples, we obtain required result.
∎

Geometric Intuition: Theorem 3 tells us that the upper bound (MEBUB) to the
objective is $\log\sigma$ of the inner product between
$\mu^{m}_{p}-\mu^{m}_{q}$ and $f$. The inner product and KL divergence
estimates have been depicted geometrically in Fig. 1. When the objective is
maximized, MEBUB may also increase which leads to an increase in the inner
product since $\log\sigma$ is monotonic. When this happens, nothing prevents
the midpoint of the mean embeddings, i.e.,
$\frac{\mu^{m}_{p}+\mu^{m}_{q}}{2}$, from going away from the origin in Fig.
1. In the next section, we show how this affects kernel complexity $S_{K}$.

Figure 1: Geometrically representing mean embeddings of two distributions, and
their relation to the maximization objective and KL divergence.

##  8 Fitting Pieces and Complexity Control

Theorem 2 shows that the error bound of the KL estimate is exponentially
controlled by the kernel complexity $S_{K}(\theta)=\underset{x,t}{sup}\hskip
2.84544ptK_{\theta}(x,t)$. Since the mean of a vector is upper bounded by
supremum,

$\displaystyle||\mu^{m}_{p}+\mu^{m}_{q}||_{\mathcal{H}_{K}}$
$\displaystyle=\sqrt{\frac{1}{m^{2}}\sum_{i,j}K_{\theta}(y_{i},y_{j})+2K_{\theta}(y_{i},z_{j})+K_{\theta}(z_{i},z_{j})}$
(13) $\displaystyle\leq
2\sqrt{\underset{x\in\\{Y,Z\\},t\in\\{Y,Z\\}}{sup}K_{\theta}(x,t)}=2\sqrt{S_{K}(\theta)}$
(14)

As the training progresses in maximizing the objective in eq.(4), the
algorithm tries to do two things: 1) align $f$ with $\mu^{m}_{p}-\mu^{m}_{q}$
and 2) increase norms $||\mu^{m}_{p}-\mu^{m}_{q}||_{\mathcal{H}_{K}}$ and
$||f||_{\mathcal{H}_{K}}$. For fixed
$\langle\mu^{m}_{p},\mu^{m}_{q}\rangle_{\mathcal{H}_{K}}/(||\mu^{m}_{p}||.||\mu^{m}_{q}||)$,
we can show that the ratio
$||(\mu^{m}_{p}-\mu^{m}_{q})||_{\mathcal{H}_{K}}/||\mu^{m}_{p}+\mu^{m}_{q}||_{\mathcal{H}_{K}}$
also remains unchanged. Under this assumption, we could say that maximizing
eq.(4) could lead to increment of
${||\mu^{m}_{p}+\mu^{m}_{q}||}_{\mathcal{H}_{K}}$, and nothing would stop the
network from going this way. When this happens, the inequality in eq.(14)
suggests that $S_{K}(\theta)$ also increases, thereby increasing the
probability of deviation-from-the-mean error in the KL estimate by Theorem 2.
In other words, as we train the neural discriminator, the neural network
parameters $\theta$ change such that the complexity of RKHS might itself keep
increasing which causes exponential growth in the sample complexity.

Figure 2: Comparison of KL divergence estimates using infinite samples
(purple) versus finite samples with (blue) and without (red) complexity
control.

To control the complexity of RKHS space, we can control $S_{K}(\theta)$ during
the training of the neural network. To do this in a scalable way compatible
with neural networks, we use gradient descent based optimization. Computation
of gradient of $S_{K}(\theta)$ w.r.t $\theta$ is straightforward using
definition of $K_{\theta}(x,t)$ and can be easily realized by using
backpropagation. Ideally, $S_{K}(\theta)$ is $max.K_{\theta}(x,t)$ over all
the data-pairs $(x,t)\in\mathcal{X}\times\mathcal{X}$, which requires passing
all the datapoints through neural network. Instead, we simply compute supremum
over the minibatch matrix which contains the $2b\times 2b$ entries
corresponding to every pair in $2b$ elements ($b$ from each distribution
$p(x)$ and $q(x)$). This is obviously a lowerbound – denoted by
$S_{mini}(\theta)$ – of $S_{K}(\theta)$. To penalize the RKHS space that are
high in complexity, we add a regularization term with parameter $\lambda$ to
maximize a modified objective:

$\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})))}-\lambda.S_{mini}^{\gamma}$ (15)

where $\gamma$ is an estimation of $\frac{n}{h}$ and treated as a
hyperparameter. Optimization of eq.(15) 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.

##  9 Experimental Results

Experimental Setup: 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 in a simple case of 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 $61.1$, respectively and use $m=5000$ samples from each
distribution to estimate KL in the finite case.

As the discriminator, we use a fully connected neural network with two hidden
layers. The number of hidden units are varied to understand the effect of the
discriminator complexity on the fluctuation of the KL estimate. The dimensions
are kept identical between the neural-net discriminator and the RKHS
discriminator, the latter being different only in that its last layer is
stochastic. We perform random estimation experiment 30 times and report the
mean, standard deviation, scatter and box plots

Finite v.s. Infinite Samples: In infinite samples experiment, we assume that
we can continuously sample from the model generating data from the two given
distributions. The results of KL estimates using infinite samples is shown in
Fig. 2 and Table.1 left, in comparison with estimates using finite samples
without controlling the complexity of the neural-net discriminator. We observe
that when we use infinite samples, we obtain an estimate with low variance and
values close to the analytical truth in KL = $1.3$ and KL = $13.8$ and an
underestimate when KL = $61.1$. In contrast, when we use finite samples
without controlling the complexity of the neural-net discriminator, the
estimates fluctuated heavily confirming our hypothesis: we need to control the
complexity of the function when the number of samples is finite, or else the
probability of estimation error increases.

Complexity Control: Fig. 2 and Table.1 left compare the estimation of KL
divergence with and without controlling the discriminator complexity. With
discriminator complexity penalized in eq. (15), the KL estimates are much more
reliable (low variance) and closer to the estimates from infinite samples.
Note in Fig.2 that, without complexity penalization, the erratic behavior of
the KL estimator worsens as the number of hidden layers increases in the
discriminator. This is consistent with our theory because increasing the
number of hidden layers increases the complexity of the discriminator. This
highlights the need of higher degree penalization of the discriminator
complexity as a neural network with increased capacity is used to estimate
higher values of KL divergence.

Table 1: Left: Comparison of KL divergence estimates using different methods
(hidden layer dim =25). Right: The effect of the regularization parameter
$\boldsymbol{\lambda}$; hidden layer dimension = 20

Method | True KL
---|---
1.3 |  13.8 |  61.1
MINE |  Unstable |  Unstable |  Unstable
VDM |  Unstable |  Unstable |  Unstable
Infinite sample |  $1.4\pm 0.05$ |  $12.6\pm 1.5$ |  $32.4\pm 7.9$
NN Disc |  $2.1\pm 0.42$ |  $20.6\pm 9.8$ |  $150\pm 65$
Complexity control |  $1.5\pm 0.15$ |  $13.6\pm 2.4$ |  $24\pm 8.2$

$\boldsymbol{\lambda}$ | True KL
---|---
1.3 |  13.8 |  61.1
5e-5 |  $1.46\pm 0.22$ |  $16.65\pm 10.4$ |  $116.7\pm 116$
1e-4 |  $1.56\pm 0.25$ |  $30.97\pm 10.5$ |  $39.17\pm 18.5$
5e-4 |  $1.47\pm 0.11$ |  $13.44\pm 2.68$ |  $18.36\pm 3.9$

Effect of Regularization Parameter: Table 1 right shows the effect of the
regularization parameter $\lambda$ that tunes the level of complexity control
in eq. (15). The fluctuation in estimates decreases as we increase the value
of $\lambda$. We also repeated these experiments varying the latent dimension
(see Fig.3 in appendix) and found that the pattern holds in all cases.
Furthermore, for a discriminator with low complexity (e.g., latent dimension =
10), a smaller value of $\lambda$ is sufficient to yield low-variance
estimate. As the size of the hidden layer increases, we need to penalize the
complexity aggressively with a higher value of $\lambda$ in order to obtain
the same level of consistency. This further supports our theory.

Underestimation for High KL Divergence We observe in Fig.2 and Table. 1 that,
for KL = $61.1$, results from both infinite samples and finite samples with
complexity control give underestimated KL divergences even though they reduce
fluctuation significantly. This is not surprising since we were focusing on
deviation-from-mean error. The total estimation error consists of two
additional errors: discriminator induced error and the bias (see eq.6). For
the small KL divergence, simply controlling complexity was sufficient to
minimize all the errors, but for higher value, it is no longer sufficient. The
underestimation might be either because of the bias or error induced by
incorrect discriminator function. High bias might be caused if we control the
function space too much such that the optimum discriminator
$f^{*}_{\mathcal{H}_{K}}$ is not close to true discriminator function, $f^{*}$
(see bias-variance trade-off [10, 6]). It would be an interesting future
direction to quantify all three error terms in eq.(6).

##  10 Conclusions & Discussion

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 penalizing
the kernel complexity in a scalable way using neural networks.

The idea of constructing a neural-net function in RKHS and complexity control
could also be useful in stabilizing GANs, or potentially in improving
generalization of neural networks. Several papers have identified issues with
the stability of GANs [14, 13, 21]. One common understanding is that, in its
raw form, we do not enforce the discriminator function to be smooth or regular
around the neighborhood of its inputs. Currently, the most successful way to
stabilize GANs is to enforce smoothness by gradient penalization. Even in
variations like Wasserstein GAN [1, 12] and MMD GAN [7], gradient penalty is
crucial to achieve stable results. On the light of the present analysis, we
believe that the gradient penalty can be thought as one way to control the
complexity of the discriminator. The objective and nature of optimization is
such that the complexity of discriminator is bound to increase and therefore
some way of decreasing complexity is a must. Similarly, generalization of
neural network classifiers and regressors could be improved with complexity
control.

## References

  * [1] Martin Arjovsky, Soumith Chintala, and Léon Bottou.  Wasserstein generative adversarial networks.  In International Conference on Machine Learning, pages 214–223, 2017.
  * [2] Francis Bach.  Breaking the curse of dimensionality with convex neural networks.  The Journal of Machine Learning Research, 18(1):629–681, 2017.
  * [3] Francis Bach.  On the equivalence between kernel quadrature rules and random feature expansions.  The Journal of Machine Learning Research, 18(1):714–751, 2017.
  * [4] Mohamed Ishmael Belghazi, Aristide Baratin, Sai Rajeshwar, Sherjil Ozair, Yoshua Bengio, Aaron Courville, and Devon Hjelm.  Mutual information neural estimation.  In International Conference on Machine Learning, pages 531–540, 2018.
  * [5] A. Berlinet and C. Thomas-Agnan.  Reproducing Kernel Hilbert Spaces in Probability and Statistics.  Springer US, 2011.
  * [6] Christopher M Bishop.  Pattern recognition and machine learning.  springer, 2006.
  * [7] Mikołaj Bińkowski, Dougal J. Sutherland, Michael Arbel, and Arthur Gretton.  Demystifying MMD GANs.  In International Conference on Learning Representations, 2018.
  * [8] Tian Qi Chen, Xuechen Li, Roger B Grosse, and David K Duvenaud.  Isolating sources of disentanglement in variational autoencoders.  In Advances in Neural Information Processing Systems, pages 2610–2620, 2018.
  * [9] Xi Chen, Yan Duan, Rein Houthooft, John Schulman, Ilya Sutskever, and Pieter Abbeel.  Infogan: Interpretable representation learning by information maximizing generative adversarial nets.  In Advances in neural information processing systems, pages 2172–2180, 2016.
  * [10] Felipe Cucker and Steve Smale.  On the mathematical foundations of learning.  Bulletin of the American mathematical society, 39(1):1–49, 2002\.
  * [11] Arthur Gretton, Karsten M Borgwardt, Malte J Rasch, Bernhard Schölkopf, and Alexander Smola.  A kernel two-sample test.  Journal of Machine Learning Research, 13(Mar):723–773, 2012.
  * [12] Ishaan Gulrajani, Faruk Ahmed, Martin Arjovsky, Vincent Dumoulin, and Aaron C Courville.  Improved training of wasserstein gans.  In I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett, editors, Advances in Neural Information Processing Systems 30, pages 5767–5777. Curran Associates, Inc., 2017.
  * [13] Naveen Kodali, Jacob Abernethy, James Hays, and Zsolt Kira.  On convergence and stability of gans.  arXiv preprint arXiv:1705.07215, 2017.
  * [14] Lars Mescheder, Andreas Geiger, and Sebastian Nowozin.  Which training methods for gans do actually converge?  In International Conference on Machine Learning, pages 3481–3490, 2018.
  * [15] Lars Mescheder, Sebastian Nowozin, and Andreas Geiger.  Adversarial variational bayes: Unifying variational autoencoders and generative adversarial networks.  In International Conference on Machine Learning (ICML), 2017.
  * [16] XuanLong Nguyen, Martin J Wainwright, and Michael I Jordan.  Estimating divergence functionals and the likelihood ratio by convex risk minimization.  IEEE Transactions on Information Theory, 56(11):5847–5861, 2010\.
  * [17] Sebastian Nowozin, Botond Cseke, and Ryota Tomioka.  f-gan: Training generative neural samplers using variational divergence minimization.  In Advances in neural information processing systems, pages 271–279, 2016.
  * [18] Casper Kaae Sønderby, Jose Caballero, Lucas Theis, Wenzhe Shi, and Ferenc Huszár.  Amortised map inference for image super-resolution.  ICLR, 2017.
  * [19] Jiaming Song and Stefano Ermon.  Understanding the limitations of variational mutual information estimators.  arXiv preprint arXiv:1910.06222, 2019.
  * [20] Bharath K Sriperumbudur, Arthur Gretton, Kenji Fukumizu, Bernhard Schölkopf, and Gert RG Lanckriet.  Hilbert space embeddings and metrics on probability measures.  Journal of Machine Learning Research, 11(Apr):1517–1561, 2010.
  * [21] Hoang Thanh-Tung, Truyen Tran, and Svetha Venkatesh.  Improving generalization and stability of generative adversarial networks.  In International Conference on Learning Representations, 2019.

##  Appendix A Bounding the Error Probability of KL Estimates

To obtain the probability of deviation-from-mean error, we first bound the
error probability in terms of the covering number in Lemma 1. Then, we use an
estimate of the covering number of RKHS due to [10] to obtain a bound of error
probability in terms of the kernel $K_{\theta}$ in Theorem 2.

######  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)}$ (16)

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}}$ (17)

Hence, we have:

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

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.(18), we obtain,

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

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

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

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}}}$
(21)
$\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}}{8M^{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.17, we have
$||f^{m}_{\mathcal{H}_{K}}||_{\infty}\leq\sqrt{S_{K}}||f^{m}_{\mathcal{H}_{K}}||_{\mathcal{H}_{K}}$.
Since by construction, $||f^{m}_{\mathcal{H}_{K}}||_{\mathcal{H}_{K}}\leq 1$,
$f^{m}_{\mathcal{H}_{K}}$ is bounded if $S_{K}$ is bounded, which is true by
assumption and seems to hold true in experiments. ∎

######  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. 17 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 [10] to obtain an error bound in
estimating KL divergence with finite samples in Theorem 2.

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

Let $K:\mathcal{X}\times\mathcal{X}\to{\rm I\\!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}}$
(22)

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_{s}\sqrt{||L_{s}||}$ (23)

where $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, the RKHS space is embedded in the Sobolev Space $H^{h/2}$
using $L_{s}$ and then covering number of Sobolev space is used. Thus the norm
of $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 this Lemma to
bound the estimation error of KL divergence.

######  Theorem 2.

Let ${KL}(f^{m}_{\mathcal{H}_{K}})$ and ${KL}_{m}(f^{m}_{\mathcal{H}_{K}})$ 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}_{K}})-{KL}(f^{m}_{\mathcal{H}_{K}})|\leq\epsilon)\geq
1-2\exp\Bigg{[}\left(\frac{4RC_{s}\sqrt{S_{K}||L_{s}||}}{\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 a RKHS space. If
we consider the hypothesis space to be a ball of radius $R$ in Lemma 1 , 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}})-{KL}(f^{m}_{\mathcal{H}_{K}})|\leq\epsilon)\geq
1-2\exp\Big{[}\left(\frac{4\sqrt{S_{K}}RC_{h}}{\epsilon}\right)^{\frac{2n}{h}}-\frac{m\epsilon^{2}}{4M^{2}}\Big{]}$

Substituting $C_{h}=C_{s}\sqrt{||L_{s}||}$ gives the required result. ∎

##  Appendix B Mean Embedding Upper Bound

Using mean embedding representation of functions in RKHS space, we derive some
geometrical insights into the maximization objective. Theorem 3 first gives an
upper bound to the optimization objective, which depends on Lemma 3 as proved
next.

######  Lemma 3.

${E_{p(x)}\log\sigma(f(x))+E_{q(y)}\log(1-\sigma(f(y)))}\leq\log\sigma[E_{p(x)}(f(x))-E_{q(y)}(f(y))]$

###### Proof.

$\displaystyle{E_{p(x)}\log\sigma(f(x))+E_{q(x)}\log(1-\sigma(f(x)))}$
$\displaystyle={E_{p(x)}\log\sigma(f(x))+E_{q(x)}\log(\frac{1}{1+\exp(f)})}$
$\displaystyle={E_{p(x)}\log\sigma(f(x))+E_{q(x)}\log(\frac{\exp(-f)}{1+\exp(-f)})}$
$\displaystyle={E_{p(x)}\log\sigma(f(x))+E_{q(x)}\log\sigma(-f(x))}$
$\displaystyle\leq\log\sigma[{E_{p(x)}(f(x))]+\log\sigma[E_{q(x)}(-f(x))}]$
$\displaystyle\leq\log\sigma[{E_{p(x)}(f(x))-E_{q(x)}(f(x))}$

where we used the fact that $\log\sigma$ is a concave function and applied
Jensen’s inequality in last two lines and linearity of expectation in the last
line. ∎

##  Appendix C Algorithm

Algorithm 1 details the algorithm to estimate KL divergence with complexity
control.

Algorithm 1 KL divergence estimation with complexity control

1: Fix minibatch size, $b$, hyperparameter $\gamma$, number of samples $m$,
$flat\\_n=100$, $idx=0$,$\ell_{min}=\infty$

2: Initialize the neural network parameters ${\theta}$, last layer
$w\sim\mathcal{N}(\bar{w},LL^{T})$, such that $\bar{w}=0,LL^{T}=\text{I}$

3: for iteration $iter$ in 1 to $iter_{max}$ do

4: $kl_{sum}=0,\ell_{adv}=0$, $n\\_batch=(m/b)$

5: for iteration k in 1 to $n\\_batch$ do

6: Sample minibatch $\\{x_{i}\\}_{i=1}^{b}$ from $p(x)$ and
$\\{y_{i}\\}_{i=1}^{b}$ from $q(x)$, and
$J=\\{\\{x_{i}\\}_{i=1}^{b},\\{y_{i}\\}_{i=1}^{b}\\}$

7: For each $x_{i},y_{i}$, sample $\epsilon\sim\mathcal{N}(0,I)$ and obtain
samples $\\{w_{j}\\}_{j=1}^{d}$ where $w_{j}=\bar{w}+L\epsilon_{j}$

8:
$f(x)_{i},=\frac{1}{d}\sum_{j}\phi_{\theta}(x_{i})^{T}w_{j}$$f(y)_{i}=\frac{1}{d}\sum_{j}\phi_{\theta}(y_{i})^{T}w_{j}$

9: $loss_{d}=-\frac{1}{b}\sum_{i}\log\sigma(f(x)_{i})+\log(1-f(y)_{i})$

10: $S_{mini}=\underset{x\in J,t\in
J}{max}\phi_{\theta}(x)^{T}(\bar{w}\bar{w}^{T}+\Sigma)\phi_{\theta}(t)$

11: Backpropagate $loss=loss_{d}+\lambda.S_{mini}^{\gamma}$ and update
$\theta,\bar{w},L$

12:
$kl_{sum}=kl_{sum}+\frac{1}{b}\sum_{i}\log\sigma(f(x)_{i})$$\ell_{adv}=\ell_{adv}+loss_{d}$

13: end for

14: $\ell=\ell_{adv}/n\\_{batch}$, $kl_{iter}=kl_{sum}/n\\_batch$

15: if $\ell<\ell_{min}$ then

16: $kl=kl_{iter}$, $idx=iter$

17: else if $iter>idx+flat\\_n$ then

18: return $kl$

19: end if

20: end for

21: return $kl$

##  Appendix D Experimental Results

###  D.1 Details of experimental setup

Neural RKHS discriminator architecture (Proposed method)
Fully connected
Leaky ReLU
Fully connected
Leaky ReLU

For the RKHS discriminator, this gives $\phi_{\theta}(x)$ for the input data
$x$. Then, $f(x)$ needs to defined as in line 8 of Algorithm 1. Similarly,
total loss with complexity penalization is computed as line 11 in Algorithm 1.

Neural network discriminator architecture
Fully connected
Leaky ReLU
Fully connected
Leaky ReLU
Fully connected

For the Neural net discriminator, this would directly give $f(x)$ for the
input data $x$. Also, loss would be defined by line 9, no penalization as in
line 11 of Algorithm 1.

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

Figure 3: The effect of the regularization parameter $\lambda$ in KL estimates
(y-axis) plotted against the varying hidden layer dimension for each KL
divergence value.

###  D.2 Effect of Regularization Parameter $\lambda$

Fig. 3 shows the effect of the regularization parameter $\lambda$ that tunes
the level of complexity control in estimating the KL divergence. As expected,
in all cases, we observe that the fluctuation in estimates decreases as we
increase the value of $\lambda$.


# A new approach to predict changes in physical condition: A new extension of
the classical Banister model

Marcos Matabuena, Rosana Rodríguez-López  Departamento de Análise Matemática,
Estatística e Optimización,
Universidade de Santiago de Compostela,
15782, Santiago de Compostela, Spain<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

In this article, a new model based on techniques of differential equations is
introduced to predict the athletic performance based training load and a data
sample of the physical form of athletes arises. This model is an extension of
the classical model of Banister but, in this case, unlike the classical
Banister model, the variation produced in the athletic performance depends,
not only on the current training load, but also on the training performed the
previous day. The model has been validated with the training data of a cyclist
taken from the reference [6], obtaining an excellent fit of the predicted data
with respect to the experimental data.

###### keywords:

Mathematical models , functional differential equations , difference
equations.

###### MSC:

[2010] 34K06 , 39A06, 39A60

††journal: Applied Mathematics and Computation

##  1 Introduction

Athletes, in their quest to optimize performance, try to use mathematical
tools (Issurin, 2010 [12]), but nowadays there are no tools to know the effect
of workout on the final performance (González Badillo & Gorostiaga, 1995
[10]). A first attempt at the problem of modeling and predicting changes
induced by fitness training appeared in 1975 with the contribution of the
physiologist Banister and his coworkers (Banister et al., 1975 [1]), whose
model has been frequently used in the literature to predict the effects of
training on physical performance, see, for example, its validation in
triathlon athletes by Banister et al. in 1999 [2], the case of a swimmer
presented in 1976 by Calvert et al. [4], studies in cross-country skiers made
in 1992 by Candau et al. [5], the convenience to exercise physiology
practitioners and researchers in the use of mathematical models (see Clarke &
Skiba, 2013 [6]), the optimal design of a training strategy (see Fitz-Clarke
et al., 1991 [9]), the performance of runners analyzed in 1990 by Morton et
al. [13], or the recent consideration of altitude training by Rodríguez et al.
[16], among others. The interest of the application of computational
intelligence in sports has also been defended (see [8]).

Since the introduction of the Banister model, there have been proposed
modified versions of the model, some of which are those introduced by Busso in
2003 [3], or Coggan in 2006 [7], among others, and some new models have been
built through regression techniques, neural networks (see Pfeiffer & Hohmann,
2012 [15]) or time series (see Pfeiffer, 2008). In this paper, we propose a
new model based on techniques of functional differential equations and
nonlinear regression that substantially improves the accuracy of the above
methods, besides overcoming some of their limitations (see Hellard et al.,
2006 [11]). The model proposed has been successfully applied to study and
predict the effect of training a cyclist in the task of optimizing a workout
plan, constituting what might be a reliable and useful tool for the future, in
their decision making and when applied to other sports.

##  2 Mathematical model

In the work [14], the authors consider the following equation:

$g^{\prime}(t)=-\frac{1}{\tau_{1}}g(t)+w(t),$

whose solution is given by

$g(t)=g(0)+\int_{0}^{t}w(s)e^{\frac{-1}{\tau_{1}}(t-s)}\,ds=\int_{0}^{t}w(s)e^{\frac{-1}{\tau_{1}}(t-s)}\,ds,$

which gives the following approximation

$g(n)=\sum_{i=0}^{n-1}w(i)e^{\frac{-1}{\tau_{1}}(n-i)}.$

If we consider $w(0)=0$, then it coincides with equation (3) in [14].

Next, we consider the delay differential model

$g^{\prime}(t)=-\frac{1}{\tau_{1}}g(t)-\frac{1}{\tau_{2}}g(t-1)+w(t).$ (1)

and consider an equivalent formulation.

######  Theorem 1.

The equation (1) is equivalent to

$g(t)e^{\frac{1}{\tau_{1}}(t-k)}-g(k)=\int_{k}^{t}\left[w(s)-\frac{1}{\tau_{2}}g(s-1)\right]e^{\frac{1}{\tau_{1}}(s-k)}\,ds,$

for $t\in(k,k+1]$, $k=0,1,\ldots$

From Theorem 1, we have

$g(t)=g(k)e^{-\frac{1}{\tau_{1}}(t-k)}+\int_{k}^{t}\left[w(s)-\frac{1}{\tau_{2}}g(s-1)\right]e^{-\frac{1}{\tau_{1}}(t-s)}\,ds,$

for $t\in(k,k+1]$, $k=0,1,\ldots$, so that we can consider the following
approximation:

$\begin{split}g(k+1)&=g(k)e^{-\frac{1}{\tau_{1}}}+\left[w(k)-\frac{1}{\tau_{2}}g(k-1)\right]e^{-\frac{1}{\tau_{1}}}\\\
&=\left[w(k)+g(k)-\frac{1}{\tau_{2}}g(k-1)\right]e^{-\frac{1}{\tau_{1}}},k=0,1,\ldots\end{split}$

Note that, for $t\in[0,1]$, we get

$g(t)=\int_{0}^{t}\left[w(s)-\frac{1}{\tau_{2}}g(s-1)\right]e^{-\frac{1}{\tau_{1}}(t-s)}\,ds=\int_{0}^{t}w(s)e^{-\frac{1}{\tau_{1}}(t-s)}\,ds,$

which is coincident with [14].

Other way to calculate the solution to (1), by integration between $0$ and
$t$, which provides a model more similar to that in [14], would be to write
the equation equivalently as

$g(t)e^{\frac{1}{\tau_{1}}t}=g(0)+\int_{0}^{t}\left[w(s)-\frac{1}{\tau_{2}}g(s-1)\right]e^{\frac{1}{\tau_{1}}s}\,ds,\quad
t\geq 0,$

or

$\begin{split}g(t)&=g(0)e^{-\frac{1}{\tau_{1}}t}+\int_{0}^{t}\left[w(s)-\frac{1}{\tau_{2}}g(s-1)\right]e^{-\frac{1}{\tau_{1}}(t-s)}\,ds\\\
&=\int_{0}^{t}\left[w(s)-\frac{1}{\tau_{2}}g(s-1)\right]e^{-\frac{1}{\tau_{1}}(t-s)}\,ds,\quad
t\geq 0,\end{split}$

therefore

$g(n)=\int_{0}^{n}\left[w(s)-\frac{1}{\tau_{2}}g(s-1)\right]e^{-\frac{1}{\tau_{1}}(n-s)}\,ds,\quad
t\geq 0,$

and we can take, similarly to the study in [14], the approximation

$g(n)=\sum_{i=0}^{n-1}\left[w(i)-\frac{1}{\tau_{2}}g(i-1)\right]e^{\frac{-1}{\tau_{1}}(n-i)}.$

If we suppose that $w(0)=0$, since $g(-1)=0$, the previous expression is equal
to

$g(n)=\sum_{i=1}^{n-1}\left[w(i)-\frac{1}{\tau_{2}}g(i-1)\right]e^{\frac{-1}{\tau_{1}}(n-i)}.$
(2)

###  2.1 Model proposed

The model proposed is based on the study of the properties of differential
equations with delay and the estimation of parameters. We first state the
problem in mathematical terms, as follows:

$p(t)=p(0)+k_{1}g(t)-k_{2}h(t),$ (3)

where $k_{1}$ and $k_{2}$ are two positive constants and $p(0)$ is the initial
physical condition. The values of the functions $g(t)$ and $h(t)$ are obtained
by solving the following system of differential equations:

$\left\\{\begin{array}[]{c}g^{\prime}(t)=w(t)-\frac{1}{\tau_{1}}g(t)-\frac{1}{\tau_{2}}g(t-1),\\\
\\\
h^{\prime}(t)=w(t)-\frac{1}{\tau_{3}}h(t)-\frac{1}{\tau_{4}}h(t-1),\end{array}\right.$

where $w(t)$ represents the load at the instant of time $t$ and $\tau_{1}$,
$\tau_{2}$, $\tau_{3}$, $\tau_{4}$ are positive constants. Since, in practice,
we work with measurements at fixed times, we consider a discrete model and
approximate the solution to the differential system by considering the
following difference equation (see (2)):

$p(n)=k_{1}\sum_{i=1}^{n-1}\left[w(i)-\frac{1}{\tau_{2}}g(i-1)\right]e^{\frac{1}{\tau_{1}}(n-i)}-k_{2}\sum_{i=1}^{n-1}\left[w(i)-\frac{1}{\tau_{4}}h(i-1)\right]e^{\frac{1}{\tau_{3}}(n-i)},$

where $n$ denotes the $n$th instant of time. The last equation represents the
explicit solution of the discretized (difference) equation of problem (3).

At this point, we can use nonlinear regression techniques, or any heuristic
algorithm, to estimate the constants
$k_{1},k_{2},\tau_{1},\tau_{2},\tau_{3},\tau_{4}$ by the method of least
squares, using a sample of physical condition $p(1),\ldots,p(s)$, at different
times. In this work, we use the experimental data of a cyclist extracted from
the 2013 article by Clarke & Skiba [6]. To optimize the constants, we use the
Solver in Microsoft Office and we represent the results with the statistical
package R.

###  2.2 Results obtained

The result of the $R^{2}$ test between the experimental data and the
mathematical model proposed is $0.9903$. In Figure 1, we find a graph where
the red points represent the data and the blue curve the solution predicted by
the model. In Figure 2, we see the athlete on his training-workload, which
determines the change in his physical condition. We observe that this graph
presents high and low values, which correspond to the effect of the
assimilation of the training and the fatigue states, respectively, what
actually happens in reality.

Figure 1: Experimental data (in red) and the solution predicted by the model
(in blue) Figure 2: Function $w(t)$

##  3 Other formulations of the model

In this section, we show other different formulations of the model.

###  3.1 Problem with three delays

Next, we consider

$g^{\prime}(t)=-\frac{1}{\tau_{1}}g(t)-\frac{1}{\tau_{2}}g(t-1)-\frac{1}{\tau_{3}}g(t-2)-\frac{1}{\tau_{4}}g(t-3)+w(t).$
(4)

The solution to (4) can be written as:

$g(t)=g(0)e^{-\frac{1}{\tau_{1}}t}+\int_{0}^{t}\left[w(s)-\sum_{k=1}^{3}\frac{1}{\tau_{k+1}}g(s-k)\right]e^{-\frac{1}{\tau_{1}}(t-s)}\,ds,\quad
t\geq 0,$

that is,

$g(t)=\int_{0}^{t}\left[w(s)-\sum_{k=1}^{3}\frac{1}{\tau_{k+1}}g(s-k)\right]e^{-\frac{1}{\tau_{1}}(t-s)}\,ds,\quad
t\geq 0,$

so that,

$g(n)=\int_{0}^{n}\left[w(s)-\sum_{k=1}^{3}\frac{1}{\tau_{k+1}}g(s-k)\right]e^{-\frac{1}{\tau_{1}}(n-s)}\,ds$

and we can take, similarly to [14], the approximation

$g(n)=\sum_{i=0}^{n-1}\left[w(i)-\sum_{k=1}^{3}\frac{1}{\tau_{k+1}}g(i-k)\right]e^{\frac{-1}{\tau_{1}}(n-i)}$

or

$g(n)=\sum_{i=0}^{n-1}\left[w(i)-\frac{1}{\tau_{2}}g(i-1)-\frac{1}{\tau_{3}}g(i-2)-\frac{1}{\tau_{4}}g(i-3)\right]e^{\frac{-1}{\tau_{1}}(n-i)}.$

If we suppose that $w(0)=0$, since $g(s)=0$, $\forall s\in[-3,0]$, we have

$\begin{split}g(n)&=\sum_{i=1}^{n-1}\left[w(i)-\frac{1}{\tau_{2}}g(i-1)-\frac{1}{\tau_{3}}g(i-2)-\frac{1}{\tau_{4}}g(i-3)\right]e^{\frac{-1}{\tau_{1}}(n-i)}\\\
&=\sum_{i=1}^{n-1}w(i)e^{\frac{-1}{\tau_{1}}(n-i)}-\frac{1}{\tau_{2}}\sum_{i=1}^{n-1}g(i-1)e^{\frac{-1}{\tau_{1}}(n-i)}-\frac{1}{\tau_{3}}\sum_{i=1}^{n-1}g(i-2)e^{\frac{-1}{\tau_{1}}(n-i)}\\\
&\hskip
30.00005pt-\frac{1}{\tau_{4}}\sum_{i=1}^{n-1}g(i-3)e^{\frac{-1}{\tau_{1}}(n-i)}\\\
&=\sum_{i=1}^{n-1}w(i)e^{\frac{-1}{\tau_{1}}(n-i)}-\frac{1}{\tau_{2}}\sum_{i=1}^{n-1}g(i-1)e^{\frac{-1}{\tau_{1}}(n-i)}-\frac{1}{\tau_{3}}\sum_{j=0}^{n-2}g(j-1)e^{\frac{-1}{\tau_{1}}(n-j-1)}\\\
&\hskip
30.00005pt-\frac{1}{\tau_{4}}\sum_{j=-1}^{n-3}g(j-1)e^{\frac{-1}{\tau_{1}}(n-j-2)}\\\
&=\sum_{i=1}^{n-1}w(i)e^{\frac{-1}{\tau_{1}}(n-i)}-\frac{1}{\tau_{2}}\sum_{i=1}^{n-1}g(i-1)e^{\frac{-1}{\tau_{1}}(n-i)}-\frac{1}{\tau_{3}}\sum_{j=1}^{n-2}g(j-1)e^{\frac{-1}{\tau_{1}}(n-j-1)}\\\
&\hskip
30.00005pt-\frac{1}{\tau_{4}}\sum_{j=1}^{n-3}g(j-1)e^{\frac{-1}{\tau_{1}}(n-j-2)}\\\
&=\sum_{i=1}^{n-1}w(i)e^{\frac{-1}{\tau_{1}}(n-i)}-\sum_{i=1}^{n-3}g(i-1)\left[\frac{1}{\tau_{2}}e^{\frac{-1}{\tau_{1}}(n-i)}+\frac{1}{\tau_{3}}e^{\frac{-1}{\tau_{1}}(n-i-1)}+\frac{1}{\tau_{4}}e^{\frac{-1}{\tau_{1}}(n-i-2)}\right]\\\
&\hskip
30.00005pt-g(n-3)\left[\frac{1}{\tau_{2}}e^{\frac{-2}{\tau_{1}}}+\frac{1}{\tau_{3}}e^{\frac{-1}{\tau_{1}}}\right]-g(n-2)\frac{1}{\tau_{2}}e^{\frac{-1}{\tau_{1}}}\\\
&=\sum_{i=1}^{n-1}w(i)e^{\frac{-1}{\tau_{1}}(n-i)}-\sum_{i=1}^{n-3}g(i-1)\left[\frac{1}{\tau_{2}}+\frac{1}{\tau_{3}}e^{\frac{1}{\tau_{1}}}+\frac{1}{\tau_{4}}e^{\frac{2}{\tau_{1}}}\right]e^{\frac{-1}{\tau_{1}}(n-i)}\\\
&\hskip
30.00005pt-g(n-3)\left[\frac{1}{\tau_{2}}e^{\frac{-2}{\tau_{1}}}+\frac{1}{\tau_{3}}e^{\frac{-1}{\tau_{1}}}\right]-g(n-2)\frac{1}{\tau_{2}}e^{\frac{-1}{\tau_{1}}}.\end{split}$

Following for problem (4) a formulation similar to the first one exposed, we
have:

$g(t)=g(k)e^{-\frac{1}{\tau_{1}}(t-k)}+\int_{k}^{t}\left[w(s)-\frac{1}{\tau_{2}}g(s-1)-\frac{1}{\tau_{3}}g(s-2)-\frac{1}{\tau_{4}}g(s-3)\right]e^{-\frac{1}{\tau_{1}}(t-s)}\,ds,$

for $t\in(k,k+1],$ $k=0,1,\ldots$, so that, we can give the following
approximation:

$\begin{split}g(k+1)&=g(k)e^{-\frac{1}{\tau_{1}}}+\left[w(k)-\frac{1}{\tau_{2}}g(k-1)-\frac{1}{\tau_{3}}g(k-2)-\frac{1}{\tau_{4}}g(k-3)\right]e^{-\frac{1}{\tau_{1}}}\\\
&=\left[w(k)+g(k)-\frac{1}{\tau_{2}}g(k-1)-\frac{1}{\tau_{3}}g(k-2)-\frac{1}{\tau_{4}}g(k-3)\right]e^{-\frac{1}{\tau_{1}}},\end{split}$
(5)

for $k=0,1,\ldots$

###  3.2 Modelling through an integral equation

We consider

$g^{\prime}(t)=-\frac{1}{\tau_{1}}g(t)+\tau_{5}\int_{0}^{t}\delta(t,s)g(s)\,ds+w(t).$
(6)

We can impose, for instance, that $\delta(t,s)=0$ for $s<t-4$. The solution to
(6) is obtained as

$\begin{split}g(t)&=g(0)e^{-\frac{1}{\tau_{1}}t}+\int_{0}^{t}\left[w(u)+\tau_{5}\int_{0}^{u}\delta(u,s)g(s)\,ds\right]e^{-\frac{1}{\tau_{1}}(t-u)}\,du\\\
&=\int_{0}^{t}\left[w(u)+\tau_{5}\int_{0}^{u}\delta(u,s)g(s)\,ds\right]e^{-\frac{1}{\tau_{1}}(t-u)}\,du,\quad
t\geq 0,\end{split}$

what gives the approximation

$\begin{split}g(n)&=\sum_{i=0}^{n-1}\left[w(i)+\tau_{5}\int_{0}^{i}\delta(i,s)g(s)\,ds\right]e^{-\frac{1}{\tau_{1}}(n-i)}\\\
&=\sum_{i=0}^{n-1}\left[w(i)+\tau_{5}\sum_{j=0}^{i-1}\delta(i,j)g(j)\right]e^{-\frac{1}{\tau_{1}}(n-i)}.\end{split}$

From this identity, taking into account that $w(0)=0$ and $g(0)=0$, we get

$\begin{split}g(n)&=\sum_{i=0}^{n-1}w(i)e^{-\frac{1}{\tau_{1}}(n-i)}+\tau_{5}\sum_{i=0}^{n-1}\sum_{j=0}^{i-1}\delta(i,j)g(j)e^{-\frac{1}{\tau_{1}}(n-i)}\\\
&=\sum_{i=1}^{n-1}w(i)e^{-\frac{1}{\tau_{1}}(n-i)}+\tau_{5}\sum_{i=0}^{n-1}\sum_{j=1}^{i-1}\delta(i,j)g(j)e^{-\frac{1}{\tau_{1}}(n-i)}\\\
\end{split}$

If we understand that $\delta(i,j)\neq 0$ only if $j=i-3,i-2,i-1$, we get

$\begin{split}&\sum_{i=0}^{n-1}\sum_{j=1}^{i-1}\delta(i,j)g(j)e^{-\frac{1}{\tau_{1}}(n-i)}\\\
&=\delta(2,1)g(1)e^{-\frac{1}{\tau_{1}}(n-2)}+\sum_{j=1}^{2}\delta(3,j)g(j)e^{-\frac{1}{\tau_{1}}(n-3)}+\sum_{j=1}^{3}\delta(4,j)g(j)e^{-\frac{1}{\tau_{1}}(n-4)}\\\
&\hskip
20.00003pt+\sum_{j=1}^{4}\delta(5,j)g(j)e^{-\frac{1}{\tau_{1}}(n-5)}+\cdots+\sum_{j=1}^{n-3}\delta(n-2,j)g(j)e^{-\frac{2}{\tau_{1}}}\\\
&\hskip 20.00003pt+\sum_{j=1}^{n-2}\delta(n-1,j)g(j)e^{-\frac{1}{\tau_{1}}}\\\
&=\delta(2,1)g(1)e^{-\frac{1}{\tau_{1}}(n-2)}+\sum_{j=1}^{2}\delta(3,j)g(j)e^{-\frac{1}{\tau_{1}}(n-3)}+\sum_{j=1}^{3}\delta(4,j)g(j)e^{-\frac{1}{\tau_{1}}(n-4)}\\\
&\hskip
20.00003pt+\sum_{j=2}^{4}\delta(5,j)g(j)e^{-\frac{1}{\tau_{1}}(n-5)}+\cdots+\sum_{j=n-5}^{n-3}\delta(n-2,j)g(j)e^{-\frac{2}{\tau_{1}}}\\\
&\hskip
20.00003pt+\sum_{j=n-4}^{n-2}\delta(n-1,j)g(j)e^{-\frac{1}{\tau_{1}}}\\\
&=\left[\delta(2,1)e^{-\frac{1}{\tau_{1}}(n-2)}+\delta(3,1)e^{-\frac{1}{\tau_{1}}(n-3)}+\delta(4,1)e^{-\frac{1}{\tau_{1}}(n-4)}\right]g(1)\\\
&\hskip
20.00003pt+\left[\delta(3,2)e^{-\frac{1}{\tau_{1}}(n-3)}+\delta(4,2)e^{-\frac{1}{\tau_{1}}(n-4)}+\delta(5,2)e^{-\frac{1}{\tau_{1}}(n-5)}\right]g(2)+\cdots\\\
&\hskip
20.00003pt+\left[\delta(n-3,n-4)e^{-\frac{3}{\tau_{1}}}+\delta(n-2,n-4)e^{-\frac{2}{\tau_{1}}}+\delta(n-1,n-4)e^{-\frac{1}{\tau_{1}}}\right]g(n-4)\\\
&\hskip
20.00003pt+\left[\delta(n-2,n-3)e^{-\frac{2}{\tau_{1}}}+\delta(n-1,n-3)e^{-\frac{1}{\tau_{1}}}\right]g(n-3)\\\
&\hskip 20.00003pt+\delta(n-1,n-2)e^{-\frac{1}{\tau_{1}}}g(n-2)\\\
&=\sum_{j=1}^{n-4}\left[\sum_{i=j+1}^{j+3}\delta(i,j)e^{-\frac{1}{\tau_{1}}(n-i)}\right]g(j)+\left[\sum_{i=n-2}^{n-1}\delta(i,n-3)e^{-\frac{1}{\tau_{1}}(n-i)}\right]g(n-3)\\\
&\hskip 20.00003pt+\delta(n-1,n-2)e^{-\frac{1}{\tau_{1}}}g(n-2).\end{split}$

We can proceed by assigning fixed values in $(0,1)$ to

$\delta(i,i-1),\>\delta(i,i-2),\>\delta(i,i-3),$

that can be always of the same type, for instance,

$\delta(i,i-1)=0.5,\>\delta(i,i-2)=0.3,\>\delta(i,i-3)=0.2.$

Following for problem (6) a formulation similar to the first one exposed, we
get:

$g(t)=g(k)e^{-\frac{1}{\tau_{1}}(t-k)}+\int_{k}^{t}\left[w(u)+\tau_{5}\int_{0}^{u}\delta(u,s)g(s)\,ds\right]e^{-\frac{1}{\tau_{1}}(t-u)}\,du,$

for $t\in(k,k+1],$ $k=0,1,\ldots$, so that, we can give the approximation:

$g(k+1)=g(k)e^{-\frac{1}{\tau_{1}}}+\left[w(k)+\tau_{5}\int_{0}^{k}\delta(k,s)g(s)\,ds\right]e^{-\frac{1}{\tau_{1}}}$

and, taking again an approximation, we get the recursive formula

$g(k+1)=\left[w(k)+g(k)+\tau_{5}\sum_{j=k-3}^{k-1}\delta(k,j)g(j)\right]e^{-\frac{1}{\tau_{1}}},k=0,1,\ldots$

If we select always

$\delta(k,k-1)=0.5,\>\delta(k,k-2)=0.3,\>\delta(k,k-3)=0.2,$

the recurrence formula is given by:

$g(k+1)=\left[w(k)+g(k)+0.5\,\tau_{5}\,g(k-1)+0.3\,\tau_{5}\,g(k-2)+0.2\,\tau_{5}\,g(k-3)\right]e^{-\frac{1}{\tau_{1}}},$
(7)

for $k=0,1,\ldots$, so that, choosing $\delta$ this way, we obtain something
similar to the retarded approach, since choosing $\tau_{2}$, $\tau_{3}$ and
$\tau_{4}$ such that: $-\frac{1}{\tau_{2}}=0.5\tau_{5}$,
$-\frac{1}{\tau_{3}}=0.3\tau_{5}$ and $-\frac{1}{\tau_{4}}=0.2\tau_{5}$, the
recurrence formula (7) coincides with (5).

##  4 Conclusions

In this work, we have shown that, using the model proposed, it is possible to
obtain an excellent fit with respect to the experimental data, calculating the
positive and negative effects of the training with the ability to predict the
effect of training loads in fitness. The improvement of the accuracy in the
predicted data leads us definitely to an important tool for planning sports
coaching, and new gadgets and software that attempt to control the effect that
sports training has on the athlete individually. Besides, this model could be
adapted for weight control or to predict the risk of injury or the emotional
state of the athlete, issues which play an important role in athletic
performance.

We have also provided different models which include several delays or an
integro-differential approach.

## References

  * [1] Banister, E., Calvert, T., Savage, M., and Bach, T. (1975), ‘A systems model of training for athletic performance’, Australian Journal of Sports Medicine 7(3), 57–61.
  * [2] Banister, E., Carter, J., and Zarkadas, P. (1999), ‘Training theory and taper: validation in triathlon athletes’, European Journal of Applied Physiology and Occupational Physiology 79(2), 182–191.
  * [3] Busso, T. (2003), ‘Variable dose-response relationship between exercise training and performance’, Medicine and Science in Sports and Exercise 35(7), 1188–1195.
  * [4] Calvert, T., Banister, E., Savage, M., and Bach, T. (1976), ‘A systems model of the effects of training on physical performance’, IEEE Transactions on Systems, Man and Cybernetics SMC-6(2), 94–102.
  * [5] Candau, R., Busso, T., and Lacour, J. (1992), ‘Effects of training on iron status in cross-country skiers’, European Journal of Applied Physiology and Occupational Physiology 64(6), 497–502.
  * [6] Clarke, D., and Skiba, P. (2013), ‘Rationale and resources for teaching the mathematical modeling of athletic training and performance’, Advances in Physiology Education 37(2), 134–152.
  * [7] Coggan, A. (2006), ‘Training and racing using a power meter: an introduction’.
  * [8] Fister, I. Jr., Ljubič, K., Suganthan, P. N., Perc, M., and Fister, I. (2015), ‘Computational intelligence in sports: Challenges and opportunities within a new research domain’, Applied Mathematics and Computation 262, 178–186.
  * [9] Fitz-Clarke, J. R., Morton, R. H., and Banister, E. W. (1991) ‘Optimizing athletic performance by influence curves’, Journal of Applied Physiology 71(3), 1151–1158.
  * [10] González Badillo, J., and Gorostiaga, E. (1995), Fundamentos del entrenamiento de la fuerza. Aplicación al alto rendimiento deportivo, INDE (in Spanish).
  * [11] Hellard, P., Avalos, M., Lacoste, L., Barale, F., Chatard, J.-C., and Millet, G. (2006), ‘Assessing the limitations of the Banister model in monitoring training’, Journal of Sports Sciences 24(5), 509–520.
  * [12] Issurin, V. (2010), ‘New horizons for the methodology and physiology of training periodization’, Sports Medicine 40(3), 189–206.
  * [13] Morton, R., Fitz-Clarke, J., and Banister, E. (1990), ‘Modeling human performance in running’, Journal of Applied Physiology 69(3), 1171–1177.
  * [14] Pfeiffer, M. (2008), ‘Modeling the relationship between training and performance - A Comparison of two antagonistic concepts’, International Journal of Computer Science in Sport, 7(2), 13–32.
  * [15] Pfeiffer, M., and Hohmann, A. (2012), ‘Applications of neural networks in training science’, Hum. Mov. Sci. 31(2), 344–359.
  * [16] Rodríguez F. A., Iglesias X., Feriche B., Calderón -Soto C., Chaverri D., Wachsmuth N. B., Schmidt W., and Levine B. D. (2015), ‘Altitude Training in Elite Swimmers for Sea Level Performance (Altitude Project)’, Med. Sci. Sports Exerc. 47(9), 1965–1978.