#  Robust VAEs via Generating Process of Noise Augmented Data

Hiroo Irobe111Equally contributed Department of Mathematical and Computing
Science, Tokyo Institute of Technology, 2-12-1 Ookayama, Meguro-ku, Tokyo,
152-8550, Japan  Wataru Aoki∗ Department of Mathematical and Computing
Science, Tokyo Institute of Technology, 2-12-1 Ookayama, Meguro-ku, Tokyo,
152-8550, Japan  Kimihiro Yamazaki  Fujitsu Limited, 4-1-1 Kamikodanaka,
Nakahara-ku, Kawasaki-shi, Kanagawa, 211-8588, Japan  Yuhui Zhang  Department
of Mathematical and Computing Science, Tokyo Institute of Technology, 2-12-1
Ookayama, Meguro-ku, Tokyo, 152-8550, Japan  Takumi Nakagawa  Department of
Mathematical and Computing Science, Tokyo Institute of Technology, 2-12-1
Ookayama, Meguro-ku, Tokyo, 152-8550, Japan  RIKEN Center for Advanced
Intelligence Project, Nihonbashi 1-chome Mitsui Building, 15th floor, 1-4-1
Nihonbashi, Chuo-ku, Tokyo, 103-0027, Japan  Hiroki Waida  Department of
Mathematical and Computing Science, Tokyo Institute of Technology, 2-12-1
Ookayama, Meguro-ku, Tokyo, 152-8550, Japan  Yuichiro Wada  Fujitsu Limited,
4-1-1 Kamikodanaka, Nakahara-ku, Kawasaki-shi, Kanagawa, 211-8588, Japan
RIKEN Center for Advanced Intelligence Project, Nihonbashi 1-chome Mitsui
Building, 15th floor, 1-4-1 Nihonbashi, Chuo-ku, Tokyo, 103-0027, Japan
Takafumi Kanamori222Corresponding author; E-mail<EMAIL_ADDRESS>Department of Mathematical and Computing Science, Tokyo Institute of
Technology, 2-12-1 Ookayama, Meguro-ku, Tokyo, 152-8550, Japan  RIKEN Center
for Advanced Intelligence Project, Nihonbashi 1-chome Mitsui Building, 15th
floor, 1-4-1 Nihonbashi, Chuo-ku, Tokyo, 103-0027, Japan

###### Abstract

Advancing defensive mechanisms against adversarial attacks in generative
models is a critical research topic in machine learning. Our study focuses on
a specific type of generative models - Variational Auto-Encoders (VAEs).
Contrary to common beliefs and existing literature which suggest that noise
injection towards training data can make models more robust, our preliminary
experiments revealed that naive usage of noise augmentation technique did not
substantially improve VAE robustness. In fact, it even degraded the quality of
learned representations, making VAEs more susceptible to adversarial
perturbations. This paper introduces a novel framework that enhances
robustness by regularizing the latent space divergence between original and
noise-augmented data. Through incorporating a paired probabilistic prior into
the standard variational lower bound, our method significantly boosts defense
against adversarial attacks. Our empirical evaluations demonstrate that this
approach, termed Robust Augmented Variational Auto-ENcoder (RAVEN), yields
superior performance in resisting adversarial inputs on widely-recognized
benchmark datasets.

##  1 Introduction

Generative models, such as Variational Auto-Encoders (VAEs) [1], generative
adversarial networks [2], normalizing flows [3], diffusion models [4], and
autoregressive models [5], have garnered significant attention in machine
learning research. These models are capable of generating novel content,
including text [6], images [1], and videos [7]. In the case of VAEs,
generation is achieved through an auto-encoder, with the training objective
being the maximization of a variational lower bound [1], hence its name.
Recently, this methodology has been applied in drug discovery, particularly
using a VAE-based model, cryoDRGN [8], for reconstructing continuous
distributions of 3D protein structures from 2D images collected by electron
microscopes.

Existing research [9] and [10] show that the original VAE [1] struggles with
reconstruction and downstream classification under adversarial inputs. Various
methods [9, 11, 10, 12, 13] have been proposed to address these issues. In
[9], it was observed that the original VAE lacked robustness against inputs
outside the empirical data distribution. Subsequently, they formulated a
modified variational lower-bound that incorporates these inputs, aiming to
construct a more robust VAE. In downstream classification tasks, their method
showed remarkable efficacy. In study [10], a Markov chain Monte Carlo sampling
approach was employed during the inference steps to counteract adversarial
inputs, by shifting the latent representation of the attacks back to more
probable regions within the posterior distribution. This technique [10]
effectively improved the model’s robustness.

Besides the above successful approaches, another strategy, presented in [14],
involves using additive random noise $\bm{\epsilon}$ on original data
$\tilde{\bm{x}}$ to fortify classifiers against adversarial inputs in
supervised learning, as supported by theoretical and experimental evidence.
Influenced by these findings, we conducted preliminary experiments on the
MNIST dataset [15], comparing the performance of the following two VAE
configurations in downstream classification tasks [9]:

(A)

a VAE trained exclusively on original dataset $\\{\tilde{\bm{x}}_{i}\\}$,

(B)

a VAE trained on both original $\\{\tilde{\bm{x}}_{i}\\}$ and corresponding
noisy dataset $\\{\bm{x}_{i}\\}$, where
$\bm{x}_{i}=\tilde{\bm{x}}_{i}+\bm{\epsilon_{i}}$, and $\bm{\epsilon_{i}}$ is
random Gaussian noise:
$\bm{\epsilon_{i}}\sim\mathcal{N}\left(\bm{0},0.05^{2}\bm{I}\right)$.

We assessed performance using the encoder’s linear classification accuracy
against adversarial inputs [9]. However, the results showed no significant
improvement in classification accuracies, suggesting that the naive usage of
the noisy data $\bm{x}$ in (B) was insufficient to enhance the original VAE’s
robustness.

Figure 1:  t-SNE visualization of latent variables by trained VAE on test
MNIST: (A) trained on original MNIST training data $\\{\tilde{\bm{x}}_{i}\\}$,
(B) trained on both original and noisy data
$\\{\tilde{\bm{x}}_{i}\\}\cup\\{\bm{x}_{i}\\}$, where
$\bm{x}_{i}=\tilde{\bm{x}}_{i}+\bm{\epsilon}_{i}$. Colors denote class labels
on the top-right.

To understand the cause of the failure, we analyzed the latent variables
$\tilde{\bm{z}}$ and $\bm{z}$, outputs of the encoder on inputs
$\tilde{\bm{x}}$ and $\bm{x}$ respectively; the t-SNE visualization [16] of
Figure 1 also shows that the representations of (B) are similar to (A). We
then discovered unexpected behavior in the VAE under condition (B): the
average squared Euclidean distance between $\tilde{\bm{z}}$ and $\bm{z}$ was
larger compared to (A), despite the inclusion of noisy data in (B)’s training
dataset; see the first and second columns in Table 1. See further details of
our preliminary experiments in Appendix C.1.

We _hypothesize_ that this failure partly stems from the distant
representation of the pair $(\tilde{\bm{z}},\bm{z})$, which hampers the
enhancement of robustness. In this study, we aim to construct a robust VAE by
bringing the variables $\tilde{\bm{z}}$ and $\bm{z}$ closer together in the
latent space. To achieve this, we integrated a probabilistic model of the
paired variables [17] into the prior of the original variational lower-bound
[1]. This model illustrates the generation mechanism of the paired variables
from original data.

Our three key contributions are summarized below:

1)

We derived a novel variational lower-bound by incorporating the generating
process for a robust VAE. The proposed bound is computationally efficient due
to its closed-form structure, and its simplicity allows seamless integration
with existing VAE extensions like $\beta$-VAE [18] and VaDE [19], as well as
methods like [10] and [13] to further enhance the robustness. We name our
method the Robust Augmented Variational auto-ENcoder (_RAVEN_).

2)

Our analysis of this bound reveals that it is a generalization of the original
bound [1], and provides an intuitive understanding of how it effectively
narrows the distance between the paired representational variables
$(\tilde{\bm{z}},\bm{z})$.

3)

We empirically prove that our hypothesis is correct by results of our
numerical experiments: RAVEN shows superior resilience against various
adversarial attacks across common datasets without compromising on
reconstruction quality. The superior performance represents a significant
advancement in the field, enhancing both the robustness and versatility of
VAEs.

Table 1:  Mean $\pm$ std on
$\\{\|\tilde{\bm{z}}_{i}-\bm{z}_{i}\|_{2}^{2}\\}_{i=1}^{10^{4}}$, where
$\tilde{\bm{z}}$ (resp. $\bm{z}$) is representation of original MNIST test
data $\tilde{\bm{x}}$ (resp. the corresponding noisy data
$\bm{x}=\tilde{\bm{x}}+\bm{\epsilon}$), obtained by an encoder. "VAE" and "VAE
(noise)" correspond to (A) and (B).

VAE | VAE (noise) | Ours
---|---|---
$0.89\pm 1.29$ | $3.24\pm 5.12$ | $0.53\pm 0.36$

##  2 Related Work

In Section 2.1, we review details of the original VAE [1], which serves as a
foundational method for our study. Subsequently, in Section 2.2, we briefly
review existing methods to build a robust VAE [12] and [13], as introduced at
the second paragraph of Section 1. At last, in Section 2.3, we describe
details of the Smooth Encoder (SE) [9], given its relative similarity to our
proposed method. The distinctions between this method and our own to ours are
clarified at the conclusion of this section.

In the following, for a random variable $\bm{z}$, let
$\mathcal{N}\left(\bm{z};\bm{\mu},\bm{\Sigma}\right)$ denote a multivariate
Gaussian distribution with mean $\bm{\mu}$ and variance matrix $\bm{\Sigma}$:

$\mathcal{N}\left(\bm{z};\bm{\mu},\bm{\Sigma}\right)=\frac{\exp\left(-\frac{1}{2}(\bm{z}-\bm{\mu})^{\top}\bm{\Sigma}^{-1}(\bm{z}-\bm{\mu})\right)}{\sqrt{\left|2\pi\bm{\Sigma}\right|}},$
(1)

and we use notation $\bm{z}\sim\mathcal{N}\left(\bm{\mu},\bm{\Sigma}\right)$.

###  2.1 Variational Auto-Encoder

The VAE stands out as a widely recognized unsupervised generative model, based
on an auto-encoder architecture. The auto-encoder is bifurcated into two
components: an encoder $\mathrm{Enc}_{\bm{\phi}}$ and the decoder
$\mathrm{Dec}_{\bm{\theta}}$, where $\bm{\phi}$ and $\bm{\theta}$ are a set of
the trainable parameters. Given a data $\bm{x}\in\mathbb{R}^{D}$, the encoder
returns
$(\bm{\mu}_{\bm{x}},\bm{\sigma}_{\bm{x}})=\mathrm{Enc}_{\bm{\phi}}(\bm{x}),\bm{\mu}_{\bm{x}}\in\mathbb{R}^{d},\bm{\sigma}_{\bm{x}}\in\mathbb{R}^{d}$
to define the latent variable $\bm{z}\in\mathbb{R}^{d}$ by
$\bm{z}=\bm{\mu}_{\bm{x}}+\bm{\sigma}_{\bm{x}}\odot\bm{\varepsilon}$, where
$\odot$ is the Hadamard product and $\bm{\varepsilon}$ is a $d$-dimensional
random vector, whose distribution is a unit Gaussian
$\bm{\varepsilon}\sim\mathcal{N}\left(\bm{0},\bm{I}\right)$. The two symbols
$\bm{0}\;\text{and}\;\bm{I}$ respectively express the zero mean vector and the
variance-covariance matrix defined by an identity matrix, respectively. By the
definition,
$\bm{z}\sim\mathcal{N}\left(\bm{\mu}_{\bm{x}},\bm{\Sigma}_{\bm{x}}\right)$,
$\bm{\Sigma}_{\bm{x}}=\mathrm{diag}\left(\bm{\sigma}_{\bm{x}}\odot\bm{\sigma}_{\bm{x}}\right)$.

The auto-encoder is trained under an objective that maximizes the log-
likelihood $\log p_{\bm{\theta}}(\bm{x})$. Since the direct maximization is
intractable, the authors of [1] proposed the variational lower-bound to
indirectly maximize $\log p_{\bm{\theta}}(\bm{x})$, where
$p_{\bm{\theta}}(\bm{x})=\int
p_{\bm{\theta}}(\bm{x}|\bm{z})p(\bm{z})\mathrm{d}\bm{z}$,
$p_{\bm{\theta}}(\bm{x}|\bm{z})$ is a likelihood defined by the decoder
$\mathrm{Dec}_{\bm{\theta}}$, and $p(\bm{z})$ is the prior distribution. The
lower-bound $\mathcal{L}\left(\bm{x};\bm{\phi},\bm{\theta}\right)$ has a form
of

$\mathcal{L}\left(\bm{x};\bm{\phi},\bm{\theta}\right):=\underbrace{\mathbb{E}_{q_{\bm{\phi}}(\bm{z}|\bm{x})}\left[\log
p_{\bm{\theta}}(\bm{x}|\bm{z})\right]}_{\text{Reconstruction
error}}-\underbrace{D_{\mathrm{KL}}\left[q_{\bm{\phi}}(\bm{z}|\bm{x})\|p(\bm{z})\right]}_{\text{Regularization
term}},$ (2)

where $q_{\bm{\phi}}(\bm{z}|\bm{x})$ is the approximated posterior
distribution defined by $\mathrm{Enc}_{\bm{\phi}}$. The symbol
$D_{\textrm{KL}}$ means the Kullback-Leibler (KL) divergence. In the original
VAE, the normal Gaussian is set to the prior, and the KL term can be expressed
as a closed form [1]. Usually, the reconstruction term is approximately
computed by a sample from the approximated posterior.

###  2.2 Existing Methods for Building Robust VAE

The researchers in [12] introduced the regularization method called _total
correlation_ to the original VAE, which serves to smoothen the embeddings in
the latent space. Additionally, by introducing a hierarchical VAE, it was able
to achieve both adversarial robustness and high reconstruction quality. The
authors of [13] suggest that VAEs do not consistently encode typical samples
generated from their own decoder. They propose an alternative construction of
the VAE with a Markov chain alternating between encoder and decoder, to
enforce their consistency. Their method is reported to perform well in
downstream classification tasks with adversarial inputs.

###  2.3 Smooth Encoder

As briefly mentioned in Section 1, the authors of [9] experimentally found
that the original VAE was not robust against an input out of the the empirical
distributional support. To solve the problem, they proposed SE. This method
assumes that the types of future attacks are known beforehand, thereby
generating adversarial data $\bm{x}$ from the original data $\tilde{\bm{x}}$
tailored to these expected attacks. The method then employs a regularization
to ensure the proximity between the representations of the original and the
adversarial data. This subsequently leads to a notable minimization of the
entropy-regularized Wasserstein distance across the data representations. In
SE’s framework, the derivation of a lower bound on the joint distribution of
$(\tilde{\bm{x}},\bm{x})$ plays a crucial role, leading to a corresponding
lower bound on the log-likelihood by integrating out $x^{\prime}$.

In contrast to SE, our method is not reliant on foreknowledge of attack types,
thus simplifying its computational process and eliminating the iterative
gradient calculations required for creating adversarial data in SE.
Additionally, our method is based on a variational lower bound directly on the
joint distribution over the augmented data pair, as we will soon discuss.

##  3 Proposed Method

Let $\bm{x}$ and $\bm{x}^{\prime}$ be two augmented data instances generated
from the original data $\tilde{\bm{x}}$. The joint distribution of this pair,
denoted as $p\left(\bm{x},\bm{x}^{\prime}\right)$, is assumed to have the
following distribution [17]:

$\begin{split}p(\bm{x},\bm{x}^{\prime})&=\mathbb{E}_{\tilde{\bm{x}}\sim
p(\tilde{\bm{x}})}\left[A(\bm{x}|\tilde{\bm{x}})A(\bm{x}^{\prime}|\tilde{\bm{x}})\right]\\\
&=\int
A(\bm{x}|\tilde{\bm{x}})A(\bm{x}^{\prime}|\tilde{\bm{x}})p(\tilde{\bm{x}})\mathrm{d}\tilde{\bm{x}},\end{split}$
(3)

where $p(\tilde{\bm{x}})$ is the original data distribution, and
$A(\cdot|\tilde{\bm{x}})$ is the augmented data distribution conditioning on
$\tilde{\bm{x}}$.

Similarly to the original VAE objective in Eq.(2), we consider a variational
lower-bound of log-likelihood of the joint distribution over
$\bar{\bm{x}}=\left(\bm{x},\bm{x}^{\prime}\right)$, as follows:

$\begin{split}\log
p_{\bm{\theta}}(\bar{\bm{x}})&\geq\mathcal{L}\left(\bar{\bm{x}};\bm{\phi},\bm{\theta}\right)\\\
&=\underbrace{\mathbb{E}_{\bar{\bm{z}}\sim
q_{\bm{\phi}}(\bar{\bm{z}}|\bar{\bm{x}})}\left[\log
p_{\bm{\theta}}(\bar{\bm{x}}|\bar{\bm{z}})\right]}_{\text{\rm\scriptsize1⃝}}\underbrace{-D_{\mathrm{KL}}\left[q_{\bm{\phi}}(\bar{\bm{z}}|\bar{\bm{x}})\|p(\bar{\bm{z}})\right]}_{\text{\rm\scriptsize2⃝}},\end{split}$
(4)

where $\bar{\bm{z}}=\left(\bm{z},\bm{z}^{\prime}\right)$ are the corresponding
latent variables to $\bar{\bm{x}}=(\bm{x},\bm{x}^{\prime})$.

Inspired by the formulation presented in Eq.(3), we reconsider the generative
process of the augmented data from the latent space perspective. We assume
that, the corresponding latent variables, $\bm{z}$ and $\bm{z}^{\prime}$,
follow some conditional distribution given the original one
$a(\cdot|\tilde{\bm{z}})$. By adopting Gaussian distributions for
$a(\cdot|\tilde{\bm{z}})$, we achieve a closed-form derivation for Eq.(4). Our
construction of a robust VAE will be based on the maximization of the novel
lower bound.

######  Definition 1.

(Variant of Example 3.8 in [17]) The prior $p(\bar{\bm{z}})$ in Eq.(4) is
given by

$p(\bar{\bm{z}})=\int
a(\bm{z}|\tilde{\bm{z}})a(\bm{z}^{\prime}|\tilde{\bm{z}})p(\tilde{\bm{z}})\mathrm{d}\tilde{\bm{z}},$
(5)

where
$p(\tilde{\bm{z}})=\mathcal{N}\left(\tilde{\bm{z}};\bm{0},\bm{I}\right)$,
$a(\bm{z}|\tilde{\bm{z}})=\mathcal{N}\left(\bm{z};\tilde{\bm{z}},\bm{\Sigma}_{\mathrm{aug}}\right)$,
and
$a(\bm{z}^{\prime}|\tilde{\bm{z}})=\mathcal{N}\left(\bm{z}^{\prime};\tilde{\bm{z}},\bm{\Sigma}_{\mathrm{aug}}\right)$.

Note that we can recover Definition 1 by setting a normal Gaussian to the
original data distribution of Example 3.8 in [17].

###  3.1 Proposed Variational Lower-Bound

######  Proposition 1.

The prior $p(\bar{\bm{z}})$ in Definition 1 has a form of

$\mathcal{N}\left(\bm{0};\bm{z}-\bm{z}^{\prime},2\bm{\Sigma}_{\mathrm{aug}}\right)\,\mathcal{N}\left(\frac{\bm{z}+\bm{z}^{\prime}}{2};\bm{0},\bm{I}+\frac{1}{2}\bm{\Sigma}_{\mathrm{aug}}\right),$
(6)

where the first and second terms are defined by Eq.(1).

The proof is given in Appendix A.2.

######  Theorem 1.

Assume the independence for $q_{\bm{\phi}}(\bar{\bm{z}}|\bar{\bm{x}})$ of
Eq.(4) as follows:
$q_{\bm{\phi}}(\bar{\bm{z}}|\bar{\bm{x}})=q_{\bm{\phi}}(\bm{z}|\bm{x})\cdot
q_{\bm{\phi}}(\bm{z}^{\prime}|\bm{x}^{\prime})$, i.e.,
$q_{\bm{\phi}}(\bar{\bm{z}}|\bar{\bm{x}})=\mathcal{N}\left(\bm{z};\bm{\mu}_{\bm{x}},\bm{\Sigma}_{\bm{x}}\right)\cdot\mathcal{N}\left(\bm{z}^{\prime};\bm{\mu}_{\bm{x}^{\prime}},\bm{\Sigma}_{\bm{x}^{\prime}}\right)$;
see the definitions of $\bm{\mu}_{\bm{x}}$ and $\bm{\Sigma}_{\bm{x}}$ in
Section 2.1. Then, \scriptsize2⃝ of Eq.(4) can be written as a closed form:

$\begin{split}&\text{\rm\scriptsize2⃝}\\!=\\!-\frac{1}{4}\Bigl{\\{}\mathrm{Tr}\left(\left(\bm{\Sigma}_{\mathrm{aug}}^{-1}+\left(2\bm{I}+\bm{\Sigma}_{\mathrm{aug}}\right)^{-1}\right)(\bm{\Sigma}_{\bm{x}}+\bm{\Sigma}_{\bm{x}^{\prime}})\right)+\text{\rm\scriptsize3⃝}\Bigr{\\}}\\\
&\;\;\;\;\;\;\;\;\;\;\;+\frac{1}{2}\Bigl{\\{}\log\left|\bm{\Sigma}_{\bm{x}}\right|+\log\left|\bm{\Sigma}_{\bm{x}^{\prime}}\right|+2d\left(1-\log
2\right)-\log\left|\bm{\Sigma}_{\mathrm{aug}}\right|-\log\left|2\bm{I}+\bm{\Sigma}_{\mathrm{aug}}\right|\Bigr{\\}},\end{split}$
(7)

where
$\text{\rm\scriptsize3⃝}=(\bm{\mu}_{\bm{x}}-\bm{\mu}_{\bm{x}^{\prime}})^{\top}\bm{\Sigma}_{\mathrm{aug}}^{-1}(\bm{\mu}_{\bm{x}}-\bm{\mu}_{\bm{x}^{\prime}})+(\bm{\mu}_{\bm{x}}+\bm{\mu}_{\bm{x}^{\prime}})^{\top}\left(2\bm{I}+\bm{\Sigma}_{\mathrm{aug}}\right)^{-1}(\bm{\mu}_{\bm{x}}+\bm{\mu}_{\bm{x}^{\prime}})$.
Moreover, assume that
$p_{\bm{\theta}}(\bar{\bm{x}}|\bar{\bm{z}})=p_{\bm{\theta}}(\bm{x}|\bm{z})\cdot
p_{\bm{\theta}}(\bm{x}^{\prime}|\bm{z}^{\prime})$. Then, based on
\scriptsize2⃝ of Eq.(7),
$\mathcal{L}\left(\bar{\bm{x}};\bm{\phi},\bm{\theta}\right)$ of Eq.(4) is
equal to

$\mathbb{E}_{\bm{z}\sim q_{\bm{\phi}}(\bm{z}|\bm{x})}\left[\log
p_{\bm{\theta}}(\bm{x}|\bm{z})\right]+\mathbb{E}_{\bm{z}^{\prime}\sim
q_{\bm{\phi}}(\bm{z}^{\prime}|\bm{x}^{\prime})}\left[\log
p_{\bm{\theta}}(\bm{x}^{\prime}|\bm{z}^{\prime})\right]+\text{\rm\scriptsize2⃝}.$
(8)

The proof sketch is given in Section 3.2. We also derived the variational
lower bound of Eq.(4) in the case that $p(\tilde{\bm{z}})$ of Definition 1 is
the Gaussian mixture model; see Appendix B.

######  Remark 1.

We remark the connection between our proposed bound in Theorem 1 and the
standard bound Eq.(2). We select a "typical" input data for $\bm{x}^{\prime}$
so that
$\bm{\mu}_{\bm{x}^{\prime}}=\bm{0},\bm{\Sigma}_{\bm{x}^{\prime}}=\bm{I}$, and
set $\bm{\Sigma}_{\mathrm{aug}}=\frac{\sqrt{5}-1}{2}\bm{I}$. Disregarding
constants that do not depend on network parameters or the input $\bm{x}$, the
KL term of Eq.(4) simplifies to
$-\frac{1}{2}\left(\bm{\mu}_{\bm{x}}^{\top}\bm{\mu}_{\bm{x}}+\mathrm{Tr}(\bm{\Sigma}_{\bm{x}})-\log\left|\bm{\Sigma}_{\bm{x}}\right|\right)$,
aligning with the KL term from the standard VAE in Eq.(2).

######  Remark 2.

Regarding \scriptsize3⃝ in Eq.(7) - the terms related to the mean vectors
$\bm{\mu}_{\bm{x}}$ and $\bm{\mu}_{\bm{x}^{\prime}}$ outputted by the encoder,
when $\bm{\Sigma}_{\mathrm{aug}}\ll\bm{I}$: $\bm{\Sigma}_{\mathrm{aug}}$ is
substantially less than $\bm{I}$ in terms of e.g., maximum eigenvalues, it
follows that
$\bm{\Sigma}_{\mathrm{aug}}^{-1}\gg(2\bm{I}+\bm{\Sigma}_{\mathrm{aug}})^{-1}$.
Consequently, the first term predominantly enforces $\bm{\mu}_{\bm{x}}$ and
$\bm{\mu}_{\bm{x}^{\prime}}$ to converge. Furthermore, as shown in the
following proposition, \scriptsize3⃝ enables the model to learn the
appropriate distance for representations from the origin.

######  Proposition 2.

The term \scriptsize3⃝ in Eq.(7) can further be derived as

$\text{\rm\scriptsize3⃝}=\bm{\mu}_{\bm{x}}^{\top}\bm{\Sigma}_{\mathrm{aug}}^{-1}\bm{\mu}_{\bm{x}}+\bm{\mu}_{\bm{x}^{\prime}}^{\top}\bm{\Sigma}_{\mathrm{aug}}^{-1}\bm{\mu}_{\bm{x}^{\prime}}\\\
-\left(\bm{\mu}_{\bm{x}}+\bm{\mu}_{\bm{x}^{\prime}}\right)^{\top}\bm{\Sigma}_{\mathrm{aug}}^{-1}\left(2^{-1}\bm{I}+\bm{\Sigma}_{\mathrm{aug}}^{-1}\right)^{-1}\bm{\Sigma}_{\mathrm{aug}}^{-1}\left(\bm{\mu}_{\bm{x}}+\bm{\mu}_{\bm{x}^{\prime}}\right).$
(9)

See the proof in Appendix A.3. In the right hand of Eq.(9), the first and
second terms force both $\bm{\mu}_{\bm{x}}$ and $\bm{\mu}_{\bm{x}^{\prime}}$
to be the zero vector, while the last term keeps these vectors away and avoid
to contract to the zero vector. Thus, the model has to learn to balance the
two.

###  3.2 Proof Sketch for Theorem 1

######  Lemma 1.

(Corollary 3.2b.1 of [20]) Consider the random variable
$\bm{x}\sim\mathcal{P}$, whose mean and the variance matrix are $\bm{\mu}$ and
$\bm{\Sigma}$, respectively. For a symmetric matrix $\bm{A}$, the following
equation holds:
$\mathbb{E}_{\bm{x}\sim\mathcal{P}}\left[\bm{x}^{\top}\bm{A}\bm{x}\right]=\mathrm{Tr}\left(\bm{A}\bm{\Sigma}\right)+\bm{\mu}^{\top}\bm{A}\bm{\mu}$.
Then, with a constant vector $\bm{a}$, the following holds:
$\mathbb{E}_{\bm{x}\sim\mathcal{N}\left(\bm{\mu},\bm{\Sigma}\right)}\left[\left(\bm{x}-\bm{a}\right)^{\top}\bm{A}\left(\bm{x}-\bm{a}\right)\right]=\mathrm{Tr}\left(\bm{A}\bm{\Sigma}\right)+\left(\bm{\mu}-\bm{a}\right)^{\top}\bm{A}\left(\bm{\mu}-\bm{a}\right)$.

Figure 2: Attack budget $\delta$ (horizontal axis) versus classification
accuracy (vertical axis) on MNIST (a, b) and Fashion-MNIST (c, d) datasets.
Our proposed method (blue) significantly outperforms baseline methods under
severe adversarial attacks, and is on par with them as the adversarial signal
vanishes. Shaded area indicates standard derivations over 5 runs.

######  Lemma 2.

For the two $d$-dimensional multivariate normal distributions
$\mathcal{N}\left(\bm{\mu}_{a},\bm{\Sigma}_{a}\right)$ and
$\mathcal{N}\left(\bm{\mu}_{b},\bm{\Sigma}_{b}\right)$, consider the cross-
entropy
$\mathbb{E}_{\bm{x}\sim\mathcal{N}\left(\bm{\mu}_{a},\bm{\Sigma}_{a}\right)}\left[-\log\mathcal{N}\left(\bm{x};\bm{\mu}_{b},\bm{\Sigma}_{b}\right)\right]$.
Then, it is written as
$\frac{1}{2}\\{\mathrm{Tr}(\bm{\Sigma}_{a}\bm{\Sigma}_{b}^{-1})+(\bm{\mu}_{a}-\bm{\mu}_{b})^{\top}\bm{\Sigma}_{b}^{-1}(\bm{\mu}_{a}-\bm{\mu}_{b})+\log\left|\bm{\Sigma}_{b}\right|+d\log
2\pi\\}$. Additionally, for the $d$-dimensional normal distribution
$\mathcal{N}\left(\bm{\mu},\bm{\Sigma}\right)$, the Shannon entropy of
$\mathcal{N}\left(\bm{\mu},\bm{\Sigma}\right)$, i.e.,
$H:=\mathbb{E}_{\bm{x}\sim\mathcal{N}\left(\bm{\mu},\bm{\Sigma}\right)}\left[-\log\mathcal{N}\left(\bm{x};\bm{\mu},\bm{\Sigma}\right)\right]$
is written as $\frac{1}{2}\left\\{d+d\log
2\pi+\log\left|\bm{\Sigma}\right|\right\\}$.

Lemma 2 can be immediately proved by using Lemma 1, as shown in Appendix A.1.

From Lemmas 1 and 2, we have

$\begin{split}&\mathbb{E}_{\bm{z}\sim\mathcal{N}\left(\bm{\mu}_{\bm{x}},\bm{\Sigma}_{\bm{x}}\right)}\left[\log\mathcal{N}\left(\bm{0};\bm{z}-\bm{z}^{\prime},2\bm{\Sigma}_{\mathrm{aug}}\right)\right]\\\
&=-\frac{1}{4}\left(\mathrm{Tr}\left(\bm{\Sigma}_{\mathrm{aug}}^{-1}\bm{\Sigma}_{\bm{x}}\right)+(\bm{\mu}_{\bm{x}}-\bm{z}^{\prime})^{\top}\bm{\Sigma}_{\mathrm{aug}}^{-1}(\bm{\mu}_{\bm{x}}-\bm{z}^{\prime})\right)-\frac{1}{2}\log\left|4\pi\bm{\Sigma}_{\mathrm{aug}}\right|.\end{split}$
(10)

Since
$q_{\bm{\phi}}(\bar{\bm{z}}|\bar{\bm{x}})=\mathcal{N}\left(\bm{z};\bm{\mu}_{\bm{x}},\bm{\Sigma}_{\bm{x}}\right)\cdot\mathcal{N}\left(\bm{z}^{\prime};\bm{\mu}_{\bm{x}^{\prime}},\bm{\Sigma}_{\bm{x}^{\prime}}\right)$,
using Eq.(10), we have

$\begin{split}&\mathbb{E}_{\bar{\bm{z}}\sim
q_{\bm{\phi}}(\bar{\bm{z}}|\bar{\bm{x}})}\left[\log\mathcal{N}\left(\bm{0};\bm{z}-\bm{z}^{\prime},2\bm{\Sigma}_{\mathrm{aug}}\right)\right]\\\
&=-\frac{1}{4}\Bigl{\\{}\mathrm{Tr}(\bm{\Sigma}_{\mathrm{aug}}^{-1}\bm{\Sigma}_{\bm{x}})+\mathbb{E}_{\bm{z}^{\prime}\sim\mathcal{N}\left(\bm{\mu}_{\bm{x}^{\prime}},\bm{\Sigma}_{\bm{x}^{\prime}}\right)}\left[(\bm{\mu}_{\bm{x}}-\bm{z}^{\prime})^{\top}\bm{\Sigma}_{\mathrm{aug}}^{-1}(\bm{\mu}_{\bm{x}}-\bm{z}^{\prime})\right]\Bigr{\\}}-\frac{1}{2}\log\left|4\pi\bm{\Sigma}_{\mathrm{aug}}\right|\\\
&=-\frac{1}{4}\mathrm{Tr}\left(\bm{\Sigma}_{\mathrm{aug}}^{-1}(\bm{\Sigma}_{\bm{x}}\\!+\\!\bm{\Sigma}_{\bm{x}^{\prime}})+(\bm{\mu}_{\bm{x}}\\!-\\!\bm{\mu}_{\bm{x}^{\prime}})^{\top}\bm{\Sigma}_{\mathrm{aug}}(\bm{\mu}_{\bm{x}}\\!-\\!\bm{\mu}_{\bm{x}^{\prime}})\right)-\frac{1}{2}\left(\log\left|\bm{\Sigma}_{\mathrm{aug}}\right|+d\log
4\pi\right).\end{split}$ (11)

Like deriving Eq.(11), we have

$\begin{split}&\mathbb{E}_{\bar{\bm{z}}\sim
q_{\bm{\phi}}(\bar{\bm{z}}|\bar{\bm{x}})}\left[\log\mathcal{N}\left(\frac{\bm{z}+\bm{z}^{\prime}}{2},\bm{0},\bm{I}+\frac{1}{2}\bm{\Sigma}_{\mathrm{aug}}\right)\right]\\\
&=-\frac{1}{4}\Bigl{\\{}\mathrm{Tr}\left(\left(\bm{\Sigma}_{\bm{x}}+\bm{\Sigma}_{\bm{x}^{\prime}}\right)(2\bm{I}+\bm{\Sigma}_{\mathrm{aug}})^{-1}\right)+\left(\bm{\mu}_{\bm{x}}+\bm{\mu}_{\bm{x}^{\prime}}\right)^{\top}(2\bm{I}+\bm{\Sigma}_{\mathrm{aug}})^{-1}(\bm{\mu}_{\bm{x}}+\bm{\mu}_{\bm{x}^{\prime}})\Bigr{\\}}\\\
&\quad\quad-\frac{1}{2}\left(\log\left|\bm{\Sigma}_{\mathrm{aug}}+2\bm{I}\right|+d\log
4\pi\right).\end{split}$ (12)

In Eq.(4), by the definition of the KL divergence, the part \scriptsize2⃝
equals to

$\underbrace{\mathbb{E}_{q_{\bm{\phi}}(\bar{\bm{z}}|\bar{\bm{x}})}\left[\log
p_{\bm{\theta}}(\bar{\bm{z}})\right]}_{\text{\rm\scriptsize4⃝}}+\underbrace{H\left(q_{\bm{\phi}}(\bar{\bm{z}}|\bar{\bm{x}})\right)}_{\text{\rm\scriptsize5⃝}},$

where $H$ means the Shannon entropy. Moreover, from Proposition 1, the part
\scriptsize4⃝ is decomposed into the summation of Eq.(11) and Eq.(12).
Finally, by applying Lemma 2 to \scriptsize5⃝, we obtain Eq.(7).

##  4 Numerical Experiments

Toward a practical evaluation of our method, we conducted experiments on MNIST
[15] and Fashion MNIST [21] (F-MNIST) datasets while keeping our focus on
model robustness.

The efficiency of our proposed model against adversarial attacks is
highlighted in those experiments. We report the corresponding details and
results in the following.

###  4.1 Experiment Setup

####  4.1.1  Adversarial attacks

An adversarial attack is some carefully crafted alternations of the input
data, designed to mislead models into yielding erroneous predictions or
classifications [22]. In our study, we focus on small perturbations on the
input data that can provoke substantial deviations in representations. Given
attackers an accessible encoder model $q_{\bm{\phi}}(\cdot|\bm{x})$, and a
specified attack budget $\|\bm{\epsilon}\|_{p}\leq\delta$, their goal is to
induce the maximum change in model representations, measured by some
dissimilarity function $\Delta$:

$\bm{\epsilon}=\arg\max_{\|\bm{\epsilon}\|_{p}\leq\delta}\Delta\left[q_{\bm{\phi}}(\cdot|\bm{x}),q_{\bm{\phi}}(\cdot|\bm{x}+\bm{\epsilon})\right].$
(13)

In Eq.(13), we adopted $p=\infty$ hence the maximum norm $\|\cdot\|_{\infty}$.
$\Delta$ takes Wasserstein distance or KL divergence.

Therefore, throughout our experiments, there are 2 types of adversarial
attacks being considered in total.

To optimize Eq.(13) w.r.t. $\bm{\epsilon}$, we apply the Projected Gradient
Descent (PGD) algorithm with 50 iterations and stepsize $\delta/25$,
maximizing either Wasserstein distance or KL divergence between the encoded
representation posteriors.

####  4.1.2 Experiment protocols

Referring the experiment protocols in [9], our experiments with adversarial
attack proceeds as follows:

As the first step, we train the VAE independently of downstream tasks until
model converges. The encoder and the decoder are symmetric 3-layer MLPs with
PReLU activations [23], without normalization layers. The model is trained
with RAdam optimizer [24] and learning rate 0.001.

We then freeze the encoder parameters, and extract representation $\bm{z}_{i}$
for each training data $\bm{x}_{i}$. The mean given by the encoder will be the
representation $\bm{z}_{i}$:

$\bm{z}_{i}=\arg\max_{\bm{z}\in\mathbb{R}^{d}}\mathcal{N}\left(\bm{z};\bm{\mu}_{\bm{x}_{i}},\bm{\Sigma}_{\bm{x}_{i}}\right)=\bm{\mu}_{\bm{x}_{i}}.$
(14)

Secondly, following the common linear evaluation method, we train a simple
linear classifier on the pairs of fixed representations
$\\{\bm{z}_{i}\\}_{i=1}^{n}$ and corresponding labels $\\{y_{i}\\}_{i=1}^{n}$
via the cross-entropy loss.

Finally, we assess model robustness by generating adversarial examples for
each test instance using Eq.(13). We report their test-set classification
accuracy under such attacks in our results.

Table 2:  Reconstruction performance results of proposed method and various
baselines. Numbers indicate mean and standard derivations over 5 runs. The
best performance is in bold font. Results with $\delta=0$ in Figure 2 are
reported in this table.

| model | Accuracy (%) $\uparrow$ | MSE $\downarrow$ | FID $\downarrow$
---|---|---|---|---
MNIST | VAE | 91.45 $\pm$ 0.26 | 862.114 $\pm$ 10.96 | 21.766 $\pm$ 0.27
VAE (noise) | 91.84 $\pm$ 0.45 | 916.383 $\pm$ 15.43 | 23.990 $\pm$ 0.86
SE | 91.44 $\pm$ 0.19 | 876.454 $\pm$ 10.82 | 22.434 $\pm$ 0.70
Ours | 92.49 $\pm$ 0.51 | 795.946 $\pm$ 16.01 | 20.636 $\pm$ 0.49
F-MNIST | VAE | 76.33 $\pm$ 0.89 | 915.934 $\pm$ 10.69 | 50.338 $\pm$ 0.85
VAE (noise) | 76.71 $\pm$ 0.58 | 953.350 $\pm$ 17.46 | 52.151 $\pm$ 0.78
SE | 77.02 $\pm$ 0.87 | 963.828 $\pm$ 19.64 | 52.896 $\pm$ 1.36
Ours | 77.20 $\pm$ 0.18 | 821.254 $\pm$ 14.84 | 45.751 $\pm$ 0.54

###  4.2 Results and Discussions

Following the protocols established in Section 4.1, we compare our method with
various baselines in practice, including vanilla and noise-augmented VAEs (for
definitions of vanilla and noise-augmented, see (A) and (B) in Section 1) and
SE [9] on MNIST and F-MNIST datasets. In our method, the augmented pair
$(\bm{x},\bm{x}^{\prime})$ in Eq.(3) is constructed as
$(\tilde{\bm{x}},\tilde{\bm{x}}+\bm{\epsilon})$, where
$\bm{\epsilon}\sim\mathcal{N}(\bm{0},0.05^{2}\bm{I})$. We set the variance
$\bm{\Sigma}_{\mathrm{aug}}$ in Definition 1 to $0.04^{2}\bm{I}$ for F-MNIST
and $0.01^{2}\bm{I}$ for MNIST. Finally, we adopt cross-entropy metric for the
first and second terms in Eq.(8) like the authors of [1]. See Appendix C.2 for
details of hype-parameter tuning. The results are shown in Figure 2 and Table
2.

As shown in the above results, for vanilla and noise-augmented VAEs, their
improvements in robustness under adversarial attacks are limited. Across
datasets and different attacks, our method consistently outperforms SE in
adversarial accuracy. Especially, the performance of SE is highly impacted by
strong adversarial attacks in F-MNIST, a relatively complex dataset, while our
method is able to sustain its performance under such circumstances, offering
robust representations.

In Table 2, we also report the reconstruction quality in pixel-wise Mean
Square Error (MSE) as in [9] and Fréchet Inception Distance (FID) [25], which
captures a deep distance between the ground-truth (clean) images and
reconstructed ones. Although SE and noise-augmented VAEs suffer a slight
decline in reconstruction fidelity—which can be attributed to their training
on altered data and subsequent testing on clean data—our method appears immune
to this trend.

Finally, we emphasize that distance of the paired representation
$(\tilde{\bm{z}},\bm{z})$ are well improved compared to VAE and VAE (noise) in
Table 1, suggesting our hypothesis introduced in Section 1 is correct.

##  5 Conclusion and Future Work

Our robust VAE, RAVEN, was proposed based on a novel variational bound via our
hypothesis. We empirically confirm that the hypothesis is correct by RAVEN’s
superior resilience against adversarial attacks. One of our future work is to
further investigate our method theoretically and experimentally.

## Acknowledgment

TK was partially supported by JSPS KAKENHI Grant Number 17H00764, 19H04071,
and 20H00576.

##  Appendix A Lemmas and Propositions

######  Lemma 3.

Let $\bm{X}$ and $\bm{Y}$ be the positive definite symmetric matrix. Then, the
following equation holds:
$\left(\bm{X}^{-1}+\bm{Y}^{-1}\right)^{-1}=\bm{X}\left(\bm{X}+\bm{Y}\right)^{-1}\bm{Y}=\bm{Y}\left(\bm{X}+\bm{Y}\right)^{-1}\bm{X}$.

###### Proof.

$\displaystyle(\text{Middle term})^{-1}$
$\displaystyle=\left(\bm{X}\left(\bm{X}+\bm{Y}\right)^{-1}\bm{Y}\right)^{-1}$
$\displaystyle=\bm{Y}^{-1}\left(\bm{X}+\bm{Y}\right)\bm{X}^{-1}$
$\displaystyle=\bm{Y}^{-1}\bm{X}\bm{X}^{-1}+\bm{Y}^{-1}\bm{Y}\bm{X}^{-1}$
$\displaystyle=\bm{Y}^{-1}+\bm{X}^{-1}$ $\displaystyle=(\text{Left term})$

∎

######  Proposition 3.

The product of the two density functions
$\mathcal{N}\left(\bm{\mu}_{a},\bm{\Sigma}_{a}\right),\mathcal{N}\left(\bm{\mu}_{b},\bm{\Sigma}_{b}\right)$
can be written as
$\mathcal{N}\left(\bm{x};\bm{\mu}_{a},\bm{\Sigma}_{a}\right)\cdot\mathcal{N}\left(\bm{x};\bm{\mu}_{b},\bm{\Sigma}_{b}\right)=\mathcal{N}\left(\bm{x};\bm{\mu}_{c},\bm{\Sigma}_{c}\right)\cdot\mathcal{N}\left(\bm{0};\bm{\mu}_{a}-\bm{\mu}_{b},\bm{\Sigma}_{a}+\bm{\Sigma}_{b}\right),$
where
$\bm{\Sigma}_{c}=\left(\bm{\Sigma}_{a}^{-1}+\bm{\Sigma}_{b}^{-1}\right)^{-1},\bm{\mu}_{c}=\bm{\Sigma}_{c}\left(\bm{\Sigma}_{a}^{-1}\bm{\mu}_{a}+\bm{\Sigma}_{b}^{-1}\bm{\mu}_{b}\right)$.

###### Proof.

Let
$P=\mathcal{N}\left(\bm{x};\bm{\mu}_{a},\bm{\Sigma}_{a}\right)\cdot\mathcal{N}\left(\bm{x};\bm{\mu}_{b},\bm{\Sigma}_{b}\right)$,
i.e.,
$P=\frac{1}{\sqrt{\left|2\pi\bm{\Sigma}_{a}\right|\left|2\pi\bm{\Sigma}_{b}\right|}}\exp\left(-\frac{1}{2}A\right),$
where
$A:=(\bm{x}-\bm{\mu}_{a})^{\top}\bm{\Sigma}_{a}^{-1}(\bm{x}-\bm{\mu}_{a})+(\bm{x}-\bm{\mu}_{b})^{\top}\bm{\Sigma}_{b}^{-1}(\bm{x}-\bm{\mu}_{b}).$
Using Lemma 3,

$\displaystyle A$
$\displaystyle=\bm{x}^{\top}\underbrace{\left(\bm{\Sigma}_{a}^{-1}+\bm{\Sigma}_{b}^{-1}\right)}_{=\bm{\Sigma}_{c}^{-1}}\bm{x}-2\left(\bm{\mu}_{a}^{\top}\bm{\Sigma}_{a}^{-1}+\bm{\mu}_{b}^{\top}\bm{\Sigma}_{b}^{-1}\right)\bm{x}+\left(\bm{\mu}_{a}^{\top}\bm{\Sigma}_{a}^{-1}\bm{\mu}_{a}+\bm{\mu}_{b}^{\top}\bm{\Sigma}_{b}^{-1}\bm{\mu}_{b}\right)$
$\displaystyle=\bm{x}^{\top}\bm{\Sigma}_{c}^{-1}\bm{x}-2\underbrace{\left(\bm{\mu}_{a}^{\top}\bm{\Sigma}_{a}^{-1}+\bm{\mu}_{b}^{\top}\bm{\Sigma}_{b}^{-1}\right)\bm{\Sigma}_{c}}_{=\bm{\mu}_{c}^{\top}}\bm{\Sigma}_{c}^{-1}\bm{x}+\left(\bm{\mu}_{a}^{\top}\bm{\Sigma}_{a}^{-1}\bm{\mu}_{a}+\bm{\mu}_{b}^{\top}\bm{\Sigma}_{b}^{-1}\bm{\mu}_{b}\right)$
$\displaystyle=\bm{x}^{\top}\bm{\Sigma}_{c}^{-1}\bm{x}-2\bm{\mu}_{c}^{\top}\bm{\Sigma}_{c}^{-1}\bm{x}+\left(\bm{\mu}_{a}^{\top}\bm{\Sigma}_{a}^{-1}\bm{\mu}_{a}+\bm{\mu}_{b}^{\top}\bm{\Sigma}_{b}^{-1}\bm{\mu}_{b}\right)$
$\displaystyle=\underbrace{(\bm{x}-\bm{\mu}_{c})^{\top}\bm{\Sigma}_{c}^{-1}(\bm{x}-\bm{\mu}_{c})}_{\text{Exponential
part
of\;}\mathcal{N}\left(\bm{x};\bm{\mu}_{c},\bm{\Sigma}_{c}\right)}+\underbrace{\bigl{(}\bm{\mu}_{a}^{\top}\bm{\Sigma}_{a}^{-1}\bm{\mu}_{a}+\bm{\mu}_{b}^{\top}\bm{\Sigma}_{b}^{-1}\bm{\mu}_{b}-\bm{\mu}_{c}^{\top}\bm{\Sigma}_{c}^{-1}\bm{\mu}_{c}\bigr{)}}_{:=B}.$

Thus,

$\displaystyle P$
$\displaystyle=\sqrt{\frac{\left|2\pi\bm{\Sigma}_{c}\right|}{\left|2\pi\bm{\Sigma}_{a}\right|\left|2\pi\bm{\Sigma}_{b}\right|}}\cdot\exp\left(-\frac{1}{2}B\right)\times\underbrace{\frac{1}{\sqrt{\left|2\pi\bm{\Sigma}_{c}\right|}}\exp\left(-\frac{1}{2}(\bm{x}-\bm{\mu}_{c})^{\top}\bm{\Sigma}_{c}^{-1}(\bm{x}-\bm{\mu}_{c})\right)}_{=\mathcal{N}\left(\bm{x};\bm{\mu}_{c},\bm{\Sigma}_{c}\right)}$
$\displaystyle=\sqrt{\frac{\left|2\pi\bm{\Sigma}_{c}\right|}{\left|2\pi\bm{\Sigma}_{a}\right|\left|2\pi\bm{\Sigma}_{b}\right|}}\exp\left(-\frac{1}{2}B\right)\cdot\mathcal{N}\left(\bm{x};\bm{\mu}_{c},\bm{\Sigma}_{c}\right).$

Moreover,

$\displaystyle B$
$\displaystyle=\bm{\mu}_{a}^{\top}\bm{\Sigma}_{a}^{-1}\bm{\mu}_{a}+\bm{\mu}_{b}^{\top}\bm{\Sigma}_{b}^{-1}\bm{\mu}_{b}-\bm{\mu}_{c}^{\top}\bm{\Sigma}_{c}^{-1}\bm{\mu}_{c}$
$\displaystyle=\bm{\mu}_{a}^{\top}\bm{\Sigma}_{a}^{-1}\bm{\mu}_{a}+\bm{\mu}_{b}^{\top}\bm{\Sigma}_{b}^{-1}\bm{\mu}_{b}-\left(\bm{\Sigma}_{a}^{-1}\bm{\mu}_{a}+\bm{\Sigma}_{b}^{-1}\bm{\mu}_{b}\right)^{\top}\bm{\Sigma}_{c}\bm{\Sigma}_{c}^{-1}\bm{\Sigma}_{c}\left(\bm{\Sigma}_{a}^{-1}\bm{\mu}_{a}+\bm{\Sigma}_{b}^{-1}\bm{\mu}_{b}\right)$
$\displaystyle=\bm{\mu}_{a}^{\top}\bm{\Sigma}_{a}^{-1}\bm{\mu}_{a}+\bm{\mu}_{b}^{\top}\bm{\Sigma}_{b}^{-1}\bm{\mu}_{b}-\left(\bm{\Sigma}_{a}^{-1}\bm{\mu}_{a}+\bm{\Sigma}_{b}^{-1}\bm{\mu}_{b}\right)^{\top}\bm{\Sigma}_{c}\left(\bm{\Sigma}_{a}^{-1}\bm{\mu}_{a}+\bm{\Sigma}_{b}^{-1}\bm{\mu}_{b}\right)$
$\displaystyle=\bm{\mu}_{a}^{\top}\underbrace{\left(\bm{\Sigma}_{a}^{-1}-\bm{\Sigma}_{a}^{-1}\bm{\Sigma}_{c}\bm{\Sigma}_{a}^{-1}\right)}_{=:C}\bm{\mu}_{a}+\bm{\mu}_{b}^{\top}\underbrace{\left(\bm{\Sigma}_{b}^{-1}-\bm{\Sigma}_{b}^{-1}\bm{\Sigma}_{c}\bm{\Sigma}_{b}^{-1}\right)}_{=:D}\bm{\mu}_{b}-2\bm{\mu}_{a}^{\top}\underbrace{\bm{\Sigma}_{a}^{-1}\bm{\Sigma}_{c}\bm{\Sigma}_{b}^{-1}}_{=:E}\bm{\mu}_{b},$
$\displaystyle C$
$\displaystyle=\bm{\Sigma}_{a}^{-1}-\bm{\Sigma}_{a}^{-1}\underline{\bm{\Sigma}_{c}}\bm{\Sigma}_{a}^{-1}$
$\displaystyle=\bm{\Sigma}_{a}^{-1}-\bm{\Sigma}_{a}^{-1}\underline{\bm{\Sigma}_{a}\left(\bm{\Sigma}_{a}+\bm{\Sigma}_{b}\right)^{-1}\bm{\Sigma}_{b}}\bm{\Sigma}_{a}^{-1}$
$\displaystyle=\bm{\Sigma}_{a}^{-1}-\left(\bm{\Sigma}_{a}+\bm{\Sigma}_{b}\right)^{-1}\bm{\Sigma}_{b}\bm{\Sigma}_{a}^{-1}$
$\displaystyle=\left(\bm{\Sigma}_{a}+\bm{\Sigma}_{b}\right)^{-1}\left(\left(\bm{\Sigma}_{a}+\bm{\Sigma}_{b}\right)\bm{\Sigma}_{a}^{-1}-\bm{\Sigma}_{b}\bm{\Sigma}_{a}^{-1}\right)$
$\displaystyle=\left(\bm{\Sigma}_{a}+\bm{\Sigma}_{b}\right)^{-1}\left(\bm{I}+\bm{\Sigma}_{b}\bm{\Sigma}_{a}^{-1}-\bm{\Sigma}_{b}\bm{\Sigma}_{a}^{-1}\right)$
$\displaystyle=\left(\bm{\Sigma}_{a}+\bm{\Sigma}_{b}\right)^{-1},$
$\displaystyle D$
$\displaystyle=\bm{\Sigma}_{b}^{-1}-\bm{\Sigma}_{b}^{-1}\bm{\Sigma}_{c}\bm{\Sigma}_{b}^{-1}$
$\displaystyle=\hskip 30.0pt\vdots$
$\displaystyle=\left(\bm{\Sigma}_{a}+\bm{\Sigma}_{b}\right)^{-1},$
$\displaystyle E$
$\displaystyle=\bm{\Sigma}_{a}^{-1}\underline{\bm{\Sigma}_{c}}\bm{\Sigma}_{b}^{-1}$
$\displaystyle=\bm{\Sigma}_{a}^{-1}\underline{\bm{\Sigma}_{a}\left(\bm{\Sigma}_{a}+\bm{\Sigma}_{b}\right)^{-1}\bm{\Sigma}_{b}}\bm{\Sigma}_{b}^{-1}$
$\displaystyle=\left(\bm{\Sigma}_{a}+\bm{\Sigma}_{b}\right)^{-1}.$

Therefore,

$\displaystyle B$
$\displaystyle=\bm{\mu}_{a}^{\top}\left(\bm{\Sigma}_{a}+\bm{\Sigma}_{b}\right)^{-1}\bm{\mu}_{a}+\bm{\mu}_{b}^{\top}\left(\bm{\Sigma}_{a}+\bm{\Sigma}_{b}\right)^{-1}\bm{\mu}_{b}-2\bm{\mu}_{a}^{\top}\left(\bm{\Sigma}_{a}+\bm{\Sigma}_{b}\right)^{-1}\bm{\mu}_{b}$
$\displaystyle=\underbrace{(\bm{\mu}_{a}-\bm{\mu}_{b})^{\top}\left(\bm{\Sigma}_{a}+\bm{\Sigma}_{b}\right)^{-1}(\bm{\mu}_{a}-\bm{\mu}_{b})}_{\text{Exponential
part
of\;}\mathcal{N}\left(\bm{0};\bm{\mu}_{a}-\bm{\mu}_{b},\bm{\Sigma}_{a}+\bm{\Sigma}_{b}\right)}.$

Since
$\left|\bm{\Sigma}_{c}\right|=\left|\bm{\Sigma}_{a}\left(\bm{\Sigma}_{a}+\bm{\Sigma}_{b}\right)^{-1}\bm{\Sigma}_{b}\right|=\frac{\left|\bm{\Sigma}_{a}\right|\left|\bm{\Sigma}_{b}\right|}{\left|\bm{\Sigma}_{a}+\bm{\Sigma}_{b}\right|},$

$\displaystyle P$
$\displaystyle=\sqrt{\frac{\left|2\pi\bm{\Sigma}_{c}\right|}{\left|2\pi\bm{\Sigma}_{a}\right|\left|2\pi\bm{\Sigma}_{b}\right|}}\exp\left(-\frac{1}{2}B\right)\cdot\mathcal{N}\left(\bm{x};\bm{\mu}_{c},\bm{\Sigma}_{c}\right)$
$\displaystyle=\underbrace{\frac{\exp\left(-\frac{1}{2}(\bm{\mu}_{a}-\bm{\mu}_{b})^{\top}\bm{\Sigma}_{a}+\bm{\Sigma}_{b}(\bm{\mu}_{a}-\bm{\mu}_{b})\right)}{\sqrt{\left|2\pi\left(\bm{\Sigma}_{a}+\bm{\Sigma}_{b}\right)\right|}}}_{=\mathcal{N}\left(\bm{0};\bm{\mu}_{a}-\bm{\mu}_{b},\bm{\Sigma}_{a}+\bm{\Sigma}_{b}\right)}\times\mathcal{N}\left(\bm{x};\bm{\mu}_{c},\bm{\Sigma}_{c}\right)$
$\displaystyle=\mathcal{N}\left(\bm{0};\bm{\mu}_{a}-\bm{\mu}_{b},\bm{\Sigma}_{a}+\bm{\Sigma}_{b}\right)\cdot\mathcal{N}\left(\bm{x};\bm{\mu}_{c},\bm{\Sigma}_{c}\right)$

∎

###  A.1 Proof of Lemma 2

###### Proof.

From Lemma 1,

$\displaystyle\mathbb{E}_{\bm{x}\sim\mathcal{N}\left(\bm{\mu}_{a},\bm{\Sigma}_{a}\right)}\left[-\log\mathcal{N}\left(\bm{x};\bm{\mu}_{b},\bm{\Sigma}_{b}\right)\right]$
$\displaystyle=\mathbb{E}_{\bm{x}\sim\mathcal{N}\left(\bm{\mu}_{a},\bm{\Sigma}_{a}\right)}\Biggl{[}\frac{1}{2}(\bm{x}-\bm{\mu}_{b})^{\top}\bm{\Sigma}_{b}^{-1}(\bm{x}-\bm{\mu}_{b})+\frac{1}{2}\log\left|2\pi\bm{\Sigma}_{b}\right|\Biggr{]}$
$\displaystyle=\frac{1}{2}\Bigl{[}\,\mathrm{Tr}\left(\bm{\Sigma}_{b}^{-1}\bm{\Sigma}_{a}\right)+(\bm{\mu}_{a}-\bm{\mu}_{b})^{\top}\bm{\Sigma}_{b}^{-1}(\bm{\mu}_{a}-\bm{\mu}_{b})+\log\left((2\pi)^{d}\cdot\left|\bm{\Sigma}_{b}\right|\right)\,\Bigr{]}$
$\displaystyle=\frac{1}{2}\Bigl{[}\,\mathrm{Tr}\left(\bm{\Sigma}_{a}\bm{\Sigma}_{b}^{-1}\right)+(\bm{\mu}_{a}-\bm{\mu}_{b})^{\top}\bm{\Sigma}_{b}^{-1}(\bm{\mu}_{a}-\bm{\mu}_{b})+\log\left|\bm{\Sigma}_{b}\right|+d\log
2\pi\,\Bigr{]}.$

Let $\bm{\mu}_{b}=\bm{\mu}_{a}$ and $\bm{\Sigma}_{b}=\bm{\Sigma}_{a}$, and
then the Shannon entropy is

$\displaystyle H$
$\displaystyle=\frac{1}{2}\Bigl{\\{}\mathrm{Tr}(\bm{\Sigma}_{a}\bm{\Sigma}_{a}^{-1})+(\bm{\mu}_{a}-\bm{\mu}_{a})^{\top}\bm{\Sigma}^{-1}(\bm{\mu}_{a}-\bm{\mu}_{a})+\log\left|\bm{\Sigma}_{a}\right|+d\log
2\pi\Bigr{\\}}$
$\displaystyle=\frac{1}{2}\left\\{\mathrm{Tr}(\bm{I})+\log\left|\bm{\Sigma}_{a}\right|+d\log
2\pi\right\\}$ $\displaystyle=\frac{1}{2}\left\\{d+d\log
2\pi+\log\left|\bm{\Sigma}_{a}\right|\right\\}.$

∎

###  A.2 Proof of Proposition 1

###### Proof.

From Proposition 3,

$\displaystyle a(\bm{z}|\tilde{\bm{z}})\cdot
a(\bm{z}^{\prime}|\tilde{\bm{z}})$
$\displaystyle=\mathcal{N}\left(\bm{z};\tilde{\bm{z}},\bm{\Sigma}_{\mathrm{aug}}\right)\cdot\mathcal{N}\left(\bm{z}^{\prime};\tilde{\bm{z}},\bm{\Sigma}_{\mathrm{aug}}\right)$
$\displaystyle=\mathcal{N}\left(\tilde{\bm{z}};\bm{z},\bm{\Sigma}_{\mathrm{aug}}\right)\cdot\mathcal{N}\left(\tilde{\bm{z}};\bm{z}^{\prime},\bm{\Sigma}_{\mathrm{aug}}\right)$
$\displaystyle=\mathcal{N}\left(\tilde{\bm{z}};\frac{\bm{z}+\bm{z}^{\prime}}{2},\frac{1}{2}\bm{\Sigma}_{\mathrm{aug}}\right)\cdot\mathcal{N}\left(\bm{0};\bm{z}-\bm{z}^{\prime},2\bm{\Sigma}_{\mathrm{aug}}\right),$

and then,

$\begin{split}&\mathcal{N}\left(\tilde{\bm{z}};\bm{0},\bm{I}\right)\cdot\mathcal{N}\left(\tilde{\bm{z}};\frac{\bm{z}+\bm{z}^{\prime}}{2},\frac{1}{2}\bm{\Sigma}_{\mathrm{aug}}\right)\\\
&=\mathcal{N}\left(\tilde{\bm{z}};\bm{\mu}_{\bar{\bm{z}}},\bm{\Sigma}_{\bar{\bm{z}}}\right)\cdot\mathcal{N}\left(\bm{0};\frac{\bm{z}+\bm{z}^{\prime}}{2},\bm{I}+\frac{1}{2}\bm{\Sigma}_{\mathrm{aug}}\right),\end{split}$

where
$\bm{\Sigma}_{\bar{\bm{z}}}=\left(\bm{I}+2\bm{\Sigma}_{\mathrm{aug}}^{-1}\right)^{-1}\;\text{and}\;\bm{\mu}_{\bar{\bm{z}}}=\bm{\Sigma}_{\bar{\bm{z}}}\bm{\Sigma}_{\mathrm{aug}}^{-1}\left(\bm{z}+\bm{z}^{\prime}\right)$.
Therefore,

$\displaystyle p_{\bm{\psi}}(\bar{\bm{z}})$
$\displaystyle=p_{\bm{\psi}}(\bm{z},\bm{z}^{\prime})$ $\displaystyle=\int
a(\bm{z}|\tilde{\bm{z}})a(\bm{z}^{\prime}|\tilde{\bm{z}})p_{\bm{\psi}}(\tilde{\bm{z}})\mathrm{d}\tilde{\bm{z}}$
$\displaystyle=\int\mathcal{N}\left(\tilde{\bm{z}};\bm{0},\bm{I}\right)a(\bm{z}|\tilde{\bm{z}})\cdot
a(\bm{z}^{\prime}|\tilde{\bm{z}})\>\mathrm{d}\tilde{\bm{z}}$
$\displaystyle=\int\mathcal{N}\left(\tilde{\bm{z}};\bm{\mu}_{\bar{\bm{z}}},\bm{\Sigma}_{\bar{\bm{z}}}\right)\mathcal{N}\left(\bm{0};\frac{\bm{z}+\bm{z}^{\prime}}{2},\bm{I}+\frac{1}{2}\>\bm{\Sigma}_{\mathrm{aug}}\right)\times\mathcal{N}\left(\bm{0};\bm{z}-\bm{z}^{\prime},2\bm{\Sigma}_{\mathrm{aug}}\right)\mathrm{d}\tilde{\bm{z}}$
$\displaystyle=\>\mathcal{N}\left(\bm{0};\frac{\bm{z}+\bm{z}^{\prime}}{2},\bm{I}+\frac{1}{2}\bm{\Sigma}_{\mathrm{aug}}\right)\mathcal{N}\left(\bm{0};\bm{z}-\bm{z}^{\prime},2\bm{\Sigma}_{\mathrm{aug}}\right)\times\underbrace{\int\mathcal{N}\left(\tilde{\bm{z}};\bm{\mu}_{\bar{\bm{z}}},\bm{\Sigma}_{\bar{\bm{z}}}\right)\mathrm{d}\tilde{\bm{z}}}_{=1}$
$\displaystyle=\mathcal{N}\left(\bm{0};\bm{z}-\bm{z}^{\prime},2\bm{\Sigma}_{\mathrm{aug}}\right)\>\mathcal{N}\left(\bm{0};\frac{\bm{z}+\bm{z}^{\prime}}{2},\bm{I}+\frac{1}{2}\bm{\Sigma}_{\mathrm{aug}}\right)$
$\displaystyle=\mathcal{N}\left(\bm{0};\bm{z}-\bm{z}^{\prime},2\bm{\Sigma}_{\mathrm{aug}}\right)\,\mathcal{N}\left(\frac{\bm{z}+\bm{z}^{\prime}}{2};\bm{0},\bm{I}+\frac{1}{2}\bm{\Sigma}_{\mathrm{aug}}\right).$

∎

###  A.3 Proof of Proposition 2

By Corollary 1.7.1 in [26], $\left(2\bm{I}+\bm{\Sigma}_{\text{aug
}}\right)^{-1}=\bm{\Sigma}_{\text{aug }}^{-1}-\bm{\Sigma}_{\text{aug
}}^{-1}\left(2^{-1}\bm{I}+\bm{\Sigma}_{\text{aug
}}^{-1}\right)^{-1}\bm{\Sigma}_{\text{aug }}^{-1}$, and note that

$\left(\bm{\mu}_{\bm{x}}-\bm{\mu}_{\bm{x}^{\prime}}\right)^{\top}\bm{\Sigma}_{\text{aug
}}^{-1}\left(\bm{\mu}_{\bm{x}}\\!-\\!\bm{\mu}_{\bm{x}^{\prime}}\right)\\!+\\!\left(\bm{\mu}_{\bm{x}}\\!+\\!\bm{\mu}_{\bm{x}^{\prime}}\right)^{\top}\bm{\Sigma}_{\text{aug
}}^{-1}\left(\bm{\mu}_{\bm{x}}+\bm{\mu}_{\bm{x}^{\prime}}\right)=2\left(\bm{\mu}_{\bm{x}}^{\top}\bm{\Sigma}_{\text{aug
}}^{-1}\bm{\mu}_{\bm{x}}+\bm{\mu}_{\bm{x}^{\prime}}^{\top}\bm{\Sigma}_{\text{aug
}}^{-1}\bm{\mu}_{\bm{x}^{\prime}}\right).$

The result follows from the above two equations.

##  Appendix B Gaussian Mixture Case

######  Proposition 4.

Let us consider the case that $p(\tilde{\bm{z}})$ in Definition 1 is defined
by a Gaussian Mixture Model (GMM):
$p_{\bm{\psi}}(\tilde{\bm{z}})=\sum_{c=1}^{C}\pi_{c}\,\mathcal{N}\left(\tilde{\bm{z}};\bm{\mu}_{c},\bm{\Sigma}_{c}\right)$.
For the definition of
$\mathcal{N}\left(\tilde{\bm{z}};\bm{\mu}_{c},\bm{\Sigma}_{c}\right)$, see
Eq.(1), and $\pi_{c}$ is the $c$-th wight satisfying
$\pi_{c}\in[0,1]\;\text{and}\;\sum_{c=1}^{C}\pi_{c}=1$. The set
$\left\\{\left(\pi_{c},\bm{\mu}_{c},\bm{\Sigma}_{c}\right)\right\\}_{c=1}^{C}$
is trainable, and it is denoted by $\bm{\psi}$. Then, the prior
$p(\bar{\bm{z}})$ in Definition 1, takes on the form of

$\begin{split}\mathcal{N}\left(\bm{0};\bm{z}-\bm{z}^{\prime},2\bm{\Sigma}_{\mathrm{aug}}\right)\sum_{c=1}^{C}\pi_{c}\,\mathcal{N}\left(\frac{\bm{z}+\bm{z}^{\prime}}{2};\bm{\mu}_{c},\bm{\Sigma}_{c}+\frac{1}{2}\bm{\Sigma}_{\mathrm{aug}}\right).\end{split}$
(15)

###### Proof.

From Proposition 3,

$\displaystyle a(\bm{z}|\tilde{\bm{z}})\cdot
a(\bm{z}^{\prime}|\tilde{\bm{z}})$
$\displaystyle=\mathcal{N}\left(\bm{z};\tilde{\bm{z}},\bm{\Sigma}_{\mathrm{aug}}\right)\cdot\mathcal{N}\left(\bm{z}^{\prime};\tilde{\bm{z}},\bm{\Sigma}_{\mathrm{aug}}\right)$
$\displaystyle=\mathcal{N}\left(\tilde{\bm{z}};\bm{z},\bm{\Sigma}_{\mathrm{aug}}\right)\cdot\mathcal{N}\left(\tilde{\bm{z}};\bm{z}^{\prime},\bm{\Sigma}_{\mathrm{aug}}\right)$
$\displaystyle=\mathcal{N}\left(\tilde{\bm{z}};\frac{\bm{z}+\bm{z}^{\prime}}{2},\frac{1}{2}\bm{\Sigma}_{\mathrm{aug}}\right)\cdot\mathcal{N}\left(\bm{0};\bm{z}-\bm{z}^{\prime},2\bm{\Sigma}_{\mathrm{aug}}\right),$

and then,

$\begin{split}&\mathcal{N}\left(\tilde{\bm{z}};\bm{\mu}_{c},\bm{\Sigma}_{c}\right)\cdot\mathcal{N}\left(\tilde{\bm{z}};\frac{\bm{z}+\bm{z}^{\prime}}{2},\frac{1}{2}\bm{\Sigma}_{\mathrm{aug}}\right)\\\
&=\mathcal{N}\left(\tilde{\bm{z}};\bm{\mu}_{\bar{\bm{z}}},\bm{\Sigma}_{\bar{\bm{z}}}\right)\cdot\mathcal{N}\left(\bm{0};\frac{\bm{z}+\bm{z}^{\prime}}{2}-\bm{\mu}_{c},\bm{\Sigma}_{c}+\frac{1}{2}\bm{\Sigma}_{\mathrm{aug}}\right),\end{split}$

where
$\bm{\Sigma}_{\bar{\bm{z}}}=\left(\bm{\Sigma}_{c}^{-1}+2\bm{\Sigma}_{\mathrm{aug}}^{-1}\right)^{-1}$
and
$\bm{\mu}_{\bar{\bm{z}}}=\bm{\Sigma}_{\bar{\bm{z}}}\left(\bm{\Sigma}_{c}^{-1}\bm{\mu}_{c}+\bm{\Sigma}_{\mathrm{aug}}^{-1}\left(\bm{z}+\bm{z}^{\prime}\right)\right)$.
Therefore,

$\displaystyle p_{\bm{\psi}}(\bar{\bm{z}})$
$\displaystyle=p_{\bm{\psi}}(\bm{z},\bm{z}^{\prime})$ $\displaystyle=\int
a(\bm{z}|\tilde{\bm{z}})a(\bm{z}^{\prime}|\tilde{\bm{z}})p_{\bm{\psi}}(\tilde{\bm{z}})\mathrm{d}\tilde{\bm{z}}$
$\displaystyle=\sum_{c=1}^{C}\pi_{c}\int\mathcal{N}\left(\tilde{\bm{z}};\bm{\mu}_{c},\bm{\Sigma}_{c}\right)a(\bm{z}|\tilde{\bm{z}})\cdot
a(\bm{z}^{\prime}|\tilde{\bm{z}})\>\mathrm{d}\tilde{\bm{z}}$
$\displaystyle=\sum_{c=1}^{C}\pi_{c}\int\mathcal{N}\left(\tilde{\bm{z}};\bm{\mu}_{\bar{\bm{z}}},\bm{\Sigma}_{\bar{\bm{z}}}\right)\times\mathcal{N}\left(\bm{0};\frac{\bm{z}+\bm{z}^{\prime}}{2}-\bm{\mu}_{c},\bm{\Sigma}_{c}+\frac{1}{2}\>\bm{\Sigma}_{\mathrm{aug}}\right)$
$\displaystyle\quad\quad\quad\quad\quad\quad\quad\quad\quad\times\mathcal{N}\left(\bm{0};\bm{z}-\bm{z}^{\prime},2\bm{\Sigma}_{\mathrm{aug}}\right)\mathrm{d}\tilde{\bm{z}}$
$\displaystyle=\sum_{c=1}^{C}\pi_{c}\>\mathcal{N}\left(\bm{0};\frac{\bm{z}+\bm{z}^{\prime}}{2}-\bm{\mu}_{c},\bm{\Sigma}_{c}+\frac{1}{2}\bm{\Sigma}_{\mathrm{aug}}\right)$
$\displaystyle\quad\quad\quad\quad\quad\quad\quad\quad\quad\times\mathcal{N}\left(\bm{0};\bm{z}-\bm{z}^{\prime},2\bm{\Sigma}_{\mathrm{aug}}\right)\times\underbrace{\int\mathcal{N}\left(\tilde{\bm{z}};\bm{\mu}_{\bar{\bm{z}}},\bm{\Sigma}_{\bar{\bm{z}}}\right)\mathrm{d}\tilde{\bm{z}}}_{=1}$
$\displaystyle=\mathcal{N}\left(\bm{0};\bm{z}-\bm{z}^{\prime},2\bm{\Sigma}_{\mathrm{aug}}\right)\times\sum_{i=1}^{C}\pi_{i}\>\mathcal{N}\left(\bm{0};\frac{\bm{z}+\bm{z}^{\prime}}{2}-\bm{\mu}_{c},\bm{\Sigma}_{c}+\frac{1}{2}\bm{\Sigma}_{\mathrm{aug}}\right)$
$\displaystyle=\mathcal{N}\left(\bm{0};\bm{z}-\bm{z}^{\prime},2\bm{\Sigma}_{\mathrm{aug}}\right)\times\sum_{c=1}^{C}\pi_{c}\,\mathcal{N}\left(\frac{\bm{z}+\bm{z}^{\prime}}{2};\bm{\mu}_{c},\bm{\Sigma}_{c}+\frac{1}{2}\bm{\Sigma}_{\mathrm{aug}}\right).$

∎

######  Theorem 2.

Set Eq.(15) to $p(\tilde{\bm{z}})$ of Eq.(4), then the following equation
holds:

$\displaystyle\mathcal{L}\left(\bar{\bm{x}};\bm{\phi},\bm{\theta}\right)$
$\displaystyle=\mathbb{E}_{q_{\bm{\phi}}(\bar{\bm{z}}|\bar{\bm{x}})}\left[\log
p_{\bm{\theta}}(\bar{\bm{x}}|\bar{\bm{z}})\right]+\mathbb{E}_{q_{\bm{\phi}}(\bar{\bm{z}}|\bar{\bm{x}})}\left[\log\sum_{c=1}^{C}\pi_{c}\>\mathcal{N}\left(\frac{\bm{z}+\bm{z}^{\prime}}{2};\bm{\mu}_{c},\bm{\Sigma}_{c}+\frac{1}{2}\bm{\Sigma}_{\mathrm{aug}}\right)\right]$
$\displaystyle\;\;\;\;\;-\frac{1}{4}\Bigl{\\{}\mathrm{Tr}\left(\bm{\Sigma}_{\mathrm{aug}}^{-1}(\bm{\Sigma}_{\bm{x}}+\bm{\Sigma}_{\bm{x}^{\prime}})\right)+(\bm{\mu}_{\bm{x}}-\bm{\mu}_{\bm{x}^{\prime}})^{\top}\bm{\Sigma}_{\mathrm{aug}}^{-1}(\bm{\mu}_{\bm{x}}-\bm{\mu}_{\bm{x}^{\prime}})\Bigr{\\}}$
$\displaystyle\;\;\;\;\;\;+\frac{1}{2}\left(\log\left|\bm{\Sigma}_{\bm{x}}\right|+\log\left|\bm{\Sigma}_{\bm{x}^{\prime}}\right|\right)+d+\frac{d}{2}\log\pi-\frac{1}{2}\log\left|\bm{\Sigma}_{\mathrm{aug}}\right|.$

###### Proof.

From Lemma 1 and 2,

$\displaystyle\mathbb{E}_{\bm{z}\sim\mathcal{N}\left(\bm{\mu}_{\bm{x}},\bm{\Sigma}_{\bm{x}}\right)}\left[\log\mathcal{N}\left(\bm{0};\bm{z}-\bm{z}^{\prime},2\bm{\Sigma}_{\mathrm{aug}}\right)\right]$
$\displaystyle=-\frac{1}{4}\left(\mathrm{Tr}\left(\bm{\Sigma}_{\mathrm{aug}}^{-1}\bm{\Sigma}_{\bm{x}}\right)+(\bm{\mu}_{\bm{x}}-\bm{z}^{\prime})^{\top}\bm{\Sigma}_{\mathrm{aug}}^{-1}(\bm{\mu}_{\bm{x}}-\bm{z}^{\prime})\right)-\frac{1}{2}\log\left|4\pi\bm{\Sigma}_{\mathrm{aug}}\right|.$

Using the above equation,

$\displaystyle\mathbb{E}_{\bar{\bm{z}}\sim
q_{\bm{\phi}}(\bar{\bm{z}}|\bar{\bm{x}})}\left[\log\mathcal{N}\left(\bm{0};\bm{z}-\bm{z}^{\prime},2\bm{\Sigma}_{\mathrm{aug}}\right)\right]$
$\displaystyle=-\frac{1}{4}\Bigl{\\{}\mathrm{Tr}(\bm{\Sigma}_{\mathrm{aug}}^{-1}\bm{\Sigma}_{\bm{x}})+\mathbb{E}_{\bm{z}^{\prime}\sim\mathcal{N}\left(\bm{\mu}_{\bm{x}^{\prime}},\bm{\Sigma}_{\bm{x}^{\prime}}\right)}\left[(\bm{\mu}_{\bm{x}}-\bm{z}^{\prime})^{\top}\bm{\Sigma}_{\mathrm{aug}}^{-1}(\bm{\mu}_{\bm{x}}-\bm{z}^{\prime})\right]\Bigr{\\}}-\frac{1}{2}\log\left|4\pi\bm{\Sigma}_{\mathrm{aug}}\right|$
$\displaystyle=-\frac{1}{4}\Bigl{\\{}\mathrm{Tr}(\bm{\Sigma}_{\mathrm{aug}}^{-1}\bm{\Sigma}_{\bm{x}})+\mathrm{Tr}(\bm{\Sigma}_{\mathrm{aug}}^{-1}\bm{\Sigma}_{\bm{x}^{\prime}})+(\bm{\mu}_{\bm{x}}-\bm{\mu}_{\bm{x}^{\prime}})^{\top}\bm{\Sigma}_{\mathrm{aug}}^{-1}(\bm{\mu}_{\bm{x}}-\bm{\mu}_{\bm{x}^{\prime}})\Bigr{\\}}$
$\displaystyle\quad\quad\quad\quad\quad\quad\quad\quad\quad-\frac{1}{2}\left(\log\left|\bm{\Sigma}_{\mathrm{aug}}\right|+d\log
4\pi\right)$
$\displaystyle=-\frac{1}{4}\mathrm{Tr}\left(\bm{\Sigma}_{\mathrm{aug}}^{-1}(\bm{\Sigma}_{\bm{x}}\\!+\\!\bm{\Sigma}_{\bm{x}^{\prime}})\\!+\\!(\bm{\mu}_{\bm{x}}\\!-\\!\bm{\mu}_{\bm{x}^{\prime}})^{\top}\bm{\Sigma}_{\mathrm{aug}}(\bm{\mu}_{\bm{x}}\\!-\\!\bm{\mu}_{\bm{x}^{\prime}})\right)-\frac{1}{2}\left(\log\left|\bm{\Sigma}_{\mathrm{aug}}\right|+d\log
4\pi\right),$

and

$\displaystyle\mathbb{E}_{q_{\bm{\phi}}(\bar{\bm{z}}|\bar{\bm{x}})}\left[\log
p_{\bm{\psi}}(\bar{\bm{z}})\right]$
$\displaystyle=\mathbb{E}_{q_{\bm{\phi}}(\bar{\bm{z}}|\bar{\bm{x}})}\Biggl{[}\log\sum_{c=1}^{C}\pi_{c}\,\mathcal{N}\left(\frac{\bm{z}+\bm{z}^{\prime}}{2};\bm{\mu}_{c},\bm{\Sigma}_{c}+\frac{1}{2}\bm{\Sigma}_{\mathrm{aug}}\right)+\log\mathcal{N}\left(\bm{0};\bm{z}-\bm{z}^{\prime},2\bm{\Sigma}_{\mathrm{aug}}\right)\Biggr{]}$
$\displaystyle=\mathbb{E}_{q_{\bm{\phi}}(\bar{\bm{z}}|\bar{\bm{x}})}\left[\log\sum_{c=1}^{C}\pi_{c}\,\mathcal{N}\left(\frac{\bm{z}+\bm{z}^{\prime}}{2};\bm{\mu}_{c},\bm{\Sigma}_{c}+\frac{1}{2}\bm{\Sigma}_{\mathrm{aug}}\right)\right]-\frac{1}{2}(\log\left|\bm{\Sigma}_{\mathrm{aug}}\right|+d\log
4\pi)$ $\displaystyle\hskip
20.0pt-\frac{1}{4}\Bigl{\\{}\mathrm{Tr}\left(\bm{\Sigma}_{\mathrm{aug}}^{-1}(\bm{\Sigma}_{\bm{x}}+\bm{\Sigma}_{\bm{x}^{\prime}})\right)+(\bm{\mu}_{\bm{x}}-\bm{\mu}_{\bm{x}^{\prime}})^{\top}\bm{\Sigma}_{\mathrm{aug}}^{-1}(\bm{\mu}_{\bm{x}}-\bm{\mu}_{\bm{x}^{\prime}})\Bigr{\\}}.$

Therefore,

$\displaystyle\mathcal{L}\left(\bar{\bm{x}};\bm{\phi},\bm{\theta}\right)$
$\displaystyle=\mathbb{E}_{q_{\bm{\phi}}(\bar{\bm{z}}|\bar{\bm{x}})}\left[\log
p_{\bm{\theta}}(\bar{\bm{x}}|\bar{\bm{z}})\right]-D_{\mathrm{KL}}\left[q_{\bm{\phi}}(\bar{\bm{z}}|\bar{\bm{x}})\|p_{\bm{\psi}}(\bar{\bm{z}})\right]$
$\displaystyle=\mathbb{E}_{q_{\bm{\phi}}(\bar{\bm{z}}|\bar{\bm{x}})}\left[\log
p_{\bm{\theta}}(\bar{\bm{x}}|\bar{\bm{z}})\right]+\mathbb{E}_{q_{\bm{\phi}}(\bar{\bm{z}}|\bar{\bm{x}})}\left[\log
p_{\bm{\psi}}(\bar{\bm{z}})\right]+H\left(q_{\bm{\phi}}(\bar{\bm{z}}|\bar{\bm{x}})\right)$
$\displaystyle=\mathbb{E}_{q_{\bm{\phi}}(\bar{\bm{z}}|\bar{\bm{x}})}\left[\log
p_{\bm{\theta}}(\bar{\bm{x}}|\bar{\bm{z}})\right]+\mathbb{E}_{q_{\bm{\phi}}(\bar{\bm{z}}|\bar{\bm{x}})}\left[\log\sum_{c=1}^{C}\pi_{c}\>\mathcal{N}\left(\frac{\bm{z}+\bm{z}^{\prime}}{2};\bm{\mu}_{c},\bm{\Sigma}_{c}+\frac{1}{2}\bm{\Sigma}_{\mathrm{aug}}\right)\right]$
$\displaystyle\hskip
20.0pt-\frac{1}{4}\Bigl{\\{}\mathrm{Tr}\left(\bm{\Sigma}_{\mathrm{aug}}^{-1}(\bm{\Sigma}_{\bm{x}}+\bm{\Sigma}_{\bm{x}^{\prime}})\right)+(\bm{\mu}_{\bm{x}}-\bm{\mu}_{\bm{x}^{\prime}})^{\top}\bm{\Sigma}_{\mathrm{aug}}^{-1}(\bm{\mu}_{\bm{x}}-\bm{\mu}_{\bm{x}^{\prime}})\Bigr{\\}}$
$\displaystyle\hskip
30.0pt+\frac{1}{2}\left(\log\left|\bm{\Sigma}_{\bm{x}}\right|+\log\left|\bm{\Sigma}_{\bm{x}^{\prime}}\right|\right)+d\left(1+\frac{1}{2}\log\pi\right)-\frac{1}{2}\log\left|\bm{\Sigma}_{\mathrm{aug}}\right|.$

∎

##  Appendix C Further Details of Numerical Experiments

Let $\tilde{\bm{x}}^{\text{tr}}_{i}$ and $\tilde{\bm{x}}^{\text{ts}}_{j}$
denote $i$-th original data in training dataset and $j$-th original data in
test dataset, respectively. In addition, let $\mathcal{D}^{\rm
tr}=\\{\tilde{\bm{x}}^{\text{tr}}_{i}\\}_{i=1}^{n}$ and $\mathcal{D}^{\rm
ts}=\\{\tilde{\bm{x}}^{\text{ts}}_{j}\\}_{j=1}^{m}$, where $n$ (resp. $m$)
denotes the number of training (resp. test) data points.

###  C.1 Preliminary Experiments

Using the original MNIST dataset [15], we firstly evaluated the performance of
the following two trained model:

(A)

a VAE trained by $\mathcal{D}^{\rm tr}$,

(B)

a VAE trained by
$\\{\tilde{\bm{x}}^{\text{tr}}_{i}\\}_{i=1}^{n}\cup\\{\bm{x}^{\text{tr}}_{i}\\}_{i=1}^{n}$,
where
$\bm{x}^{\text{tr}}_{i}=\tilde{\bm{x}}^{\text{tr}}_{i}+\bm{\epsilon}_{i}$, and
$\bm{\epsilon}_{i}$ is $i$-th random Gaussian noise.

The dimension of the latent space was ten for both VAEs, and both
architectures (encoder and decoder) were the exactly same as described at the
beginning of Section 4.1.2. Additionally, $(n,m)=(60000,10000)$ with the
original MNIST. Let $\textrm{Enc}^{\rm(A)}_{\bm{\phi}^{\ast}}$ (resp.
$\textrm{Enc}^{\rm(B)}_{\bm{\phi}^{\ast}}$) denote the trained encoder of (A)
(resp. (B)).

As shown in a) and b) of Figure 2, the accuracies of (A) and (B) were almost
the same. We then investigated the following two sets of the latent variables:
$\mathcal{Z}_{\rm(A)}=\left\\{\tilde{\bm{z}}^{\rm
ts}_{j,{\rm(A)}}\right\\}_{j=1}^{m}\cup\left\\{\bm{z}^{\rm
ts}_{j,{\rm(A)}}\right\\}_{j=1}^{m}$ and
$\mathcal{Z}_{\rm(B)}=\left\\{\tilde{\bm{z}}^{\rm
ts}_{j,{\rm(B)}}\right\\}_{j=1}^{m}\cup\left\\{\bm{z}^{\rm
ts}_{j,{\rm(B)}}\right\\}_{j=1}^{m}$, where $\tilde{\bm{z}}^{\rm
ts}_{j,{\rm(\cdot)}}:=\textrm{Enc}^{\rm(\cdot)}_{\bm{\phi}^{\ast}}\left(\tilde{\bm{x}}^{\text{ts}}_{j}\right)$
and $\bm{z}^{\rm
ts}_{j,{\rm(\cdot)}}:=\textrm{Enc}^{\rm(\cdot)}_{\bm{\phi}^{\ast}}\left(\tilde{\bm{x}}^{\text{ts}}_{j}+\bm{\epsilon}_{j}\right)$.
In Figure 1, the t-SNE visualization of $\mathcal{Z}_{\rm(A)}$ and
$\mathcal{Z}_{\rm(B)}$ are shown in (A) and (B), respectively. Consider the
set $S:=\left\\{\left\|\tilde{\bm{z}}^{\rm ts}_{j}-\bm{z}^{\rm
ts}_{j}\right\|_{2}^{2}:j=1,...,m\right\\}$ on $\mathcal{Z}$. Then, the
average and std over $S$ for (A) and (B) are shown in the first and second
column in Table 1.

###  C.2 Hyper-Parameter Tuning and Model Architecture

In Section 4, we fixed the dimension of the latent space to $10$ for all
methods. In addition, the hyper-parameters in all the existing methods were
sufficiently tuned in our preliminary experiments using both MNIST and
F-MNIST. For our method, Raven, the primary hyperparameter
$\bm{\Sigma}_{\mathrm{aug}}$ is set to be $\sigma^{2}I$. We explored a range
for $\sigma\in\\{0.1,0.2,0.4,0.1,0.2,0.4,1\\}$ and found the best performing
$\sigma=0.04$ for F-MNIST and $\sigma=0.01$ for MNIST, based on the metrics
described at the end of Section 4.1.2.

In our experiments, we used a fully connected feed-forward neural network with
layers arranged as $784-500-250-10$ (encoder), $10-250-500-784$ (decoder) with
PReLU activation functions [23].

For SE, the key hyperparameter is $\eta$, which controls the coupling strength
for the regularized Wasserstein distance between the latent representations of
the original data points and adversarial data points. From
$\eta\in\\{50,100,250,500,1000,2500,5000,10000\\}$, we found $\eta=500$
performs the best for MNIST and $\eta=5000$ for F-MNIST.

## References

  * [1] D. P. Kingma and M. Welling, “Auto-encoding variational bayes,” in _International Conference on Learning Representations_ , 2014. [Online]. Available: http://arxiv.org/abs/1312.6114
  * [2] I. Goodfellow, J. Pouget-Abadie, M. Mirza, B. Xu, D. Warde-Farley, S. Ozair, A. Courville, and Y. Bengio, “Generative adversarial nets,” in _Advances in Neural Information Processing Systems_ , Z. Ghahramani, M. Welling, C. Cortes, N. Lawrence, and K. Weinberger, Eds., vol. 27. Curran Associates, Inc., 2014. [Online]. Available: https://proceedings.neurips.cc/paper_files/paper/2014/file/5ca3e9b122f61f8f06494c97b1afccf3-Paper.pdf
  * [3] D. Rezende and S. Mohamed, “Variational inference with normalizing flows,” in _Proceedings of the 32nd International Conference on Machine Learning_ , ser. Proceedings of Machine Learning Research, F. Bach and D. Blei, Eds., vol. 37. Lille, France: PMLR, 07–09 Jul 2015, pp. 1530–1538. [Online]. Available: https://proceedings.mlr.press/v37/rezende15.html
  * [4] Y. Song, J. Sohl-Dickstein, D. P. Kingma, A. Kumar, S. Ermon, and B. Poole, “Score-based generative modeling through stochastic differential equations,” in _International Conference on Learning Representations_ , 2021\. [Online]. Available: https://openreview.net/forum?id=PxTIG12RRHS
  * [5] A. van den Oord, N. Kalchbrenner, and K. Kavukcuoglu, “Pixel recurrent neural networks,” in _Proceedings of The 33rd International Conference on Machine Learning_ , ser. Proceedings of Machine Learning Research, M. F. Balcan and K. Q. Weinberger, Eds., vol. 48. New York, New York, USA: PMLR, 20–22 Jun 2016, pp. 1747–1756. [Online]. Available: https://proceedings.mlr.press/v48/oord16.html
  * [6] H. Touvron, T. Lavril, G. Izacard, X. Martinet, M.-A. Lachaux, T. Lacroix, B. RoziÃ¨re, N. Goyal, E. Hambro, F. Azhar, A. Rodriguez, A. Joulin, E. Grave, and G. Lample, “Llama: Open and efficient foundation language models,” in _arXiv_ , 2023.
  * [7] U. Singer, A. Polyak, T. Hayes, X. Yin, J. An, S. Zhang, Q. Hu, H. Yang, O. Ashual, O. Gafni, D. Parikh, S. Gupta, and Y. Taigman, “Make-a-video: Text-to-video generation without text-video data,” in _The Eleventh International Conference on Learning Representations_ , 2023. [Online]. Available: https://openreview.net/forum?id=nJfylDvgzlq
  * [8] E. D. Zhong, T. Bepler, B. Berger, and J. H. Davis, “CryoDRGN: reconstruction of heterogeneous cryo-EM structures using neural networks,” _Nature Methods_ , vol. 18, no. 2, pp. 176–185, 2021.
  * [9] T. Cemgil, S. Ghaisas, K. D. Dvijotham, and P. Kohli, “Adversarially robust representations with smooth encoders,” in _International Conference on Learning Representations_ , 2020. [Online]. Available: https://openreview.net/forum?id=H1gfFaEYDS
  * [10] A. Kuzina, M. Welling, and J. Tomczak, “Alleviating adversarial attacks on variational autoencoders with mcmc,” in _Advances in Neural Information Processing Systems_ , S. Koyejo, S. Mohamed, A. Agarwal, D. Belgrave, K. Cho, and A. Oh, Eds., vol. 35. Curran Associates, Inc., 2022, pp. 8811–8823. [Online]. Available: https://proceedings.neurips.cc/paper_files/paper/2022/file/39e9c5913c970e3e49c2df629daff636-Paper-Conference.pdf
  * [11] S. Eduardo, A. Nazabal, C. K. I. Williams, and C. Sutton, “Robust variational autoencoders for outlier detection and repair of mixed-type data,” in _Proceedings of the Twenty Third International Conference on Artificial Intelligence and Statistics_ , ser. Proceedings of Machine Learning Research, S. Chiappa and R. Calandra, Eds., vol. 108. PMLR, 26–28 Aug 2020, pp. 4056–4066. [Online]. Available: https://proceedings.mlr.press/v108/eduardo20a.html
  * [12] M. J. Willetts, A. Camuto, T. Rainforth, S. Roberts, and C. C. Holmes, “Improving VAEs’ robustness to adversarial attack,” in _International Conference on Learning Representations_ , 2021. [Online]. Available: https://openreview.net/forum?id=-Hs_otp2RB
  * [13] T. Cemgil, S. Ghaisas, K. Dvijotham, S. Gowal, and P. Kohli, “The autoencoding variational autoencoder,” in _Advances in Neural Information Processing Systems_ , H. Larochelle, M. Ranzato, R. Hadsell, M. Balcan, and H. Lin, Eds., vol. 33. Curran Associates, Inc., 2020, pp. 15 077–15 087. [Online]. Available: https://proceedings.neurips.cc/paper_files/paper/2020/file/ac10ff1941c540cd87c107330996f4f6-Paper.pdf
  * [14] B. Li, C. Chen, W. Wang, and L. Carin, “Certified adversarial robustness with additive noise,” in _Advances in Neural Information Processing Systems_ , H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett, Eds., vol. 32. Curran Associates, Inc., 2019. [Online]. Available: https://proceedings.neurips.cc/paper_files/paper/2019/file/335cd1b90bfa4ee70b39d08a4ae0cf2d-Paper.pdf
  * [15] Y. Lecun, L. Bottou, Y. Bengio, and P. Haffner, “Gradient-based learning applied to document recognition,” _Proceedings of the IEEE_ , vol. 86, no. 11, pp. 2278–2324, 1998.
  * [16] L. van der Maaten and G. Hinton, “Visualizing data using t-sne,” _Journal of Machine Learning Research_ , vol. 9, no. 86, pp. 2579–2605, 2008\. [Online]. Available: http://jmlr.org/papers/v9/vandermaaten08a.html
  * [17] J. Z. HaoChen, C. Wei, A. Gaidon, and T. Ma, “Provable guarantees for self-supervised deep learning with spectral contrastive loss,” in _Advances in Neural Information Processing Systems_ , M. Ranzato, A. Beygelzimer, Y. Dauphin, P. Liang, and J. W. Vaughan, Eds., vol. 34. Curran Associates, Inc., 2021, pp. 5000–5011. [Online]. Available: https://proceedings.neurips.cc/paper_files/paper/2021/file/27debb435021eb68b3965290b5e24c49-Paper.pdf
  * [18] I. Higgins, L. Matthey, A. Pal, C. Burgess, X. Glorot, M. Botvinick, S. Mohamed, and A. Lerchner, “beta-VAE: Learning basic visual concepts with a constrained variational framework,” in _International Conference on Learning Representations_ , 2017. [Online]. Available: https://openreview.net/forum?id=Sy2fzU9gl
  * [19] Z. Jiang, Y. Zheng, H. Tan, B. Tang, and H. Zhou, “Variational deep embedding: An unsupervised and generative approach to clustering,” in _Proceedings of the Twenty-Sixth International Joint Conference on Artificial Intelligence, IJCAI-17_ , 2017, pp. 1965–1972. [Online]. Available: https://doi.org/10.24963/ijcai.2017/273
  * [20] A. M. Mathai and S. B. Provost, _Quadratic Forms in Random Variables_. CRC Press, 1992.
  * [21] H. Xiao, K. Rasul, and R. Vollgraf, “Fashion-mnist: a novel image dataset for benchmarking machine learning algorithms,” in _arXiv_ , 2017.
  * [22] I. Goodfellow, J. Shlens, and C. Szegedy, “Explaining and harnessing adversarial examples,” in _International Conference on Learning Representations_ , 2015. [Online]. Available: http://arxiv.org/abs/1412.6572
  * [23] K. He, X. Zhang, S. Ren, and J. Sun, “Delving deep into rectifiers: Surpassing human-level performance on imagenet classification,” in _2015 IEEE International Conference on Computer Vision (ICCV)_ , 2015, pp. 1026–1034.
  * [24] L. Liu, H. Jiang, P. He, W. Chen, X. Liu, J. Gao, and J. Han, “On the variance of the adaptive learning rate and beyond,” in _International Conference on Learning Representations_ , 2020. [Online]. Available: https://openreview.net/forum?id=rkgz2aEKDr
  * [25] M. Heusel, H. Ramsauer, T. Unterthiner, B. Nessler, and S. Hochreiter, “Gans trained by a two time-scale update rule converge to a local nash equilibrium,” in _Advances in Neural Information Processing Systems_ , I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett, Eds., vol. 30. Curran Associates, Inc., 2017. [Online]. Available: https://proceedings.neurips.cc/paper_files/paper/2017/file/8a1d694707eb0fefe65871369074926d-Paper.pdf
  * [26] J. R. Schott, _Matrix analysis for statistics_ , 2nd ed. John Wiley & Sons, 2005.