# A moment-matching metric for latent variable generative models

Cédric Beaulac

###### Abstract

It is difficult to assess the quality of a fitted model when facing
unsupervised learning problems. Latent variable models, such as variational
autoencoders and Gaussian mixture models, are often trained with likelihood-
based approaches. In the scope of Goodhart’s law, when a metric becomes a
target it ceases to be a good metric and therefore we should not use
likelihood to assess the quality of the fit of these models. The solution we
propose is a new metric for model comparison or regularization that relies on
moments. The key idea to study the difference between the data moments and the
model moments using a matrix norm, such as the Frobenius norm. We show how to
use this new metric for model comparison and then for regularization. We show
that our proposed metric is faster to compute and has a smaller variance than
the commonly used procedure of drawing samples from the fitted distribution.
We conclude this article with a proof of concept for both applications and we
discuss future work.

Keywords : Moment estimators, Latent Variable models, Gaussian Mixture Models,
Variational Auto-Encoders, Frobenius norm

##  1 Introduction

When fitting supervised models, statisticians and computer scientists alike
have come up with a variety of metrics in order to evaluate the quality of
their predictions, from simple mean-squared error to general loss functions.
However, in the context of unsupervised learning there is no direct measure of
success and it can be difficult to assess the validity of the fitted model
[12].

Assume an unsupervised learning context where we have a data set
$S=\\{\mathbf{x}_{1},...,\mathbf{x}_{N}\\}$ and the abstract goal of capturing
the distribution $p(\mathbf{x})$. One way to train such model is to assume a
distribution and then maximize the likelihood of the data set with respect to
the distribution parameters. Then the trained models are assessed and compared
using the likelihood as well. Goodhart’s law [10, 27] states that when a
measure becomes a target, it ceases to be a good measure and thus we should
not strictly rely on the likelihood to evaluate models that were trained using
the likelihood.

In this article we propose a new way to assess the quality of the fit of a
large family of unsupervised models with respect to our abstract goal of
capturing the distribution $p(\mathbf{x})$. In other words, we propose a way
to measure if estimated distribution $\hat{p}(x)$ resembles the observed
distribution $p_{S}(x)$. More precisely, we offer a diagnostic technique for
parametric latent variable models such as Variational AutoEncoders (VAEs) [15,
16] and Gaussian Mixture Models (GMMs). Our proposed metric evaluates the
quality of the fitted model; we compare the learned parameters of the
unsupervised model with the observed data distribution the model is trying to
capture. We do so by building two distinct moment estimators. The main purpose
of such metric is to provide a way to compare the fit of multiple models from
different families by assessing how well these models captured the first two
moments of the data. Though capturing the first and second moment of a data
set is not a sufficient condition to claim the trained model has captured the
data distribution it certainly is a necessary condition under some assumption
we discuss later.

In statistics and machine learning, Goodhart’s law is often compared with the
concept of overfitting. One popular way to circumvent model overfitting has
been regularization. Consequently, we offer a second perspective on our new
metric; it can be used for regularization. Our metric favours simple models
such as a simple Gaussian distribution, thus it can be easily integrated in
the optimization procedure as a regularize.

The technique we propose is fast to compute, works for a wide range of models
and is built upon a rigorous mathematical formulation. It provides a new way
to compare multiple models or regularize them and behaves similarly to
previously used heuristic techniques.

In the next section, we establish the family of latent variable models suited
for this metric. After we discuss related work in section 3, we discuss the
moments estimators used in section 4. In section 5, we present our metric and
its implementation and next we demonstrate how it performs on simple examples
in section 6. We then introduce the framework for the application of our
metric for regularization in section 7 and we demonstrate how it behaves as a
regularizer in section 8. Finally, we discuss the limitations of our approach
in section 9 before some concluding remarks in section 10.

##  2 Latent variable generative models

In this section we define the family of Latent variable generative models
(LVGMs).

Assume we have a data set $S=\\{\mathbf{x}_{1},...,\mathbf{x}_{N}\\}$
consisting of $N$ observations of a $D$-dimensional variable $\mathbf{x}$.
Assume $\mathbf{x}$ $\in$ $\mathcal{X}$ which is $D$-dimensional. We want to
estimate the distribution of the random variable $\mathbf{x}$ but it is too
complicated to be captured by a simple distribution. Latent variable models
suppose there exist an unobserved latent variable, say $\mathbf{z}$, that has
a direct influence on the distribution of $\mathbf{x}$

$\displaystyle p(\mathbf{x})=\int_{z}p(\mathbf{x}|z)p(\mathbf{z})dz,$ (1)

where we assume $\mathbf{z}$ $\in$ $\mathcal{Z}$ which is $M$-dimensional. The
model proposed by equation (1) is quite general but allows relatively complex
marginal distributions over observed variables $\mathbf{x}$ to be expressed in
terms of more tractable conditional distributions $p(\mathbf{x}|z)$ [7].
Similarly, it leads to a tractable joint distribution as well

$\displaystyle p(\mathbf{x},\mathbf{z})=p(\mathbf{z})p(\mathbf{x}|z),$ (2)

and this is quite often represented using a simple graph as seen in Figure 1.

$\mathbf{z}$$\mathbf{x}$ Figure 1: Graphical representation of latent
variables models with joint distribution
$p(\mathbf{x},\mathbf{z})=p(\mathbf{z})p(\mathbf{x}|z)$.

These models are generative models because learning
$p(\mathbf{x},\mathbf{z})=p(\mathbf{z})p(\mathbf{x}|z)$ allows us to generate
new samples of $\mathbf{x}$ using ancestral sampling. What makes this model
probabilistic is that the mapping from $\mathbf{z}$ to $\mathbf{x}$ is not a
deterministic function $f:\mathcal{Z}\rightarrow\mathcal{X}$ but instead a
probabilistic mapping from $\mathcal{Z}$ to $\Theta$, where $\Theta$ is the
parameter space of $p_{\theta}(\mathbf{x}|z)$; $\theta\in\Theta$. We call
$p_{\theta}(\mathbf{x}|z)$ the emission distribution or observation
distribution interchangeably.

When training or fitting such models, we train the function
$f:\mathcal{Z}\rightarrow\Theta$ to maximize the likelihood of the data set
$S$ under the model of equation (1). This mapping $f$ explains the effect of
$\mathbf{z}$ on $\mathbf{x}$ and is at the centre of latent variable models.
Learning this function $f$ is the main challenge of training latent variable
models and the Expectation-Maximization algorithm (EM) or variational
inference are commonly used strategies to learn this function. In most cases,
$p(\mathbf{z})$ is assumed to be known and fixed but in some cases the
parameters of $p(\mathbf{z})$ are estimated as well.

Usually $p_{\theta}(\mathbf{x}|z)$ is a simple parametric distribution and the
latent variable increases the complexity of $p_{\theta}(\mathbf{x})$.
Additionally, the function $f$ can take many forms, from simple linear
combination to neural network functions. We use $f(\mathbf{z})$
interchangeably with $f$ or the distribution parameters it outputs directly.

Let us introduce a simple example. Assume the emission distribution is
Poisson: $p_{\theta}(\mathbf{x}|z)=\text{Poisson}(\lambda)$ then
$f:\mathcal{Z}\rightarrow\mathbb{R}^{+}$ because
$\theta=\lambda\in\mathbb{R}^{+}$ and we use $f(\mathbf{z})$ and
$\lambda(\mathbf{z})$ interchangeably. If there exist a simple mapping from
the parameters of the distribution to its expectation and its variance, we
also use them interchangeably. For the Poisson example,
$\textbf{E}_{x}[\mathbf{x}|\mathbf{z}]=\lambda(\mathbf{z})$ and
$\textbf{Var}_{x}[\mathbf{x}|\mathbf{z}]=\lambda(\mathbf{z})$.

One important detail to bring up is that the moments are only meaningful for a
certain family of emission distributions. For the application of our metric,
we will consider the family of emission distribution for which the moment
generating function exists.

###  2.1 Probabilistic Principal Component Analysis

The Probabilistic Principal Component Analysis (pPCA) [29, 7] is a member of
the LVGM family we just described where $p(\mathbf{z})$ is assumed to be a
normal distribution $N(0,I)$. We also assume the emission distribution to be
Normal: $p(\mathbf{x}|\mathbf{z})=N(W\mathbf{z}+b,I\sigma^{2})$. In this
formulation, we see that $f:\mathbb{R}_{M}\rightarrow\mathbb{R}_{D}$ is a
linear function that maps the latent variable $\mathbf{z}$ to
$\textbf{E}_{x}[\mathbf{x}|\mathbf{z}]$:
$\textbf{E}_{x}[\mathbf{x}|\mathbf{z}]=\mu(\mathbf{z})=W\mathbf{z}+b$ .
Outside of estimating $W$ and $b$ as part of $f$, the model also estimates the
parameter $\sigma^{2}$ though it is not a function of $z$. However it is a
function of $D$ and $M$, the dimension of $\mathcal{X}$ and $\mathcal{Z}$.

The parameters of pPCA can be obtained analytically as the solution of a
direct maximization of the likelihood or with the EM algorithm.

###  2.2 Variational AutoEncoders

The VAE is also a member of the LVGM family. It is assumed that $\mathbf{z}$
is a continuous variable where $p(\mathbf{z})$ is assumed to be $N(0,I)$ in
the introductory papers [15, 16]. $p(\mathbf{x}|\mathbf{z})$ can be any
parametric distribution where $f(\mathbf{z})$ outputs the parameters of this
distribution. For instance, if $p(\mathbf{x}|\mathbf{z})$ is normal then
$f(\mathbf{z})$ will output a mean and a variance parameter,
$f:\mathbb{R}_{M}\rightarrow\mathbb{R}_{D}\times\mathbb{R}_{(D,D)}^{+}$.

One novelty of VAEs is that the function $f$ proposed is much more flexible
than a linear combination; it is a neural network. In turn, this makes the
posterior distribution $p(\mathbf{z}|\mathbf{x})$ intractable and prevents the
model from being fitted by the EM algorithm. The solution proposed is to
assume a variational distribution $q(\mathbf{z}|\mathbf{x})$ and optimize the
likelihood by maximizing the Evidence Lower BOund (ELBO) [7, 15], a lower
bound of the observed-data log-likelihood.

###  2.3 Gaussian Mixture Models

For a GMM with $K$-components we define $\mathbf{z}$ as a k-class categorical
variable, $\mathbf{z}\in\\{1,...,K\\}=\mathcal{Z}$, $p(\mathbf{z})$ is a
categorical distribution where $\pi_{j}=p(\mathbf{z}=j)$, and
$\sum_{j=1}^{K}\pi_{j}=1$. Finally, setting
$p(\mathbf{x}|z=j)=N(\mu_{j},\Sigma_{j})$ leads to a GMM:

$\displaystyle p(\mathbf{x})=\sum_{j=1}^{K}\pi_{j}p(\mathbf{x}|z=j).$ (3)

In this situation, $f$ maps the latent variable to a pair of distribution
parameters, $\mu$ and $\Sigma$,
$f:\\{1,...,K\\}\rightarrow\mathbb{R}_{D}\times\mathbb{R}_{D\times D}^{+}$. In
this particular case $\textbf{E}_{x}[\mathbf{x}|\mathbf{z}]=f_{1}(\mathbf{z})$
where $f_{1}$ is the first output of $f(\mathbf{z})$ or simply
$\mu(\mathbf{z})$ and
$\textbf{Var}_{x}[\mathbf{x}|\mathbf{z}]=f_{2}(\mathbf{z})=\Sigma(\mathbf{z})$.

The GMM is a special case of LVGM where we also estimate the parameters
$\\{\pi_{j}:j\in 1,...K\\}$ of $p(\mathbf{z})$, identifiable up to a
permutation. A GMM is usually trained with the EM algorithm.

##  3 Related Works

In this article, we propose a metric to evaluate the goodness-of-fit of the
family of latent variable models defined in section 2 and in this section we
discuss some common alternatives. When proposing new LVGMs researchers rely
either on the likelihood of the data under the fitted model or a heuristic
analysis of generated data points [15, 19, 14, 31, 30]. To evaluate the
performance of the models, both those techniques have their fair share of
problems which we discuss in this section. The metric we propose is an
alternative to those techniques.

A problem with evaluating models with the likelihood is that a high likelihood
does not necessarily mean that the proposed model captured the distribution of
the observed data. For instance, Bishop et al. [7] demonstrate that for GMMs
it is possible to have a likelihood be infinite ($\infty$) by setting the mean
of one component to be exactly one of the observed points, say
$\mathbf{x}_{i}$, and then pushing the variance of that component to $0$, thus
the likelihood of $\mathbf{x}_{i}$ under that particular component will be
infinite. Similarly, Zhao et al. [31] built a toy example where the ELBO (a
lower bound of the log-likelihood) would converge to infinity.

Another commonly employed strategy is to generate new observations and try to
determine if they look like real data with a simple visual inspection, this is
a technique used by many authors [15, 19, 14, 16, 31, 30]. A weakness of
visual inspection is that it is subjective, it incentivizes cherry-picking of
results and it is not a rigorous means of comparing models.

Comparing data moments with model moments have been proposed in the past but
only for training purposes [3, 4, 8, 19, 23]. Anandkumar [3, 4] proposes
efficient ways to code and optimize latent variable models by comparing the
true data set with a sample generated from the model. Podosinnikova [23]
approaches this topic very thoroughly in their thesis where they also discuss
the use of moment-generating function estimators. What we propose here is
different; we propose a metric. Even though we discuss the possibility of
using our metric for optimization in later sections we believe there already
exists a rich literature that discusses new ways of optimizing latent variable
models but few publications that address the lack of evaluation metrics.

##  4 Moment estimators

In this section we define two different moment estimators for the first and
the second moment. Our goal is to build different estimators that contain
different information. To begin, we define moment estimators for the data set,
we call those Data Estimators (DE). Then, we define another set of moment
estimators that represent the distribution captured by the LVGM, we call those
Forward Model Estimators (FME).

###  4.1 Second moment

We build two different estimators of the same quantity, the second moment. One
uses observed data while the other uses the proposed generative model. What
makes the proposed FME different is that we do not sample new data points from
the LVGM but instead rely on a simple probability identity to build the FME.

To do so, let us introduce the well-know Law of Total Variance

$\displaystyle\textbf{Var}_{x}(\mathbf{x})=\textbf{E}_{z}[\textbf{Var}_{x}(\mathbf{x}|\mathbf{z})]+\textbf{Var}_{z}[\textbf{E}_{x}(\mathbf{x}|\mathbf{z})],$
(4)

and notice the second term is

$\displaystyle\textbf{Var}_{z}[\textbf{E}_{x}(\mathbf{x}|\mathbf{z})]$
$\displaystyle=\textbf{E}_{z}[\textbf{E}_{x}(\mathbf{x}|\mathbf{z})^{2}]-(\textbf{E}_{z}[\textbf{E}_{x}(\mathbf{x}|\mathbf{z})])^{2}$
(5)
$\displaystyle=\textbf{E}_{z}[\textbf{E}_{x}(\mathbf{x}|\mathbf{z})^{2}]-(\textbf{E}_{x}[\mathbf{x}])^{2}.$
(6)

We combine and reorganize both equations

$\displaystyle\textbf{Var}_{x}(\mathbf{x})+(\textbf{E}_{x}[\mathbf{x}])^{2}=\textbf{E}_{z}[\textbf{Var}_{x}(\mathbf{x}|\mathbf{z})]+\textbf{E}_{z}[\textbf{E}_{x}(\mathbf{x}|\mathbf{z})^{2}].$
(7)

We have reorganized both terms in this particular way because the left-hand
side of equation (7) is independent of the latent variables and can be
estimated from $S$ independently from the choice of model while the right-hand
side contains information about both the expectation and the variance of the
generative model. Additionally, notice the left-hand side is actually
$\textbf{E}_{x}[\mathbf{x}^{2}]$, the second moment of $\mathbf{x}$ and thus
our work here consists of comparing two different estimators of
$\textbf{E}_{x}[\mathbf{x}^{2}]$ which we introduce next. The left-hand side
of equation (7) can be estimated using the observed data

$\displaystyle\textbf{Var}_{x}(\mathbf{x})+(\textbf{E}_{x}[\mathbf{x}])^{2}\approx\frac{\sum_{i=1}^{n}(\mathbf{x}_{i}-\bar{\mathbf{x}})^{T}(\mathbf{x}_{i}-\bar{x})}{n-1}+\bar{\mathbf{x}}^{T}\bar{\mathbf{x}}:=\text{DE},$
(8)

where $\bar{\mathbf{x}}$ is the mean vector. The right-hand side of equation
(7) can be estimated using the proposed generative model through a Monte Carlo
sample of $p_{\theta}(z)$ and using both
$\textbf{Var}_{x}(\mathbf{x}|\mathbf{z})$ and
$\textbf{E}_{x}(\mathbf{x}|\mathbf{z})$.

$\displaystyle\textbf{E}_{z}[\textbf{Var}_{x}(\mathbf{x}|\mathbf{z})+\textbf{E}_{x}(\mathbf{x}|\mathbf{z})^{2}]$
$\displaystyle=\int_{z}(\textbf{Var}_{x}(X|Z=z)+\textbf{E}_{x}(\mathbf{z}|\mathbf{z}=z)^{2})p(z)dz$
(9)
$\displaystyle\approx\frac{1}{m}\sum_{i=1}^{m}\left[\textbf{Var}_{x}(\mathbf{x}|\mathbf{z}=z_{i})+\textbf{E}_{x}(\mathbf{x}|\mathbf{z}=z_{i})^{T}\textbf{E}_{x}(\mathbf{x}|\mathbf{z}=z_{i})\right]:=\text{FME},$
(10)

where $z_{i}\sim p(\mathbf{z})$, and
$\textbf{E}_{x}(\mathbf{x}|\mathbf{z}=z_{i})$ and
$\textbf{Var}_{x}(\mathbf{x}|\mathbf{z}=z_{i})$ are expressed as functions of
$f(z_{i})$. Thus this estimator relies on both components of the fitted LVGM:
$p(\mathbf{z})$ and $f(\mathbf{z})$. This is the forward model estimates
(FME). Notice that this estimator does not require we sample from
$p_{\theta}(x|z)$ and directly uses the estimated parameters of the emission
distribution. It is faster to sample a large amount of $z$ than sample a large
amount of $x$ because traditionally $M<<D$. Additionally, this is a rather
simple Monte Carlo sample and it is unbiased [25].

Thus it follows from equation (7) that

$\displaystyle(\textbf{Var}_{x}(\mathbf{x})+(\textbf{E}_{x}[\mathbf{x}])^{2})-(\textbf{E}_{z}[\textbf{Var}_{x}(\mathbf{x}|\mathbf{z})]+\textbf{E}_{z}[\textbf{E}_{x}(\mathbf{x}|\mathbf{z})^{2}])=0,$
(11)

and since both DE and FME are unbiased estimators of the second moment then
consequently, the gap between those two estimators reflects how the LVGM
captured the second moment of the data set; the bigger the gap is, the poorer
the fit is. Thus, we propose to analyse the following moment estimator gap

$\displaystyle\text{DE}-\text{FME},$ (12)

which is a matrix of dimension $D\times D$.

###  4.2 First moment

We have a similar result for the first moment

$\displaystyle\textbf{E}_{x}[\mathbf{x}]=\textbf{E}_{z}[\textbf{E}_{x}(\mathbf{x}|\mathbf{z})],$
(13)

where we estimate the left-hand side with $\bar{x}$ (DE) and the right-hand
side with
$\frac{1}{m}\sum_{i=1}^{m}\textbf{E}_{x}(\mathbf{x}|\mathbf{z}=z_{i})$ where
$z_{i}\sim p_{\theta}(z)$ (FME). However, to get a complete picture of the
fitted LGVM proficiency at estimating the distribution of the observed-data we
proposed looking at the gap between moment estimators for both the first and
the second moment.

###  4.3 Additional justifications

For any generative model, it is always possible to generate a new sample of
$G$ points say $S_{LVGM}=\\{\tilde{x_{1}},...,\tilde{x_{G}}\\}$ and a simple
model estimator for the first and second moment would be to simply compute
$\frac{1}{G}\sum_{i=1}^{G}(\tilde{x_{i}})$ for the first moment and similarly
$\frac{1}{G}\sum_{i=1}^{G}(\tilde{x_{i}}\tilde{x_{i}}^{T})$, let us call those
sample estimators (SE). These estimators could replace both FMEs defined
previously as they also reflect the distribution learned by the LVGM.
Conceptually, the SEs are simple and easy to use and thus we want to justify
why our estimators (FMEs) are better.

To begin, FMEs are faster to compute by a slight margin. In order to draw the
same sample size for the Monte Carlo estimates, say $m$ sample, we need $m$
sample for the FMEs. However, for the SEs have to draw twice as many samples
($2m$), since we must sample from $p(\mathbf{z})$ $m$ times and then from
$p(\mathbf{x}|\mathbf{z})$ an additional $m$ times. However, more importantly,
$M<<D$, this means that not only we FME requires half as many samples, but
they are samples from a much lower dimension distribution which further
increases the difference in computational cost between the FME and the SE.

Next, our estimator has a smaller variance for both the first and the second
moment. A complete proof is located in the appendices.

Finally, we proceeded with a simulation that demonstrates which estimators is
closer to the true LVGM moments. For simple models, such as GMMs, we can
compute analytically $\textbf{E}_{z}[\textbf{Var}_{x}(\mathbf{x}|\mathbf{z})+$
$\textbf{E}_{x}(\mathbf{x}|\mathbf{z})^{2}]$, the LVGM second moment. We
compared the gap between the true LVGM second moment and the FME to the gap
between the true LVGM second moment and the SE and plot both against $m$ the
number of Monte Carlo samples. We can see in figure 2 that the estimator we
proposed (FME) is overall closer to the true LVGM second moment than the SE
while being faster to compute.

Figure 2: The evolution of both gap plotted against $m$. This gap is computed
using the Frobenius norm as justified in section 5.

##  5 MEGA: a new metric for comparing models

For latent variable models as defined in section 2 we assess the ability of
the model to capture the moment of the data set $S$ by comparing the DE with
the FME. Since we are looking at the difference between two different moment
estimators of the same value, we named this metric the Moment Estimators Gap
(MEGA). We study the difference between our two second moment estimators and
we refer to this matrix as 2MEGA. Similarly, we refer to the vector
representing the first moment estimator gap as 1MEGA.

###  5.1 Selecting a matrix norm

In order to make the MEGA tangible and comparable we propose to use a matrix
norm of the MEGA as our metric. There exist a wide range of possible
candidates, let us introduce a few and justify our finale choice.

We want to use a norm that looks at the global properties of the matrix and
fortunately Rigollet [24] introduces and studies the behaviour of well-
established matrix norms. To begin we use norms inspired by vector norms.
Given the vector $v$, assume $|v|_{q}$ is the following vector norm

$\displaystyle|v|_{q}=(\sum_{i}|v_{i}|^{q})^{(1/q)},$ (14)

and given the matrix $M$, its matrix equivalent is $|M|_{q}$

$\displaystyle|M|_{q}=(\sum_{ij}|M_{ij}|^{q})^{(1/q)}.$ (15)

When $q=2$, this is a special case call the Frobenius norm

$\displaystyle|M|_{2}=|M|_{F}=(\sum_{ij}|M_{ij}|^{2})^{(1/2)}=\sqrt{\text{Tr}(M^{T}M)},$
(16)

where Tr is the trace operator that sums the elements of the diagonal of the
input matrix.

This norm is also a member of the Schatten $q$-norms (for $q=2$) which is a
family of matrix norms defined as the vector norm of Equation 14 for the
singular values of the matrix. Since we work with second moment estimators,
our matrix $M$ is a squared matrix (dimension $m\times m$) and symmetric and
thus the singular values of $M$ are equal to its eigenvalues. We identify the
vector of eigenvalues as $\lambda$. Consequently, the Schatten $q$-norm for
matrix $M$ is $||M||_{q}=|\lambda|_{q}$.

Another member of this family is considered, when $q=\infty$ we define
$||M||_{\infty}=\lambda_{max}=||M||_{op}$ and this is referred to as the
operator norm.

In random matrix theory and in matrix estimation, these two norms appear
frequently and are consequently well known in the statistical research
community. For instance, in covariance matrix estimation it is possible, under
mild assumptions, to bound the operator norm of the difference between the
true covariance matrix and simple estimators [24]. Because of the popularity
and the known properties of both the Frobenius norm and the operator norm,
they are both legitimate options to measure the MEGA.

The bigger the 2MEGA is, the further away our model second moment is from the
data second moment. However, one can perceive the Frobenius norm as the length
of the hypotenuse of a multidimensional triangle whose sides are given by the
eigenvectors of the matrix while the operator norm is the length of the
longest cathetus of the same multidimensional triangle. For that reason and
for computational reasons expressed in the next section, we prefer MEGA-F and
this will serve as our metric in the sections to come. In other words, the
metric we propose to evaluate the quality of the fit for the second moment is

2MEGA-F $\displaystyle=|\text{2MEGA}|_{F}=|\text{DE-FME}|_{F}.$ (17)

Additionally, if we also consider 1MEGA, we can then again simply use the
Frobenius norm (the vector q-norm with $q=2$) on the vector 1MEGA

$\displaystyle\text{1MEGA-F}=|\text{1MEGA}|_{2}.$ (18)

###  5.2 Implementation of the selected norm

Implementing the Frobenius norm is quite straight forward. To compute to
Frobenius norm, we need square of the MEGA , $m^{2}$ operations, and then sum
its diagonal which results in $m^{2}+m$ operations. The largest eigenvalue of
a matrix can be estimated with the von Mises iteration [21], where each of the
$p$ iterations will require $m^{2}$ operations resulting in $pm^{2}$
operations in total. Consequently, because we can implement an exact
computation of the Frobenius norm and because it is computationally faster (as
soon as $p\geq 2$) we strictly consider the Frobenius norm for the rest of
this article. However, there could be some merit to exploring the operator
norm as well in future work.

We implemented the vector Frobenius norm and the matrix Frobenius norm in
Python using the NumPy library [11] and using the Pytoch library [22]. We also
implemented a MEGA function that takes as input a sample of
$\textbf{Var}_{x}(\mathbf{x}|\mathbf{z})$ and
$\textbf{E}_{x}(\mathbf{x}|\mathbf{z})$ alongside the data set $S$ and returns
1MEGA-F and 2MEGA-F. These implementations are publicly available on the
author’s GitHub [5].

##  6 Experiments : comparing models

Now that we established the metric, let us use it to compare two distributions
$\hat{p}(x)$ learned from two different models.

For this demonstration we use simple 2-dimensional observations. We
demonstrated that large MEGAs are associated with poor fit when visualizing
generated data from LVGMs. In other words, our metric is concordant with the
currently used techniques but it provides more precise quantification of the
mismatch.

(a) Model 1

(b) Model 2

Figure 3: Three-cluster data set | Model 1 | Model 2
---|---|---
1MEGA-F | 0.08100187 | 0.02810075
2MEGA-F | 1.9527672 | 0.320804925
Table 1: MEGA for both model when trained on the three-cluster data set.

Figure 3 contains both the training observations (in blue) and generated data
points (in red) from both tested LVGM models, Model 1 and Model 2. By looking
at Figure 3, Model 2 seems to have a better fit, mostly because it seems like
it captures the true variance much better.

We computed 1MEGA-F and 2MEGA-F for both trained distributions, the results
are containing in Table 1. The results agree with our intuitive visualisation
of the situation. Though the distribution learned by fitting Model 2 has a
first moment that matches better than training observations than the
distribution learning by fitting Model 1, this difference is much more drastic
when comparing the second moment of both trained distributions according to
our metric.

We have repeated this small demonstration with a second data set as illustrate
in figure 4.

(a) Model 1

(b) Model 2

Figure 4: Moon data sets | Model 1 | Model 2
---|---|---
1MEGA-F | 0.01490262 | 0.00728456
2MEGA-F | 0.07236395 | 0.00862108
Table 2: MEGA for both model when trained on the moon data set.

We observe similar results using the moon data set as illustrated in Figure 4
and table 2; the MEGA metric favours the model that visually looks better.

##  7 MEGA for regularization

As we previously mentioned, based on the principle of Goodhart’s Law, we
proposed to use the proposed metrics 1MEGA-F and 2MEGA-F to compare complex
LVGMs such as IAF-VAEs [17] or NVAEs[30]. In section 3 we discussed some
issues with the likelihood as an evaluation metric in some special cases.
Additionally, it is unfair to use likelihood as a metric to compare two models
when one is trained using a likelihood approach and the other is not.

However, if this metric has some merit when comparing model, why could we not
simply included in the optimization process ? In this section we address the
effect of incorporating MEGA as part of the optimization process.

Because the likelihood and the moments reflect different aspects of the
distribution $p(\mathbf{x})$ then incorporating MEGA as part of the
optimization process for models trained by maximum likelihood acts as
regularization.

This is easy to see for models with Gaussian emission distribution such as
GMMs and VAEs. To simplify the explanation we will focus on the GMM case for
now. When fitting a single Gaussian distribution with maximum likelihood, we
set its parameters $\mu$ and $\sigma$ to the data mean and the data standard
error and thus a single Gaussian performs very well according to MEGAs. As we
increase the number of components $k$ in a GMM, we increase its likelihood but
we also increase its MEGA. This means that we can use the MEGA as a model
selection (or regularization metric) to fit GMM in place of alternatives such
as AIC and BIC in order to get a model that balances likelihood and moment
matching.

One way to do so is to select the model that maximizes:

$\displaystyle\text{ll}(S)-\alpha(\text{1MEGA-F}+\sqrt{\text{2MEGA-F}}),$ (19)

where $\text{ll}(S)$ is the log-likelihood of the data under the fitted model
and $\alpha$ is a hyper-parameter that we manually fixed in the experiments of
section 8. Basically in this particular case, the $\alpha$ parameters
interpolated between maximum-likelihood and moment estimator. Additionally, we
can draw an entire regularization path for different values of $\alpha$,
similar to what is produced when using Lasso [28, 26, 9].

Similarly, we can use a MEGA term in the objective function while training a
VAE. Our proposed metric could again serve as regularization in order to
ensure the VAE model captures those moments by maximizing:

$\displaystyle\textbf{E}_{q}[\ln
p(\mathbf{x}|\mathbf{z})]-KL(q(\mathbf{z}|\mathbf{x})|p(\mathbf{z}))-\alpha(\text{1MEGA-F}+\sqrt{\text{2MEGA-F}}).$
(20)

Arguments have been made in the past [16, 14] that the KL divergence term
serves as a regularizer and that adding a $\beta$ parameter can provide a way
to adjust the strength of the regularization. Adding a MEGA term to the ELBO
should provide a different type of regularization. The KL divergence term
provides regularization for the distribution $q(\mathbf{z}|\mathbf{x})$ and
the MEGA term provides regularization for $p(\mathbf{x})$:

$\displaystyle\underset{\text{Reconstruction error}}{\textbf{E}_{q}[\ln
p(\mathbf{x}|\mathbf{z})]}-\underset{\text{Regularization
for}q(\mathbf{z}|\mathbf{x})}{\beta
KL(q(\mathbf{z}|\mathbf{x})|p(\mathbf{z}))}-\underset{\text{Regularization
for}p(\mathbf{x})}{\alpha(\text{1MEGA-F}+\sqrt{\text{2MEGA-F}})}.$ (21)

As previously demonstrated [6], tunning a VAE to obtain both good
reconstruction and generative performances is difficult. For instance, when
using the $\beta$-VAE [14, 13] we can make the hyper-parameter $\beta$
arbitrary small and obtain close to perfect reconstruction but this leads to
poor generative performance due to the distribution $q(z)$ used for training
being widely different from $p(z)$ used for generation. The additional term we
add to the VAE objective function should help with generative perspective of
VAEs.

Thus, to understand the benefits of integrating this additional term to the
objective function we have to understand how different it is from the first
term. The first term measure the reconstruction error because the likelihood
$p(x|z)$ is computed on latent variable $z$ sampled from $q(z|x)$. In
contrast, the MEGA penalty term evaluate the first and second moment of $p(x)$
independently from $q(z|x)$ but rather based on the generative model where
$p(x)=\int_{z}p(\mathbf{x}|z)p(\mathbf{z})dz$. We believe this additional
constraint in the objective function should lead to samples that matches the
true distribution of $p(x)$ more closely.

##  8 Experiments: regularization

As discussed in section 7 another way to integrate such metric in a workflow
is to regularize LVGMs with Gaussian observation distributions. In the
following experiments, we explore the use of our proposed metric as a
regularizer for GMMs and VAEs.

###  8.1 Regularizer for GMMs

For GMMs, increasing the number of components increases the likelihood and
thus it is necessary to regularize this model; we cannot simply select the
ideal number of components based on the likelihood alone. We have generated a
simple data set from a 3-component GMM. Figure 5 is a scatter of the simulated
data, by looking at the figure we would like the model selection procedure to
settle on three components.

Figure 5: Scatter plot of the simulated data set.

The Akaike Information Criteria (AIC) [2, 1] is a well-establish penalized
likelihood function that can be used to select the number of components in a
GMM. The AIC can be expressed as

$\displaystyle\text{AIC}=2p-2\ln(L),$ (22)

where $L$ is the likelihood of the model, and $p$ is the number of parameters
($p=2k$ in the GMM case). In equation 22 we see a natural tradeoff; $-2\ln(L)$
goes down as we increase the number of components but $2p$ goes up. To select
the right GMM model, we can fit multiple GMM with various number of components
and we select the model that minimizes the AIC.

Figure 6: AIC plotted against the number of components for trained GMM models
(lower is better).

Based on Figure 6, the AIC is minimized at $k=3$ which is concordant the model
used to generate the data.

Next, we use our MEGA-penalized likelihood function described in equation 19
to select the number of components.

Figure 7: MEGA-penalized likelihood plotted against the number of components
for trained GMM models (higher is better)

As observed above, when using our proposed metric for GMM regularization we
also obtain a function that favours 3 components, which is concordant with
both the model used to generate the data and the AIC. However, in this
instance the flexibility given by the hyper-parameter $\alpha$ is both a
benefit and a drawback. Using the hyper-parameter, we can easily adjust how
strong we want the penalty to be but on the flip side there are no automated
ways to adjust it. As previously mentioned, we can build a regularization path
for different values of $\alpha$ leading to a sequence of possible models from
which the user can select one.

###  8.2 Regularizer for VAEs

Finally we experiment with using MEGA as a regularization for VAEs. Now that
we used MEGA as part of the objective function we can no longer use it to
assess the quality of the samples so we are forced to visually inspect samples
and their parameters.

We run this demonstration on a subset of the MNIST data set [18] that contains
only the digit four. We have experimented with a wide range of parameters
$\alpha$ and $\beta$, as defined in equation 21 and fixed the hyper-parameters
to the values leading to the more realistic-looking images.

Figures 8,9 and 10 illustrate our results. The images on in the left column
were produced by a VAE trained without MEGA and the right one with MEGA. We
have included images of a sample, its mean and its standard deviation.

(a) Model train without MEGA

(b) Model train with MEGA

Figure 8:  A sample of 64 images from
$p_{\theta}(\mathbf{x}\|\mathbf{z})=N(\mu(\mathbf{z}),\sigma(\mathbf{z}))$
where $\mathbf{z}\sim N(0,1)$.

(a) Model train without MEGA

(b) Model train with MEGA

Figure 9:  The 64 sampled mean of the images: $\mu(\mathbf{z})$ where
$\mathbf{z}\sim N(0,1)$.

(a) Model train without MEGA

(b) Model train with MEGA

Figure 10:  The 64 sampled standard deviation for each pixel of the images:
$\sigma(\mathbf{z})$ where $\mathbf{z}\sim N(0,1)$. For those images, the
whiter the pixel is the larger the standard deviation of that pixel is.

The samples are still grainy, and this is due to the pixel-independence
assumption of the simple VAE which only learns individual pixel variance and
not a full-covariance matrix. However, the images generated by a VAE trained
with the MEGA penalty seems to have sharper edges; the contrast between the
white of the digit and the black of the exterior is bigger. This is noticeable
in all three figures but mostly in figure 10. It seems like the MEGA penalty
changes rather drastically $\sigma(\mathbf{z})$, the variance of the observed
distribution. Both the exterior of the digit and the inner part of the digit
are darker than the model trained without MEGA. The model trained with MEGA
perceives the contour and extremities of the digit as the location with higher
variability which is not as prevalent for the model trained without MEGA.

We also notice a bigger variability in the shapes of fours that we are able to
generate, mostly visible in figures 8 and 9. The lack of variability in the
left column is most likely caused by posterior collapse [20, 6], a common
problem when using $\beta$-VAEs that the MEGA penalty seems to counteract.

##  9 Discussion

As a comparative metric, MEGA is fast to compute and easy to interpret. The
larger 1MEGA is, the larger is the gap between the first moment of trained
distribution $p(x)$ and the empirical first moment of the data set $S$.
Similarly, the larger 2MEGA is, the larger the gap is for the second moment.
This can give us a quick and easy to compare and evaluate the quality of
fitted models.

However, there are still some limitations to this approach. The most obvious
is that the formulation we propose only allow us to quickly evaluate the gap
for the first two moments. This leads to an incomplete comparison of the
learned and empirical distribution which can create some problems in niche
cases. A example of this problem is the high performance of the simple
Gaussian distribution. Usually, when fitting a Gaussian distribution to a data
set we set the parameters $\mu$ and $\sigma$ to be the empirical mean and the
empirical, thus the Gaussian mixture with a single component shows very good
results (low MEGA).

Fortunately, we have designed this metric for complex latent variable model
and there are no reasons to use it when assessing the fit of a single
Gaussian. Additionally, other comparative strategies discussed in section 3
can still be used in parallel of our proposed one.

We also want to point out that if a trained distribution has low MEGA but bad-
looking generated samples this still provides us with insightful information.
It indicates that the problem is in the LVGM distribution’s higher moments and
our metric was able to provide us with that information very quickly.

For regularization applications, we were able to successfully used MEGA with
two different LVGMs, GMMs and VAEs. Though in both these cases we managed to
achieve good results, selecting the appropriate constraint using the hyper-
parameters was not an easy task and had to be done heuristically. This is
certainly a weakness that should be addressed in future work.

Let us now discuss some of ongoing work. We are currently working on a
generalization of MEGA. We want to extend our metric not only to higher
moments of $\mathbf{x}$ but to any functions $g(\mathbf{x})$. This would not
only allow us to compare the skewness of the trained model compared to the
data but also more complicated properties of distributions such as
multimodality. This generalized MEGA would provide a more complete evaluation
of models than the currently proposed metric.

Because we believe our theoretical arguments to be quite solid, the
experiments above are a quick showcase of possible applications of the MEGA
metric, they serve as a proof of concept rather than exhaustive tests.
However, some components of our experiments were done heuristically and more
thorough test should be done to broaden our understanding of the possible
applications of our metric.

##  10 Conclusion

In this article we have introduced a fairly simple, easy and computationally
fast way to check the generative model’s distribution of a large class of
LVGMs. This metric, MEGA, evaluates the gap between the LVGM distribution
first and second moment and the training (or validation) data set first and
second moment.

The premise of the proposed metric is theoretically simple and quite
intuitive. Both the DE and the FME are unbiased estimators of both the first
and the second moment and if a gap exists between them then the LVGM
distribution does not match key aspects of the data distribution.

To support our theoretical arguments, we have demonstrated how to use this
metric for two different purposes. First, as an evaluation metric that can
replace the more heuristics approaches that rely on eyeballing generated
samples. Second, since this metric is currently available for the first two
moments, it favours a simple model, such a single Gaussian, and thus can be
used as regularization for models such as GMMs and VAEs.

However, we believe we have only scratched the surface of all of the
applications and ways to incorporate these moment-gap-based metrics in model
fitting and model selection pipelines. We hope to make further progress in
this direction in future work. Another future work direction is to extend
these moment gap estimators to sequential LVGMs, such as Hidden Markov Models
and State-Space Models. Finally, the biggest improvement we could work on is
to extend the moment estimators to higher moments; this would make the
evaluation metrics much more valuable.

Nonetheless, we believe this work is a first step in data-driven automated
model selection and we hope it inspires similar contributions.

## Ackowledgement

I would like to thank Yanbo Tang for its help when selecting the matrix norm
and providing sketches for the proofs that FMEs have smaller variance than
SEs. I would also like to thank David Duvenaud, Michael Lalancette, Jeffrey S.
Rosenthal, Anthony Coache and Renaud Alie for their insightful comments.

I would like to acknowledge the Ontario Graduate Scholarship (OGS) program and
the University of Toronto School of Graduate Studies for their financial
support during the early parts of this project and the Canadian Statistical
Sciences Institute (CANSSI) for the financial support during the later parts.

## Appendices

### Proof that the FME has small variance than the SE

#### First moment

For the SE, we first sample $z_{1},...z_{m}$, a set of $m$ latent variables
where $z\sim p(\mathbf{z})$ then sample $x_{1},...,x_{m}$ with $x_{i}\sim
p(\mathbf{x}|\mathbf{z}_{i})$; this implies that $x\sim
p(\mathbf{x}|\mathbf{z})p(\mathbf{z})=p(\mathbf{x})$. Then, the SE for the
first moment is $\frac{1}{m}\sum_{i=1}^{m}\mathbf{x}_{i}$.

$\displaystyle\textbf{Var}_{x}(\text{SE})=\textbf{Var}_{x}\left(\frac{\sum_{i=1}^{m}x_{i}}{m}\right)=\frac{1}{m^{2}}\sum_{i=1}^{m}\textbf{Var}_{x}(x_{i})=\frac{\textbf{Var}_{x}(\mathbf{x})}{m}.$
(23)

When building our FME, we sample $z_{1},...z_{m}$, a set of $m$ latent
variables where $z\sim p(\mathbf{z})$ and then compute
$\frac{1}{m}\sum_{i=1}^{m}\textbf{E}_{x}(\mathbf{x}|\mathbf{z}=z_{i})$.

$\displaystyle\textbf{Var}_{z}(\text{FME})=\textbf{Var}_{z}\left(\frac{\sum_{i=1}^{m}\textbf{E}_{x}(\mathbf{x}|z_{i})}{m}\right)=\frac{1}{m^{2}}\sum_{i=1}^{m}\textbf{Var}_{z}(\textbf{E}_{x}(\mathbf{x}|z_{i}))=\frac{\textbf{Var}_{z}[\textbf{E}_{x}(\mathbf{x}|z_{i})]}{m}.$
(24)

Then using the Law of Total Variance we have that:

$\displaystyle\textbf{Var}_{x}(\mathbf{x})$
$\displaystyle=\textbf{E}_{z}[\textbf{Var}_{x}(\mathbf{x}|\mathbf{z})]+\textbf{Var}_{z}[\textbf{E}_{x}(\mathbf{x}|\mathbf{z})]$
$\displaystyle\geq\textbf{Var}_{z}[\textbf{E}_{x}(\mathbf{x}|\mathbf{z})]$
$\displaystyle\Rightarrow\frac{\textbf{Var}_{x}(\mathbf{x})}{m}$
$\displaystyle\geq\frac{\textbf{Var}_{z}[\textbf{E}_{x}(\mathbf{x}|z_{i})]}{m}$
(25) $\displaystyle\Rightarrow\textbf{Var}(\text{SE})$
$\displaystyle\geq\textbf{Var}(\text{FME}),$

and thus our FME has lower variance than the commonly used alternative.

#### Second moment

For the SE, we first sample $z_{1},...z_{m}$, a set of $m$ latent variables
where $z\sim p(\mathbf{z})$ then sample $x_{1},...,x_{m}$ with $x_{i}\sim
p(\mathbf{x}|\mathbf{z}_{i})$; this implies that $x\sim
p(\mathbf{x}|\mathbf{z})p(\mathbf{z})=p(\mathbf{x})$. Then, the SE for the
second moment is $\frac{1}{m}\sum_{i=1}^{m}\mathbf{x}_{i}^{2}$.

$\displaystyle\textbf{Var}_{x}(\text{SE})=\textbf{Var}_{x}\left(\frac{\sum_{i=1}^{m}x_{i}^{2}}{m}\right)=\frac{\textbf{Var}_{x}(\mathbf{x}^{2})}{m}.$
(26)

When building our FME, we sample $z_{1},...z_{m}$, a set of $m$ latent
variables where $z\sim p(\mathbf{z})$ and then compute
$\frac{1}{m}\sum_{i=1}^{m}\left(\textbf{Var}_{x}(\mathbf{x}|\mathbf{z}=z_{i})+[\textbf{E}_{x}(\mathbf{x}|\mathbf{z}=z_{i})]^{2}\right)$.

$\displaystyle\textbf{Var}_{z}(\text{FME})$
$\displaystyle=\textbf{Var}_{z}\left(\frac{\sum_{i=1}^{m}\left(\textbf{Var}_{x}(\mathbf{x}|\mathbf{z}=z_{i})+[\textbf{E}_{x}(\mathbf{x}|\mathbf{z}=z_{i})]^{2}\right)}{m}\right)$
$\displaystyle=\frac{\textbf{Var}_{z}\left(\textbf{Var}_{x}(\mathbf{x}|\mathbf{z}=z_{i})+[\textbf{E}_{x}(\mathbf{x}|\mathbf{z}=z_{i})]^{2}\right)}{m}.$
(27)

Now let’s take a closer look at the numerator of equation (27).

$\displaystyle\textbf{Var}_{z}\left(\textbf{Var}_{x}(\mathbf{x}|\mathbf{z}=z_{i})+[\textbf{E}_{x}(\mathbf{x}|\mathbf{z}=z_{i})]^{2}\right)$
$\displaystyle=\textbf{Var}_{z}\left(\textbf{E}_{x}[\mathbf{x}^{2}|\mathbf{z}]-\textbf{E}_{x}[\mathbf{x}|\mathbf{z}]^{2}+\textbf{E}_{x}[\mathbf{x}|\mathbf{z}]^{2}\right)$
$\displaystyle=\textbf{Var}_{z}\left(\textbf{E}_{x}[\mathbf{x}^{2}|\mathbf{z}]-\textbf{E}_{x}[\mathbf{x}|\mathbf{z}]^{2}+\textbf{E}_{x}[\mathbf{x}|\mathbf{z}]^{2}\right)$
$\displaystyle=\textbf{Var}_{z}\left(\textbf{E}_{x}[\mathbf{x}^{2}|\mathbf{z}]\right)$
$\displaystyle=\textbf{E}_{z}\left(\textbf{E}_{x}[\mathbf{x}^{2}|\mathbf{z}]^{2}\right)-\textbf{E}_{z}\left(\textbf{E}_{x}[\mathbf{x}^{2}|\mathbf{z}]\right)^{2}$
(28)
$\displaystyle=\textbf{E}_{z}\left(\textbf{E}_{x}[\mathbf{x}^{2}|\mathbf{z}]^{2}\right)-\textbf{E}_{x}\left(\mathbf{x}^{2}\right)^{2}$
$\displaystyle\Rightarrow\textbf{Var}(\text{FEM})$
$\displaystyle=\frac{\textbf{E}_{z}\left(\textbf{E}_{x}[\mathbf{x}^{2}|\mathbf{z}]^{2}\right)-\textbf{E}_{x}\left(\mathbf{x}^{2}\right)^{2}}{m}$

Finally, let us apply the Law of Total Variance to $\mathbf{x}^{2}$:

$\displaystyle\textbf{Var}_{x}(\mathbf{x}^{2})$
$\displaystyle=\textbf{E}_{z}(\textbf{Var}_{x}[\mathbf{x}^{2}|\mathbf{z}])+\textbf{Var}_{z}(\textbf{E}_{x}[\mathbf{x}^{2}|\mathbf{z}])$
$\displaystyle=\textbf{E}_{z}(\textbf{Var}_{x}[\mathbf{x}^{2}|\mathbf{z}])+\textbf{E}_{z}\left(\textbf{E}_{x}[\mathbf{x}^{2}|\mathbf{z}]^{2}\right)-\textbf{E}_{z}(\textbf{E}_{x}[\mathbf{x}^{2}|\mathbf{z}])^{2}$
$\displaystyle=\textbf{E}_{z}(\textbf{Var}_{x}[\mathbf{x}^{2}|\mathbf{z}])+\textbf{E}_{z}\left(\textbf{E}_{x}[\mathbf{x}^{2}|\mathbf{z}]^{2}\right)-\textbf{E}_{x}[\mathbf{x}^{2}]^{2}$
(29) $\displaystyle\Rightarrow\textbf{Var}_{x}(\mathbf{x}^{2})$
$\displaystyle\geq\textbf{E}_{z}\left(\textbf{E}_{x}[\mathbf{x}^{2}|\mathbf{z}]^{2}\right)-\textbf{E}_{x}[\mathbf{x}^{2}]^{2}$
$\displaystyle\Rightarrow\frac{\textbf{Var}_{x}(\mathbf{x}^{2})}{m}$
$\displaystyle\geq\frac{\textbf{E}_{z}\left(\textbf{E}_{x}[\mathbf{x}^{2}|\mathbf{z}]^{2}\right)-\textbf{E}_{x}[\mathbf{x}^{2}]^{2}}{m}$
$\displaystyle\Rightarrow\textbf{Var}(\text{SE})$
$\displaystyle\geq\textbf{Var}(\text{FME})$

## References

  * [1] Hirotogu Akaike.  Information theory and an extension of the maximum likelihood principle.  In Selected papers of hirotugu akaike, pages 199–213. Springer, 1998.
  * [2] Hirotugu Akaike.  A new look at the statistical model identification.  IEEE transactions on automatic control, 19(6):716–723, 1974.
  * [3] Anima Anandkumar, Dean P Foster, Daniel J Hsu, Sham M Kakade, and Yi-kai Liu.  A spectral algorithm for latent dirichlet allocation.  In F. Pereira, C. J. C. Burges, L. Bottou, and K. Q. Weinberger, editors, Advances in Neural Information Processing Systems, volume 25. Curran Associates, Inc., 2012.
  * [4] Animashree Anandkumar, Rong Ge, Daniel Hsu, Sham M. Kakade, and Matus Telgarsky.  Tensor decompositions for learning latent variable models.  Journal of Machine Learning Research, 15(80):2773–2832, 2014.
  * [5] Cedric Beaulac.  Mega for lvgm.  https://github.com/CedricBeaulac/LVM_Analysis, 2021.
  * [6] Cédric Beaulac.  Performance and accessibility of statistical learning algorithms for applied data analysis.  PhD thesis, University of Toronto, 2021.
  * [7] C. M. Bishop.  Pattern Recognition and Machine Learning.  Springer, 2007.
  * [8] Arun Tejasvi Chaganty and Percy Liang.  Estimating latent-variable graphical models using moments and likelihoods.  In International Conference on Machine Learning, pages 1872–1880. PMLR, 2014.
  * [9] Jerome Friedman, Trevor Hastie, and Robert Tibshirani.  Regularization paths for generalized linear models via coordinate descent.  Journal of Statistical Software, 33(1):1–22, 2010.
  * [10] Charles AE Goodhart.  Problems of monetary management: the uk experience.  In Monetary theory and practice, pages 91–121. Springer, 1984.
  * [11] Charles R. Harris, K. Jarrod Millman, St’efan J.van der Walt, Ralf Gommers, Pauli Virtanen, David Cournapeau, Eric Wieser, Julian Taylor, Sebastian Berg, Nathaniel J. Smith, Robert Kern, Matti Picus, Stephan Hoyer, Marten H. van Kerkwijk, Matthew Brett, Allan Haldane, Jaime Fern’andez del R’ıo, Mark Wiebe, Pearu Peterson, Pierre G’erard-Marchant, Kevin Sheppard, Tyler Reddy, Warren Weckesser, Hameer Abbasi, Christoph Gohlke, and Travis E. Oliphant.  Array programming with NumPy.  Nature, 585(7825):357–362, September 2020.
  * [12] T. Hastie, R. Tibshirani, and J. Friedman.  The Elements of Statistical Learning.  Springer, 2 edition, 2009.
  * [13] Irina Higgins, David Amos, David Pfau, Sebastien Racaniere, Loic Matthey, Danilo Rezende, and Alexander Lerchner.  Towards a definition of disentangled representations.  arXiv preprint arXiv:1812.02230, 2018.
  * [14] Irina Higgins, Loic Matthey, Arka Pal, Christopher Burgess, Xavier Glorot, Matthew Botvinick, Shakir Mohamed, and Alexander Lerchner.  beta-vae: Learning basic visual concepts with a constrained variational framework.  ICLR, 2(5):6, 2017.
  * [15] D. P Kingma and M. Welling.  Auto-Encoding Variational Bayes.  ArXiv e-prints, December 2013.
  * [16] Diederik P. Kingma.  Variational Inference & Deep Learning : A New Synthesis.  PhD thesis, Universiteit van Armsterdam, 10 2017.
  * [17] Durk P Kingma, Tim Salimans, Rafal Jozefowicz, Xi Chen, Ilya Sutskever, and Max Welling.  Improved variational inference with inverse autoregressive flow.  Advances in neural information processing systems, 29:4743–4751, 2016.
  * [18] Y. Lecun, L. Bottou, Y. Bengio, and P. Haffner.  Gradient-based learning applied to document recognition.  Proceedings of the IEEE, 86(11):2278–2324, 1998.
  * [19] Yujia Li, Kevin Swersky, and Rich Zemel.  Generative moment matching networks.  In International Conference on Machine Learning, pages 1718–1727. PMLR, 2015.
  * [20] James Lucas, George Tucker, Roger B Grosse, and Mohammad Norouzi.  Don’t blame the elbo! a linear vae perspective on posterior collapse.  In Advances in Neural Information Processing Systems, pages 9408–9418, 2019.
  * [21] RV Mises and Hilda Pollaczek-Geiringer.  Praktische verfahren der gleichungsauflösung.  ZAMM-Journal of Applied Mathematics and Mechanics/Zeitschrift für Angewandte Mathematik und Mechanik, 9(1):58–77, 1929.
  * [22] Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Kopf, Edward Yang, Zachary DeVito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala.  Pytorch: An imperative style, high-performance deep learning library.  In H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett, editors, Advances in Neural Information Processing Systems 32, pages 8024–8035. Curran Associates, Inc., 2019.
  * [23] Anastasia Podosinnikova.  On the method of moments for estimation in latent linear models.  PhD thesis, PSL Research University, 2016.
  * [24] Phillippe Rigollet and Jan-Christian Hütter.  High dimensional statistics.  Lecture notes for course 18S997, 813:814, 2015.
  * [25] Jeffrey S. Rosenthal.  Monte carlo methods.  Lecture notes for course STA2431, 2019.
  * [26] Noah Simon, Jerome Friedman, Trevor Hastie, and Rob Tibshirani.  Regularization paths for cox’s proportional hazards model via coordinate descent.  Journal of Statistical Software, 39(5):1–13, 2011.
  * [27] Marilyn Strathern.  ‘improving ratings’: audit in the british university system.  European review, 5(3):305–321, 1997.
  * [28] Robert Tibshirani.  Regression shrinkage and selection via the lasso.  Journal of the Royal Statistical Society: Series B (Methodological), 58(1):267–288, 1996.
  * [29] Michael E Tipping and Christopher M Bishop.  Probabilistic principal component analysis.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 61(3):611–622, 1999.
  * [30] Arash Vahdat and Jan Kautz.  Nvae: A deep hierarchical variational autoencoder.  arXiv preprint arXiv:2007.03898, 2020.
  * [31] Shengjia Zhao, Jiaming Song, and Stefano Ermon.  Infovae: Balancing learning and inference in variational autoencoders.  In Proceedings of the aaai conference on artificial intelligence, volume 33, pages 5885–5892, 2019.


# Effects of Neutral Hydrogen on Cosmic Ray Precursors in Supernova Remnant
Shock Waves

John C. Raymond,11affiliation: Harvard-Smithsonian Center for Astrophysics, 60
Garden St., Cambridge, MA 02138, USA<EMAIL_ADDRESS>J.
Vink,22affiliation: Astronomical Institute, Utrecht University, P.O. Box
80000, 3508TA Utrecht, The Netherlands E.A. Helder,22affiliation: Astronomical
Institute, Utrecht University, P.O. Box 80000, 3508TA Utrecht, The Netherlands
& A. de Laat22affiliation: Astronomical Institute, Utrecht University, P.O.
Box 80000, 3508TA Utrecht, The Netherlands

###### Abstract

Many fast supernova remnant shocks show spectra dominated by Balmer lines. The
H$\alpha$ profiles have a narrow component explained by direct excitations and
a thermally Doppler broadened component due to atoms that undergo charge
exchange in the post-shock region. However, the standard model does not take
into account the cosmic-ray shock precursor, which compresses and accelerates
plasma ahead of the shock. In strong precursors with sufficiently high
densities, the processes of charge exchange, excitation and ionization will
affect the widths of both narrow and broad line components. Moreover, the
difference in velocity between the neutrals and the precursor plasma gives
rise to frictional heating due to charge exchange and ionization in the
precursor. In extreme cases, all neutrals can be ionized by the precursor.

In this paper we compute the ion and electron heating for a wide range of
shock parameters, along with the velocity distribution of the neutrals that
reach the shock. Our calculations predict very large narrow component widths
for some shocks with efficient acceleration, along with changes in the broad-
to-narrow intensity ratio used as a diagnostic for the electron-ion
temperature ratio. Balmer lines may therefore provide a unique diagnostic of
precursor properties. We show that heating by neutrals in the precursor can
account for the observed H$\alpha$ narrow component widths, and that the
acceleration efficiency is modest in most Balmer line shocks observed thus
far.

shock waves — acceleration of particles — ISM: supernova remnants

##  1 Introduction

Cosmic rays are widely believed to originate in supernova remnant (SNR) shock
waves, because the cosmic-ray energy spectrum agrees with model predictions,
because power-law distributions of energetic electrons are seen in SNRs, and
because the power required to maintain the cosmic ray population could be
supplied by about 10% of kinetic energy of Galactic supernovae. The standard
theory for the process is Diffusive Shock Acceleration (DSA), which is a first
order Fermi process requiring that particles scatter between a gasdynamic
subshock and plasma turbulence in a shock precursor. Evidence for non-linear
DSA comes from curved synchrotron spectra (Reynolds & Ellison, 1992; Vink et
al., 2006; Allen et al., 2008), evidence for high compression factors (Warren
et al., 2005; Cassam-Chenaï et al., 2008) and evidence for lower than expected
downstream temperatures (Hughes et al., 2000; Helder et al., 2009). However,
all this evidence is based on observation of downstream properties. The
effects of precursor physics on the H$\alpha$ emission described here offer a
direct probe of the properties of the precursor.

A crucial parameter for these models is the diffusion coefficient $\kappa$, as
it determines the precursor scale length, which is typically $\kappa$ divided
by the shock speed $V_{S}$. Gas is compressed in the precursor and accelerated
to a fraction of the shock speed, and this compression is related to $V_{S}$,
to the efficiency of particle acceleration and to the escape of energetic
particles from the region (Bykov, 2005; Vink et al., 2010). Neutrals can
impede the acceleration process by damping the turbulence needed to scatter
particles back to the shock. However, Drury et al. (1996) found that the
acceleration efficiency can be high as long as the density and neutral
fraction are not too large, though the maximum particle energy is reduced.

One set of diagnostics for the physics of collisionless shocks is based on the
emission from particles in the narrow ionization zone just behind a
nonradiative shock (Raymond, 1991; Heng, 2010). In particular, H$\alpha$
photons from a nonradiative shock in partly neutral gas originate very close
to the shock, and Coulomb collisions do not have time to erase such signatures
as unequal electron and ion temperatures or non-Maxwellian velocity
distributions (Laming et al., 1996; Ghavamian et al., 2001; Raymond et al.,
2008, 2010). In the optical these shocks are seen as pure Balmer line
filaments whose profiles show a narrow component characteristic of the pre-
shock kinetic temperature and a broad component closely related to the post-
shock proton temperature (Chevalier & Raymond, 1978; Heng, 2010; van Adelsberg
et al., 2008). The intensity ratio of the broad and narrow components is
determined by the electron to ion temperature ratio at the shock (Ghavamian et
al., 2001; van Adelsberg et al., 2008; Helder et al., 2010).

The Balmer line profiles also contain signatures of shock precursors. In
general, the narrow component line widths are 40 to 50 $\rm km~{}s^{-1}$,
indicating temperatures around 40,000 K. If that were the ambient ISM
temperature, there would be no neutrals to create the Balmer line filament, so
the width is interpreted as an indication of heating in a narrow precursor too
thin to completely ionize the hydrogen (Smith et al., 1994; Hester et al.,
1994; Lee et al., 2007; Sollerman et al., 2003). Faint emission ahead of the
sharp filament is interpreted as emission from the compressed and heated
precursor gas (Hester et al., 1994; Lee et al., 2007, 2010).

This paper considers the role of neutrals in heating the precursor plasma and
computes the properties of precursor H$\alpha$ emission. While cosmic ray
pressure in the precursor can compress, heat and accelerate ions and electrons
by means of plasma turbulence and magnetic fields, the neutrals only interact
with the precursor by means of collisions with protons and electrons. If the
density is very high, neutrals and protons are tightly coupled by charge
transfer. In that case, the neutrals are compressed along with the protons and
adiabatically heated. They also share in any other heating of the protons,
such as dissipation of Alfvén waves generated by cosmic-ray streaming. On the
other hand, if the density is very low, neutrals pass through the precursor
and the shock without interacting at all, preserving their pre-shock velocity
distribution.

The intermediate case is more complex. A shock that efficiently accelerates
cosmic rays is strongly modified, and gas reaches a significant fraction of
the shock speed in the precursor (Vladimirov et al., 2008; Wagner et al.,
2009). If neutrals and ions are fairly well coupled, they can be described as
fluids whose relative speed gives a frictional heating similar to that in
C-shocks (Draine & McKee, 1993). If a neutral encounters this high speed
compressed plasma without having been brought gradually up to speed by many
previous charge transfers, it can be ionized and become a pickup ion (Raymond
et al., 2008; Ohira & Takahara, 2010) like those observed in the solar wind
(Moebius et al., 1985). It can then have an energy on the order of 1 keV,
which it can share with the other protons. Electron heating is more uncertain,
but it can occur by means of Lower Hybrid waves (Cairns & Zank, 2002). If the
electrons are heated they can excite and ionize H atoms, changing the
H$\alpha$ profile and the broad-to-narrow line ratio used as an electron
temperature diagnostic (Ghavamian et al., 2001).

In this paper we compute the proton, neutral and electron temperatures in the
precursors for a variety of parameters, along with the ionization and
excitation of H atoms. We consider the effects of these processes on Balmer
line diagnostics currently in use. Ohira & Takahara (2010) considered the
effects of neutrals on the velocity structure of the precursor, the
compression ratio and the acceleration process. They found that the pickup
ions can reduce the compression by the subshock and enhance proton injection
into the acceleration process. Morlino et al. (2010) self-consistently
computed the particle acceleration and heating due to neutrals, but within the
fluid approximation for both neutrals and ions. In this paper we emphasize the
effects on the H$\alpha$ line profile.

##  2 Model Calculations

We parameterize the precursor structure in a relatively simple manner. We
assume that the precursor accelerates and compresses the interstellar gas over
a length scale $\kappa$/$V_{S}$, where $\kappa$ is the diffusion coefficient
for cosmic rays near the cutoff. Effective cosmic-ray acceleration requires
$\kappa$ on the order of $10^{24}~{}\rm cm^{2}~{}s^{-1}$, and estimates based
on the scales of H$\alpha$ precursors are 2 to 4$\times 10^{24}~{}\rm
cm^{2}~{}s^{-1}$ (Lee et al., 2007, 2010). We do not consider the second order
effects of momentum and energy deposition by the neutrals on the precursor
length scale.

We assume an exponential form, so that the compression is given by

$\chi~{}=~{}1+(\chi_{1}-1)~{}~{}e^{(xV_{S}/\kappa)}$ (1)

where x is negative ahead of the shock and $\chi_{1}$ is the compression ratio
just upstream of the subshock. It is related to the fractional pressure of
cosmic rays behind the shock, w=$P_{CR}/(P_{G}+P_{CR})$, by equation 9 of Vink
et al. (2010). We simplify this equation with the assumption that for w$<$0.8
the compression in the gas subshock equals 4, so that

$\chi_{1}~{}=~{}(1-w/4)/(1-w)$ (2)

Mass conservation implies velocities $V=V_{S}/\chi$ in the frame of the shock.

To compute the proton and electron temperatures we include adiabatic
compression, Coulomb energy transfer between protons and electrons, energy
losses due to ionization and excitation of Hydrogen, and heating terms.

We assume that any neutral that interacts with the plasma at position $x_{i}$
joins the proton flow at that position. If the interaction was charge
transfer, a new neutral is formed with the bulk speed and thermal speed of the
protons at $x_{i}$. Thus the neutrals arriving at $x_{i}$ are those that last
went through charge transfer at all upstream positions $x_{j}$, and they have
the speeds, $v_{j}$, of the plasma at $x_{j}$. Each ionization of a neutral
from $x_{j}$ at $x_{i}$ deposits energy $0.5m_{p}(v_{j}-v_{i})^{2}$. We assume
that the energy is quickly thermalized among the protons, unlike Ohira &
Takahara (2010), who assumed a pickup ion velocity distribution. The
thermalization time scale is very uncertain, because full kinetic calculations
have not been carried out. However, in the highly turbulent precursor there
are many wave modes besides Alfvén waves that can thermalize the protons, in
particular those associated with bump-on-tail, mirror and firehose-like
instabilities (Winske et al., 1985; Gary, 1978; Sagdeev et al., 1986). We also
ignore heating due to Alfvén wave damping or shocks excited by the cosmic-ray
pressure gradient in order to isolate the effect of the neutrals. Therefore,
we compute a lower limit to the heating.

For electrons, we follow Cairns & Zank (2002), who found that ionization of
fast neutrals forms a ring beam, in which all the particles gyrate around the
magnetic field with the same speed but different phases. The ring beam is
unstable, and provided that the beam velocity (in this case the relative
velocity of bins i and j) is less than 5 times the Alfvén speed, it transfers
a significant amount of heat to electrons via Lower Hybrid waves. We follow
Cairns & Zank (2002) in taking this fraction to be 10%. Again, to isolate the
effects of neutrals we ignore any heating of electrons by Lower Hybrid waves
generated by cosmic-ray streaming (Ghavamian et al., 2007; Rakowski et al.,
2008).

Charge transfer rates are taken from (Schultz et al., 2008) using the
quadrature sum of the thermal speed and the ion-neutral relative speed.
Ionization and excitation rates are computed from cross sections from Janev &
Smith (1993) by integrating over the electron velocity distribution including
the relative electron-neutral flow speed.

##  3 Results

Figure 1 shows a set of models for a shock speed of 2000 $\rm km~{}s^{-1}$
with $\kappa=2.0\times 10^{24}~{}\rm cm^{2}~{}s^{-1}$, a pre-shock density of
0.2 $\rm cm^{-3}$ and a neutral fraction of 0.2. The four models have ratios
of cosmic-ray partial pressure to total pressure behind the subshock, w, of
0.1, 0.3, 0.5 and 0.7. The compression ratios just ahead of the subshock are
1.0833, 1.3214, 1.75 and 2.75. In the high $V_{S}$, high efficiency models the
neutrals are not compressed to this level, because of collisional ionization
and because some pass through without charge transfer. The protons and
electrons are strongly heated in the more efficient models, but the electrons
are much cooler than the protons. The drop in heating just before the subshock
in the 70% efficient model results from the reduced number of neutrals.

Figure 2 shows the velocity distributions of the neutrals perpendicular to the
shock just before the subshock. Note that the w=50% model shows a narrow
component due to neutrals that last experience charge transfer far upstream,
along with a broader component of particles that undergo charge transfer close
to the subshock.

Figure 3 shows a grid of models of the neutral velocity distribution at the
shock for a range of shock speeds and cosmic-ray partial pressures. The panels
show the FWHM measured directly from the computed velocity distribution and
the kurtosis, which would be 3.0 for a Gaussian distribution. Kurtosis is a
problematic statistical moment for real data because it is sensitive to noise
far from the line center and the choice of background level. However, for the
theoretical profiles computed here it highlights cases in which some neutrals
undergo charge transfer close to the subshock and others do not. We also show
the fraction of incident neutrals that survive up to the subshock and the
average number of excitations to the n=3 level per incident H atom. Not all of
these excitations will result in H$\alpha$ photons because some Ly$\beta$
photons escape, but this is a convenient comparison to the 0.2 to 0.25
H$\alpha$ photons per H atom produced in the post-shock region.

##  4 Discussion

The interaction of neutral hydrogen with the ionized plasma in cosmic-ray
precursors described above offers an important tool to measure the properties
of cosmic-ray precursors. The outcome of DSA is very much influenced by
physical processes in the precursor, which are not well determined. For
example, non-adiabatic heating and magnetic field amplification due the
presence of cosmic rays tend to decrease the overall compression factor from
$\chi_{12}>>20$ (e.g. Berezhko & Ellison, 1999; Blasi et al., 2005) to
$7\lesssim\chi_{12}\lesssim 15$ (Vladimirov et al., 2008; Caprioli et al.,
2008). In addition, if the Alfvén waves in the precursor have some drift
velocity this will affect the cosmic-ray pressure profile (Zirakashvili &
Ptuskin, 2008), which limits the escape of energy from the shock region. A
lower energy escape automatically implies a lower downstream cosmic-ray
pressure (Vink et al., 2010).

As shown here, neutrals will influence the physics of the precursor. Morlino
et al. (2010) treated the neutrals as a fluid coupled to the ions by charge
transfer for a unified treatment of the heating and dynamics of the precursor,
but the fluid approximation is only appropriate if neutrals and ions are
coupled fairly well. They obtained a FWHM of 46 $\rm km~{}s^{-1}$ for the
H$\alpha$ line in a 2000 $\rm km~{}s^{-1}$ shock with modest efficiency, a
pre-shock density of 1 $\rm cm^{-3}$ and 50% neutral fraction. For similar
parameters we find a non-Maxwellian profile with smaller FWHM and broader
wings.

Neutrals can also damp plasma waves, which limits the efficiency of cosmic-ray
acceleration. This damping is caused by the central processes described above:
charge exchange and ionization. Drury et al. (1996) found that the maximum
particle energy, and therefore the maximum acceleration efficiency, is
considerably higher than suggested by Draine & McKee (1993). In addition, the
heating due to neutrals penetrating the precursor is a form of non-adiabatic
heating. Energy dissipated in the precursor limits the amount of free energy
available for shock acceleration. If the neutrals ionized in the precursor
behave as pickup ions rather than thermalizing with the protons, the injection
efficiency and particle spectrum will be affected (Ohira & Takahara, 2010). In
any case, the heating of electrons in the precursor is poorly know, and that
will stongly affect the ionization and excitation of H atoms, which in turn
will affect the intensity ratio of the broad and narrow components as well as
the narrow component line width.

The physics of neutral-ion coupling means these processes are not only
sensitive to cosmic-ray pressure and the structure of the precursor, but also
to the pre-shock density and neutral fraction. For example, the protons and
neutrals in the Cygnus Loop nonradiative shocks (Salvesen et al., 2009) are
tightly coupled, and they behave nearly adiabatically, while the pre-shock
density in SN1006 is so low (Acero et al., 2007) that neutrals pass straight
through it. This may explain the narrow line width seen by Sollerman et al.
(2003) in SN 1006 in comparison to broader narrow lines observed for other
young SNRs. The number of charge transfer events for an average neutral in the
precursor can be estimated roughly as

$N_{\rm CT}=nL\sigma\chi_{1}/V_{S},$ (3)

where L is the precursor length scale, and the charge transfer cross section,
$\sigma$, declines slowly with velocity below about 2000 km s-1, then very
rapidly.

Comparing our results to observations, it is obvious that the observed narrow-
line H$\alpha$ widths are in general smaller than predicted by our
calculations for efficient shock acceleration (i.e., $w>0.5$). This may
indicate that none of the shocks investigated so far accelerate particles
efficiently. However, more work is needed before such a conclusion can be
drawn, as the line width depends also on pre-shock density and shock velocity.
For very high shock velocities combined with low densities the neutrals hardly
interact in the precursor, leading to narrow line widths. Such may be the case
for the northeastern region of RCW 86, for which Helder et al. (2009) reported
a high cosmic-ray acceleration efficiency ($w\gtrsim 0.5$, see also Vink et
al., 2010). For this region the pre-shock density may be as low as $n\lesssim
0.1$ (Vink et al., 2006), which, combined with the high velocity
($V_{s}\gtrsim 3000$ km/s), gives few interactions in the precursor and widths
$\lesssim 100$ km/s. Note that this is smaller than could be measured given
the moderate spectral resolution of the measurement.

Perhaps the most striking result from these calculations is that for efficient
shocks near 1000 $\rm km~{}s^{-1}$ a substantial number of H$\alpha$ photons
will be produced in the precursor. These will usually be included in the
narrow component, potentially affecting the electron temperature estimate
based on the broad-to-narrow intensity ratio (Ghavamian et al., 2001). Narrow
component emission from the precursor could explain the broad-to-narrow
intensity ratios that cannot be fit by models of post-shock emission (van
Adelsberg et al., 2008; Rakowski et al., 2009). In extreme cases, emission
from the precursor might also contribute to the broad component, possibly
accounting for the non-Maxwellian profile seen in Tycho’s SNR (Raymond et al.,
2010) and generally leading to an underestimate of the shock speed. Both of
these conclusions depend on the diffusion coefficient and the electron
heating, however.

Other H$\alpha$ line measurements show that the narrow lines are broader than
one might expect for temperatures of typical HII regions, but smaller than 50
km/s (Sollerman et al., 2003). (Not all of the narrow line emission comes from
the precursor, but the narrow line emission downstream is determined by the
velocity distribution in the precursor.) Another effect of charge exchange in
the precursor is that neutrals enter the downstream shock region with a
velocity offset with respect to the local interstellar medium, as seen in
Tycho’s supernova remnant (Lee et al., 2007). For shocks observed face on this
should produce a narrow line offset, which for the combined front and back
side of the remnant should lead to two narrow lines. The spectra of several
LMC remnants (Smith et al., 1994) do not show such an effect. For one of the
remnants in this set, SNR 0509-67.5, the cosmic-ray acceleration efficiency
was estimated to be $w\approx 0.2$ (Helder et al., 2010).

We note that several improvements should be made to the calculations presented
here. Additional heating due to wave dissipation can heat the protons,
resulting in larger narrow component line widths, or it can heat electrons,
increasing the H$\alpha$ narrow component intensity and reducing the number of
neutrals that reach the shock, especially if the electron velocity
distribution is non-Maxwellian (Laming & Lepri, 2007). In addition, ionization
and excitation by proton and helium ion impact are important at high relative
velocities (Laming et al., 1996), and at high shock speeds the velocity
distributions of particles are anisotropic (Heng & McCray, 2007; Heng et al.,
2007; van Adelsberg et al., 2008). Amplification of the magnetic field may
also be important, and radiative transfer calculations in the Ly$\beta$ line
must be done to compute the H$\alpha$ emission. We plan to address these
issues in future work.

This work was carried out while JCR was visiting the Astronomical Institute
Utrecht as Minnaert Professor. It was supported by NASA grant GO-11184.01-A-R
to the Smithsonian Astrophysical Observatory. JV and EH are supported by the
VIDI grant awarded to JV by the Netherlands Science Foundation (NWO).

## References

  * Acero et al. (2007) Acero, F., Ballet, J. & Decourchelle, A. 2007, A&A, 475, 883
  * Allen et al. (2008) Allen. G.E., Houck, J.C. & Sturner, S.J. 2008, ApJ, 683, 773
  * Berezhko & Ellison (1999) Berezhko, E.G. & Ellison, D.G. 1999, ApJ, 526, 385
  * Blasi et al. (2005) Blasi, P., Gabici, S., & Vannoni, G. 2005, MNRAS, 361, 907
  * Bykov (2005) Bykov, A. 2005, Ad. Sp. Res., 36, 738
  * Cairns & Zank (2002) Cairns, I.H. & Zank, G.P. 2002, Geophys. Res. Lett., 29, 47
  * Caprioli et al. (2008) Caprioli, D., Blasi, P., Amato, E. & Vietri, M. 2008, ApJ, 679, L139
  * Cassam-Chenaï et al. (2008) Cassam-Chenaï, G., Hughes, J.P., Reynoso, E.M., Badenes, C. & Moffett, D. 2008, ApJ, 680, 1180
  * Chevalier & Raymond (1978) Chevalier, R.A., & Raymond, J.C. 1978, ApJL, 225, L27
  * Draine & McKee (1993) Draine, B.T. & McKee, C.F. 1993, ARA$A, 31, 373
  * Drury et al. (1996) Drury, L.O., Duff, P., & Kirk, J.G. 1996, A&A, 309, 1002
  * Gary (1978) Gary, S.P. 1978, JGR, 85, 2304
  * Ghavamian et al. (2001) Ghavamian, P., Raymond, J.C., Smith, R.C., & Hartigan, P. 2001, ApJ, 547, 995
  * Ghavamian et al. (2007) Ghavamian, P., Laming, J.M. & Rakowski, C.E. 2007, ApJ, 645, L69
  * Helder et al. (2009) Helder, E.A., Vink, J., Bassa, C.G., Bamba, A., Bleeker, J.A.M., Funk, S., Ghavamian, P., van der Heyden, K.J., Verbunt, F. & Yamazaki, R. 2009, Science, 325, 719
  * Helder et al. (2010) Helder, E.A., Kosenko, D. & Vink, J. 2010, ApJ, 719, L14
  * Helder et al. (2010) Helder, E.A. 2010, PhD thesis, Utrecht University
  * Heng (2010) Heng, K. 2010, PASA, 27, 23
  * Heng & McCray (2007) Heng, K., & McCray, R. 2007, ApJ, 654, 923
  * Heng et al. (2007) Heng, K., van Adelsberg, M., McCray, R., & Raymond, J.C. 2007, ApJ, 668, 275
  * Hester et al. (1994) Hester, J.J., Raymond, J.C., & Blair, W.P. 1994, ApJ, 420, 721
  * Hughes et al. (2000) Hughes, J.P., Rakowski, C.E., Decourchelle, A. 2000, ApJ, 543, L61
  * Janev & Smith (1993) Janev, R.K. & Smith, J.J. 1993, Cross Sections for Collisional Processes of Hydrogen Atoms with Electrons, Protons and Multiply Charged Ions (Vienna: Int. At. Energy Agency)
  * Laming et al. (1996) Laming, J.M., Raymond, J.C., McLaughlin, B.M. & Blair, W.P. 1996, ApJ, 472, 267
  * Laming & Lepri (2007) Laming, J.M. & Lepri, S.T. 2007, ApJ, 660, 1642
  * Lee et al. (2007) Lee, J.-J., Koo, B.-C., Raymond, J.C., Ghavamian, P., Pyo, T.-S., Tajitsu, A., & Hayashi, M. 2007, ApJ, 659, L133
  * Lee et al. (2010) Lee, J.-J., Raymond, J.C., Park, S., Blair, W.P., Ghavamian, P., Winkler, P.-F. & Korreck, K. 2010, ApJ, 715, L146
  * Moebius et al. (1985) Moebius, E., Hovestadt, D., Klecker, B., Scholer, M, & Gloeckler, G. 1985, Nature, 318, 426
  * Morlino et al. (2010) Morlino, G., Amato, E., Blasi, P. & Caprioli, D. 2010, astro-ph 1012:2966
  * Ohira & Takahara (2010) Ohira, Y. & Takahara, F. 2010, ApJ, 721, L43
  * Rakowski et al. (2008) Rakowski, C.E., Laming, J.M., & Ghavamian, P. 2008, ApJ, 684, 408
  * Rakowski et al. (2009) Rakowski, C.E., Ghavamian, P. & Laming, J.M. 2009, ApJ, 696, 2195
  * Raymond (1991) Raymond, J.C. 1991, PASP, 103, 781
  * Raymond et al. (2008) Raymond, J.C., Isenberg, P.A., & Laming, J.M. 2008, ApJ, 682, 408
  * Raymond et al. (2010) Raymond, J.C., Winkler, P.F., Blair, W.P., Lee, J.-J. & Park, S. 2010, ApJ, 712, 901
  * Reynolds & Ellison (1992) Reynolds S.P. & Ellison, D.G. 1992, ApJ, 399, L75
  * Sagdeev et al. (1986) Sagdeev, R.Z., Shapiro, V.D., Shevchenko, V.I. & Szego, K. 1986, GRL, 13, 85
  * Salvesen et al. (2009) Salvesen, G.R., Raymond, J.C., & Edgar, R.J. 2009, ApJ, 702, 327
  * Schultz et al. (2008) Schultz, D.R., Krstic, P.S., Lee, T.G. & Raymond, J.C. 2008, ApJ, 678, 950
  * Smith et al. (1994) Smith, R.C., Kirshner, R.P., Blair, W.P., & Winkler, P.F. 1991, ApJ, 375, 652
  * Sollerman et al. (2003) Sollerman, J., Ghavamian, P., Lundqvist, P. & Smith, R.C. 2003, A&A, 407, 249
  * van Adelsberg et al. (2008) van Adelsberg, M., Heng, K., McCray, R., & Raymond, J.C. 2008, ApJ, 689, 1089
  * Vink et al. (2006) Vink, J., Bleeker, J., van der Heyden, K., Bykov, A., Bamba, A., & Yamazaki, R. 2006, ApJ, 648, 33
  * Vink et al. (2010) Vink, J., Yamazaki, R., Helder, E.A. & Schure, K.M. 2010, ApJ, 722, 1727
  * Vladimirov et al. (2008) Vladimirov, A.E., Bykov, A.M., & Ellison, D.C. 2008, ApJ, 688, 1084
  * Wagner et al. (2009) Wagner, A.Y., Lee, J.-J., Raymond, J.C., Hartquist, T.W., & Falle, S.A.E.G. 2009, ApJ, 690, 1412
  * Warren et al. (2005) Warren, J.S., et al. 2005, ApJ, 634, 376
  * Winske et al. (1985) Winske, D., Wu, C.S., Li, Y.Y., Mou, Z.Z. & Guo, S.Y. 1985, JGR, 90, 2713
  * Zirakashvili & Ptuskin (2008) Zirakashvili, V.N. & Ptuskin, V.S. 2008, in AIP Conference Series, Vol. 1085, (American Institute of Physics), p. 336

Figure 1: Plots of a) proton temperature, b) electron temperature, c) neutral
density and d) heating rate for for models having $\kappa=2.0\times
10^{24}~{}\rm cm^{2}~{}s^{-1}$, pre-shock density = 0.2 $\rm cm^{-3}$ and
neutral fraction 0.2. The models assume a post-shock cosmic-ray pressure of
10%, 30%, 50% and 70% of the total pressure, with the 10% curves at the bottom
and the 70% curves at the top in all four plots.  Figure 2: Velocity
distributions perpendicular to the flow direction at the gasdynamic subshock
for the four models shown in Figure 1. Note the narrow central component and
broader wings in the model with 50% cosmic-ray pressure and the large width
predicted by the 70% cosmic-ray pressure model.  Figure 3: Model grid for
$\kappa=2.0\times 10^{24}$, a pre-shock density of 0.2 $\rm cm^{-3}$ and
neutral fraction of 0.2. Models are shown for $P_{CR}/(P_{CR}+P_{G})$ = 0.1 to
0.8 at the shock front. The panels show the FWHM and kurtosis of neutrals that
reach the subshock, the fraction of neutrals that reach the subshock and the
number of excitations to n=3 in the precursor per incident neutral hydrogen
atom. As in Figure 1, the 80% models are the extreme cases.


#  Proton-induced activation cross sections in the energy range below 1 GeV

I. V. Moskalenko1, A. A. Andrianov2, A. V. Bytenko3, T. A. Frolova2, R. S.
Khalikov4, A. Yu. Konobeev5, Yu. A. Korovin2, T. V. Kulevoy4, I. S. Kuptsov2,
K. V. Pavlov4, A. Yu. Stankovskiy6, E. M. Syresin3, Yu. E. Titarenko4, V. D.
Vu2

(August 27, 2024)

hmargin=2.5cm, [] 1 HEPL & KIPAC, Stanford University, Stanford, CA 94305,
USA; 2 Obninsk Institute for Nuclear Power Engineering, NRNU MEPhI, Obninsk
249040, Russia; 3 Joint Institute for Nuclear Research, Dubna, Moscow Oblast,
141980 Russia; 4 NRC “Kurchatov Institute” – ITEP, Moscow, 117218 Russia; 5
Institute for Neutron Physics and Reactor Technology, KIT, 76021 Karlsruhe,
Germany; 6 SCK CEN, Boeretang 200, B-2400 Mol, Belgium

3

##  1 Introduction

The last decade brought spectacular advances in astrophysics of cosmic rays
(CRs) and $\gamma$-ray astronomy. Launches of missions that employ forefront
detector technologies were followed by a series of remarkable discoveries even
in the energy range that deemed as well-studied. Among those missions are the
Payload for Antimatter Matter Exploration and Light-nuclei Astrophysics
(PAMELA), the Fermi Large Area Telescope (Fermi-LAT), the Alpha Magnetic
Spectrometer-02 (AMS-02), NUCLEON experiment, CALorimetric Electron Telescope
(CALET), DArk Matter Particle Explorer mission (DAMPE), and Cosmic-Ray
Energetics and Mass investigation (ISS-CREAM). Outstanding results have been
also delivered by mature missions, such as the Cosmic Ray Isotope Spectrometer
onboard of the Advanced Composition Explorer (ACE-CRIS) operating for more
than two decades, and by Voyager 1, 2 spacecraft, currently at 151 AU/126 AU
from the Sun, respectively.

The modern technologies employed by many of these missions have enabled
measurements with unmatched precision, enabling searches for the subtle
signatures of dark matter (DM) and new physics in CR and $\gamma$-ray data.
The reached level of precision demonstrates that we are on the verge of major
discoveries.

Interpretation of CR data, observations of the diffuse photon emissions, and
search of new physics are critically dependent on the accuracy of the isotopic
production cross sections. The latter are the centerpiece of any propagation
code, but their accuracy has always been the Achilles’ heel of propagation
models. Handling this long-standing issue requires combined efforts of
professional nuclear physicists, both theorists and experimentalists, and
astrophysics community. Our goal is a reliable library of nuclear cross
section that can be used in many applications.

##  2 Current state

Interactions of nucleons and particles with nuclei in the energy range from
several dozens of MeV to several GeV/nucleon are usually described by a
combination of Intranuclear cascade (INC), pre-equilibrium, and equilibrium
models. More sophisticated Quantum Molecular Dynamics (QMD) models, due to
their complexity, did not spread much to really compete with INC models,
although they generally show better results at low energies.

The first stage of the reaction assumes that particle-nucleus interactions can
be represented by free particle-particle collisions inside the nucleus. Free
particle collisions form a sequence often called INC [1]. Nowadays, three
variants of INC models having statistical nature are recognized as most
developed: Bertini, CEM, INCL.

Bertini [2] considers nucleus as three-region approximation to the
continuously changing density distribution of the matter in the nucleus.
Nucleons are considered as Fermi gas and the cascade initiated by incident
particle is followed either until secondary particles escape the nucleus or
when their energy falls below a certain cut-off. Most sophisticated
developments of this model are realized in FLUKA [3] and GEANT4 [4] transport
codes. In FLUKA, in its nuclear interaction model called PEANUT, the
Generalized INC (GINC) model considers more radial zones, more quantum
effects, smooth transition to pre-equilibrium stage. Binary cascade model of
GEANT4 takes some advantages of QMD models to propagate particles in nuclear
potential.

The INC Liège (INCL) model [5] is based on classical mechanics, the only
quantum effect considered is dynamic Pauli blocking. All collisions are
considered as instantaneous at a point because the duration of events is short
compared to the time between successive collisions. The most recent version,
INCL4.6, has replaced binary cascade as default model choice in GEANT4.

The Cascade-Exciton Model (CEM) [6] describes the INC using the time-
independent Dubna cascade model [7] that makes use of 3D geometry for all
cascades.

All the INC models are followed by pre-equilibrium and equilibrium models,
describing intermediate and slow processes of nucleus de-excitation,
respectively. However, INCL4.6 does not consider pre-equilibrium stage and is
directly coupled to equilibrium models GEM [8] in PHITS transport code [9] or
ABLA [10] in GEANT4.

PEANUT of FLUKA and CEM consider their own sequences of models. PEANUT is
using modified hybrid exciton model of Blann [11] for pre-equilibrium part and
Weisskopf-Eving statistical model [12] for evaporation. For light nuclei Fermi
break-up model is implemented. In CEM, Modified Exciton Model [13] and GEM
equilibrium model are employed, Fermi break-up is also considered.

Figure 1: Example cross sections for reactions 56Fe($p,X$)54Mn (top),
56Fe($p,X$)53Fe (middle), and 28Si($p,X$)19F (bottom). The two top plots show
individual cross sections calculated using various nuclear codes and available
data. The bottom plot with 28Si target shows cumulative cross section with
contributions from the decay chains originating in secondary nuclides: 19C,
19,20N, 19O, 19Ne. In the absence of data, calculations show an order of
magnitude discrepancies at low energies, which translates into large
uncertainties in the interpretation of astrophysical data.

A restriction for all INC models at low energies is that De Broglie wavelength
of particles is becoming larger than characteristic nuclear size. Despite
their significant improvement at $\leq$100 MeV/nucleon, deterministic models
are generally recommended for this energy range.

Most complete description of nuclear reactions in the range of 10s MeV is
provided by TALYS [14], although it has constraints for atomic masses $<$10
and for energies $<$200 MeV. A comparison of evaporation mechanisms manifests
in favor of Hauser-Feshbach approach in the intermediate energy range $<$200
MeV/nucleon. This led to attempts of combining INC stage with Hauser-Feshbach
statistical emission. Such codes as CASCADEX (combining CASCADE/INPE realizing
Dubna INC model with TALYS 1.0) [18], cascadeXO (CASCADE/INPE + TALYS 1.95),
cemcas (CEM03.03 + TALYS 1.95) demonstrate better accuracy for certain
reactions. However, there is no unique model that can describe the cross
section data for all nuclei up to several GeV/nucleon. Therefore, it is
reasonable to represent the excitation functions as best estimate value
accompanied by an uncertainty band.

Table 1: Summary of INC models | Pre-equi- | De-exci- | Validity | Transport
---|---|---|---|---
INC | librium | tation | range | code
CEM03 | own | GEM |  $\leq$5 GeV | MCNP6
INCL4.6 | — | GEM |  $\leq$3 GeV | PHITS
Bertini/INCL++ | own | GEM/ABLA |  $\leq$10 GeV | GEANT4
PEANUT | own | own |  $\leq$5 GeV | FLUKA

##  3 New Facility

The major obstacle in improving an accuracy of the nuclear codes is the
scarcity of the available measurements (Fig. 1). We are aiming at using the
NICA facility (Nuclotron-based Ion Collider fAcility, https://nica.jinr.ru/) –
a new accelerator complex at the Joint Institute for Nuclear Research (JINR,
Dubna, Russia). NICA will provide variety of beam species from $p$ to Au ions.
Heavy ions will be accelerated up to kinetic energy of 4.5 GeV/nucleon,
protons – up to 12.6 GeV. Measurements of the isotopic cross sections will be
conducted using proton and ion beams and thin targets. Besides, NICA can be
used for studies of secondary neutron and $\gamma$-ray emission – of interest
for simulations of the radiation environment at the rocky surfaces of planets
and asteroids. New data will be used to tune models of nuclear reactions and
for astrophysical applications.

###  3.1 References

[1] D. Filges & F. Goldenbaum, Handbook of Spallation Research, Wiley-VCH
Verlag GmbH & Co (2009); [2] H. W. Bertini, Low-Energy Intranuclear Cascade
Calculation, Phys. Rev. 131 (1963) 1801; [3] G. Battistoni et al., Overview of
the FLUKA code, Ann. Nucl. En. 82 (2015) 10; [4] J. Allison et al., Recent
developments in GEANT4, NIMA 835 (2016) 186; [5] A. Boudard et al., New
potentialities of the Liège intranuclear cascade model for reactions induced
by nucleons and light charged particles, PRC 87 (2013) 014606; [6] S. G.
Mashnik, A. J. Sierk, CEM03.03 User Manual, LANL report LA-UR-12-01364, 2012;
[7] V. S. Barashenkov and V. D. Toneev, Interaction of High Energy Particle
and Nuclei with Atomic Nuclei, Atomizdat, Moscow, 1972; [8] S. Furihata,
Statistical Analysis of Light Fragment Production from Medium Energy Proton-
Induced Reactions, NIMB 171 (2000) 252; [9] T. Sato, et al., Features of
Particle and Heavy Ion Transport Code System PHITS Version 3.02, J. Nucl. Sci.
& Technology, 55 (6) (2018) 684; [10] A. Kelic, et al., Proc. Joint ICTP-IAEA
Adv. Workshop on Model Codes for Spallation Reactions, ICTP Trieste, Italy,
4-8 February 2008, IAEA INDC(NDS)-530 (Vienna 2008, p181); [11] M. Blann,
Hybrid model for pre-equilibrium decay in nuclear reactions. PRL 27 (1971)
337; [12] V. E. Weisskopf, D. H. Ewing, On the Yield of Nuclear Reactions with
Heavy Elements, Phys. Rev. 57 (1940) 472; [13] K. K. Gudima, G. A. Ososkov and
V. D. Toneev, Model for the Pre-Equilibrium Decay of Excited Nuclei, Sov. J.
Nucl. Phys. 21 (1975) 138; [14] A. J. Koning, D. Rochman, Modern nuclear data
evaluation with the TALYS code system, Nucl. Data Sheets, 113 (2012) 2841;
[15] J. Raynal, “Notes on ECIS94,” tech. rep., CEA Saclay Report CEA-N-2772,
1994; [16] A. Koning and M. Duijvestijn, A global pre-equilibrium model from 7
to 200 MeV based on the optical model potential, NPA 744 (2004) 15; [17] W.
Hauser, H. Feshbach, The Inelastic Scattering of Neutrons, Phys. Rev. 87
(1952) 366; [18] A. Stankovskiy and A. Konobeyev, CASCADEX – a combination of
intranuclear cascade model from CASCADE/INPE with the Hauser-Feshbach
evaporation/fission calculations from TALYS, NIMA 594 (2008) 420.

Acknowledgements: IVM acknowledges support from NASA grant No. NNX17AB48G.