# Conditional COT-GAN for Video Prediction with Kernel Smoothing

Tianlin Xu  Beatrice Acciaio

###### Abstract

Causal Optimal Transport (COT) results from imposing a temporal causality
constraint on classic optimal transport problems, which naturally generates a
new concept of distances between distributions on path spaces. The first
application of the COT theory for sequential learning was given in Xu et al.
[2020], where COT-GAN was introduced as an adversarial algorithm to train
implicit generative models optimized for producing sequential data. Relying on
[Xu et al., 2020], the contribution of the present paper is twofold. First, we
develop a conditional version of COT-GAN suitable for sequence prediction.
This means that the dataset is now used in order to learn how a sequence will
evolve given the observation of its past evolution. Second, we improve on the
convergence results by working with modifications of the empirical measures
via kernel smoothing due to [Pflug and Pichler, 2016]. The resulting _kernel
conditional COT-GAN_ algorithm is illustrated with an application for video
prediction.

Machine Learning, ICML

##  1 Introduction

Time series prediction is a challenging task. Given past observations, a
desirable model should not only capture the distribution of features at each
time step, but also predict its complex evolution over time. Autoregressive
models which predict one time step after another seem to be a natural choice
for learning such a task, see e.g. [Denton and Fergus, 2018, Kalchbrenner et
al., 2017, Oh et al., 2015, Weissenborn et al., 2020]. However, the drawbacks
of autoregressive models are the compounding error due to multi-step sampling
and their high computational cost, see e.g. [Kalchbrenner et al., 2017, Reed
et al., 2017]. Most existing models for time series prediction tend to ignore
the temporal dependencies in sequences in the loss function, merely relying on
certain specific network architectures, such as recurrent neural network (RNN)
and 1D and 3D convolutional neural network (CNN), to capture the underlying
dynamics, see e.g. [Srivastava et al., 2015, Aigner and Körner, 2018, Saito et
al., 2017, Vondrick et al., 2016, Tulyakov et al., 2018]. For this learning
task, the loss function used to compare prediction and real evolution plays a
crucial role. However, a loss function that is blind to the sequential nature
of data will almost certainly disappoint.

Yoon et al. [2019] proposed TimeGAN to tackle this problem by introducing an
auxiliary step-wise loss function to the original GAN objective, which indeed
leads to more coherent and accurate predictions. More recently, the advances
in the field of causal optimal transport (COT) have shown a promising
direction for sequential modeling, see e.g. [Backhoff et al., 2017, 2020,
Pflug and Pichler, 2012, Xu et al., 2020]. This type of transport constrains
the transport plans to respect temporal causality, in that the arrival
sequence at any time $t$ depends on the starting sequence only up to time $t$.
In this way, at every time we only use information available up to that time,
which is a natural request in sequential learning. This is the foundation of
COT-GAN [Xu et al., 2020], where the training objective is tailored to
sequential data. This proved to be an efficient tool, leading to generation of
high-quality video sequences. Although the sharpness of single frames remains
a challenge in video modeling, COT-GAN demonstrates that the evolution of
motions can be reproduced in a smooth manner without further regularization.

While COT-GAN is trained to produce sequences, the algorithm we propose here
is learning _conditional sequences_ , that is, how a sequence is likely to
evolve given the observation of its past evolution. For this task, we employ a
modification of the empirical measure that was introduced by Backhoff et al.
[2020] in the framework of adapted Wasserstein ($\mathcal{A}\mathcal{W}$)
distance. $\mathcal{A}\mathcal{W}$-distance is the result of an optimal
transport problem where the plans are constrained to be causal in both
direction (so-called _bicausal optimal transport_); see [Pflug and Pichler,
2012, 2016]. This turns out to be the appropriate distance to measure how much
two processes differ, when we want to give importance to the evolution of
information, see e.g. [Backhoff-Veraguas et al., 2020]. As noted in [Pflug and
Pichler, 2016] and [Backhoff et al., 2020], the
$\mathcal{A}\mathcal{W}$-distance between a distribution and the empirical
measure of a sample from it may not vanish while the size of the sample goes
to infinity. To correct for this, Pflug and Pichler [2016] proposed a
convoluted empirical measure with a scaled smoothing kernel, while Backhoff et
al. [2020] suggested an adapted empirical measure obtained by quantization -
both aiming to smooth the empirical measure in some way in order to yield a
better convergence. In this paper, we follow the approach of adapting the
empirical measure by kernel smoothing as done in [Pflug and Pichler, 2016],
and show that this smoothed empirical measure improves the performance of
conditional COT-GAN.

The process described above gives rise to _kernel conditional COT-GAN_. The
main contributions of the current paper can then be summarized as follows:

  * •

we extend the COT-GAN to a conditional framework, powered by an encoder-
decoder style generator structure;

  * •

we employ a new kernel empirical measure in the learning structure, which is a
strongly consistent estimator with respect to COT;

  * •

we show that our kernel conditional COT-GAN algorithm achieves state-of-the-
art results for video prediction.

##  2 Framework

We are given a dataset consisting of $n$ i.i.d. $d$-dimensional sequences
$(x_{1}^{i},\ldots,x^{i}_{T})_{i=1}^{n}$ where $T\in\mathbb{N}$ is the number
of time steps and $d\in\mathbb{N}$ is the dimensionality at each time. This is
thought of as a random sample from an underlying distribution $\mu$ on
$\mathbb{R}^{d\times T}$, from which we want to extract other sequences. More
precisely, we want to learn the conditional distribution of
$(x_{k+1},\ldots,x_{T})$ given $(x_{1},\ldots,x_{k})$ under $\mu$, for any
fixed $k\in\\{1,\ldots,T-1\\}$. In the application of video prediction, an
entire video contains $T$ frames, each of which has resolution $d$. The first
$k$ frames of the video are taken as an input sequence, and later frames from
time $k+1$ to $T$ are the target sequence. We will use the notation
$x_{s:t}=(x_{s},...,x_{t})$, for $1\leq s\leq t\leq T$.

The conditional learning will be done via a conditional generative adversarial
structure, based on a specific type of optimal transport tailored for
distributions on path spaces, as introduced in the next section, in the wake
of what is done in [Xu et al., 2020].

##  3 Optimal Transport and Causal Optimal Transport

Given two probability measures $\mu,\nu$ defined on $\mathbb{R}^{D}$,
$D\in\mathbb{N}$, and a cost function
$c:\mathbb{R}^{D}\times\mathbb{R}^{D}\to\mathbb{R}$, the classical
(Kantorovich) optimal transport of $\mu$ into $\nu$ is formulated as

$\mathcal{W}_{c}(\mu,\nu):=\inf_{\pi\in\Pi(\mu,\nu)}\mathbb{E}^{\pi}[c(x,y)],$
(1)

where $\Pi(\mu,\nu)$ is the set of probability measures on
$\mathbb{R}^{D}\times\mathbb{R}^{D}$ with marginals $\mu,\nu$, which are
called transport plans between $\mu$ and $\nu$. Here $c(x,y)$ is interpreted
as the cost of transporting a unit of mass from $x$ to $y$.
$\mathcal{W}_{c}(\mu,\nu)$ is thus the minimal total cost to transport the
mass $\mu$ to $\nu$. When $c(x,y)$ is a distance function between $x$ and $y$
(usually $\|x-y\|_{p}$ for some $1\leq p<\infty$), $\mathcal{W}_{c}(\mu,\nu)$
is known as Wasserstein distance or Earth mover distance.

We are interested in transports between path spaces, that is, $D=d\times T$ in
the above notations. Since now there is a time component intrinsic in the
space $\mathbb{R}^{D}$, we are adopting a particular kind of transport which
is tailored for path spaces. We denote by $x=(x_{1},...,x_{T})$ and
$y=(y_{1},...,y_{T})$ the first and second half of the coordinates on
$\mathbb{R}^{d\times T}\times\mathbb{R}^{d\times T}$, respectively. A
probability measure $\pi$ on $\mathbb{R}^{d\times T}\times\mathbb{R}^{d\times
T}$ is called _causal transport plan_ if it satisfies the constraint

$\pi(dy_{t}|dx_{1:T})=\pi(dy_{t}|dx_{1:t})\qquad\text{for all\,
$t=1,\cdots,T-1$}.$ (2)

Intuitively, the probability mass moved to the arrival sequence at time $t$
only depends on the starting sequence up to time $t$. The set of causal plans
between $\mu$ and $\nu$ is denoted by $\Pi^{\mathcal{K}}(\mu,\nu)$, and
restricting the space of transport plans in (1) to such a set gives rise to
the causal optimal transport problem:

$\mathcal{W}^{\mathcal{K}}_{c}(\mu,\nu):=\inf_{\pi\in\Pi^{\mathcal{K}}(\mu,\nu)}\mathbb{E}^{\pi}[c(x,y)].$
(3)

COT has already found wide application in dynamic problems in stochastic
calculus and mathematical finance, see e.g. [Acciaio et al., 2019b, a, 2020,
Backhoff-Veraguas et al., 2020, Backhoff et al., 2020], and first numerical
results are given in [Acciaio et al., 2020, Xu et al., 2020].

##  4 COT-GAN and CCOT-GAN

In this section we will recall the main steps that led to the COT-GAN
algorithm for sequential learning in Xu et al. [2020], and refer to Appendix A
for the details. We then introduce a conditional version, called _conditional_
COT-GAN (CCOT-GAN), suited for sequential prediction.

Solving (causal) optimal transport problems is typically computational costly
for large datasets. One way to circumvent this challenge is to resort to
approximations of transport problems by means of efficiently solvable
auxiliary problems. Notably, Genevay et al. [2018] proposed the _Sinkhorn
divergence_ , which allows for the use of the Sinkhorn algorithm [Cuturi,
2013]. The first observation is that (1) is the limit for $\varepsilon\to 0$
of the entropy-regularized transport problems

$\mathcal{P}_{c,\varepsilon}(\mu,\nu):=\inf_{\pi\in\Pi(\mu,\nu)}\\{\mathbb{E}^{\pi}[c(x,y)]-\varepsilon
H(\pi)\\},\quad\varepsilon>0,$ (4)

where $H(\pi)$ is the Shannon entropy of $\pi$. Denoting by
$\pi_{c,\varepsilon}(\mu,\nu)$ the optimizer in (4), and by
$\mathcal{W}_{c,\varepsilon}(\mu,\nu):=\mathbb{E}^{\pi_{c,\varepsilon}(\mu,\nu)}[c(x,y)]$
the resulting total cost, the Sinkhorn divergence is defined as

$\widehat{\mathcal{W}}_{c,\varepsilon}(\mu,\nu):=2\mathcal{W}_{c,\varepsilon}(\mu,\nu)-\mathcal{W}_{c,\varepsilon}(\mu,\mu)-\mathcal{W}_{c,\varepsilon}(\nu,\nu).$
(5)

Similarly, in a causal setting, we consider the entropy-regularized COT
problems

$\mathcal{P}^{\mathcal{K}}_{c,\varepsilon}(\mu,\nu):=\inf_{\pi\in\Pi^{\mathcal{K}}(\mu,\nu)}\\{\mathbb{E}^{\pi}[c(x,y)]-\varepsilon
H(\pi)\\},\quad\varepsilon>0,$ (6)

approximating (3). By using an equivalent characterization of causality (see
Appendix A), this can be reformulated as a maximization over regularized
transport problems with respect to a specific family of cost functions:

$\mathcal{P}^{\mathcal{K}}_{c,\varepsilon}(\mu,\nu)=\sup_{c^{\mathcal{K}}\in\mathcal{C}^{\mathcal{K}}(\mu,c)}\mathcal{P}_{c^{\mathcal{K}},\varepsilon}(\mu,\nu).$
(7)

The family of costs $\mathcal{C}^{\mathcal{K}}(\mu,c)$ is given by

$\displaystyle\mathcal{C}^{\mathcal{K}}(\mu,c):=$
$\displaystyle\Bigg{\\{}c(x,y)+\sum_{j=1}^{J}\sum_{t=1}^{T-1}h^{j}_{t}(y)\Delta_{t+1}M^{j}(x):$
$\displaystyle J\in\mathbb{N},(h^{j},M^{j})\in\mathcal{H}(\mu)\Bigg{\\}},$ (8)

where $\Delta_{t+1}M(x):=M_{t+1}(x_{1:t+1})-M_{t}(x_{1:t})$ and
$\mathcal{H}(\mu)$ is a set of functions depicting causality:

$\displaystyle\mathcal{H}(\mu):=$
$\displaystyle\\{(h,M):h=(h_{t})_{t=1}^{T-1},h_{t}\in\mathcal{C}_{b}(\mathbb{R}^{d\times
t}),$ $\displaystyle M=(M_{t})_{t=1}^{T}\in\mathcal{M}(\mu),\
M_{t}\in\mathcal{C}_{b}(\mathbb{R}^{d\times t})\\},$

with $\mathcal{M}(\mu)$ being the set of martingales on $\mathbb{R}^{d\times
T}$ w.r.t. the canonical filtration and the measure $\mu$, and
$\mathcal{C}_{b}(\mathbb{R}^{d\times t})$ the space of continuous, bounded
functions on $\mathbb{R}^{d\times t}$. This suggests the following as a robust
version of the Sinkhorn divergence from (5) that takes into account causality:

$\sup_{c^{\mathcal{K}}\in\mathcal{C}^{\mathcal{K}}(\mu,c)}\widehat{\mathcal{W}}_{c^{\mathcal{K}},\varepsilon}(\mu,\nu).$

This is the distance used by the discriminator in COT-GAN [Xu et al., 2020] in
order to evaluate the discrepancy between real data and generated one (up to a
slightly different definition of Sinkhorn divergence, see Appendix A), and it
is the one we will use in the current paper for sequential prediction.

Furthermore, [Xu et al., 2020] makes the two following adjustments needed to
make computations feasible. First, rather than considering the whole set of
costs in (4), in (7) we optimize over a subset
$\mathcal{C}^{\mathcal{K}}(\mu,c)$, by considering
$\mathbf{h}:=(h^{j})_{j=1}^{J}$ and $\mathbf{M}:=(M^{j})_{j=1}^{J}$ of
dimension bounded by a fixed $J\in\mathbb{N}$. Second, instead of requiring
$\mathbf{M}$ to be a martingale, we consider all continuous bounded functions
and introduce a regularization term which penalizes deviations from being a
martingale. For a mini-batch of size $m$, $\\{x^{i}_{1:T}\\}_{i=1}^{m}$,
sampled from the dataset, the martingale penalization for $\mathbf{M}$ is
defined as

${p}_{\mathbf{M}}(\widehat{\mu}):=\frac{1}{mT}\sum_{j=1}^{J}\sum_{t=1}^{T-1}\Bigg{|}\sum_{i=1}^{m}\frac{M^{j}_{t+1}(x^{i}_{1:t+1})-M^{j}_{t}(x^{i}_{1:t})}{\sqrt{\text{Var}[M^{j}]}+\eta}\Bigg{|},$

where $\widehat{\mu}$ is the empirical measure corresponding to the mini-batch
sampled from the dataset, $\text{Var}[M]$ is the empirical variance of $M$
over time and batch, and $\eta>0$ is a small constant. This leads to the
following objective function for COT-GAN in [Xu et al., 2020]:

$\widehat{\mathcal{W}}_{c^{\mathcal{K}}_{\varphi},\varepsilon}(\widehat{\mu},\widehat{\nu}_{\theta})-{\lambda}p_{{\bf
M}_{\varphi_{2}}}(\widehat{\mu}),$ (9)

where $\widehat{\nu}_{\theta}$ is the empirical measure corresponding to the
mini-batch produced by the generator, parameterized by $\theta$,
$\mathbf{h}_{\varphi_{1}}$ and $\mathbf{M}_{\varphi_{2}}$ represent the
discriminator who learns the worst-case cost $c^{\mathcal{K}}_{\varphi}$,
parameterized by $\varphi:=(\varphi_{1},\varphi_{2})$, and $\lambda$ is a
positive constant (see Appendix A for details).

We now extend the analysis developed in [Xu et al., 2020] to a conditional
framework for sequence prediction. Given the past history of a sequence up to
time step $k$, the aim of CCOT-GAN is learning to predict the evolution from
time step $k+1$ to $T$. The learning is done by stochastic gradient descent
(SGD) on mini-batches. Given a sample $\\{{x}^{i}_{1:T}\\}_{i=1}^{m}$ from the
dataset and a sample $\\{{z}^{i}_{k+1:T}\\}^{m}$ from a distribution $\zeta$
(noise) on some latent space $\mathcal{Z}$, we define the generator as a
conditional model $g_{\theta}$, parameterized by $\theta$, which predicts the
future evolution
$\hat{x}_{k+1:T}^{i}=g_{\theta}({x}^{i}_{1:k},{z}^{i}_{k+1:T})$. The
prediction $\hat{x}_{k+1:T}^{i}$ is then concatenated with the corresponding
input sequence ${x}^{i}_{1:k}$ over the time dimension in order to be compared
with the training sequence ${x}^{i}_{1:T}$ by the discriminator. We denote the
empirical distributions of real and concatenated data by

$\widehat{\mu}:=\frac{1}{m}\sum_{i=1}^{m}\delta_{{x}^{i}_{1:T}},\quad\widehat{\nu}^{c}_{\theta}:=\frac{1}{m}\sum_{i=1}^{m}\delta_{\text{concat}({x}^{i}_{1:k},\hat{x}_{k+1:T}^{i})},$

where $\widehat{\nu}^{c}_{\theta}$ incorporates the parameterization of
$g_{\theta}$ through $\\{\hat{x}_{k+1:T}^{i}\\}_{i=1}^{m}$. Following COT-
GAN’s formulation of adversarial training, we arrive at the parameterized
objective function for CCOT-GAN:

$\widehat{\mathcal{W}}_{c_{\varphi}^{\mathcal{K}},\varepsilon}(\widehat{\mu},\widehat{\nu}^{c}_{\theta})-{\lambda}p_{{\bf
M}_{\varphi_{2}}}(\widehat{\mu}).$ (10)

In the implementation of CCOT-GAN, the generator $g_{\theta}$ is broken down
into two components: an encoder that learns the features of input sequences
$\\{{x}^{i}_{1:k}\\}_{i=1}^{m}$ and a decoder that predicts future evolutions
given the features of inputs and noise $\\{{z}^{i}_{k+1:T}\\}_{i=1}^{m}$. The
discriminator role is played by $\mathbf{h}_{\varphi_{1}}$ and
$\mathbf{M}_{\varphi_{2}}$, which are parameterized separately by two neural
networks that respect temporal causality. These can take the shape of RNNs or
1D or 3D CNNs that are constrained to causal connections only, see Appendix B
for details. We maximize the objective function (10) over $\varphi$ to search
for a robust (worst-case) distance between the two empirical measures
$\widehat{\mu}$ and $\widehat{\nu}^{c}_{\theta}$, and minimize it over
$\theta$ to learn a conditional model that produces sequential prediction.

##  5 Adapted Empirical Measure and KCCOT-GAN

It was noted by Backhoff et al. [2020] and Pflug and Pichler [2016] that the
(classical) empirical measures are not necessarily consistent estimators with
respect to distances originating from transport problems where transports
plans respect causality constraints. The _nested distance_ [Pflug and Pichler,
2012] or _adapted Wasserstein_ distance [Backhoff et al., 2020] is the result
of an optimal transport problem where plans are required to satisfy the
causality constraint (2) as well as its symmetric counterpart, when inverting
the role of $x$ and $y$:

$\displaystyle\mathcal{A}\mathcal{W}_{c}(\mu,\nu):=\inf\\{\mathbb{E}^{\pi}[c(x,y)]:$
$\displaystyle\pi\in\Pi^{\mathcal{K}}(\mu,\nu),$
$\displaystyle\pi^{\prime}\in\Pi^{\mathcal{K}}(\nu,\mu)\\},$ (11)

where $\pi^{\prime}(dx,dy)=\pi(dy,dx)$.

Now, for any measure $\mu$, and for the empirical measures $\widehat{\mu}_{m}$
relative to a random sample of size $m$ from it, it is known (see e.g.
[Fournier and Guillin, 2015]) that

$\mathcal{W}_{c}(\mu,\widehat{\mu}_{m})\rightarrow 0\quad\text{as
$m\to\infty$},$

whereas [Backhoff et al., 2020, Pflug and Pichler, 2016] observe that this is
not necessarily true when substituting the Wasserstein distance
$\mathcal{W}_{c}$ with the adapted Wasserstein distance
$\mathcal{A}\mathcal{W}_{c}$. This is of course undesirable, in particular
thinking of the fact that the discriminator will evaluate discrepancies
between real and generated measures by relying on empirical measures of the
corresponding minibatches, see Section 4 and [Xu et al., 2020].

In [Backhoff et al., 2020] and [Pflug and Pichler, 2016], two different ways
of adapting the empirical measure are suggested: by smoothing using a scaled
kernel and by a quantization technique, respectively. The quantization
technique[Backhoff et al., 2020] divides the data space into sub-cubes, and
maps every value to the center of the sub-cube to which it belongs. We did not
adopt this approach for two reasons: first, the convergence property proved in
Theorem 1.3 in [Backhoff et al., 2020] only holds when the number of sub-cubes
is extremely small if the dimensionality of the data is large (typically a few
hundreds). To see why too few sub-cubes can be problematic, consider this
technique with two sub-cubes. This will map all data into only two possible
values, which discards substantial information from the original data. Second,
the quantization technique is non-differentiable, requiring an approximation
so the gradients can flow back via back-propagation in the stage of learning.
We therefore adopt the kernel smoothing approach which we describe in detail
in the remainder of this section.

For a probability measure $\mu$ with density $f$, and a density function
$k_{h}(x):=\frac{1}{h}k(\frac{x}{h})$ where $h$ is the bandwidth parameter,
the density estimator $\hat{f}$ is defined as

$\displaystyle\hat{f}(x)=\int k_{h}(x-y)f(y)dy=f*k_{h}(x),$ (12)

where $*$ denotes the convolution of densities.

Denoting the measure induced by density $k_{h}$ as $K^{f}$, we can write the
convoluted measures with density $k_{h}$ as the weighted empirical measures of
$\widehat{\mu}$ and $\widehat{\nu}^{c}_{\theta}$:

$\displaystyle\widehat{\mu}^{f}$
$\displaystyle:=\widehat{\mu}*K^{f}=\sum_{i=1}^{m}w_{i}\delta_{{x}^{i}_{1:T}},$
(13) $\displaystyle\widehat{\nu}^{c,f}_{\theta}$
$\displaystyle:=\widehat{\nu}^{c}_{\theta}*K^{f}=\sum_{i=1}^{m}w_{i}\delta_{\text{concat}({x}^{i}_{1:k},\hat{x}_{k+1:T}^{i})},$
(14)

where the weight $w_{i}$ is determined by $k_{h}$. Intuitively, this smooths
the observations by taking a weighted average of all observations, typically
with more influence from neighboring points.

Pflug and Pichler [2016] proved that the adapted Wasserstein distance of the
convoluted measures converges, i.e.,

$P(\mathcal{A}\mathcal{W}_{c}(\widehat{\mu}^{f},\widehat{\nu}^{c,f}_{\theta})>\varepsilon)\rightarrow
0\quad\text{as }m\rightarrow\infty,$

provided that

  1. 1.

the kernel $k_{h}$ is nonnegative and compactly supported on $\mathbb{R}^{D}$,

  2. 2.

the density f is bounded and uniformly continuous,

  3. 3.

the bandwidth $h$ is a function of the sample size $m$ that satisfies

$\displaystyle h_{m}\rightarrow 0,$ $\displaystyle\quad\frac{mh_{m}}{|\log
h_{m}|}\rightarrow\infty,\quad\frac{|\log h_{m}|}{\log\log
m}\rightarrow\infty,$ $\displaystyle\textit{and}\quad
mh_{m}\rightarrow\infty,\quad\text{as}\quad m\to\infty,$ (15)

  4. 4.

the measures $\mu$ and $\nu$ are conditionally Lipschitz.

For proofs and detailed discussions, please see Theorem 2 and 4 in [Pflug and
Pichler, 2016].

Note that convergence result above is derived for the adapted Wasserstein
distance $\mathcal{A}\mathcal{W}_{c}$. In order to deduce the results on
$\mathcal{W}^{\mathcal{K}}_{c}$, notice that

$\mathcal{W}^{\mathcal{K}}_{c}(\mu,\nu)\leq\mathcal{A}\mathcal{W}_{c}(\mu,\nu)$
(16)

for any probability measures $\mu,\nu$ and any cost function $c$, given that
the set of transports over which minimization is done for causal optimal
transport is bigger than that for $\mathcal{A}\mathcal{W}$-distance, cf. (3)
and (11).

Relying on this convergence result, we now introduce the CCOT-GAN with kernel
smoothing (KCCOT-GAN). The objective function of KCCOT-GAN at the level of
minibatches is computed on the adapted empirical measures:

$\widehat{\mathcal{W}}_{c_{\varphi}^{\mathcal{K}},\varepsilon}(\widehat{\mu}^{f},\widehat{\nu}^{c,f}_{\theta})-{\lambda}p_{{\bf
M}_{\varphi_{2}}}(\widehat{\mu}^{f}).$ (17)

We maximize the objective function over $\varphi$ to search for a worst-case
distance between the two adapted empirical measures, and minimize it over
$\theta$ to learn a conditional distribution that is as close as possible to
the real distribution. The algorithm is summarized in Algorithm 1. Its time
complexity scales as $\mathcal{O}((J+2d)2LTm^{2})$ in each iteration. The
distance
$\widehat{\mathcal{W}}_{c_{\varphi}^{\mathcal{K}},\varepsilon}(\widehat{\mu}^{f},\widehat{\nu}^{c,f}_{\theta})$
is approximated by the means of the Sinkhorn algorithm iteratively with a
fixed number of iterations, see Appendix A.

Algorithm 1 training KCCOT-GAN by SGD

Input: $\\{{x}^{i}_{1:T}\\}_{i=1}^{n}$(data), $\zeta$(distribution on latent
space)

Parameters: $\theta_{0}$, $\varphi_{0}$(initialization of parameters),
$m$(batch size), $\varepsilon$(regularization parameter), $\alpha$(learning
rate), $\lambda$(martingale penalty coefficient), $h$(bandwidth parameter)

repeat

(1) Sample $\\{{x}^{i}_{1:T}\\}_{i=1}^{m}$ from real data;

(2) Learn features from input sequences:
$\qquad\quad\\{{e}^{i}_{1:T}\\}_{i=1}^{m}\leftarrow
f_{\theta_{e}}(\\{{x}^{i}_{1:T}\\}_{i=1}^{m})$;

(3) Sample $\\{{z}^{i}_{k:T-1}\\}_{i=1}^{m}$ from $\zeta$;

(4) Predict conditioned on features and inputs:
$\\{\hat{x}_{k+1:T}^{i}\\}_{i=1}^{m}$ $\qquad\quad\leftarrow
f_{\theta_{d}}(\\{{e}^{i}_{1:T}\\}_{i=1}^{m},\\{{x}^{i}_{k:T-1}\\}_{i=1}^{m},\\{{z}^{i}_{k:T-1}\\}_{i=1}^{m})$;

(5) Obtain smoothed measures: $\widehat{\mu}^{f}$ and
$\widehat{\nu}^{c,f}_{\theta}$;

(6) Compute
$\widehat{\mathcal{W}}_{c_{\varphi}^{\mathcal{K}},\varepsilon}(\widehat{\mu}^{f},\widehat{\nu}_{\theta}^{c,f})$
by the Sinkhorn algorithm;

(7) Update discriminator parameter:
${\varphi}\leftarrow{\varphi}+\alpha\nabla_{\varphi}\Big{(}\widehat{\mathcal{W}}_{c_{\varphi}^{\mathcal{K}},\varepsilon}(\widehat{\mu}^{f},\widehat{\nu}_{\theta}^{c,f})-{\lambda}p_{{\bf
M}_{\varphi_{2}}}(\widehat{\mu}^{f})\Big{)}$;

(8) Repeat step (2) - (6);

(9) Update generator parameter:
$\qquad\theta\leftarrow\theta-\alpha\nabla_{\theta}\left(\widehat{\mathcal{W}}_{c_{\varphi}^{\mathcal{K}},\varepsilon}(\widehat{\mu}^{f},\widehat{\nu}_{\theta}^{c,f})\right)$;

until convergence

##  6 Implementation of KCCOT-GAN

The generator of KCCOT-GAN consists of an encoder that learns features from
the input sequences, and a decoder that generates predictions conditioned on
the input features and noise, supported by convolutional LSTM (convLSTM)[Shi
et al., 2015]. The decoder was trained using a hierarchical version of the
Teacher Forcing algorithm [Williams and Zipser, 1989] which feeds the real
values from observations as inputs during the training stage, in order to
reduce the compounding error from multi-step predictions. To make it concrete,
we proceed to formulate the implementation of KCCOT-GAN.

To avoid confusion, we refer to the entire input $x_{1:T}$ as the input
sequence, and to the sequence $x_{1:k}$ upon which the prediction $x_{k+1:T}$
is made as the context sequence. Since the full input sequence is available to
us at the stage of training, we first learn the hierarchical features of it
through an encoder with $n$ layers,

$\displaystyle e^{1}_{1:T}$ $\displaystyle=f_{\theta^{1}_{e}}(x_{1:T}),$
$\displaystyle e^{2}_{1:T}$ $\displaystyle=f_{\theta^{2}_{e}}(e^{1}_{1:T}),$
$\displaystyle\quad\vdots$ $\displaystyle e^{n}_{1:T}$
$\displaystyle=f_{\theta^{n}_{e}}(e^{n-1}_{1:T}).$

From here on, we denote the encoder as $f_{\theta_{e}}$ parametrized by
$\theta_{e}:=\\{\theta^{1}_{e},\theta^{2}_{e},...,\theta^{n}_{e}\\}$, and the
features extracted by the encoder as
$e_{1:T}:=\\{e_{1:T}^{1},...,e_{1:T}^{n}\\}$.

To deploy the teacher forcing algorithm, we make use of the hierarchical
features as well as the input sequence. At time step $k+1$, we predict
$\hat{x}_{k+1}$ conditioned on $(e_{k},x_{k})$, under the assumption that the
feature $e_{k}$ contains all the information about the context sequence.
Instead of feeding the prediction $\hat{x}_{k+1}$ back to the model to make
next prediction, we continue to predict $\hat{x}_{k+2}$ conditioned on
$(e_{k+1},x_{k+1})$ in an effort to prevent the model to derail from the truth
by making a mistake in an intermediate step. As a result, we train the model
to predict $\hat{x}_{k+1:T}$ conditioned on ($e_{k:T-1}$, $x_{k:T-1}$). In the
inference stage, however, we do not have the information beyond the context
sequence. The prediction is therefore completed in an auto-regressive manner.

Given Gaussian noise $z_{k:T-1}$, the decoder $f_{\theta_{d}}$ with $l$ layers
for $l\geq n+1$ learns to predict the future steps by

$\displaystyle d^{1}_{k+1:T}$
$\displaystyle=f_{\theta^{1}_{d}}(e^{n}_{k:T-1},z_{k:T-1}),$
$\displaystyle\quad\vdots$ $\displaystyle d^{l-1}_{k+1:T}$
$\displaystyle=f_{\theta^{l-1}_{d}}(e^{1}_{k:T-1},d^{l-2}_{k+1:T})$
$\displaystyle\hat{x}_{k+1:T}$
$\displaystyle=f_{\theta^{l}_{d}}(x_{k:T-1},d^{l-1}_{k+1:T}).$

As usual, the generator parameters $\theta:=\\{\theta_{e},\theta_{d}\\}$ and
discriminator parameters $\varphi$ are learned on the level of mini-batches
via Stochastic Gradient Descent (SGD). To yield better convergence property,
we smooth the mini-batches in each iteration using a scaled Gaussian kernel
with zero mean,

$k_{h}(x)=\frac{1}{h}e^{-\frac{x^{2}}{2h^{2}}}.$

Differently from the technique of Gaussian blur widely used in image
processing, see e.g. [Haddad et al., 1991, Reinhard et al., 2010, Nixon and
Aguado, 2019, Getreuer, 2013], we apply a 3D scaled Gaussian kernel to both
spatio and temporal dimensions. In another line of work, Zhang et al. [2020]
show that convoluting measures with a kernel density estimator is also a valid
approach to tackle the problem of disjoint supports in divergence
minimization.

The choices of the bandwidth parameter $h$ are restricted by the conditions in
Eq. (3). In the implementation, we relax this assumption by deploying a
decaying bandwidth as a function of the number of the training iterations,
rather than a function of sample size $m$. We realize that this simplification
may lead to inferior theoretical guarantee of convergence. However, we will
leave the exploration of a more appropriate approach to satisfy the
theoretical assumptions to future research.

##  7 Related Work

Video prediction is an active area of research. Methods relying on Variational
inference[Blei et al., 2017] and VAE [Kingma and Welling, 2013], e.g. SV2P
[Babaeizadeh et al., 2017], SVP-LP [Denton and Fergus, 2018], VTA [Kim et al.,
2019], and VRNN [Castrejon et al., 2019], have shown promising results. The
majority of adversarial models adopted in this domain were trained on the
original GAN objective [Goodfellow et al., 2014] or the Wasserstein GAN
objective [Arjovsky et al., 2017], both of which provide step-wise comparison
of sequences. SAVP [Lee et al., 2018] combined the objective function of the
original GAN and VAE to achieve the state of the art performance.

Substantial efforts have been devoted to designing specific architectures that
tackle the spatio-temporal dependencies, e.g. [Vondrick et al., 2016, Saito et
al., 2017, Tulyakov et al., 2018, Clark et al., 2019, Mathieu et al., 2016,
Villegas et al., 2017], and training schemes that facilitate learning, e.g.
[Mathieu et al., 2016, Villegas et al., 2017, Aigner and Körner, 2018]. Whilst
some works such as TGAN [Saito et al., 2017] and VGAN [Vondrick et al., 2016]
combined a static content generator with a motion generator, others, e.g.
[Tulyakov et al., 2018, Clark et al., 2019], designed two discriminators to
evaluate the spatial and temporal components separately. Mathieu et al. [2016]
explored a loss that measures gradient difference at frame level on top of an
adversarial loss trained with a multi-scale architecture. As a result, better
performance was achieved in comparison to a simple mean square error loss
commonly used in the literature. MCnet [Villegas et al., 2017] extended
[Mathieu et al., 2016] by adopting _convolutional long short-term memory
(ConvLSTM)_ [Shi et al., 2015] in the networks. Alternatively, 3D CNN with
progressively growing training scheme [Karras et al., 2018] was also shown to
be successful by FutureGAN [Aigner and Körner, 2018].

However, it may not be sufficient to rely solely on the network architecture
to capture the temporal structure of data. An important development in time
series synthesis and prediction is the identification of more suitable loss
functions. TimeGAN [Yoon et al., 2019] combined the original GAN loss with a
step-wise loss that computes the distance between the conditional
distributions in a supervised manner. By matching a conditional model to the
real conditional probability $p(x_{t}|x_{1:t-1})$ at every time step, it
explicitly encouraged the model to consider the temporal dependencies in the
sequence. In comparison, COT-GAN [Xu et al., 2020] explored a more natural
formulation for sequential generation which leads to convincing results.

##  8 Experiments

Figure 1: GQN Mazes results on the test set. Only the last 2 frames from the
context sequence are shown.

We compare KCCOT-GAN to CCOT-GAN without kernel smoothing as an ablation
study, to SVP-LP (Denton and Fergus [2018]), to SAVP (Lee et al. [2018]), and
to VRNN (Castrejon et al. [2019]), on three well-established video prediction
datasets. The source code and video results are available at
https://github.com/neuripss2020/kccotgan. In all our experiments, the choice
of cost function is $c(x,y)=\sum_{t}\|x_{t}-y_{t}\|_{2}^{2}$, and initial
bandwidth $h$ is 1.5 and is gradually decayed to 0.1 as training progresses.
We select the first 15 frames and downsample them to a resolution of $64\times
64$. We use the first 5 frames as the context sequence and the rest 10 frames
as the target sequence. All results are evaluated on test sets. Note that the
maximum number of hidden units used for the layers in the generator and
discriminator networks is 256 for the GQN Mazes and BAIR Push Small datasets
and 128 for the Moving MNIST dataset, due to the constraint of available
computation power. This is at most half of the baseline model sizes. Although
a compromised model capacity is expected, KCCOT-GAN still produces excellent
results on various tasks. Network architectures and more training details are
given in Appendix B.

#### GQN Mazes.

The GQN Mazes was first introduced by [Eslami et al., 2018] for training
agents to learn their surroundings by moving around. The dataset contains
random mazes generated by a game engine. A camera traverses one or two rooms
with multiple connecting corridors in each maze. The dataset comes with a
training set that contains 900 sequences and a test set with a size of 120.
The original sequences have a length of 300 and resolution of 84 $\times$ 84.

Figure 1 demonstrates that all models successfully captured the spatial
structure in the frames well. However, predictions produced by SVG-LP lack of
the evolution of motions, which is observed in many reproduced results of the
model across various dataset. This could be attributed to the fact that SVG-LP
is conditioned on a single frame from the previous time step, which makes it
impossible for the model to pick up any information about past evolution.
Visually, KCCOT-GAN and VRNN produced the sharpest frames out of all. Whilst
samples from VRNN show more variations, those from KCCOT-GAN tend to be closer
to the ground truth which may contribute to the better numerical evaluations
in Table 1.

#### BAIR Push Small.

Due to computation and storage constraint, we opted for this smaller version
of the original BAIR Push dataset. The BAIR Push Small contains about 44,000
example with a resolution of $64\times 64$. Each example shows a sequence of
motions of robot arm pushing objects on a table.

For this dataset, the results from SVG-LP and VRNN are extremely good in terms
of both the image quality and the variation in samples, see Figure 2. It is
clearly a very difficult task to outperform these two baselines. On the other
hand, SAVP has failed in producing high quality predictions.

On this dataset, although KCCOT-GAN underperforms the SVG-LP and VRNN
baselines, we observe a clear improvement in sharpness from CCOT-GAN to KCCOT-
GAN. As these two models share the same network structure and hyper-parameter
settings, we can confirm that this improvement solely comes from the adaption
of empirical measures via kernel smoothing.

Figure 2: BAIR Push Small results on the test set. Only the last 2 frames from
the context sequence are shown.

#### Moving MNIST Dataset.

Moving MINST[Srivastava et al., 2015] contains two digits that move with
velocities sampled uniformly in the range of 2 to 6 pixels per frame and
bounce within the edges of each frame. The dataset has 10000 sequences, of
which we use 8000 for training and the rest for testing. Each of the original
sequence contains 20 frames with resolution $64\times 64$. Results are given
in Table 1 and Appendix C.

#### Evaluation.

We evaluate the video predictions using three metrics: Structural Similarity
index[Wang et al., 2004] (SSIM, higher is better), Learned Perceptual Image
Patch Similarity[Zhang et al., 2018] (LPIPS, lower is better), Fréchet Video
Distance[Unterthiner et al., 2018] (FVD, lower is better).

The evaluation scores are reported in Table 1. We can see that KCCOT-GAN
outperforms the baseline models on GQN Mazes dataset based on the three
metrics. However, VRNN are well ahead other models in BAIR Push Small dataset.
The performances of VRNN and KCCOT-GAN on the Moving MMNIST dataset is
reasonably close with KCCOT-GAN leading in SSIM and LPIPS but VRNN having
better FVD score.

Table 1: Evaluations for video datasets. Lower values in the metrics indicate better sample quality for LPIPS and FVD, whereas higher values in SSIM are better. GQN Mazes | SSIM | LPIPS | FVD
---|---|---|---
SAVP | 0.49 | 0.077 | 488.35
VRNN | 0.56 | 0.062 | 345.51
SVG-LP | 0.43 | 0.094 | 575.22
CCOT-GAN | 0.60 | 0.061 | 323.28
KCCOT-GAN | 0.64 | 0.060 | 267.90
BAIR Push Small |  |  |
SAVP | 0.502 | 0.090 | 280.32
VRNN | 0.825 | 0.054 | 148.51
SVG-LP | 0.822 | 0.059 | 158.80
CCOT-GAN | 0.723 | 0.063 | 201.72
KCCOT-GAN | 0.765 | 0.060 | 167.94
Moving MMNIST |  |  |
SAVP | 0.571 | 0.123 | 129.33
VRNN | 0.770 | 0.116 | 59.14
SVG-LP | 0.668 | 0.160 | 101.39
CCOT-GAN | 0.661 | 0.139 | 74.20
KCCOT-GAN | 0.788 | 0.975 | 60.33

##  9 Discussion

In the present paper we introduce KCCOT-GAN, the first algorithm for sequence
prediction that is based on recently developed modifications of optimal
transport specifically tailored for path spaces. For this we build on the
results by Xu et al. [2020], where COT was first applied for the task of
sequential generation. Our experiments show the ability of KCCOT-GAN to not
only capture the spatial structure in the frames, but also learn the complex
dynamics evolving over time.

A limitation of the KCCOT-GAN algorithm is the restricted sample variations in
comparison to the baseline models that emphasize stochastic components in the
model design. An improvement on KCCOT-GAN could be achieved by encoding more
stochasiticity. Another direction for future work is to explore alternative
choices of the kernel function convoluted over the empirical measures as well
as a bandwidth parameter that better satisfies the conditions required for the
convergence guarantee. One may also construct a learned kernel in a similar
manner as done in MMD-GAN [Li et al., 2017], whose parameters are updated
along with those in the generator and discriminator.

## References

  * Acciaio et al. [2019a] B. Acciaio, J. Backhoff-Veraguas, and R. Carmona.  Extended mean field control problems: stochastic maximum principle and transport perspective.  _SIAM Journal on Control and Optimization_ , 57(6), 2019a.
  * Acciaio et al. [2019b] B. Acciaio, J. Backhoff-Veraguas, and A. Zalashko.  Causal optimal transport and its links to enlargement of filtrations and continuous-time stochastic optimization.  _Stochastic Processes and their Applications_ , 2019b.
  * Acciaio et al. [2020] B. Acciaio, J. Backhoff-Veraguas, and J. Jia.  Cournot-nash equilibrium and optimal transport in a dynamic setting.  _arXiv preprint arXiv:2002.08786_ , 2020.
  * Aigner and Körner [2018] S. Aigner and M. Körner.  Futuregan: Anticipating the future frames of video sequences using spatio-temporal 3d convolutions in progressively growing gans.  _arXiv preprint arXiv:1810.01325_ , 2018.
  * Arjovsky et al. [2017] M. Arjovsky, S. Chintala, and L. Bottou.  Wasserstein generative adversarial networks.  In _International conference on machine learning_ , pages 214–223. PMLR, 2017.
  * Babaeizadeh et al. [2017] M. Babaeizadeh, C. Finn, D. Erhan, R. H. Campbell, and S. Levine.  Stochastic variational video prediction.  _ICLR_ , 2017.
  * Backhoff et al. [2017] J. Backhoff, M. Beiglbock, Y. Lin, and A. Zalashko.  Causal transport in discrete time and applications.  _SIAM Journal on Optimization_ , 27(4):2528–2562, 2017.
  * Backhoff et al. [2020] J. Backhoff, D. Bartl, M. Beiglböck, and J. Wiesel.  Estimating processes in adapted Wasserstein distance.  _arXiv preprint arXiv:2002.07261_ , 2020.
  * Backhoff-Veraguas et al. [2020] J. Backhoff-Veraguas, D. Bartl, M. Beiglböck, and M. Eder.  Adapted Wasserstein distances and stability in mathematical finance.  _Finance and Stochastics_ , 24(3):601–632, 2020\.
  * Blei et al. [2017] D. M. Blei, A. Kucukelbir, and J. D. McAuliffe.  Variational inference: A review for statisticians.  _Journal of the American statistical Association_ , 112(518):859–877, 2017.
  * Castrejon et al. [2019] L. Castrejon, N. Ballas, and A. Courville.  Improved conditional vrnns for video prediction.  In _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , pages 7608–7617, 2019.
  * Clark et al. [2019] A. Clark, J. Donahue, and K. Simonyan.  Adversarial video generation on complex datasets.  _arXiv preprint arXiv:1907.06571_ , 2019.
  * Cuturi [2013] M. Cuturi.  Sinkhorn distances: Lightspeed computation of optimal transport.  In _NeurIPS_ , 2013.
  * Denton and Fergus [2018] E. Denton and R. Fergus.  Stochastic video generation with a learned prior.  In _International Conference on Machine Learning_ , pages 1174–1183. PMLR, 2018.
  * Eslami et al. [2018] S. A. Eslami, D. J. Rezende, F. Besse, F. Viola, A. S. Morcos, M. Garnelo, A. Ruderman, A. A. Rusu, I. Danihelka, K. Gregor, et al.  Neural scene representation and rendering.  _Science_ , 360(6394):1204–1210, 2018.
  * Fournier and Guillin [2015] N. Fournier and A. Guillin.  On the rate of convergence in Wasserstein distance of the empirical measure.  _Probability Theory and Related Fields_ , 162(3):707–738, 2015.
  * Genevay et al. [2018] A. Genevay, G. Peyre, and M. Cuturi.  Learning generative models with sinkhorn divergences.  In _AISTATS_ , 2018.
  * Getreuer [2013] P. Getreuer.  A survey of gaussian convolution algorithms.  _Image Processing On Line_ , 2013:286–310, 2013.
  * Goodfellow et al. [2014] I. J. Goodfellow, J. Pouget-Abadie, M. Mirza, B. Xu, D. Warde-Farley, S. Ozair, A. Courville, and Y. Bengio.  Generative adversarial networks.  _NIPS_ , 2014.
  * Haddad et al. [1991] R. A. Haddad, A. N. Akansu, et al.  A class of fast gaussian binomial filters for speech and image processing.  _IEEE Transactions on Signal Processing_ , 39(3):723–727, 1991.
  * Kalchbrenner et al. [2017] N. Kalchbrenner, A. Oord, K. Simonyan, I. Danihelka, O. Vinyals, A. Graves, and K. Kavukcuoglu.  Video pixel networks.  In _International Conference on Machine Learning_ , pages 1771–1779. PMLR, 2017.
  * Karras et al. [2018] T. Karras, T. Aila, S. Laine, and J. Lehtinen.  Progressive growing of GANs for improved quality, stability, and variation.  _ICLR_ , 2018.
  * Kim et al. [2019] T. Kim, S. Ahn, and Y. Bengio.  Variational temporal abstraction.  _NeurIPS_ , 2019.
  * Kingma and Welling [2013] D. P. Kingma and M. Welling.  Auto-encoding variational bayes.  _arXiv preprint arXiv:1312.6114_ , 2013.
  * Lee et al. [2018] A. X. Lee, R. Zhang, F. Ebert, P. Abbeel, C. Finn, and S. Levine.  Stochastic adversarial video prediction.  _arXiv preprint arXiv:1804.01523_ , 2018.
  * Li et al. [2017] C.-L. Li, W.-C. Chang, Y. Cheng, Y. Yang, and B. Póczos.  Mmd gan: Towards deeper understanding of moment matching network.  _arXiv preprint arXiv:1705.08584_ , 2017.
  * Mathieu et al. [2016] M. Mathieu, C. Couprie, and Y. LeCun.  Deep multi-scale video prediction beyond mean square error.  _ICLR_ , 2016.
  * Nixon and Aguado [2019] M. Nixon and A. Aguado.  _Feature extraction and image processing for computer vision_.  Academic press, 2019.
  * Oh et al. [2015] J. Oh, X. Guo, H. Lee, R. Lewis, and S. Singh.  Action-conditional video prediction using deep networks in atari games.  _NIPS_ , 2015.
  * Pflug and Pichler [2012] G. C. Pflug and A. Pichler.  A distance for multistage stochastic optimization models.  _SIAM Journal on Optimization_ , 22(1):1–23, 2012\.
  * Pflug and Pichler [2016] G. C. Pflug and A. Pichler.  From empirical observations to tree models for stochastic optimization: convergence properties.  _SIAM Journal on Optimization_ , 26(3):1715–1740, 2016.
  * Reed et al. [2017] S. Reed, A. Oord, N. Kalchbrenner, S. G. Colmenarejo, Z. Wang, Y. Chen, D. Belov, and N. Freitas.  Parallel multiscale autoregressive density estimation.  In _International Conference on Machine Learning_ , pages 2912–2921. PMLR, 2017.
  * Reinhard et al. [2010] E. Reinhard, W. Heidrich, P. Debevec, S. Pattanaik, G. Ward, and K. Myszkowski.  _High dynamic range imaging: acquisition, display, and image-based lighting_.  Morgan Kaufmann, 2010.
  * Saito et al. [2017] M. Saito, E. Matsumoto, and S. Saito.  Temporal generative adversarial nets with singular value clipping.  In _ICCV_ , 2017.
  * Shi et al. [2015] X. Shi, Z. Chen, H. Wang, D.-Y. Yeung, W.-K. Wong, and W.-c. Woo.  Convolutional LSTM network: A machine learning approach for precipitation nowcasting.  _arXiv preprint arXiv:1506.04214_ , 2015.
  * Srivastava et al. [2015] N. Srivastava, E. Mansimov, and R. Salakhudinov.  Unsupervised learning of video representations using LSTMs.  In _International conference on machine learning_ , pages 843–852. PMLR, 2015.
  * Tulyakov et al. [2018] S. Tulyakov, M.-Y. Liu, X. Yang, and J. Kautz.  Mocogan: Decomposing motion and content for video generation.  In _CVPR_ , 2018.
  * Unterthiner et al. [2018] T. Unterthiner, S. van Steenkiste, K. Kurach, R. Marinier, M. Michalski, and S. Gelly.  Towards accurate generative models of video: A new metric & challenges.  _arXiv preprint arXiv:1812.01717_ , 2018.
  * Villegas et al. [2017] R. Villegas, J. Yang, S. Hong, X. Lin, and H. Lee.  Decomposing motion and content for natural video sequence prediction.  _ICLR_ , 2017.
  * Vondrick et al. [2016] C. Vondrick, H. Pirsiavash, and A. Torralba.  Generating videos with scene dynamics.  In _NeurIPS_ , 2016.
  * Wang et al. [2004] Z. Wang, A. C. Bovik, H. R. Sheikh, and E. P. Simoncelli.  Image quality assessment: from error visibility to structural similarity.  _IEEE transactions on image processing_ , 13(4):600–612, 2004.
  * Weissenborn et al. [2020] D. Weissenborn, O. Täckström, and J. Uszkoreit.  Scaling autoregressive video models.  _ICLR_ , 2020.
  * Williams and Zipser [1989] R. J. Williams and D. Zipser.  A learning algorithm for continually running fully recurrent neural networks.  _Neural computation_ , 1(2):270–280, 1989.
  * Xu et al. [2020] T. Xu, L. K. Wenliang, M. Munn, and B. Acciaio.  COT-GAN: Generating Sequential Data via Causal Optimal Transport.  In _NeurIPS_ , 2020.
  * Yoon et al. [2019] J. Yoon, D. Jarrett, and M. van der Schaar.  Time-series generative adversarial networks.  In _NeurIPS_. 2019.
  * Zhang et al. [2020] M. Zhang, P. Hayes, T. Bird, R. Habib, and D. Barber.  Spread divergence.  In _International Conference on Machine Learning_ , pages 11106–11116. PMLR, 2020.
  * Zhang et al. [2018] R. Zhang, P. Isola, A. A. Efros, E. Shechtman, and O. Wang.  The unreasonable effectiveness of deep features as a perceptual metric.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pages 586–595, 2018.

Conditional COT-GAN for Video Prediction with Kernel Smoothing:
Supplementary material

##  Appendix A Details on regularized Causal Optimal Transport

###  A.1 Sinkhorn algorithm

The entropy-regularized transport problems (4) is obtained by considering an
entropic constraint. For transport plans with marginals $\mu$ supported on a
finite set $\\{x^{i}\\}_{i}$ and $\nu$ on a finite set $\\{y^{j}\\}_{j}$, any
$\pi\in\Pi(\mu,\nu)$ is also discrete with support on the set of all possible
pairs $\\{(x^{i},y^{j})\\}_{i,j}$. Denoting $\pi_{ij}=\pi(x^{i},y^{j})$, the
Shannon entropy of $\pi$ is given by
$\textstyle{H(\pi):=-\sum_{i,j}\pi_{ij}\log(\pi_{ij})}.$ A transport plan in
the discrete case can be considered as a table identified with a joint
distribution. The intuition of imposing such a regularization is to restrict
the search of couplings to tables with sufficient smoothness in order to
improve efficiency.

When the measures are discrete, such a regularized optimal transport problem
becomes easily solvable by using the Sinkhorn algorithm for a given number of
iterations, say $L$, in order to approximate a solution to the Sinkhorn
divergence (5), see [Genevay et al., 2018] for detail. Generally speaking, the
stronger the regularization is (that is, the bigger the parameter
$\varepsilon$ is), the fewer number of iterations $L$ is needed in order to
yield a good approximation.

###  A.2 Sinkhorn divergence at the level of mini-batches

To correct the fact that $\mathcal{W}_{c,\varepsilon}(\alpha,\alpha)\neq 0$,
the Sinkhorn divergence proposed by Genevay et al. [2018] at the mini-batch
level is written as

$\widehat{\mathcal{W}}_{c,\epsilon}(\widehat{\mu},\widehat{\nu}_{\theta}):=\mathcal{W}_{c,\varepsilon}(\widehat{\mu},\widehat{\nu}_{\theta})-\mathcal{W}_{c,\varepsilon}(\widehat{\mu},\widehat{\mu})-\mathcal{W}_{c,\varepsilon}(\widehat{\nu}_{\theta},\widehat{\nu}_{\theta}),$
(18)

where the empirical measures $\widehat{\mu}$ and $\widehat{\nu}_{\theta}$
correspond to mini-batch sampled from the dataset and that produced by the
model, respectively.

This is an attempt to correct the bias introduced by the entropic
regularization via eliminating the differences brought by the variations in
both mini-batches of the real and generated samples. However, an experiment in
[Xu et al., 2020] shows that the above formulation (18) failed to reduce the
bias and recover the optimizer set up as a known quantity. Therefore, the
authors propose the mixed Sinkhorn divergence,

$\displaystyle\widehat{\mathcal{W}}_{c,\epsilon}^{\text{mix}}(\widehat{\mu},\widehat{\mu}^{\prime},\widehat{\nu}_{\theta},\widehat{\nu}^{\prime}_{\theta}):=$
$\displaystyle\mathcal{W}_{c,\varepsilon}(\widehat{\mu},\widehat{\nu}_{\theta})+\mathcal{W}_{c,\varepsilon}(\widehat{\mu}^{\prime},\widehat{\nu}^{\prime}_{\theta})$
$\displaystyle-\mathcal{W}_{c,\varepsilon}(\widehat{\mu},\widehat{\mu}^{\prime})-\mathcal{W}_{c,\varepsilon}(\widehat{\nu}_{\theta},\widehat{\nu}^{\prime}_{\theta}),$

where $\widehat{\mu}$ and $\widehat{\mu}^{\prime}$ correspond to different
mini-batches from the dataset, and $\widehat{\nu}$ and
$\widehat{\nu}^{\prime}$ from generated samples. Instead of considering the
variations within a batch, the mixed Sinkhorn divergence reduces the bias by
excluding the variations in different mini-batches from the same underlying
distribution.

Alternative mini-batch Sinkhorn divergences are also investigated in [Xu et
al., 2020], for example,

$\displaystyle\widehat{\mathcal{W}}_{c,\epsilon}^{6}(\widehat{\mu},,\widehat{\mu}^{\prime},\widehat{\nu}_{\theta},\widehat{\nu}^{\prime}_{\theta})$
$\displaystyle=\mathcal{W}_{c,\varepsilon}(\widehat{\mu},\widehat{\nu}_{\theta})+\mathcal{W}_{c,\varepsilon}(\widehat{\mu}^{\prime},\widehat{\nu}_{\theta})$
$\displaystyle+\mathcal{W}_{c,\varepsilon}(\widehat{\mu},\widehat{\nu}^{\prime}_{\theta})\quad+\mathcal{W}_{c,\varepsilon}(\widehat{\mu}^{\prime},\widehat{\nu}^{\prime}_{\theta})$
$\displaystyle-2\mathcal{W}_{c,\varepsilon}(\widehat{\mu}^{\prime},\widehat{\mu}^{\prime})-2\mathcal{W}_{c,\varepsilon}(\widehat{\nu},\widehat{\nu}^{\prime}_{\theta}).$

In sequential generation (without conditioning), the results in [Xu et al.,
2020] suggest that $\widehat{\mathcal{W}}_{c,\epsilon}^{\text{mix}}$ and
$\widehat{\mathcal{W}}_{c,\epsilon}^{6}$ outperform all other formulations of
mini-batch Sinkhorn divergence in both the low-dimensional experiments and
video generation. Although $\widehat{\mathcal{W}}_{c,\epsilon}^{\text{mix}}$
and $\widehat{\mathcal{W}}_{c,\epsilon}^{6}$ produce equally good results,
$\widehat{\mathcal{W}}_{c,\epsilon}^{6}$ is computationally more expensive
because it requires two more terms in the computation.

In the case of sequential prediction,
$\widehat{\mathcal{W}}_{c,\epsilon}(\widehat{\mu},\widehat{\nu}_{\theta})$ is
employed in the KCCCOT-GAN algorithm. Recall that $\widehat{\nu}_{\theta}$
denotes the empirical measure of the concatenated sequences which share the
input sequences with the real sequences up to time step $k$. As a result, it
is not sensible to account for the variations in two batches from the same
distribution that do not coincide before time step $k$ as $\widehat{\mu}$ and
$\widehat{\nu}_{\theta}$ do. Hence, we consider
$\widehat{\mathcal{W}}_{c,\epsilon}(\widehat{\mu},\widehat{\nu}_{\theta})$ a
more appropriate objective function for prediction under the setting of KCCOT-
GAN.

###  A.3 An equivalent characterization of causality

The expression (7) obtained in Section 4 relies on the following
characterization of causality, proved in [Backhoff et al., 2017]: a transport
plan $\pi\in\Pi(\mu,\nu)$ is causal if and only if

$\textstyle{\mathbb{E}^{\pi}\left[\sum_{t=1}^{T-1}h_{t}(y)\Delta_{t+1}M(x)\right]=0\;\;\text{for
all $(h,M)\in\mathcal{H}(\mu)$}}.$ (19)

With an abuse of notation we write $h_{t}(y)$, $M_{t}(x)$, $\Delta_{t+1}M(x)$
rather than $h_{t}(y_{1:t})$, $M_{t}(x_{1:t})$, $\Delta_{t+1}M(x_{1:t+1})$.

###  A.4 Details about COT-GAN

Adopting the mixed Sinkhorn divergence, COT-GAN is trained on the following
objective function

$\widehat{\mathcal{W}}_{c,\epsilon}^{\text{mix},L}(\widehat{\mu},\widehat{\mu}^{\prime},\widehat{\nu}_{\theta},\widehat{\nu}^{\prime}_{\theta})-\lambda
p_{\mathbf{M}_{\varphi_{2}}}(\widehat{\mu}),$ (20)

where $L$ indicates the number of iterations required for approaching a
solution to the mixed Sinkhorn divergence.

To formulate an adversarial training algorithm for implicit generative models,
COT-GAN approximates the set of functions (4) by truncating the sums at a
fixed $J$, and parameterizes
$\mathbf{h}_{\varphi_{1}}:=(h_{\varphi_{1}}^{j})_{j=1}^{J}$ and
$\mathbf{M}_{\varphi_{2}}:=(M_{\varphi_{2}}^{j})_{j=1}^{J}$ as two separate
neural networks, and let $\varphi:=(\varphi_{1},\varphi_{2})$. To capture the
characteristics of those processes, the choices of network architecture are
restricted to those with causal connections only. The mixed Sinkhorn
divergence is then calculated with respect to a parameterized cost function

$c^{\mathcal{K}}_{\varphi}(x,y):=c(x,y)+\sum_{j=1}^{J}\sum_{t=1}^{T-1}h^{j}_{\varphi_{1},t}(y)\Delta_{t+1}M^{j}_{\varphi_{2}}(x),$
(21)

where the cost function is chosen to be $c(x,y)=\|x-y\|^{2}_{2}$ in COT-GAN.

While the generator $g_{\theta}:\mathcal{Z}\to\mathcal{X}$ is incorporated in
$\widehat{\nu}_{\theta}$, the discriminator role in COT-GAN is played by
$\mathbf{h}_{\varphi_{1}}$ and $\mathbf{M}_{\varphi_{2}}$. COT-GAN learns a
robust (worst-case) distance between the real data distribution and the
generated distribution by maximizing the objective (20) over $\varphi$, and a
strong generator to fool the discriminator by minimizing the mixed divergence
over $\theta$.

##  Appendix B Experiment details

###  B.1 Network architectures and training details

Table 2: Encoder and decoder architecture. | Encoder Configuration
---|---
Input |  $x_{1:T}$ with shape $T\times 64\times 64\times 3$
1 | convLSTM2D(N32, K6, S2, P=SAME), LN
2 | convLSTM2D(N64, K6, S2, P=SAME), LN
3 | convLSTM2D(N128, K5, S2, P=SAME), LN
4 | convLSTM2D(N256, K5, S2, P=SAME), LN
5 | output features $e_{1:T}$ with shape $T\times 4\times 4\times 256$
| Decoder Configuration
Input |  $z_{k:T-1}$, $e_{k:T-1}$, $x_{k:T-1}$
1 | DCONV(N256, K2, S2, P=SAME), LN
2 | convLSTM2D(N128, K4, S1, P=SAME), LN
3 | DCONV(N128, K4, S2, P=SAME), LN
4 | convLSTM2D(N64, K6, S1, P=SAME), LN
5 | DCONV(N64, K6, S2, P=SAME), LN
6 | convLSTM2D(N32, K6, S1, P=SAME), LN
4 | DCONV(N16, K6, S1, P=SAME), LN
5 | convLSTM2D(N8, K8, S1, P=SAME), LN
7 | DCONV(N3, K8, S1, P=SAME), Sigmoid

All experiments on the three datasets share the same GAN architectures. The
generator is split into an encoder and a decoder, supported by convolutional
LSTM (convLSTM). The encoder learns both the spatial and temporal features of
the input sequences, whereas the decoder predicts the future evolution
conditioned on the learned features and a latent variable.

The features from the last encoding layer has a shape of $4\times 4$ (height
$\times$ width) per time step. A latent variable $z$ is sampled from a
multivariate standard normal distribution with the same shape as the features
(same number of channels too depending on the model size). We then concatenate
the features, input sequence, and latent variables over the channel dimension
as input for the decoder. The encoder and decoder structures are detailed in
Table 2. As the discriminator, the process $\mathbf{h}$ and $\mathbf{M}$ are
parameterized with two separate networks that share the same structure, shown
in Table 3. In all tables, we use DCONV to represent a de-convolutional
(convolutional transpose) layer. The layers may have N filter size, K kernel
size, S strides and P padding option. We adopt both batch-normalization(BN)
and layer-normalization(LN), and the LeakyReLU activation function. All
hyperparameter setting are the same for all three datasets except that the
filter size is halved for the Moving MNIST dataset.

Table 3: Discriminator architecture. Discriminator | Configuration
---|---
Input | 64x64x3
0 | CONV(N32, K5, S2, P=SAME), BN
1 | CONV(N64, K5, S2, P=SAME), BN
2 | CONV(N128, K5, S2, P=SAME), BN
3 | reshape 3D array for LSTM
4 | LSTM(state size = 128), LN
5 | LSTM(state size = 64), LN
6 | LSTM(state size = 32), LN

During training, we apply exponential decay to the learning rate by
$\eta_{t}=\eta_{0}r^{s/c}$ where $\eta_{0}$ is the initial learning rate, $r$
is decay rate, $s$ is the current number of training steps and $c$ is the
decaying frequency. The bandwidth parameter $h$ are also annealed from $1.5$
to $0.1$ in a similar manner. In all experiments, the initial learning rate is
$0.0005$, decay rate $0.985$, decaying frequency $10000$, and batch size
$m=8$. The settings of hyper-parameters in the Sinkhorn algorithm are also
shared across the three datasets with $\lambda=1.0$, $\varepsilon=0.8$ and the
Sinkhorn iterations $L=100$. We train KCCOT-GAN and CCOT-GAN on a single
NVIDIA GTX 1080 Ti GPU. Each iteration takes roughly 3.5 seconds. Each
experiment is run for around 100000 iterations.

###  B.2 Results on Moving MNIST

Predictions from KCCOT-GAN conditioned on the first 5 context frames from the
test set of the Moving MNIST dataset are presented in Figure 3.

Figure 3: Moving MNIST results on test set. The first 5 frames are context
sequence and last 10 frames are predictions from KCCOT-GAN, separated by the
yellow vertical line.


A type inference algorithm finding a possible typisation of this FJ program
also solves the boolean expression $B$. Figure 18 shows an example of this.
The classes True, False and Operations always stay the same. Here we assume
that the boolean expression only consists out of $\neg\land$ (NAND) operators.
Now any boolean expression $B_{\text{in}}=v_{1}\land\neg(v_{2}\land
v_{3})\land\ldots$ can be expressed as a Java method. The example in figure 18
represents the problem $B_{\text{in}}=\neg(v_{1}\land\neg(v_{2}\land v_{3}))$.
Additionally we force the return type of the sat method to have the type True
by instancing the SATExample class, which requires the type True. When using
the FJTypeInference algorithm on the generated FJ code it will assign each
parameter of the sat method with either the type True or False. This
represents a valid assignment for the expression $B_{\text{in}}$. If
FJTypeInference fails to compute a solution the $B_{\text{in}}$ has no
possible solution. A correct solution for the sat method in figure 18 would
be:
True sat(False v1, True v2, True v3, Nand4 o1, Nand1 o2)

Any SAT problem can be transferred in polynomial time to a typeless FJ
program. Every literal $v$ in the SAT problem becomes a method parameter of
the sat method, as well as every instance of a NAND operator used.

This reduction of SAT to our type inference algorithm proofs that its
complexity is at least NP-Hard. $\square$

##  Appendix F NP-Complete Complexity Proof

We know the algorithm is NP-hard (see 6.4). To proof NP-Completeness we have
to show that it is possible to verify a solution in polynomial time. The
verification of a type solution is the FJ typecheck.

It is easy to see that the expression typing rules can be checked in
polynomial time as long as subtyping between two types is verifiable in
polynomial time.

Subtyping is also solvable in polynomial time in FGJ-GT. Assume
$\mathtt{C\texttt{<}\overline{\tt
X}\texttt{>}}\leq\mathtt{D\texttt{<}\overline{\tt Y}\texttt{>}}$ with the
number of generics $\overline{\tt X}$ and $\overline{\tt Y}$ less or equal
$n$. Also the number of classes in the subtyperelation is less or equal to
$n$. With $n$ classes the S-TRANS rule can be applied a maximum of $n$ times.
Each time the S-CLASS rules is applied which sets in the variables
$\overline{\tt X}$ into the supertype. This operations also runs in polynomial
time, so the subtyping relation is decidable in polynomial time.

This shows that the time complexity of the GFJ type check is at least
polynomial or better. $\square$


#  2-loop short-distance constraints for the $g-2$ HLbL ∗

Johan Bijnens  Department of Astronomy and Theoretical Physics, Lund
University, Sölvegatan 14 A, 223 62, Lund, Sweden<EMAIL_ADDRESS>Nils
Hermansson–Truedsson  Albert Einstein Center for Fundamental Physics,
Institute for Theoretical Physics, Universität Bern, Sidlerstrasse 5, 3012
Bern, Switzerland<EMAIL_ADDRESS>Laetitia Laub<EMAIL_ADDRESS>Antonio
Rodríguez–Sánchez  Université Paris-Saclay, CNRS/IN2P3, IJCLab, 91405 Orsay,
France<EMAIL_ADDRESS>
###### Abstract

The recent experimental measurement of the muon $g-2$ at Fermilab National
Laboratory, at a $4.2\sigma$ tension with the Standard Model prediction,
highlights the need for further improvements on the theoretical uncertainties
associated to the hadronic sector. In the framework of the operator product
expansion in the presence of a background field, the short-distance behaviour
of the hadronic light-by-light contribution was recently studied. The leading
term in this expansion is given by the massless quark-loop, which is
numerically dominant compared to non-perturbative corrections. Here, we
present the perturbative QCD correction to the massless quark-loop and
estimate its size numerically. In particular, we find that for scales above 1
GeV it is relatively small, in general roughly $-10\%$ the size of the
massless quark-loop. The knowledge of these short-distance constraints will in
the future allow to reduce the systematic uncertainties in the Standard Model
prediction of the hadronic light-by-light contribution to the $g-2$.

###### keywords:

Muon Anomalous Magnetic Moment, Short-distance Constraints, Operator Product
Expansion, Higher Order Calculations, Standard Model Tests.

##  1 Introduction

The most recently updated Standard Model (SM) prediction of the muon anomalous
magnetic moment, or, $a_{\mu}=(g-2)_{\mu}/2$, is [1]

$a_{\mu}^{\mathrm{SM}}=116\,591\,810(43)\times 10^{-11}\,.$ (1)

Earlier this year, the Fermilab National Laboratory released its first
experimental value for $a_{\mu}$ [2], which when combined with the previous
Brookhaven National Laboratory result [3] gives

$\displaystyle a_{\mu}^{\textrm{exp}}=116\,592\,061(41)\times 10^{-11}\,.$ (2)

The difference between experiment and theory is thus

$\displaystyle\Delta
a_{\mu}=a_{\mu}^{\textrm{exp}}-a_{\mu}^{\mathrm{SM}}=251(59)\times
10^{-11}\,,$ (3)

which corresponds to a $4.2\sigma$ tension. Further effort is required to
improve the precision on both theoretical calculations and experimental
measurements, this in order to finally understand whether this is a sign of
new physics beyond the SM.

On the theoretical side, the uncertainty on $a_{\mu}^{\mathrm{SM}}$ is
dominated by contributions from the hadronic sector [1]. One of the two
hadronic contributions responsible for this is the hadronic light-by-light
(HLbL), $a_{\mu}^{\mathrm{HLbL}}$. It is $a_{\mu}^{\mathrm{HLbL}}$ that we
will be concerned with here, diagrammatically represented as in Fig. 1. This
contribution can be calculated either in a dispersive data-driven approach or
on the lattice [1, 4, 5], and the work presented here is relevant for the
former method.

Figure 1:  The HLbL contribution $a_{\mu}^{\mathrm{HLbL}}$, where the blob
contains all possible hadrons.

In Fig. 1 we have denoted the three virtual loop-momenta by $q_{1,2,3}$, and
the external photon momentum by $q_{4}$. The external photon is soft, i.e. in
the static limit $q_{4}\rightarrow 0$. In the HLbL loop-integration several
momentum scales are therefore mixed. Defining the Euclidean virtualities
$Q_{i}^{2}=-q_{i}^{2}$ for $i=1,\,2,\,3$ we can single out two of the
kinematical regions of particular interest to us: $(a)$ The mixed region where
$Q_{i}^{2},\,Q_{j}^{2}\gg Q_{k}^{2},\,\Lambda_{\textrm{QCD}}^{2}$. Note that
there is no inherent ordering between $Q_{k}^{2}$ and
$\Lambda_{\textrm{QCD}}^{2}$. $(b)$ The purely short-distance (SD) region
$Q_{i}^{2}\gg\Lambda_{\textrm{QCD}}^{2}$. In these two regions one can perform
operator product expansions (OPEs) between currents in the definition of the
HLbL correlation function, giving rise to so-called SD constraints (SDCs) on
the HLbL amplitude [6, 7, 8, 9]. These are relevant to constrain model
calculations of some of the intermediate states arising in the dispersive
treatment of the HLbL amplitude, see e.g. Refs. [1, 10, 11, 12] for recent
applications. Alternatively, they may be used to complement model calculations
in the given kinematical limits. Below, we will discuss SDCs derived in the
purely SD region, $(b)$ above, after a proper definition of the HLbL tensor
has been given in Sec. 2.

##  2 The HLbL tensor

The HLbL tensor $\Pi^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}(q_{1},q_{2},q_{3})$ is
defined as the correlation function of four light-quark electromagnetic
currents $J^{\mu_{i}}(x_{i})=\bar{q}\,Q_{q}\gamma^{\mu}q$ for $q=(u,d,s)$ and
$Q_{q}=\textrm{diag}(e_{q})=\textrm{diag}(2/3,-1/3,-1/3)$ being the charge-
matrix of the quarks in question. This can be written [7]

$\displaystyle\Pi^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}(q_{1},q_{2},q_{3})=-i\int\frac{d^{4}q_{4}}{(2\pi)^{4}}\left(\prod_{i=1}^{4}\int
d^{4}x_{i}\,e^{-iq_{i}x_{i}}\right)$ $\displaystyle\times\langle
0|T\left(\prod_{j=1}^{4}J^{\mu_{j}}(x_{j})\right)|0\rangle\,.$ (4)

The tensor $\Pi^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}$ satisfies the Ward identities
$q_{i,\,\mu_{i}}\,\Pi^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}=0$, for any $q_{i}$ .
This in turn implies that one can rewrite the tensor in terms of its own
derivative

$\Pi^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}=-q_{4,\,\nu_{4}}\frac{\partial\Pi^{\mu_{1}\mu_{2}\mu_{3}\nu_{4}}}{\partial
q_{4,\,\mu_{4}}}\,.$ (5)

From this one sees that to calculate $a_{\mu}^{\mathrm{HLbL}}$ it is
sufficient to consider the derivative above. Moreover, since the external
photon is static, a Lorentz decomposition of the derivative yields only 19
independent structures. By defining 19 projectors
$P^{\tilde{\Pi}_{i}}_{\mu_{1}\mu_{2}\mu_{3}\mu_{4}\nu_{4}}$ one can then find
19 functions

$\tilde{\Pi}_{i}=P^{\tilde{\Pi}_{i}}_{\mu_{1}\mu_{2}\mu_{3}\mu_{4}\nu_{4}}\lim_{q_{4}\to
0}\frac{\partial\Pi^{\mu_{1}\mu_{2}\mu_{3}\nu_{4}}}{\partial
q_{4}^{\mu_{4}}}\,.$ (6)

The projectors can be found in Ref. [9], and are here left out for brevity.
The scalar functions $\tilde{\Pi}_{i}$ can in turn be related to six functions
$\hat{\Pi}_{1,4,7,17,39,54}$ that entirely determine
$a_{\mu}^{\mathrm{HLbL}}$. The HLbL loop-integral can be written as [1]

$\displaystyle
a_{\mu}^{\mathrm{HLbL}}=\frac{2\alpha^{3}}{3\pi^{2}}\int_{0}^{\infty}dQ_{1}\int_{0}^{\infty}dQ_{2}\int_{-1}^{1}d\tau\,\sqrt{1-\tau^{2}}\,$
$\displaystyle\times\,Q_{1}^{3}Q_{2}^{3}\sum_{i=1}^{12}T_{i}(Q_{1},Q_{2},\tau)\,\overline{\Pi}_{i}(Q_{1},Q_{2},\tau)\,,$
(7)

where $T_{i}$ are known kernels and $\overline{\Pi}_{i}$ are linear
combinations of $\hat{\Pi}_{1,4,7,17,39,54}$. The relation between the
$\tilde{\Pi}_{i}$ and $\hat{\Pi}_{j}$ is straight-forward to obtain, see Ref.
[9]. In other words, by calculating the derivative of the HLbL tensor and
projecting onto the $\tilde{\Pi}_{i}$, one can obtain the $\hat{\Pi}_{j}$ and
thus $a_{\mu}^{\mathrm{HLbL}}$ can be determined from (2).

The integration in (2) is for the full HLbL contribution, but as discussed in
the previous section one may divide the integration domain into different
regions. In the purely SD limit where
$Q_{i}^{2}\gg\Lambda_{\textrm{QCD}}^{2}$, one can use OPE techniques to
calculate the SD contribution to $a_{\mu}^{\mathrm{HLbL}}$ [7]. The onset of
this asymptotic domain is not clear, so one defines a $Q_{\textrm{min}}$ as a
variable lower cut-off in the integration (2) via
$Q_{i}^{2}>Q_{\textrm{min}}^{2}$ which should lie above the onset. In Refs.
[7, 8] SDCs for the HLbL were derived for arbitrary $Q_{\textrm{min}}$. Since
the external photon for the $g-2$ kinematics is soft, it is not possible to
simply do an OPE of the tensor in (2). Instead, as was shown in Ref. [7], one
may treat the external photon as a static background field and consider a
3-point function in the presence of this field, i.e.

$\displaystyle\Pi^{\mu_{1}\mu_{2}\mu_{3}}(q_{1},q_{2})=-\frac{1}{e}\int\frac{d^{4}q_{3}}{(2\pi)^{4}}\left(\prod_{i=1}^{3}\int
d^{4}x_{i}\,e^{-iq_{i}x_{i}}\right)$ $\displaystyle\times\langle
0|T\left(\prod_{j=1}^{3}J^{\mu_{j}}(x_{j})\right)|\gamma(q_{4})\rangle\,.$ (8)

The reason this works is that $\Pi^{\mu_{1}\mu_{2}\mu_{3}}$ can be connected
to the derivative of $\Pi^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}$ [7]. The leading
term in this OPE is the massless perturbative quark-loop (see Fig. 2(a)) which
scales as $1/Q^{2}$ [7]. Massive corrections to this as well as leading and
next-to-leading non-perturbative corrections (see Fig. 2(b)–(d)), scaling up
to $1/Q^{6}$, were calculated and numerically estimated in Refs. [7, 8]. For
$Q_{\textrm{min}}\gtrsim 1$ GeV, it was found that these corrections were
numerically suppressed with respect to the massless perturbative quark-loop by
at least one to two orders of magnitude, as reported in this conference last
year [13]. This lies beyond the current precision goal on the HLbL. Note that
2-loop gluonic perturbative corrections to the massless quark-loop still can
be sizeable. Below, we present our calculation of these contributions. Further
details can be found in Ref. [9].

Figure 2:  The previously calculated terms in the background field OPE. The
insertion of the soft external field on hard lines has been indicated as a
crossed vertex. In diagrams without a crossed vertex the external field
induces the non-perturbative condensates present.

##  3 2-loop calculation

Diagrammatic examples of the 2-loop gluonic corrections are shown in Fig. 3.
These are given from two insertions from the Dyson series into (2), i.e.

$\displaystyle\Pi^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}$
$\displaystyle=-i\int\frac{d^{4}q_{3}}{(2\pi)^{4}}\left(\prod_{i=1}^{4}\int
d^{4}x_{i}\,e^{-iq_{i}x_{i}}\right)$ $\displaystyle\times\langle
0|T\Bigg{\\{}\prod_{j=1}^{4}J^{\mu_{j}}(x_{j})\,\frac{1}{2}\prod_{i=1}^{2}\int
d^{4}z_{i}\,i\mathcal{L}_{\textrm{int}}^{\textrm{qgq}}(z_{i})\Bigg{\\}}|0\rangle\,,$
(9)

where

$\displaystyle\mathcal{L}_{\textrm{int}}^{\textrm{qgq}}(z_{i})=g_{S}\frac{\lambda^{a_{i}}_{\bar{\gamma}_{i}\bar{\delta}_{i}}}{2}B^{a_{i}}_{\nu_{i}}(z_{i})\,\bar{q}^{\bar{\gamma}_{i}}(z_{i})\gamma^{\nu_{i}}q^{\bar{\delta}_{i}}(z_{i})\,,$
(10)

for the strong coupling $g_{S}$, gluon field $B^{a_{i}}_{\nu_{i}}$ and Gell-
Mann matrices $\lambda^{a_{i}}$ with barred colour-space indices.

Figure 3:  Examples of the 2-loop topologies. The insertion of the soft
external field on the hard quark-lines has been indicated as a crossed vertex.

Taking the derivative of (3) with respect to $q_{4}$ and the static limit, one
can project onto the $\tilde{\Pi}_{i}$ as in (6). The $\tilde{\Pi}_{i}$ are
then given as linear combinations of roughly 6000 scalar 2-loop integrals with
three external scales $q_{1,2,3}^{2}$, defined in $d=4-2\epsilon$ dimensions.
To reduce this large set of integrals we use KIRA [14], and a minimal set of
six master-integrals to describe them is shown diagrammatically in Fig. 4. The
exact definitions of these are given in Ref. [9]. Note that although the
$\tilde{\Pi}_{i}$ are all separately finite without need of renormalisation,
the six master-integrals are separately divergent in $\epsilon$. These
divergences have to cancel. The needed expansions in $\epsilon$ are provided
in Refs. [15, 16]111As explained in Ref. [9], we found minor typos in the
expansions in Ref. [16]. . The finite expressions for the $\tilde{\Pi}_{i}$
are of the form

$\displaystyle\tilde{\Pi}_{m}$
$\displaystyle=f_{m,ijk}^{pqr}F_{ijk}(2)Q_{1}^{2p}Q_{2}^{2q}Q_{3}^{2r}+w_{m,ijk}^{pqr}W_{ijk}(0)Q_{1}^{2p}Q_{2}^{2q}Q_{3}^{2r}$
$\displaystyle+c_{m,ijk}^{pqr}C_{ijk}(0)Q_{1}^{2p}Q_{2}^{2q}Q_{3}^{2r}+n_{m,1}^{pqr}Q_{1}^{2p}Q_{2}^{2q}Q_{3}^{2r}\log\frac{Q_{1}^{2}}{Q_{3}^{2}}$
$\displaystyle+n_{m,2}^{pqr}Q_{1}^{2p}Q_{2}^{2q}Q_{3}^{2r}C_{ijk}(0)\log\frac{Q_{2}^{2}}{Q_{3}^{2}}$
$\displaystyle+l_{m,ijk1}^{pqr}Q_{1}^{2p}Q_{2}^{2q}Q_{3}^{2r}C_{ijk}(0)\log\frac{Q_{1}^{2}}{Q_{3}^{2}}$
$\displaystyle+l_{m,ijk2}Q_{1}^{2p}Q_{2}^{2q}Q_{3}^{2r}C_{ijk}(0)\log\frac{Q_{2}^{2}}{Q_{3}^{2}}\,,$
(11)

where $F_{ijk}(2)$, $W_{ijk}(0)$ and $C_{ijk}(0)$ are finite coefficients from
the $\epsilon$–expansions. These can used in linear combinations to yield the
$\tilde{\Pi}_{j}$ needed for $a_{\mu}^{\mathrm{HLbL}}$. In the supplementary
material of Ref. [9] we provide all of these analytical results. The numerical
evaluation of $a_{\mu}^{\mathrm{HLbL}}$ is discussed in the next subsection.

Figure 4:  The six different master integrals contributing.

##  4 Numerical evaluation

Combining the $\tilde{\Pi}_{i}$ into $\hat{\Pi}_{j}$ analytically is very
simple, but unfortunately gives rise to numerical problems in the loop-
integration (2). The reason is the appearance of negative powers of the Källén
function $\lambda=(Q_{1}^{2}+Q_{2}^{2}-Q_{3}^{2})^{2}-4Q_{1}^{2}Q_{2}^{2}$ in
the $\hat{\Pi}_{j}$, which at certain points in the integration domain are
close to zero. These spurious divergences of course have to cancel in the end
with zeros in numerators. This corresponds to points where $Q_{i}\approx
Q_{j}+Q_{k}$ for some $i$, $j$ and $k$, a region in the integration domain
which we denote the side region. Also the so-called corner regions where
$Q_{i}\ll Q_{j},Q_{k}$ at the same time as $Q_{j}\approx Q_{k}$ prove to be
problematic. We can thus divide the integration domain in the
$(Q_{1},Q_{2})$–plane as in Fig. 5, for a fixed $\Lambda=Q_{1}+Q_{2}+Q_{3}$
and $\mu=Q_{\min}$. We have performed analytic Taylor expansions in the small
parameters for each of the problematic side and corner regions, which when
implemented into a numerical integration routine give a finite
$a_{\mu}^{\mathrm{HLbL}}$. These analytical expressions are provided in the
supplementary material of Ref. [9]. We note in the passing that the corner
expansions in fact correspond to a subset of the mixed region, namely
$Q_{i}^{2},\,Q_{j}^{2}\gg Q_{k}^{2}\gg\Lambda_{\textrm{QCD}}^{2}$. Extending
the work in Refs. [6, 10] for $Q_{i}^{2},\,Q_{j}^{2}\gg
Q_{k}^{2},\,\Lambda_{\textrm{QCD}}^{2}$ by going to higher orders should
therefore reproduce our results by choosing
$Q_{k}^{2}\gg\Lambda_{\textrm{QCD}}^{2}$.

Figure 5:  The integration domain of $a_{\mu}^{\textrm{HLbL}}$ divided into
different regions.

Having taken care of the numerical issues described above we can perform the
loop-integration in (2). First of all, in Table 1 we separately present for
$Q_{\textrm{min}}=1$ GeV the values of the twelve contributions to
$a_{\mu}^{\textrm{HLbL}}$ from the sum in (2), both for the massless quark-
loop and the gluonic corrections themselves. As can be seen, the latter are,
with $\alpha_{s}\approx 0.33$ at $1$ GeV, roughly $-10\%$ of the massless
quark-loop.

| Quark-loop | Gluon-loop $\times\left(\frac{\alpha_{s}}{\pi}\right)^{-1}$
---|---|---
$\bar{\Pi}_{1}$ | $0.0084$ | $-0.0077$
$\bar{\Pi}_{2}$ | $13.28$ | $-12.30$
$\bar{\Pi}_{3}$ | $0.78$ | $-0.87$
$\bar{\Pi}_{4}$ | $-2.25$ | $0.62$
$\bar{\Pi}_{5}$ | $0.00$ | $0.20$
$\bar{\Pi}_{6}$ | $2.34$ | $-1.43$
$\bar{\Pi}_{7}$ | $-0.097$ | $0.056$
$\bar{\Pi}_{8}$ | $0.035$ | $0.41$
$\bar{\Pi}_{9}$ | $0.623$ | $-0.87$
$\bar{\Pi}_{10}$ | $1.72$ | $-1.61$
$\bar{\Pi}_{11}$ | $0.696$ | $-1.04$
$\bar{\Pi}_{12}$ | $0.165$ | $-0.16$
Total | $17.3$ | $-17.0$
Table 1: The twelve contributions to $a_{\mu}^{\textrm{HLbL}}$ from the sum in
(2) with $Q_{\mathrm{min}}=1\,\mathrm{GeV}$. Units are in $10^{-11}$.

In Fig. 6 we plot the massless quark-loop (LO) as well as the massless quark-
loop together with the gluonic corrections (NLO) as functions of
$Q_{\textrm{min}}$. The error band on the NLO correction comes from the choice
of scale $\mu$ where $\alpha_{s}(\mu)$ is defined. In order to estimate this
uncertainty, we first of all vary $\mu^{2}\in(1/2,2)Q_{\textrm{min}}^{2}$. In
addition, we take into the uncertainty the running of
$\alpha_{s}^{N_{f}=5}(M_{Z})$, by running it at five loops to both
$\alpha_{s}^{N_{f}=3}(m_{\tau})$ and $\alpha_{s}^{N_{f}=3}(Q_{\textrm{min}})$.
From this we clearly see the perturbative breakdown of QCD at low energies
below $1$ GeV. Again, we see the general $-10\%$ correction from the gluon
corrections above $1$ GeV.

Figure 6:  Contributions to $a_{\mu}^{\textrm{HLbL}}$ from the massless
perturbative quark-loop (LO), as well as the gluonic correction calculated
herein (NLO).

##  5 Conclusions

We have in a series of papers derived SDCs for the HLbL tensor in the purely
SD region where $Q_{i}^{2}\gg\Lambda_{\textrm{QCD}}^{2}$ are all large [7, 8,
9]. We have shown that the leading term in the background field OPE used to
derive these SDCs is the massless perturbative quark-loop. With non-
perturbative and massive corrections to the leading term being numerically
insignificant for the current precision goal on $a_{\mu}^{\textrm{HLbL}}$,
only the 2-loop gluonic corrections to the massless quark-loop are relevant,
which in general are $-10\%$ of the size of the leading term. In the future,
we will use OPE techniques to study SDCs for the mixed region
$Q_{i}^{2},\,Q_{j}^{2}\gg Q_{k}^{2},\,\Lambda_{\textrm{QCD}}^{2}$, which will
go beyond the current results in Refs. [6, 10].

## Acknowledgements

N. H.–T. and L. L. are funded by the Albert Einstein Center for Fundamental
Physics at Universität Bern and the Swiss National Science Foundation
respectively. J. B. is supported in part by the Swedish Research Council
grants contract numbers 2016-05996 and 2019-03779. A. R.–S. is partially
supported by the Agence Nationale de la Recherche (ANR) under grant
ANR-19-CE31-0012 (project MORA).

## References

  * [1] T. Aoyama, et al., The anomalous magnetic moment of the muon in the Standard Model, Phys. Rept. 887 (2020) 1–166.  arXiv:2006.04822, doi:10.1016/j.physrep.2020.07.006.
  * [2] B. Abi, et al., Measurement of the Positive Muon Anomalous Magnetic Moment to 0.46 ppm, Phys. Rev. Lett. 126 (14) (2021) 141801.  arXiv:2104.03281, doi:10.1103/PhysRevLett.126.141801.
  * [3] G. W. Bennett, et al., Final Report of the Muon E821 Anomalous Magnetic Moment Measurement at BNL, Phys. Rev. D73 (2006) 072003.  arXiv:hep-ex/0602035, doi:10.1103/PhysRevD.73.072003.
  * [4] M. Knecht, Present status of the muon g-2, QCD 21 (2021).
  * [5] A. Gérardin, Muon g-2 from the lattice, QCD 21 (2021).
  * [6] K. Melnikov, A. Vainshtein, Hadronic light-by-light scattering contribution to the muon anomalous magnetic moment revisited, Phys. Rev. D70 (2004) 113006.  arXiv:hep-ph/0312226, doi:10.1103/PhysRevD.70.113006.
  * [7] J. Bijnens, N. Hermansson-Truedsson, A. Rodríguez-Sánchez, Short-distance constraints for the HLbL contribution to the muon anomalous magnetic moment, Phys. Lett. B798 (2019) 134994.  arXiv:1908.03331, doi:10.1016/j.physletb.2019.134994.
  * [8] J. Bijnens, N. Hermansson-Truedsson, L. Laub, A. Rodríguez-Sánchez, Short-distance HLbL contributions to the muon anomalous magnetic moment beyond perturbation theory, JHEP 10 (2020) 203.  arXiv:2008.13487, doi:10.1007/JHEP10(2020)203.
  * [9] J. Bijnens, N. Hermansson-Truedsson, L. Laub, A. Rodríguez-Sánchez, The two-loop perturbative correction to the $(g\textendash 2)_{\mu}$ HLbL at short distances, JHEP 04 (2021) 240.  arXiv:2101.09169, doi:10.1007/JHEP04(2021)240.
  * [10] G. Colangelo, F. Hagelstein, M. Hoferichter, L. Laub, P. Stoffer, Longitudinal short-distance constraints for the hadronic light-by-light contribution to $(g-2)_{\mu}$ with large-$N_{c}$ Regge models, JHEP 03 (2020) 101.  arXiv:1910.13432, doi:10.1007/JHEP03(2020)101.
  * [11] G. Colangelo, F. Hagelstein, M. Hoferichter, L. Laub, P. Stoffer, Short-distance constraints for the longitudinal component of the hadronic light-by-light amplitude: an update (6 2021).  arXiv:2106.13222.
  * [12] K. Melnikov, A. Vainshtein, On dispersion relations and hadronic light-by-light scattering contribution to the muon anomalous magnetic moment (11 2019).  arXiv:1911.05874.
  * [13] J. Bijnens, N. Hermansson-Truedsson, L. Laub, A. Rodriguez-Sanchez, Short-distance HLbL contributions to the muon g-2, in: 23rd High-Energy Physics International Conference in Quantum Chromodynamics: QCD20 - 35 years later, 2020.  arXiv:2011.12123.
  * [14] P. Maierhöfer, J. Usovitsch, P. Uwer, Kira—A Feynman integral reduction program, Comput. Phys. Commun. 230 (2018) 99–112.  arXiv:1705.05610, doi:10.1016/j.cpc.2018.04.012.
  * [15] T. Birthwright, E. Glover, P. Marquard, Master integrals for massless two-loop vertex diagrams with three offshell legs, JHEP 09 (2004) 042.  arXiv:hep-ph/0407343, doi:10.1088/1126-6708/2004/09/042.
  * [16] F. Chavez, C. Duhr, Three-mass triangle integrals and single-valued polylogarithms, JHEP 11 (2012) 114.  arXiv:1209.2722, doi:10.1007/JHEP11(2012)114.