\section{Introduction}

A flow-based model~\cite{DinhKB14,dinh2016density,KingmaD18}, also known as normalizing flow, represents a data distribution by a non-linear invertible transformation of a simple distribution, e.g., Gaussian distribution. Thanks to the special design of the computation layers~of~the transformation, the flow-based model defines a
normalized probability distribution explicitly, making it convenient for learning, sampling and inference. Specifically, its learning can be handily accomplished by the maximum likelihood estimation with the tractable density function, its sampling can be efficiently achieved by a direct mapping of samples drawn from a Gaussian distribution, and its inference can be easily realized by the inverse transformation of data. Therefore, flow-based~models~have been widely used for direct modeling of various types of data, such as images~\cite{KingmaD18}, videos~\cite{KumarBEFLDK20}, point clouds~\cite{yang2019pointflow}, audio~\cite{PingPZ020} etc.

In spite of the computational convenience, the expressivity of a normalizing flow is constrained by the unnatural design of the invertible transformation. Increasing the model expressivity typically requires significantly enlarging the network structure, thus leading to difficulty in constructing the network architecture and optimizing the network parameters. In order to foster the strengths and circumvent the shortcomings of the normalizing flow, we propose to learn a normalizing flow in the low-dimensional latent space of observed data. The resulting generative model becomes a normalizing flow standing on top of a top-down generator network. The latent space normalizing flow model can not only reserve the computational merit, but also focus on capturing regularities and rules of the data in the continuous latent space without having to care about the data type and data complexity.
Also, the learning and the design of a normalizing flow in the latent space is much easier than those in the data space due to the low dimensionality of the latent space. Typically, it can be parameterized by a small invertible multi-layer perceptron that consists of several layers.

Besides, the proposed latent space normalizing flow can be considered an informative prior distribution of the latent variables in the top-down generator model. Compared to the commonly used isotropic Gaussian (or uniform) non-informative prior in the conventional generator model, the informative prior represented by a trainable normalizing flow can significantly increase the expressive power of the entire generator.
The normalizing flow model is a good fit for the role of prior distribution in the deep latent variable framework in the sense that (i) the normalizing-flow-based prior can be jointly trained with the top-down network due to the tractability of the prior probability density and (ii) the data generation from the learned model can be explicit and fast because of the ease of the prior sampling. We call our model the latent space flow-based prior model (LFBM).

In this work, we propose a novel and principled way to jointly train the flow-based prior model and the generator via maximum likelihood estimation (MLE). The training process alternates the inference step and the learning step. Each inference step involves Markov chain Monte Carlo (MCMC) sampling of the posterior distribution, which is guided by both the flow-based prior distribution and the top-down generator, for the sake of inferring latent variables of each observed example, while, in each learning step, the parameters of the flow-based prior is updated by maximizing the likelihood of the inferred latent variables and the parameters of the generator are updated by minimizing the regression error between the inferred latent variables and the corresponding observed examples.
As to the MCMC inference, we adopt a short-run Langevin dynamics, i.e., a finite-step Langevin dynamics initialized with a fixed Gaussian distribution. Even though the short-run MCMC is an approximate inference due to the non-convergence, the normalizing flow prior can be adjusted based on the approximately inferred latent variables to mend the discrepancy between the aggregate posterior distribution and the prior distribution. We show that the learning algorithm based on short-run MCMC inference is a perturbation of MLE. In contrast to the latent space energy-based prior model (LEBM)~\cite{NijkampP0ZZW20}, in which short-run MCMCs are applied to not only the posterior sampling but also the EBM prior sampling, our model is able to perform an exact prior sampling via ancestral sampling and thus less biased from MLE than the LEBM.
Compared to the baseline that trains the same LFBM along with an extra bottom-up inference model (with additional parameters) for variational auto-encoders (VAEs) ~\cite{KingmaW13}, our framework with a Langevin flow for approximate inference exhibits superiorities, e.g., better empirical performance, less model parameters, and free of inference network design.
Table \ref{tab:propertyanalysis} presents an analytical comparison of different latent variable models using various combinations of prior and posterior components. This comparison can clearly show the advantage of our proposed model and the motivation of this work. Specifically, this paper is motivated by proposing a practical, design-friendly and lightweighted latent variable framework with informative prior.

We summarize the main contributions of our paper as follows: (1) We propose a simple but novel deep generative model where a latent space normalizing flow, which serves as a prior, stands on a single top-down generator network. (2) We propose a natural and principled maximum likelihood learning algorithm to jointly train the latent space normalizing flow and the top-down network with MCMC-based inference over latent variables. (3) We also propose to use short-run Langevin flow as an approximate inference for efficient training. (4) We provide theoretical understanding of the proposed learning framework.  (5) We provide extensive and strong experimental results on different aspects, including image synthesis, inference, reconstruction, inpainting and recovery, to validate the effectiveness of the proposed models and learning algorithms.

\begin{table*}[ht]
\centering
\begin{tabular}{l|ccccccc}
\multirow{2}{*}{Methods (prior + posterior)}   & easy   & informative & fast      & fast      & less       & fast   & practical  \\
                          & design & prior       & sample & inference & parameters & training \\ \hline
Normalizing Flow + short-run MCMC       &    \circletfill    &      \circletfill       &     \circletfill      &     \circletfillhl      &     \circletfill       &   \circletfillhl     & \circletfill \\
Normalizing Flow + long-run MCMC        &   \circletfill    & \circletfill            &    \circletfill       &      \circlet     &     \circletfill        &  \circlet      &  \circlet  \\
Normalizing Flow + inference net        &   \circlet   &   \circletfill          &     \circletfill      &       \circletfill     &       \circlet      &    \circletfill   &  \circletfill \\
EBM + short-run MCMC      &   \circletfill     &  \circletfill           &   \circlet        &         \circletfillhl   &        \circletfill     &  \circletfillhl    &  \circletfill   \\
EBM + long-run MCMC       &     \circletfill    &  \circletfill           &    \circlet       &   \circlet        &      \circletfill       &    \circlet  &   \circlet  \\
EBM + inference net       &    \circlet     &    \circletfill         &    \circlet       &    \circletfill       &     \circlet        &  \circletfill   &  \circletfill   \\
Gaussian + short-run MCMC &    \circletfill     &  \circlet           &     \circletfill      &  \circletfillhl         &     \circletfill        &      \circletfillhl  &  \circletfill \\
Gaussian + long-run MCMC  &   \circletfill      &  \circlet           &    \circletfill       &    \circlet       &     \circletfill        &   \circlet   &    \circlet \\
Gaussian + inference net  &    \circlet     &   \circlet          &     \circletfill      &  \circletfill         &     \circlet        &   \circletfill   &   \circletfill  \\
\end{tabular}
\caption{A comparison of latent variable models with different priors and inference process. ($\circletfill$ is good, $\circletfillhl$ is OK, and $\circlet$ is bad.) }
\label{tab:propertyanalysis}
\end{table*}

\section{Related Work}

\paragraph{Data Space Flow-based Models and Latent Space Flow-based Models.} This work builds upon the techniques of normalizing flows, which have been widely applied to data space for representing high-dimensional data distributions directly~\cite{DinhKB14,dinh2016density,KingmaD18,yang2019pointflow,PingPZ020,KumarBEFLDK20,ShiXZZZT20}. This work moves the flow-based models from data space to latent space for representing low-dimensional latent manifolds from which the observed data are generated. Some prior works, such as~\citet{0022KSDDSSA17,huang2017learnable,MaZH19,xiao2019generative,DingG21}, study flow-based priors in VAE scenarios. Our paper mainly focus on maximum likelihood learning of flow-based priors with short-run Langevin MCMC~inference.

\paragraph{Top-down Models with Implicit Density and Explicit Density.} A top-down model is a nonlinear generalization of the factor analysis model~\cite{rubin1982algorithms} and its density is generally implicit and intractable. Likelihood-based training of a top-down generator model typically relies on either the Markov chain Monte Carlo sampling that evaluates the posterior distribution, such as~\citet{HanLZW17,XieGZZW19} or the variational inference that approximates the posterior distribution, such as~\citet{KingmaW13}. A normalizing flow is a special top-down model that defines an explicit and normalized probability density and it can be learned by maximizing the exact likelihood of the observed data through a gradient ascent algorithm. Our model also belongs to a top-down model that stacks a top-down model with explicit density on top of the other one with implicit density. Similar to other generic top-down generator models, sampling from our model can also be achieved via direct ancestral sampling. However, the flow-based model at the top serves as a flexible and tractable prior distribution for the top-down model at the~bottom.

\paragraph{Latent Space Energy-based Prior.} Recently,~\citet{Pang0NZW20,PangZ0W21,PangW21,nipsZhang21} propose to use deep energy-based models (EBM) \cite{XieLZW16} to serve as flexible priors for top-down generation models. Because an EBM prior defines an unnormalized density, drawing samples from such a prior requires iterative MCMC sampling, thus bringing in inconveniences during the training and sampling of the whole model. Given the inferred latent variables, the update of the latent space EBM prior involves a non-convergent short-run MCMC that samples from the prior, which leads to a biased estimation of the EBM prior distribution \cite{NijkampHZW19,coopflow}. In contrast, our normalizing-flow-based prior model doesn't have such an issue due to its properties of tractable density and fast generation.

\paragraph{MCMC-based Inference.}~Our~LFBM~is~related~to MLE of top-down generator model with MCMC-based inference, which is also called the alternating back-propagation (ABP) learning in~\citet{HanLZW17,XieGZZW19,XieGZZW20,ZhuXLE19,ZhangXB20,AnXL21,XingGHZW22}. The original ABP algorithm is designed to train a generator network with a known and simple Gaussian prior. The LEBM~\cite{Pang0NZW20} generalizes the ABP algorithm to a framework with a trainable and unnormalized energy-based prior model, while our paper generalizes the ABP algorithm to a framework with a trainable and normalized flow-based prior model.

\paragraph{Variational Inference.}
Most regular variational auto-encoder (VAE) frameworks employ a fixed Gaussian distribution as the prior and use reparameterization trick~\cite{KingmaW13} to design the amortized inference network that outputs a Gaussian posterior. Some prior works, such as~\citet{RezendeM15,kingma2016improved}, improve the amortized inference network by using a normalizing flow that outputs a flexible posterior. Our LFBM parameterizes the prior by a normalizing flow, treats the built-in short-run Langevin MCMC as the inference model, and learns them simultaneously from the observed data by following the empirical Bayes principle.

\section{Proposed Framework}
\subsection{Normalizing Flow in Latent Space}
Let $x\in \mathbb{R}^D$ be an observed signal such as an image, and let $z \in \mathbb{R}^d$ be the latent variables of $x$. The joint distribution of the signal and the latent variables $(x,z)$ is given by $p_{\theta}(x,z)=p_{\alpha}(z)p_{\beta}(x|z)$, where the distribution of $z$, i.e., $p_{\alpha}(z)$, is the prior model parameterized by $\alpha$, and the conditional distribution of $x$ given $z$, i.e., $p_{\beta}(x|z)$, is the top-down generation model parameterized by $\beta$. For notational convenience, let $\theta=(\alpha,\beta)$.

The top-down generaton model is a non-linear transformation of the latent variables $z$ to generate the signal $x$, in which the transformation is parameterized by a neural network $g_{\beta}: \mathbb{R}^d \rightarrow \mathbb{R}^D$, which is similar to the decoder of a variational auto-encoder (VAE)~\cite{KingmaW13}. To be specific, $x=g_{\beta}(z)+\epsilon$, where $\epsilon \sim \mathcal{N}(0,\sigma^2I_D)$ is an observation residual. Thus, $p_{\beta}(x|z)=\mathcal{N}(g_{\beta}(z), \sigma^2I_D)$, where the standard deviation $\sigma$ is a hyper-parameter and assumed to be given.

For most existing top-down generative models, the prior models are assumed to be non-informative isotropic Gaussian distributions. In this paper, we propose to formulate the prior $p_{\alpha}(z)$ as a flow-based model~\cite{dinh2016density, DinhKB14, KingmaD18}, which is of the form
\begin{eqnarray}
z_0 \sim q_0(z_0), \hspace{1mm} z =f_{\alpha}(z_0),
\label{eq:prior_gen}
\end{eqnarray}
where $q_{0}(z_0)$ is a base distribution and has a simple and tractable density, such as a Gaussian white noise distribution: $q_0(z_0)=\mathcal{N}(0,I_d)$. $f_{\alpha}: \mathbb{R}^d \rightarrow \mathbb{R}^d  $ is an invertible or bijective function, which is a composition of a sequence of invertible transformations, i.e., $f_{\alpha}(z_0) = f_{\alpha_L} \circ  \cdots \circ f_{\alpha_2} \circ f_{\alpha_1}(z_0)$, whose inverse and logarithm of the determinants of the Jacobians can be explicitly obtained in closed form. Examples of such architectures include NICE~\cite{dinh2016density}, RealNVP~\cite{DinhKB14}, and Glow~\cite{KingmaD18}.

For $l=1,...,L$, let $z_{l}=f_{\alpha_l}(z_{l-1})$ be a sequence of random variables transformed from $z_{0}$, and define $z:=z_L$. According to the change-of-variable law of probabilities, $q_{0}(z_0)dz_0 = p_{\alpha}(z)dz$, the density of the flow-based prior model can be written as
\begin{eqnarray}
p_{\alpha}(z)&=q_{0}(z_0) \frac{dz_0}{dz}= q_{0}(f_{\alpha}^{-1}(z))\left|\det \left(\frac{\partial f_{\alpha}^{-1}(z)}{\partial z}\right)\right| \nonumber\\
&= q_{0}(f_{\alpha}^{-1}(z)) \prod_{l=1}^L\left|\det \left(\frac{\partial z_{l-1}}{\partial z_l}\right)\right|,
\label{eq:prior_den}
\end{eqnarray}
where $f_{\alpha}^{-1}(z) =f_{\alpha_1}^{-1}  \circ \cdots \circ f_{\alpha_{L-1}}^{-1} \circ f_{\alpha_L}^{-1}(z) $, and the determinant of the Jacobian matrix $(\partial z_{l-1}/\partial z_l)$ can be easy to compute with well-designed transformation functions in the flow-based models. Typically, choosing transformations that can lead to a triangle Jacobian matrix will simplify the computation of the determinant: $\left|\det  ( \partial z_{l-1} /\partial z_l  )\right|=\prod |\text{diag}(\partial z_{l-1} /\partial z_l)|$, where $\text{diag}()$ takes the diagonal elements of the Jacobian matrix. Thus, the flow-based prior model has two nice properties: (i) analytically tractable normalized density (i.e., Eq. (\ref{eq:prior_den})),  and (ii) easy to draw samples from by using ancestral sampling (i.e., Eq. (\ref{eq:prior_gen})).

The flow-based prior model in Eq.~(\ref{eq:prior_gen}) can be regarded as a flow-based refinement or calibration of the original Gaussian prior distribution $q_{0}$, which is a widely used prior distribution in most top-down generative frameworks. There are two advantages of using flow-based prior models for generative learning: (i) This enables us to learn a flexible prior from data to capture more meaningful latent components compared to those using a simple Gaussian prior model; (ii) The original generator learns $g_{\beta}$ to map from a fixed unimodal prior $q_0$ to the highly multi-modal data distribution, while in the proposed framework, the flow-based prior model corrects $q_0$ such that $g_{\beta}$ can be easier to transform the calibrated distribution $p_{\alpha}$ to the data distribution.

\subsection{Maximum Likelihood Learning with Langevin Inference}
The top-down generative model with a flow-based prior can be trained via maximum likelihood estimation. For the training examples $\{x_i,i=1,...,N\}$, the observed-data log-likelihood function is given by
\begin{equation}
    \mathcal{L}(\theta)=\frac{1}{N}\sum_{i=1}^{N} \log p_{\theta}(x_i), \label{eq:mle}
\end{equation}
where the marginal distribution is obtained by integrating out the latent variables $z$: $p_{\theta}(x)=\int p_{\theta}(x,z)dz =\int p_{\alpha}(z)p_{\beta}(x|z)dz$. Maximizing the log-likelihood function $\mathcal{L}(\theta)$ is equivalent to minimizing the Kullback-Leibler (KL) divergence between the model $p_{\theta}(x)$ and the data distribution $p_{\text{data}}(x)$. The gradient of $\mathcal{L}(\theta)$ is computed according~to
\begin{eqnarray}
\begin{aligned}
    \nabla_{\theta}& \log p_{\theta}(x) = \mathbb{E}_{p_{\theta}(z|x)} \left[\nabla_{\theta} \log p_{\theta}(x,z)\right]  \\     &=\mathbb{E}_{p_{\theta}(z|x)}[\nabla_{\theta} (\log p_{\alpha}(z)+  \log p_{\beta}(x|z))], \label{eq:gradient}
    \end{aligned}
\end{eqnarray}
where the posterior distribution of $z$ is given by $p_{\theta}(z|x)=p_{\theta}(x,z)/p_{\theta}(x) \propto p_{\alpha}(z)p_{\beta}(x|z)$. The inference distribution $p_{\theta}(z|x)$ is dependent on both the prior model $\alpha$ and the generation model $\beta$.

For the flow-based prior model,
\begin{eqnarray}
\begin{aligned}
    \log p_{\alpha}(z) &= \log q_{0}(f_{\alpha}^{-1}(z)) + \sum_{l=1}^{L} \log |\text{det} (\frac{\partial z_{l-1}}{\partial z_l})|  \\
    &= \log q_{0}(z_0) + \sum_{l=1}^{L} \text{sum}(\log |\text{diag} (\frac{\partial z_{l-1}}{\partial z_l})|)
    \end{aligned}
\end{eqnarray}
where $\log()$ takes element-wise logarithm, and $\text{sum}()$ takes the sum over all elements in a vector. Given a datapoint $z$, computing its log-likelihood only need one  pass of the inverse function $f_{\alpha}^{-1}$. We define $l_{\alpha}(z)=\log p_{\alpha}(z)$ to explicitly indicate that computing the log-likelihood of $z$ under the flow-based model is computationally tractable and can be regarded as a function of $z$. The learning gradient of $\alpha$ for a datapoint $x$ is
\begin{eqnarray}
\begin{aligned}
    \nabla_{\alpha} \log p_{\theta}(x)&= \mathbb{E}_{p_{\theta}(z|x)}[\nabla_{\alpha} \log p_{\alpha}(z)] \\
    &=\mathbb{E}_{p_{\theta}(z|x)}[\nabla_{\alpha} l_{\alpha}(z)]. \label{eq:gradient_prior}
    \end{aligned}
\end{eqnarray}
The update of the latent space flow-based prior model depends on the observed example $x$, but different from the original data space flow-based model, it treats the latent variables inferred from $x$ as a training example and seeks to maximize the log-likelihood of the inferred latent variables. The updated $p_{\alpha}(z)$ will further influence the inference accuracy of $p_{\theta}(z|x)$.

As to the generation model, the learning gradient of $\beta$ for a datapoint $x$ is
\begin{eqnarray}
    \nabla_{\beta} \log p_{\theta}(x)= \mathbb{E}_{p_{\theta}(z|x)}[\nabla_{\beta} \log p_{\beta}(x|z)]. \label{eq:gradient_generation}
\end{eqnarray}
Since $p_{\beta}(x|z)$ is in the form of a Gaussian distribution with a mean of $g_{\beta}(z)$ and a standard deviation of $\sigma$, $\nabla_{\beta} \log p_{\beta}(x|z)= \nabla_{\beta} (-\frac{1}{2\sigma^2}||x-g_{\beta}(z)||^2+\text{const})= \frac{1}{\sigma^2}(x-g_{\beta}(z))\nabla_{\beta}g_{\beta}(z)$.

The learning gradient of $\theta$ in Eq.(\ref{eq:gradient}) is decomposed into the learning gradient of $\alpha$ in Eq.(\ref{eq:gradient_prior}) and the learning gradient of $\beta$ in Eq.(\ref{eq:gradient_generation}), both of which involve the expectation with respect to the intractable posterior $p_{\theta}(z|x)$, which can be approximated by drawing samples from $p_{\theta}(z|x)$ and then computing the Monte Carlo average. Sampling from $p_{\theta}(z|x)$ can be achieved by Langevin dynamics that iterates
\begin{eqnarray}
    z_{(k+1)} = z_{(k)} + \xi  \nabla_{z} \log p_{\theta}(z_{(k)}|x)+\sqrt{2 \xi} \epsilon_{(k)}; \nonumber\\
     \hspace{2mm} z_{(0)}\sim q_{0}(z), \epsilon_{(k)} \sim \mathcal{N}(0,I_d),
    \label{eq:MCMC_post}
\end{eqnarray}
where $k$ indexes the Langevin time step, $\nabla_{z} \log p_{\theta}(z|x)=\nabla_{z} \log p_{\alpha}(z)+ \nabla_{z} \log p_{\beta}(x|z) = \nabla_{z} l_{\alpha}(z)+ \frac{1}{\sigma^2}(x-g_{\beta}(z))\nabla_{z}g_{\beta}(z)$, and $\xi$ is the Langevin step size. Algorithm~\ref{alg1} presents the learning algorithm of~LFBM. We use the Adam optimizer \cite{KingmaB14} to update the parameters.

\subsection{Perturbation of Maximum Likelihood}

The Langevin inference in Eq.(\ref{eq:MCMC_post}) is practically  non-mixing and non-convergent because MCMC chains from different stating points can get trapped in local modes of the posterior distribution or MCMC chains are not sufficiently long to converge. Using a short-run non-convergent MCMC, e.g., a finite-step Langevin dynamics, as a flow-like approximate inference is more computationally efficient than using a long-run MCMC and more implementationally convenient than using an amortized inference network. We use $\hat{p}_{\theta}(z|x)$ to denote the short-run Langevin flow distribution, which is obtained by running a $K$-step Langevin dynamics starting from a fixed initial distribution $q_0(z)$ toward the true posterior ${p}_{\theta}(z|x)$.
Strictly speaking, Algorithm~\ref{alg1} with a short-run MCMC inference is a perturbation of maximum likelihood. That is, given $\theta_t$ at iteration $t$, the learning gradient of $\theta$ is
\begin{eqnarray}
\mathcal{\hat{L}}(\theta) = \mathcal{L}(\theta)- \frac{1}{N}\sum_{i=1}^{N} \mathbb{D}_{\text{KL}}(\hat{p}_{\theta_t}(z_i|x_i)||p_{\theta}(z_i|x_i)),
\label{eq:short_run_loss0}
\end{eqnarray}
which is a lower bound of  $\mathcal{L}(\theta)$. Theoretically, if $K~\rightarrow~\infty$, $\xi \rightarrow 0$, then $\mathbb{D}_{\text{KL}}(\hat{p}_{\theta}(z|x)||p_{\theta}(z|x)) \rightarrow 0$. The estimate of $\theta$ is an MLE solution. However, $\mathbb{D}_{\text{KL}}(\hat{p}_{\theta}(z|x)||p_{\theta}(z|x))\neq0$ in practise, the estimate of $\theta$ is biased from~MLE.

In Eq.(\ref{eq:short_run_loss0}), the first term $\mathcal{L}(\theta)$ corresponds to MLE, i.e., we want $p_{\theta}(x)$ to be close to $p_{\text{data}}$, while the second KL-divergence term means that we want to bias $\theta$ from MLE so that $p_{\theta}(z|x)$ is close to $\hat{p}_{\theta_t}(z|x)$. In fact, this is not a bad thing because it means we intend to bias the model so that short-run MCMC is close to the true posterior.

The algorithm converges when learning gradients equal to zeros. The resulting estimators of $\theta=(\alpha,\beta)$ solve the following estimating equations:
\begin{align}
& \frac{1}{N}\sum_{i=1}^{N} \mathbb{E}_{\hat{p}_{\theta}(z_i|x_i)}[\nabla_{ \alpha} \log p_{\alpha}(z_i)]=0, \\
& \frac{1}{N}\sum_{i=1}^{N} \mathbb{E}_{\hat{p}_{\theta}(z_i|x_i)}[\nabla_{ \beta} \log p_{\beta}(x_i|z_i)]=0.
\end{align}

In contrast to LFBM, the latent space energy-based prior model (LEBM)~\cite{Pang0NZW20} adopts an unnormalized energy-based density $p^{\text{ebm}}_{\alpha}(z)\propto \exp(f_{\alpha}(z))q_0(z)$ as the prior. The MLE of LEBM requires not only MCMC inference from $p_{\beta}(z|x)$ but also MCMC sampling from $p^{\text{ebm}}_{\alpha}(z)$. The short-run MCMC inference and sampling cause the learning algorithm of the LEBM to be a perturbation of MLE. Given $\theta_t$, the update of $\theta$ is based on the gradient of
\begin{eqnarray}
\mathcal{\hat{L}}_{\text{LEBM}}(\theta) &= \mathcal{L}(\theta)- \frac{1}{N}\sum_{i=1}^{N} \mathbb{D}_{\text{KL}}(\hat{p}_{\theta_t}(z_i|x_i)||p_{\theta}(z_i|x_i)) \nonumber\\
&+ \mathbb{D}_{\text{KL}}(\hat{p}^{\text{ebm}}_{\alpha_t}(z)||p^{\text{ebm}}_{\alpha}(z)), \label{eq:short_run_ebm_loss0}
\end{eqnarray}
where $\hat{p}^{\text{ebm}}_{\alpha}(z)$ represents the short-run MCMC distribution. The second KL-divergence term in Eq.(\ref{eq:short_run_ebm_loss0}) is due to imperfect sampling from the EBM prior. Our flow-based prior $p_{\alpha}(z)$ is capable of exact sampling, therefore the second KL-divergence term disappears in Eq.(\ref{eq:short_run_loss0}). Thus, $\mathcal{\hat{L}}_{\text{LEBM}}(\theta)$ is more biased from MLE than $\mathcal{\hat{L}}(\theta)$.

\begin{algorithm}
\caption{Maximum likelihood learning of latent space normalizing flow model}\label{alg1}
\textbf{Input}: (1) Observed signals for training $\{x_i\}_{i}^{N}$; (2) Maximal number of learning iterations $T$; (3) Numbers of Langevin steps for posterior $K$; (4) Langevin step size $\xi$ for the posterior; (5) Learning rates for flow-based prior model and the generation model $\{\eta_\alpha,\eta_\beta\}$.\\
\textbf{Output}:
Parameters $\beta$ for the generation model and $\alpha$ for the flow-based prior model
\begin{algorithmic}[1]
\State Randomly initialize $\alpha$ and $\beta$
\For{$t \leftarrow  1$ to $T$}
\State Sample a batch of observed examples $\{x_i\}_i^n$
\State For each $x_i$, sample the posterior $z_i \sim p_{\theta}(z|x_i)$ using $K$ Langevin steps in  Eq.(\ref{eq:MCMC_post}) with a step size $\xi$.
\State Update flow-based prior by Adam optimizer  with the gradient $\nabla{\alpha}$ in Eq.(\ref{eq:gradient_prior}) and a learning rate $\gamma_\alpha$.
\State Update generation model by Adam optimizer with the gradient $\nabla{\beta}$ in Eq.(\ref{eq:gradient_generation}) and a learning rate $\gamma_\beta$.
\EndFor
\end{algorithmic}
\end{algorithm}

\subsection{Matching Normalizing Flow Prior and Aggregated Langevin Flow Posterior}

We further reveal the interaction between the normalizing flow prior and the Langevin flow
posterior in latent space during maximum likelihood estimation. The maximum likelihood learning minimizes the KL-divergence between the aggregated posterior and the prior.
Let $\tilde{p}(x,z)=p_{\text{data}}(x)p_{\theta}(z|x)$, and then the aggregated posterior is $\tilde{p}(z)=\int \tilde{p}(x,z) dx =\mathbb{E}_{p_{\text{data}}(x)}[p_{\theta}(z|x)]$. Also, we have $\tilde{p}(x,z)=\tilde{p}(z)\tilde{p}(x|z)$.
We show that
\begin{eqnarray}
\begin{aligned}
&\mathbb{D}_{\text{KL}}(p_{\text{data}}(x)||p_{\theta}(x))\\
=&\mathbb{D}_{\text{KL}}(p_{\text{data}}(x)p_{\theta}(z|x)||p_{\theta}(x)p_{\theta}(z|x))\\
=&\mathbb{D}_{\text{KL}}(\tilde{p}(z)\tilde{p}(x|z)||p_{\alpha}(z)p_{\beta}(x|z))\\
=&\mathbb{D}_{\text{KL}}(\tilde{p}(z)||p_{\alpha}(z)) + \mathbb{D}_{\text{KL}}(\tilde{p}(x|z)||p_{\beta}(x|z)),
\end{aligned}
\end{eqnarray}
which means that the maximum likelihood learning of the whole latent variable models involves a behavior that minimizes the KL divergence between the aggregated posterior and the normalizing flow in the latent space.

The above analysis is based on the ideal scenario with a long-run convergent Langevin inference. In our framework, we assume and allow short-run Langevin sampling, and then our objective turns into maximizing the perturbation of maximum likelihood, i.e., $\log p_{\theta}(x) - \mathbb{D}_{\text{KL}}(\hat{p}_{\theta_t} (z|x)  ||  p_{\theta}(z|x) )$, or equivalently minimizing $\mathbb{D}_{\text{KL}}(p_{\text{data}}(x)||p_{\theta}(x)) + \mathbb{D}_{\text{KL}}(\hat{p}_{\theta_t} (z|x)  ||  p_{\theta}(z|x) )$.
We now generalize the above analysis to the short-run non-convergent Langevin flow scenario. We show that
\begin{eqnarray}
\begin{aligned}
&\mathbb{D}_{\text{KL}}(p_{\text{data}}(x) || p_\theta(x)) +\mathbb{D}_{\text{KL}}(\hat{p}_{\theta_t}(z|x) || p_\theta(z|x))\\
= &\mathbb{D}_{\text{KL}}(p_{\text{data}}(x) \hat{p}_{\theta_t}(z|x) || p_\theta(x) p_\theta(z|x))\\
=&\mathbb{D}_{\text{KL}}(p_{\text{data}}(x) \hat{p}_{\theta_t}(z|x) || p_\alpha(z) p_\beta(x|z)) \\
= &\mathbb{D}_{\text{KL}}(p_{\text{data}}(z) \hat{p}_{\theta_t}(x|z) || p_\alpha(z)p_\beta(x|z))\\
= &\mathbb{D}_{\text{KL}}(p_{\text{data}}(z) || p_\alpha(z)) + \mathbb{D}_{\text{KL}}(\hat{p}_{\theta_t}(x|z) || p_\beta(x|z)),
\end{aligned}
\end{eqnarray}
where $p_{\text{data}}(z) = \int p_{\text{data}}(x) \hat{p}_{\theta_t}(z|x) dx$ is the aggregated posterior.  $\hat{p}_{\theta_t}(x|z) = p_{\text{data}}(x)\hat{p}_{\theta_t}(z|x)/p_{\text{data}}(z) $ is the generation model based on data distribution and the inference model (i.e., short-run Langevin flow). KL-divergence between conditional distributions is understood to be averaged over the variable $z$ being conditioned upon. In the above, we update $p_\alpha(z)$ to be close to $p_{\text{data}}(z)$. We update $p_{\beta}(x|z)$ to be close to $\hat{p}_{\theta_t}(x|z)$. Both $p_{\text{data}}(z)$ and $\hat{p}_{\theta_t}(x|z)$ depends on the short-run Langevin flow $\hat{p}_{\theta_t}(z|x)$.

Therefore, the learning of our model accomplishes two things: (1) Let the flow-based prior match the aggregated posterior produced by the short-run Langevin flow inference model. (2) Bias the model from MLE so that the short-run Langevin flow inference is close to the true posterior.

\begin{figure*}[ht]
\centering
\begin{subfigure}{.3\linewidth}
    \centering
    \includegraphics[width=0.95\textwidth]{figs/cifar10_mcmc_synthesis.png}
    \caption{CIFAR10}
\end{subfigure}
    \hspace{-0.5em}
\begin{subfigure}{.3\linewidth}
    \centering
    \includegraphics[width=0.95\textwidth]{figs/svhn_mcmc_synthesis.png}
    \caption{SVHN}
\end{subfigure}
  \hspace{-0.5em}
\begin{subfigure}{.3\linewidth}
    \centering
    \includegraphics[width=0.95\textwidth]{figs/celeba_mcmc_synthesis.png}
    \caption{CelebA}
\end{subfigure}
\caption{Generated examples from the LFBM-MCMC models traind on the CIFAR10 (32$\times$ 32), SVHN, (32$\times$ 32) and CelebA (64$\times$ 64) datasets. The LFBM-MCMC model is trained with short-run Langevin flow as approximate inference.}
\label{fig:sample_mcmc}
\end{figure*}

\section{Experiments}
\subsection{Experiment Settings} We conduct a set of experiments to examine the effectiveness of the proposed latent space flow-based prior model in image modeling. In particular, we evaluate the model from the perspectives of synthesis, reconstruction and inference.
We use SVHN~\cite{netzer2011reading}, CelebA~\cite{LiuLWT15}, and CIFAR-10~\cite{krizhevsky2009learning} image datasets for generation and reconstruction tasks, and use MNIST~\cite{mnist2010} dataset for anomaly detection. All the inpainting and recovery experiments are conducted on CelebA face dataset. The network architectures and the hyperparameters
are presented in Supplementary Material. The experiments are mainly for the sake of proof of concepts rather than state-of-the-art performance. We will compare our framework with two important baseline methods, which are (1) LFBM-VAE baseline, which trains LFBM using variational inference instead of MCMC inference, (2) LEBM baseline, which uses an energy-based model in the latent space. For fair comparison, the architecture designs of the top-down generators $g_{\beta}$ in LEBM, LFBM-VAE and ours are the same, and the  flow-based prior models $f_{\alpha}$ in LFBM-VAE and ours are also identical. The design of the inference network of LFBM-VAE follows \citet{kingma2016improved}, who parameterize the inference network by an encoder followed by an inverse autoregressive flow.
The designed network architectures are expected to vary in different datasets due to various image sizes and pattern complexities. We use LFBM-MCMC to denote our LFBM using Langevin inference in order to distinguish from the baseline LFBM-VAE.

\subsection{Image Synthesis and Reconstruction}
A well-trained latent variable generative model can be useful for generation and reconstruction. We train models on training images of SVHN, CelebA and CIFAR-10 respectively, and generate synthesized examples by first sampling latent vectors from the learned latent space normalizing flow and then transforming the vectors to image space. We calculate the Fréchet inception distance (FID)~\cite{HeuselRUNH17} to measure the quality of the synthesized images in Table~\ref{tab:compare}.
We show generated samples of the learned LFBM-MCMC models in Figure~\ref{fig:sample_mcmc}.
The models can reconstruct images by first inferring the latent vectors from the images, and then mapping the inferred latent vectors back to data space. The inference of latent variables can be achieved by the MCMC in our LFBM-MCMC. The quality of reconstruction images are measured by mean squared error (MSE) in Table~\ref{tab:compare}. We mainly compare our methods with likelihood-based top-down models rather than adversarial frameworks. Specifically, we compare with traditional variational inference baselines, such as VAE~\cite{KingmaW13} and SRI~\cite{NijkampP0ZZW20} which use Gaussian priors and two VAE variants, 2sVAE and RAE, which learn their priors from posterior samples in the second stage. We also compare with MCMC inference frameworks, such as ABP~\cite{HanLZW17} and LEBM~\cite{Pang0NZW20} whose prior distributions are fixed Gaussian distribution and energy-based model, respectively. We also compare our LFBM-MCMC with the VAE variant, i.e., the LFBM-VAE. Our models outperform other baselines in terms of MSE and FID.

\begin{table*}[]
\begin{center}
\resizebox{.9\textwidth}{!}{
\begin{small}
\begin{tabular}{cc|ccccccc|cc}
\hline
\multicolumn{2}{c|}{\multirow{2}{*}{Models}}   & \multirow{2}{*}{VAE}    & \multirow{2}{*}{2sVAE}  & \multirow{2}{*}{RAE}  &  \multirow{2}{*}{SRI}   & \multirow{2}{*}{SRI (L=5)} & \multirow{2}{*}{ABP} & \multirow{2}{*}{LEBM}   & \multicolumn{2}{c}{LFBM}   \\
  & &     &   &   &    &  &  &   & VAE  & MCMC \\ \hline
\multirow{2}{*}{SVHN}    & MSE & 0.019  & 0.019  & 0.014 &  0.018 & 0.011  &  - & 0.008 &     0.005&\textbf{0.005}  \\
                         & FID & 46.78  & 42.81  & 40.02 &  44.86 & 35.23   & 49.71 & 29.44 & 24.96     &\textbf{23.64} \\ \hline
 \multirow{2}{*}{Cifar10} & MSE & 0.057  & 0.056  & 0.027  &-     & -      & 0.018  & 0.020 & 0.020 & \textbf{0.016}
                          \\
                          & FID & 106.37 & 72.90 & 74.16 &  -     & -       & 90.30 & 70.15 & 69.70  & \textbf{66.41}  \\ \hline
\multirow{2}{*}{CelebA}  & MSE & 0.021  & 0.021  & 0.018  & 0.020 & 0.015  &  - & 0.013 &   0.014  & \textbf{0.011} \\
                         & FID & 65.75  & 44.40  & 40.95  & 61.03 & 47.95   & 51.50 & 37.87 & 33.64     & \textbf{33.64} \\ \hline
\end{tabular}
\end{small}}
\end{center}
\caption{Quantitative results of image reconstruction and generation on different datasets.
}
\label{tab:compare}
\vspace{-0.1cm}
\end{table*}

\subsection{Supervised Image Inpainting}
\label{sec:image_inpainting}
We can train an LFBM-MCMC from fully-observed training images, and then use the learned model to complete the missing pixels of testing images. Let $m$ be a matrix, which has the same number of dimension as that of an image $x$, with values ones indicating the visible pixels and zeros indicating the invisible ones (corrupted or occluded) of the image $x$, respectively. Suppose there is an incomplete image $x_{m}$ with missing pixels indicated by a mask $m$. With the learned flow-based prior model $p_{\alpha}(z)$ and the generation model $p_{\beta}(x|z)$, we can restore the missing pixels by first inferring the latent vectors $z$ from $x_{m}$, and then using $g_{\beta}(z_m)$ to generate an image $x$, which is a complete version of $x_m$. The inference can be easily performed by Langevin dynamics that follows Eq.(\ref{eq:MCMC_post}) with the gradient of a modified log-posterior $\nabla_{z} \log p_{\theta}(z|x_m,m)=\nabla_{z} \log p_{\alpha}(z)+ \nabla_{z} \log p_{\beta}(x_m|z,m) = \nabla_{z} l_{\alpha}(z)+ \frac{1}{\sigma^2}(m \odot (x_m-g_{\beta}(z)))\nabla_{z}g_{\beta}(z)$, where $\odot$ is the element-wise multiplication operator. That is, we evaluate the posterior over only the visible pixels of the images. Each testing image is occluded by a $40 \times 40$ region mask at the center. The resolution of the images are $64 \times 64$ pixels. For each testing image, our LFBM-MCMC model can generate diverse and meaningful inpainting results by using different $z_0$ sampled from $q_0(z)$ to initialize the Langevin chains for inference.

\begin{figure*}[t!]
\centering
\includegraphics[width=0.955\textwidth]{figs/impant_4rows2.png}
\caption{Supervised image inpainting results on the CelebA dataset. Images in the first column are the original images. Images in the second column are the masked images to be inpainted. Images  in column 3 to column 12 (yellow panel) are inpainting results using the learned LFBM-MCMC model. Images in column 13 to column 22 (green panel) are inpainting results using the trained LFBM-VAE. For each panel, different columns correspond to different initializations of the inference process.}
\label{fig:pretrained_inpaint}
\end{figure*}

Figure \ref{fig:pretrained_inpaint} displays some qualitative results, where we compare our LFBM-MCMC model with the baseline LFBM-VAE. Both LFBM-MCMC and LFBM-VAE are trained on fully observed training images of CelebA dataset. For task of inpainting, given an incomplete image with a mask occluding pixels at the center, our LFBM-MCMC first infers the latent variables of the incomplete image by the Langevin flow, and then generate the occluded region of the image by the top-down generator taking the inferred variables as input. As to the LFBM-VAE, the learned inference model in the LFBM-VAE is hard to employ for inferring latent variables from an incomplete data, which means that even though the variational inference is computationally efficient, it is not suitable to infer latent variables from a partially observed image because the encoder can not map a portion of the image into the latent space. In order to use the LFBM-VAE for image inpainting, we abandon the learned inference network, and derive the Langevin dynamics from the learned prior and generator models in LFBM-VAE. We adopt the derived Langevin flow as in the LFBM-MCMC for inference and inpainting. As shown in Figure \ref{fig:pretrained_inpaint}, each row illustrates one inpainting task. The first column displays the original images, and the second column shows the testing images that need to be inpainted. The yellow panel shows inpainted images by the LFBM-MCMC, while the green panel shows inpainted images by the LFBM-VAE. Different columns in each panel show different inpainting results due to the different initialization (i.e., randomness) of the inference process. Although the LFBM-VAE adopts the MCMC inference for inpainting, its results are still not as good as those by our LFBM-MCMC. Some obvious artifacts are observed in the inpainting results of the LFBM-VAE. Further, we quantitatively evaluate the image inpainting performance of both LFBM-MCMC and LFBM-VAE in Table~\ref{tab:image_inpainting}. We use two metrics, i.e., FID and MSE, to measure the quality of the inpainting performance. We test the models in 10,000 images. For the MSE, we compute the per pixel difference between the inpianting result and the ground truth image within the masked region. Our model outperforms the baseline.

\begin{table}[h!]
\centering

\begin{tabular}{ccc}
\hline
Method & FID & MSE \\ \hline
LFBM-MCMC (ours) & \textbf{0.0557} & \textbf{0.1937}\\
LFBM-VAE  & 0.0562 & 0.2074 \\
\hline
\end{tabular}
\caption{Quantitative results of image inpainting.}
\label{tab:image_inpainting}
\end{table}

\subsection{Anomaly Detection} We evaluate our generative model on anomaly detection of MNIST~\cite{mnist2010} data. Given a latent variable model that is well-trained on normal examples, we can perform anomaly detection on a testing image $x$ by firstly inferring its latent variables $z$ and then computing the logarithm of the joint probability $\log p_{\theta}(x, z)=\log p_{\alpha}(z) + \log p_{\beta}(x|z)=l_{\alpha}(z)-\frac{1}{2\sigma^2}||x-g_{\beta}(z)||^2- \log \sigma \sqrt{2 \pi}$ as a decision score. The score should be high for a normal example and low for an anomalous one. We can see that a normal example can be well-reconstructed by the learned generator (i.e., small reconstruction error $||x-g_{\beta}(z)||^2$) with a correctly inferred latent variables (i.e., high log likelihood in the learned prior distribution $l_{\alpha}(z)$).

Following the same experiment setup as~\citet{zenati2018efficient,kumar2019maximum, Pang0NZW20}, we treat each class of digit images in MNIST dataset as anomaly examples and the remaining 9 classes of digit images as normal examples. We train models only on the normal examples, and then test the models on both the normal and anomalous examples. To quantify the model performance, we compute the area under the precision-recall curve (AUPRC) based on the decision function $\log p_{\theta}(x, z)$. We report the mean and variance of AUPRC scores over 10 runs of each experiment.We show the results of our LFBM in Table~\ref{tab:anormaly}, and compare with the related models, including the VAE~\cite{KingmaW13}, MEG~\cite{kumar2019maximum}, BiGAN-$\sigma$~\cite{zenati2018anomaldetection}, EBM-VAE~\cite{Han2020VAE-EBM}, LEBM~\cite{Pang0NZW20} and ABP model~\cite{HanLZW17}. From Table~\ref{tab:anormaly}, we can find that the proposed LFBM can obtain much better results than those of other~methods.

\begin{table*}[]
\begin{center}
\resizebox{1.02\textwidth}{!}{
\begin{small}
\begin{tabular}{c|ccccc}
\hline
Heldout Digit   & 1 & 4 & 5 & 7 & 9 \\  \hline
VAE~\cite{KingmaW13}  & 0.063 & 0.337 & 0.325 & 0.148 & 0.104 \\
MEG ~\cite{kumar2019maximum} & 0.281 $\pm$ 0.035 & 0.401 $\pm$ 0.061 & 0.402 $\pm$ 0.062 & 0.290 $\pm$ 0.040 & 0.342 $\pm$ 0.034 \\
BiGAN-$\sigma$~\cite{zenati2018anomaldetection}   & 0.287 $\pm$ 0.023 & 0.443 $\pm$ 0.029 & 0.514 $\pm$ 0.029 & 0.347 $\pm$ 0.017 & 0.307 $\pm$ 0.028 \\

  EBM-VAE~\cite{Han2020VAE-EBM} & 0.297 $\pm$ 0.033 & 0.723 $\pm$ 0.042 & 0.676 $\pm$ 0.041 & 0.490 $\pm$ 0.041 & 0.383 $\pm$ 0.025 \\

LEBM ~\cite{Pang0NZW20} & 0.336 $\pm$ 0.008 & 0.630 $\pm$ 0.017 & 0.619 $\pm$ 0.013 & 0.463 $\pm$ 0.009 & 0.413 $\pm$ 0.010 \\
ABP~\cite{HanLZW17}  & 0.095 $\pm$ 0.028 & 0.138 $\pm$ 0.037 & 0.147 $\pm$ 0.026 & 0.138 $\pm$ 0.021 & 0.102 $\pm$ 0.033 \\ \hline
LFBM (ours)  & \bf{0.349 $\pm$ 0.002} & \bf{0.812 $\pm$ 0.007} & \bf{0.823 $\pm$ 0.009} & \bf{0.682 $\pm$ 0.004} & \bf{0.514 $\pm$ 0.008} \\
\hline
\end{tabular}
\end{small}}
\end{center}
\caption{AUPRC scores (larger is better) for unsupervised anomaly detection.
Results are averaged over 10 experiments.}
\label{tab:anormaly}
\end{table*}

\subsection{Analysis of Hyperparameters}
In Figure~\ref{ablation_study}(a), we show the effects of using different Langevin step sizes and numbers of steps for LFBM. We carry out experiments on SVHN dataset. We can see that the best FID score is achieved when the number of Langevin steps is around 40 and the model with 20 steps can already give us good results. On the other hand, the optimal step sizes for Langevin flows using different numbers of steps are different. Figure~\ref{ablation_study}(b) displays the effect of tuning latent size. We can see the best result can be achieved in the range from 50 to 100. Figure~\ref{ablation_study}(c) shows FID scores using different depths of the flow-based prior on Cifar10. We find that increasing the depth of the prior can improve the performance.

\begin{figure*}[h]
	\centering
	\begin{subfigure}{.28\linewidth}
	\includegraphics[width=1\linewidth]{figs/stepsize_n.png}
	\caption{number of steps and step size}
	\end{subfigure}
	\begin{subfigure}{.28\linewidth}
	\includegraphics[width=1\linewidth]{figs/latent_size_n.png}
	\caption{number of latent dimension}
	\end{subfigure}
	\begin{subfigure}{.31\linewidth}
	\includegraphics[width=1\linewidth]{figs/f_depth_cifar.png}
	\caption{depth of normalizing flow}
	\end{subfigure}
\caption{FIDs of LFBMs using different (a)  numbers of Langevin steps and step sizes, (b)  latent sizes, and (c) depths of the normalizing flow prior model.
}
\label{ablation_study}
\end{figure*}

\subsection{Unsupervised Image Recovery}

The LFBM can be learned from incomplete training data, e.g., images with occluded pixels. The learning algorithm updates the model parameters by maximizing the likelihood of the visible pixels in training images. This also belongs to unsupervised image inpainting, which is essentially different from the task of supervised image inpainting shown in Section~\ref{sec:image_inpainting}. Traditional VAE-based frameworks, such as LFBM-VAE, are incapable of learning from incomplete data. To demonstrate this ability of the LFBM, we experiment on 10,000 occluded images that we create. They are selected from the CelebA dataset
 and occluded with different types of masks. A mask is randomly placed in each image. Two types of masks are designed, including single region mask and salt-and-pepper mask. We specify three different sizes of single regions, e.g., $20 \times 20$, $30 \times 30$, and $40 \times 40$, and three different occlusion percentages of salt-and-pepper
masks, e.g., $30\%$, $50\%$, and $70\%$. We compare our LFBM with a fixed Gaussian prior and an energy-based prior learned by MCMC inference. We can evaluate the performance from two aspects: (i) Recovery quality: The incomplete training images are gradually recovered as the learning algorithm proceeds. MSE between recovered images and the corresponding ground truth images at masked regions are calculated; (ii) Generation quality: We use FID to measure the visual quality of the synthesized images, which are generated by models learned from incomplete data. Table~\ref{tab:MSE_mask_recovery} and Table~\ref{tab:FID_mask_recovery} show comparisons of our method with the baselines using EBM prior and Gaussian prior in terms of MSE and FID, respectively.
 Figure~\ref{fig:rec_recovery} shows the recovery results for learning from incomplete training images of the CelebA dataset. We further evaluate the generation capacities of the models trained in this scenario. We present qualitative results in Figure~\ref{fig:recovery_gen}, which shows the randomly generated images by the models learned from incomplete images with different occlusion levels. We can see that frameworks using Gaussian priors have the worst generated results. Our LFBM model using Langevin flow as an inference process can generate more realistic images than other baselines in the task of unsupervised learning.

\begin{figure*}[th]
    \centering
     \includegraphics[width=0.935\textwidth]{figs/inpainting.jpg}
     \caption{A comparison of unsupervised image recovery results by different methods on training images with different levels of occlusions. In each panel, the first row shows some original images that are used in the training process, the second row shows the corresponding occluded images with a certain occlusion level, and the third, fourth and fifth rows show the recovered images by models using Gaussian prior, EBM prior and normalizing flow prior, respectively.}
    \label{fig:rec_recovery}
\end{figure*}

\begin{figure*}[h!]
    \centering
     \includegraphics[width=0.93\textwidth]{figs/inpainting_gen.jpg}
    \caption{Image synthesis by models learned from incomplete images. Each panel represents a different level of occlusion.}
    \label{fig:recovery_gen}
\end{figure*}

\begin{table}[h!]
\centering
\begin{tabular}{cccc}
\hline
\multicolumn{4}{c}{Salt and pepper mask}                     \\ \hline
Occ \%    & 30\%           & 50\%           & 70\%           \\ \hline
flow (ours)      & \textbf{0.0244}  & \textbf{0.0317} & \textbf{0.0464} \\
EBM       & 0.0256          & 0.0319          & 0.0465          \\
Gaussian  & 0.0259         & 0.0326         & 0.0472         \\ \hline
\multicolumn{4}{c}{Single region mask}                       \\ \hline
mask size & $20 \times20$          & $30\times 30$          & $40 \times40$          \\ \hline
flow (ours)      & 0.0420 & 0.0587 & \textbf{0.0864} \\
EBM       & 0.0429         & 0.0684          & 0.0957          \\
Gaussian  & \textbf{0.0404} & \textbf{0.0572} & 0.0918         \\ \hline
\end{tabular}
\caption{MSEs of methods with different priors (e.g., flow, EBM and Gaussian) in unsupervised image recovery.}
\label{tab:MSE_mask_recovery}
\end{table}
\begin{table}[h!]
\centering
\begin{tabular}{cccc}
\hline
\multicolumn{4}{c}{Salt and pepper mask}                     \\ \hline
Occ \%    & 30\%           & 50\%           & 70\%           \\ \hline
flow (ours)     & \textbf{46.2}  & \textbf{59.14} & \textbf{86.77} \\
EBM       & 52.78          & 61.91          & 88.27          \\
Gaussian  & 153.01         & 156.71         & 172.77         \\ \hline
\multicolumn{4}{c}{Single region mask}                       \\ \hline
mask size & $20\times20$          & $30\times30$          & $40\times40$          \\ \hline
flow (ours)      & \textbf{42.39} & \textbf{47.52} & \textbf{72.47} \\
EBM       & 49.16          & 51.59          & 77.39          \\
Gaussian  & 150.95         & 146.41         & 184.53         \\ \hline
\end{tabular}
\caption{FIDs of methods with different priors (e.g., flow, EBM and Gaussian) in unsupervised image recovery.}
\label{tab:FID_mask_recovery}
\end{table}

\section{Conclusion}

\paragraph{Summary} In this paper, we study modeling the latent space of data by the normalizing flow model and follow the philosophy of empirical Bayes to learn the latent space flow-based model (LFBM) from observed data. Specifically, we propose a novel top-down generative model with a latent space flow-based prior model. We propose a novel learning framework, in which the LFBM built on top of a top-down generative network serves as the prior model of the latent space, and is trained simultaneously with the top-down network in an MCMC-based MLE algorithm. We show that the learning algorithm with a short-run Langevin flow is a perturbation of MLE. The LFBM is more flexible and informative than Gaussian prior and less biased to MLE than the unnormalized energy-based prior (LEBM). We test our LFBM on a variety of tasks to validate its effectiveness.
\paragraph{Limitation and Future Work} In contrast to EBM prior, the normalizing flow prior allows for efficient and unbiased sampling and training of the prior distribution. However, these good properties do not come for free, and the usage of the flow-based prior implies the assumption that the latent variables are computationally efficient to normalize and can be generated by a finite sequence of invertible transformations. This assumption might potentially limit the expressive power of the prior distribution for the latent space. Given the fact that EBM prior is an unnormalized probability density and doesn’t constrain the prior distribution by invertible functions, we may consider combining the representational flexibility of the EBM and the computational tractability of the normalizing flow model to tackle the above limitations. To be specific, we specify a latent space EBM (LEBM) and a latent space flow-based model (LFBM) together to represent the prior distribution. We  generalize the data space flow contrastive estimation (FCE) \cite{FlowEBM2020} to latent space for jointly training the LEBM and the LFBM without relying on MCMC. Flow contrastive estimation \cite{FlowEBM2020} can be considered as an improved version of noise contrastive estimation (NCE) \cite{GutmannH10} for training EBM, where the original Gaussian noise is replaced by a normalizing flow to make it closer to the target distribution and a stronger contrast to the EBM. Meanwhile, the flow model is updated by approximately minimizing the Jensen-Shannon divergence between the flow model and the target distribution. In our future work, we will learn both LFBM and LEBM together in the latent space via FCE.

\section*{Acknowledgments}
 The authors sincerely thank Dr. Ying Nian Wu at statistics department of  the University of California, Los Angeles (UCLA) for the helpful discussion on theoretical understanding of the learning algorithm. The authors would also like to thank the anonymous reviewers for providing constructive comments and suggestions to improve the work. Our work is also supported by XSEDE grant CIS210052.


\section{\label{introduction}Introduction}

The discovery of $X(3872)$ opened a new era of exotic
hadrons called charmonium-like states~\cite{X3872-Belle}.
Understanding the nature of this state and of other charmonium-like
states, in general, provides an opportunity to study the nonperturbative
behaviors of quantum chromodynamics.
In searching for other charmonium-like states,
$X(3915)$ was found by the Belle experiment~\cite{X3915-0,X3915-1}
and confirmed by the BaBar experiment~\cite{X3915-2,X3915-3},
initially in the study of
$B^- \rightarrow J/\psi\omega K^-$~\footnote{Charge-conjugate process
will be implicitly included.}
and later in the no-tag two-photon interactions.
This state, $X(3915)$, was classified as $\chi_{c0}(3915)$ in the latest
listing by the Particle Data Group~\cite{PDG-2},
but the assignment is not firmly established.
Spin-parity of $X(3915)$ is consistent with $J^P = 0^+$ based on
the experimental analysis~\cite{X3915-3}; it has also a small
possibility of $2^+$~\cite{PDG,Zhou}.
If $X(3915)$ is a conventional $c\bar{c}$ state, it should
also decay to $D^{(*)}\bar{D}^{(*)}$ or its charge conjugate.
In an amplitude analysis of $B^-\rightarrow K^-D^+D^-$ by the LHCb
experiment, $0^{++}$ and $2^{++}$ states near 3.930~GeV/$c^2$
are reported~\cite{LHCb-DD}; they are assigned
as $\chi_{c0}(3930)$ and $\chi_{c2}(3930)$, respectively.
However, no peaks in $M(D\bar{D}^{(*)})$ have been seen in the studies of
$B^- \rightarrow D\bar{D}K^-$ and $B^- \rightarrow D\bar{D}^*K^-$
performed by the B-factories~\footnote{Kinematically
allowed decay is only $D\bar{D}$ if the mass is 3.920~GeV/$c^2$
and $J=0$. Also in $D\bar{D}^*$, no peak is seen.}
\cite{Belle-Ds0D0,DDbar-1,DDbar-3,DDbar-2}.
Non-$c\bar{c}$ models such as $c\bar{c}s\bar{s}$ models
or $D_s\bar{D}_s$ molecule models can predict this
result~\cite{Exotics1,Exotics2,Exotics3}.
The $0^{++}$ state(s) reported by LHCb and B-factories can
be different states: $\chi_{c0}(2P)$ and non-$c\bar{c}$ state,
respectively.

In this paper, we report on a study of the production of $X(3915)$
by highly virtual photons, $\gamma^*$. The reaction used is
$\gamma^*\gamma\rightarrow X(3915) \rightarrow J/\psi\omega$,
where $\omega$ decays to $\pi^+\pi^-\pi^0$ and $J/\psi$ decays to
either $e^+e^-$ or $\mu^+\mu^-$, shown in Fig.~\ref{fig:feyn}.
The highly virtual photon is identified by tagging either $e^-$ or $e^+$
in the final state where its partner, $e^+$ or $e^-$, respectively, is missed
going into the beam pipe. This type of interaction is referred to as a
``single-tag" two-photon interaction.
If $X(3915)$ is a non-$c\bar{c}$ state,
naively it should have a larger
spatial size than $c\bar{c}$. This larger size is predicted
for charm-molecule models~\cite{Exotics3,Exotics4}.
In such case, its production rate should
decrease steeply at high virtuality.
To test a deviation from a pure $c\bar{c}$,
we use a reference $c\bar{c}$ model calculated
by Schuler, Berends, and van Gulik (SBG)~\cite{SBG}.
In this test, we use the parameter $Q^2$, appearing in its production,
where $Q^2 (= - q^2)$ is the negative mass-squared
of the virtual photon;
$q$ is the four-momentum of the virtual photon.

\begin{figure}[ht]
\centering
  \includegraphics*[width=0.32\textwidth]{feyn.pdf}
  \caption{\label{fig:feyn}Single-tag two-photon $X(3915)$ production.
Virtual photon, $\gamma^*$, is produced in the tagging side;
$q$ is the four-momentum of the $\gamma^*$.
$W$ is the energy of the two-photon
system in its rest frame which corresponds to the invariant
mass of $J/\psi\omega$, $M(J/\psi\omega)$, in this case.
Tagging is either $e^-$ or $e^+$.
  }
\end{figure}

We will use the term ``electron'' for either electron or positron.
Quantities calculated in the initial-state $e^+e^-$
center-of-mass (c.m.) system are indicated by an asterisk(*).

\section{\label{detector}Detector and Data}

The analysis is based on 825~fb\textsuperscript{-1} of data collected by the
Belle detector operated at the KEKB $e^+e^-$ asymmetric
collider~\cite{KEKB,PTEP-2}. The data were taken at the
$\Upsilon(nS)$ resonances ($n\le 5$) and nearby energies,
9.42~GeV $< \sqrt{s} <$ 11.03~GeV.
The Belle detector was a general-purpose magnetic spectrometer
asymmetrically enclosing the interaction point (IP) in
4$\pi$~\cite{Detector,PTEP-1}.
Charged-particle momenta are measured by a silicon vertex detector and a
cylindrical drift chamber (CDC).
Electron and charged-pion identification relies on a combination
of the drift chamber,
time-of-flight scintillation counters (TOF),
aerogel Cherenkov counters (ACC),
and electromagnetic calorimeter
(ECL) made of CsI(Tl) crystals.
Muon identification relies on resistive plate chambers
(RPC) in the iron return yoke.
Photon detection and energy measurement utilize ECL.

We use Monte Carlo (MC) simulations to set selection
criteria and to derive the reconstruction efficiency.
Signal events, $e^+e^-\rightarrow e^{\pm}(e^{\mp})(\gamma^*\gamma
\rightarrow J/\psi\omega)$, are generated using
TREPSBSS~\cite{TREPSBSS-1,TREPSBSS-2} with a mass
distribution, centered at $M=3.918$~GeV/$c^2$ and
width $\Gamma=0.020$~GeV/$c^2$~\cite{PDG}, with constant
transition form factor, $F(Q^2)$=const.
For $\omega$ decays, the Dalitz
method~\cprotect\footnote{
\url{https://evtgen.hepforge.org/doc/models.html}, entry EvtOmegaDalitz}
is used. Radiative $J/\psi$ decays are
simulated by PHOTOS~\cite{PHOTOS-1,PHOTOS-2}.
Detector response is simulated employing
GEANT3~\cite{Geant}.

\section{\label{PID}Particle Identification}

Final-state particles in this reaction are
$\ell^+\ell^-\pi^+\pi^-\gamma\gamma$ where $\ell^+\ell^-$ is either
an electron pair or a muon pair.

Electrons are identified using a combination of five
discriminants:
$E/p$, where $E$ is the energy measured by ECL
and $p$ is the momentum of the particle, then, transverse shower shape
in ECL, position matchings between the energy cluster and the extrapolated
track at ECL, ionization loss in CDC, and light yield in ACC.
For these, probability density functions are derived
and likelihoods, $L_i$'s, are calculated, where $i$'s stand
for the discriminants. Electron likelihood ratio, ${\cal L}_e$, is
obtained by $\Pi_i L_i^{\rm electron}/(\Pi_i L_i^{\rm electron} +
\Pi_i L_i^{\rm non-electrons})$~\cite{eID}.

Muons are identified using a combination of two measurements:
penetration depth in RPC, and deviations of hit-positions in RPC
from the extrapolated track.
From these, muon likelihood ratio, ${\cal L}_{\mu}$, is
obtained by $P_{\mu}/(P_{\mu}+P_{\pi}+P_K)$, where $P_{\mu}$,
$P_{\pi}$, and $P_K$ are probabilities for muon, pion, and kaon,
respectively~\cite{muonID}.

Charged pions and kaons are identified using the combination of
three measurements: ionization loss in CDC, time-of-flight by TOF,
light-yield in ACC. From these, pion likelihood ratio,
${\cal L}_{\pi}$, is calculated by $P_{\pi}/(P_K + P_{\pi})$ where $P_{\pi}$
and $P_K$ are pion and kaon probabilities, respectively~\cite{PID}.

Photons are identified by position isolations
between the energy cluster and the extrapolated track at ECL.

\section{\label{selection}Event Selection}

Event-selection criteria share the ones in our previous
publication~\cite{X3872-our}.
We select events with five charged tracks coming from the IP
since one final-state electron goes into the beam pipe
and stays undetected.
Each track has to have $p_T > 0.1$~GeV/$c$, with two or more
having $p_T > 0.4$~GeV/$c$, where $p_T$ is the transverse
momentum with respect to the $e^+$ direction.
Total charge has to be $\pm 1$.

$J/\psi$ candidates are reconstructed by their decays to lepton pairs:
$e^+e^-$ or $\mu^+\mu^-$. Electron is identified by requiring
${\cal L}_e$ to be greater than 0.66 having 90\% efficiency.
Similarly, muon is identified by requiring ${\cal L}_{\mu}$
to be greater than 0.66 having 80\% efficiency.
We require the invariant mass of the lepton
pair to be in the range [3.047 GeV/$c^2$; 3.147 GeV/$c^2$].
In the calculation of the invariant mass of an $e^+e^-$ pair, we
include the four-momenta of radiated photons if the photons have
energies less than 0.2~GeV and polar angles, relative to
the electron direction at the IP, less than 0.04~rad.

For the tagging electron, a charged track has to satisfy
${\cal L}_e$ greater than 0.95 or $E/p$ greater
than 0.87.
In addition, we require $p > 1.0$~GeV/$c$ and $p_T > 0.4$~GeV/$c$.
In the calculation of $p$, four-momenta of radiated photons are
included using the same requirements as for the electrons
from $J/\psi$ decays.

Charged pions are identified by satisfying its
${\cal L}_{\pi}$ be greater than 0.2,
${\cal L}_{\mu}$ less than 0.9,
${\cal L}_e$ less than 0.6
and $E/p$ less than 0.8, having 90\% efficiency.

Neutral pions are reconstructed from their decay photons,
where the photons are identified as energy clusters
in the electromagnetic calorimeter and isolated from charged tracks.
These photons have to fulfill the requirements
$E_{\gamma{\rm H}}<-7 E_{\gamma{\rm L}}
+ 0.54$~GeV and $E_{\gamma{\rm H}}>0.12$~GeV, where $E_{\gamma{\rm H}}$ is
the energy of the higher-energy photon, and
$E_{\gamma{\rm L}}$ is the energy of the lower-energy photon, both in GeV.
Neutral-pion candidates have to satisfy $\chi^2 < 9$
for their mass-constraint fit.
If there is only one $\pi^0$ candidate with $p_T > 0.1$~GeV/$c$,
we accept the one as $\pi^0$. If there is no such $\pi^0$,
but there are one or more
$\pi^0$ candidates with $p_T < 0.1$~GeV/$c$, we calculate the invariant mass,
$M(\pi^+\pi^-\pi^0)$, for each $\pi^0$ candidate. If there is only one
candidate having its $M(\pi^+\pi^-\pi^0)$ in the $\omega$-mass region
[0.7326~GeV/$c^2$; 0.8226~GeV/$c^2$], we accept the one as $\pi^0$.
If more than one candidate satisfy the $\omega$-mass condition,
we accept the one with the smallest mass-constraint fit $\chi^2$ as $\pi^0$.
If there are more than one $\pi^0$ candidate with $p_T > 0.1$~GeV/$c$,
we test the $\omega$-mass condition for each $\pi^0$ candidate.
If there is only one candidate that satisfies the $\omega$-mass condition,
we accept it as $\pi^0$. If more than one such candidate exist,
we accept the one with the smallest mass-constraint fit $\chi^2$ as $\pi^0$.

Events should not have $e^+e^-$ pairs from $\gamma\rightarrow e^+e^-$.
Therefore, we discard the event if the invariant mass of the pair of
any oppositely charged tracks is less than 0.18~GeV/$c^2$,
calculated assuming them as electrons.
We require that the event has no photon with energy above 0.4~GeV.
Events must have one $\omega$ identified by the $\omega$-mass condition.

The tagging electron and the rest of the particles should
be back-to-back, projected in the plane perpendicular to the
$e^+$ beam axis. For this, we require $||\phi({\rm tag})
- \phi({\rm rest \ combined})| - \pi | < 0.15$~rad, where $\phi$
is the azimuthal angle about the $e^+$ beam axis.

A missing momentum arises from the momentum of the final-state
electron that goes undetected into the beam pipe. We require the
missing-momentum projection in the $e^-$ beam direction in the c.m. system
to be less than $-0.2$~GeV/$c$ for $e^-$-tagging events and greater than
0.2~GeV/$c$ for $e^+$-tagging events.

The total visible transverse momentum of the event, $p^*_T$, should
be less than 0.2~GeV/$c$. Measured energy of
the $J/\psi\pi^+\pi^-\pi^0$ system, $E^*_{\rm obs}$,
should be equal to the expected energy, $E^*_{\rm exp}$,
calculated from the momentum of the tagging electron and the direction
and invariant mass of the $J/\psi\pi^+\pi^-\pi^0$ system.
Since energy and $p^*_T$ are correlated, we impose a
two-dimensional criterion
\begin{equation}\label{eq:hyper}
\begin{array}{l}
{\displaystyle
  (p^*_T + 0.04~{\rm GeV}/c)\left(\frac{|E^*_{\rm obs}-E^*_{\rm exp}|}
{E^*_{\rm exp}}+0.003\right)
} \\
\multicolumn{1}{r}{<0.012~{\rm GeV}/c.}
\end{array}
\end{equation}
Figure~\ref{fig:hyper} shows the $p^*_T$ vs.~$E^*_{\rm obs}/E^*_{\rm exp}$
distribution from MC events with the selection criteria.

\begin{figure}[ht]
\centering
  \includegraphics*[width=0.48\textwidth]{pub_hyper.pdf}
  \caption{\label{fig:hyper}$p^*_T$ vs.~$E^*_{\rm obs}/E^*_{\rm exp}$
distribution (MC events). The (black) line shows the selection criteria
applied to $p^*_T$ and $E^*_{\rm obs}/E^*_{\rm exp}$; events below
the line are accepted.
  }
\end{figure}

A non-signal event imitates $X(3915)$ if a $\psi(2S)$ is produced by
a virtual photon from internal bremsstrahlung
and if it accompanies either a $\pi^0$ or a fake $\pi^0$
and also the $\pi^+\pi^-\pi^0$ combination satisfies the
$\omega$-mass condition.
To suppress this background, we reject the event having the invariant
mass of $J/\psi\pi^+\pi^-$ in the $\psi(2S)$
window [3.6806~GeV/$c^2$; 3.6914~GeV/$c^2$]. This window is defined as
$\pm 2\sigma$ of the $\psi(2S)$ mass resolution.
The mass resolutions of $\psi(2S)$(= $J/\psi\pi^+\pi^-$)
and $X(3915)$(= $J/\psi\pi^+\pi^-\pi^0$)
are approximately 2.7~MeV/$c^2$.

\section{\label{result}Results}
\subsection{\label{signal}Signals and Backgrounds}

Figure~\ref{fig:open-2D} shows the $Q^2$ vs.~$M(J/\psi\omega)$
distribution from the selected data. Here, $Q^2$ is calculated by
$Q^2 = 2(p_{\rm beam}p_{\rm tag} - m_e^2)$, where $p_{\rm beam}$ and
$p_{\rm tag}$ are the four-momenta of the beam $e^{\pm}$ and tagging $e^{\pm}$,
respectively, and $m_e$ is the electron mass.

\begin{figure}[ht]
\centering
  \includegraphics*[width=0.47\textwidth]{pub_Open-2D.pdf}
  \caption{\label{fig:open-2D}$Q^2$ vs.~$M(J/\psi\omega)$
distribution from data. Dashed (green) line indicates the kinematical
limit: 3.8795~GeV/$c^2$.
  }
\end{figure}

To derive the numbers of signal and background events,
we fit a combination of the threshold-corrected
relativistic Breit--Wigner (BW) function and
a constant to the  $M(J/\psi\omega)$ distribution.
The threshold-corrected BW function, $f_{\rm BW}(W)$, is

\begin{equation}\label{eq:std-BW}
{\displaystyle
f_{\rm BW}(W) =
\frac{\alpha M^2\Gamma'}{(W^2-M^2)^2+M^2\Gamma'^2}
}
\end{equation}
where $M$ is the resonance mass,
$\alpha$ is a dimensionless normalization factor,
and $\Gamma'$ is the threshold-corrected resonance width
defined by
\begin{equation}
{\displaystyle
\Gamma' = \Gamma \cdot \frac{\rho(W)}{\rho(M)}
}
\end{equation}
where $\Gamma$ is the resonance width,
$\rho(W)$ is the phase space factor for $W$, which is
\begin{equation}\label{eq:phase-space}
{\displaystyle
\rho(W) = \frac{1}{16\pi}\frac{\lambda^{1/2}(W^2, m_J^2, m_{\omega}^2)}{W^2}
}
\end{equation}
and $\lambda$ is the K{\"a}ll{\'e}n function~\cite{PDG-2,Kallen}.
It is defined as
\begin{equation}
\begin{array}{l}
{\displaystyle
\lambda^{1/2}(W^2, m_{J/\psi}^2, m_{\omega}^2)
} \\
{\displaystyle
  = \sqrt{(m_{J/\psi}^2 + m_{\omega}^2 - W^2)^2
    - 4m_{J/\psi}^2m_{\omega}^2}
}
\end{array}
\end{equation}
where $m_{J/\psi}$ is the mass of $J/\psi$(= 3.0969~GeV/$c^2$) and
$m_{\omega}$ that of $\omega$(= 0.78265~GeV/$c^2$)~\cite{PDG}.
In the fit, we set $M=3.918$~GeV/$c^2$,
$\Gamma = 0.020$~GeV/$c^2$~\cite{PDG}, and $\alpha = 2/\pi$,
with the fit function (modified BW combined with a flat component)
\begin{equation}\label{eq:BW-fit}
{\displaystyle
f_{\rm BW + flat} = a_{\rm BW}\cdot f_{\rm BW} + a_{\rm flat},
}
\end{equation}
where the fit parameters $a_{\rm BW}$ and $a_{\rm flat}$ are
the magnitudes of the BW and the flat component, respectively.
We ignore a possible distortion of the fit-distribution
due to the energy dependence of the detection sensitivity,
because the effect is small.
Energy dependence of the detection sensitivity for
$J/\psi\omega$, which is defined by the production of
detection efficiency times luminosity function, is estimated
as $0.1\Delta W$~\%, where $\Delta W$ is in the MeV unit.
We use the ROOT/MINUIT implementation of the binned maximum-likelihood
method with a 5~MeV/$c^2$ bin width and perform the fit in the
$M(J/\psi\omega)$ range of [3.880 GeV/$c^2$; 4.100 GeV/$c^2$].
The units of $f_{\rm BW + flat}$ and $f_{\rm BW}$ are events/(5~MeV/$c^2$)
and (GeV/$c^2$)$^{-1}$, respectively.
The result of the fit is shown in Fig.~\ref{fig:BW_lin}.

\begin{figure}[ht]
  \centering
  \includegraphics*[width=0.47\textwidth]{pub_BW.pdf}
  \caption{\label{fig:BW_lin}$M(J/\psi\omega)$ distribution
with the Breit--Wigner + flat function fit.
Abscissa is $M(J/\psi\omega)$ in GeV/$c^2$. Ordinate is the
number of events per 5~MeV/$c^2$.
Solid (magenta) curve shows the result of the fit.
Horizontal dashed (magenta) line shows the flat component.
Vertical dashed (blue) lines indicate the fit
region.}
\end{figure}

The obtained parameters are
$a_{\rm BW} = 0.049\pm 0.018$~GeV/$c^2$/(5~MeV/$c^2$)
and $a_{\rm flat} = 0.022\pm 0.035$~/(5~MeV/$c^2$).
The number of signal events is $n_{\rm sig} = 9.0\pm 3.2$,
obtained by integrating $f_{\rm BW}$ with $a_{\rm BW}$ over the fit
region [3.8795~GeV/$c^2$; 4.1000~GeV/$c^2$].
The number of background events is $n_{\rm bg}^{\rm fit} = 0.3 \pm 0.4$,
calculated for the $X(3915)$ band, which we define 60~MeV/$c^2$.
It is obtained by multiplying $a_{\rm flat}$ by
the ratio of the $X(3915)$ band width, 60~MeV/$c^2$,
to the bin width, 5~MeV/$c^2$.

To confirm the number of background events, it is
also derived using the number of events in the $\omega$ sidebands.
Figure~\ref{fig:Mx-Mhad} shows the $M(\pi^+\pi^-\pi^0)$
vs.~$M(J/\psi\pi^+\pi^-\pi^0)$ distribution.
The sideband regions are set as two rectangles with heights
$M(\pi^+\pi^-\pi^0)$ in [0.60,0.70]~GeV/$c^2$ and
[0.83,0.93]~GeV/$c^2$ and the width $M(J/\psi\pi^+\pi^-\pi^0)$
in [3.88,4.10]~GeV/$c^2$. There are in total four events in the
$\omega$ sideband rectangles.
For the signal region, a rectangle of 0.080~GeV/$c^2$ high in
$M(\pi^+\pi^-\pi^0)$ and 0.060~GeV/$c^2$ wide in
$M(J/\psi\pi^+\pi^-\pi^0)$ is used.
From this, the obtained number of background events is
$n_{\rm bg}^{\omega}= 0.4 \pm 0.3$.
As $n_{\rm bg}^{\omega}$ is calculated using non-$\omega$
events while $n_{\rm bg}^{\rm fit}$ is obtained from identified $\omega$
events, the contents in the samples are different. Nevertheless,
the results from the two methods are approximately the same.
We use a conservative number: $n_{\rm bg}=0.4 \pm 0.4$.
The resulting signal significance for 9 observed events is then 5.6$\sigma$.

\begin{figure}[ht]
  \centering
  \includegraphics*[width=0.47\textwidth]{pub_Mx-Mhad.pdf}
    \caption{\label{fig:Mx-Mhad}$M(\pi^+\pi^-\pi^0)$
vs.~$M(J/\psi\pi^+\pi^-\pi^0)$ distribution to estimate the background rate
using $\omega$ sideband events.
Dashed (green) rectangles show the $\omega$ sidebands. Horizontal
(red) line-pair shows the $\omega$ signal band. Vertical (red) line-pair
shows the $X(3915)$ signal band.
    }
\end{figure}

The expected number of signal events, $n_{\rm sig}^{\rm exp} = 4.1\pm 0.7$,
is obtained from measurements by no-tag two-photon interactions,
$\Gamma_{\gamma\gamma}(0){\cal B}(X\rightarrow J/\psi\omega)=(54\pm 9)$~eV/$c^2$~\cite{PDG,PDG-2},
extrapolated to higher $Q^2$
using Eqs.~(\ref{eq:SBG-2}) and (\ref{eq:SBG-0+}), where
$\Gamma_{\gamma\gamma}(0)$
is the $\gamma\gamma$ decay width of $X(3915)$ at
$Q^2=0$, and ${\cal B}(X\rightarrow J/\psi\omega)$ is the
branching fraction of $X(3915)$ decaying to $J/\psi\omega$.
For a different prediction, if we assume the spin-parity of
$X(3915)$ as $J^P=2^+$, the expectation is $7.5\pm 1.3$
events using
$\Gamma_{\gamma\gamma}(0){\cal B}(X\rightarrow J/\psi\omega)=16$~eV
in Ref.~\cite{X3915-1} with the $J=2$ SBG model, Eq.~(\ref{eq:SBG-2+}),
assuming $\epsilon=1.0$.

\subsection{\textbf{\textit{Q}\textsuperscript{2}} Distribution}

Looking at the $Q^2$ distribution in Fig.~\ref{fig:open-2D}, the events
are classified into two groups: a cluster of events
below $Q^2=10$~(GeV/$c$)$^2$ and one event
at high $Q^2$ [$\approx 30$~(GeV/$c$)$^2$].
Since the high-$Q^2$ event is isolated, we discuss some of the
characteristics of this event in the following.

Figure~\ref{fig:psi2S-pT-omega}(a) shows the
$M(\pi^+\pi^-\pi^0)$ vs.~$M(J/\psi\pi^+\pi^-)$ distribution,
applying neither the $\omega$ selection
nor the $\psi(2S)$ veto.
The high-$Q^2$ event locates at 0.9~MeV/$c^2$
above the upper boundary of the $\psi(2S)$ veto. There are six
events in the $\psi(2S)$ veto. Of the six events, two pass the
$\omega$ selection.
Figure~\ref{fig:psi2S-pT-omega}(b) shows the $Q^2$
vs.~$M(J/\psi\pi^+\pi^-)$ distribution.
The two $\psi(2S)$-vetoed events, in addition to the high-$Q^2$ event,
have high $Q^2$s: $Q^2 > 10$~(GeV/$c$)$^2$.
All the other events that pass the $\psi(2S)$ veto have a lower $Q^2$, i.e.,
$Q^2 < 10$~(GeV/$c$)$^2$.
From this we conclude that $\psi(2S)$-vetoed events have
significantly higher $Q^2$ than the $X(3915)$ events.

As for a possibility of the high-$Q^2$ event being
an $X(3915)$ signal, the Belle experiment had little sensitivity
to measure single-tag two-photon events with $Q^2$ around 30~(GeV/$c$)$^2$
as detailed in the Appendix (see, e.g., Fig.~\ref{fig:eff-lum}).
Hence, it is improbable for the high-$Q^2$ event to be a
single-tag two-photon event.

\begin{figure}[ht]
  \centering
  \includegraphics*[width=0.47\textwidth]{pub_psi2S_pT_omega.pdf}
    \caption{\label{fig:psi2S-pT-omega}(a) $M(\pi^+\pi^-\pi^0)$
vs.~$M(J/\psi\pi^+\pi^-)$ and (b) $Q^2$
vs.~$M(J/\psi\pi^+\pi^-)$ distributions.
Star (red): high-$Q^2$ event. Open circles (blue): $\psi(2S)$-vetoed
events. Closed circles (blue): $X(3915)$ candidates.
(a) neither  $\omega$ selection nor $\psi(2S)$ veto are applied.
(b) events pass the $\omega$ selection but no $\psi(2S)$ veto is applied.
Vertical red lines indicate the $\psi(2S)$-veto window;
the horizontal red lines in panel (a) are the $\omega$ signal band.
}
\end{figure}

To estimate the probability of having one $\psi(2S)\pi^0$ event in the
region adjacent to the $\psi(2S)$ veto window, where the high-$Q^2$ event
is located, we estimate the probability of $\psi(2S)$ events escaping the
veto and having a $\pi^0$.
For this, we employ the data sample used in the
$X(3872)$ search and examine the $M(J/\psi\pi^+\pi^-)$
distribution~\cite{X3872-our}. There are 231 events
in the $\psi(2S)$-veto window of $\pm 5.4$~MeV/$c^2$ used
in the current study. There are 12 events in the 2.7~MeV bin,
adjacent to the upper boundary of the veto, where the high-$Q^2$
event locates. If we normalize the number of events in the veto window
to six that we observe as $\psi(2S)\pi^0$s in this study,
those 12 events correspond to 0.31 events/bin, or 0.11 events/MeV.
As seen in Fig.~\ref{fig:psi2S-pT-omega}(a), two out of six
events are inside the $\omega$ region. Hence, the expected number
of veto-leaks is 0.04~events/MeV. Then, by assuming the width
of the leak region as 2~MeV and the uncertainty in the number of events
as 0.1~events, the number of expected events is estimated
to be $0.1\pm 0.1$~events.
Significance of that number exceeding one event is 1.5~$\sigma$,
or 7\%.

A possible way of producing $\psi(2S)\pi^0$ is by a virtual photon,
radiated by internal bremsstrahlung from $e^-$ or $e^+$,
similar to the case of $\psi(2S)$ production.
However, there are suppressions to the $\psi(2S)\pi^0$ production compared
to $\psi(2S)$. The $\psi(2S)$s are produced as resonances,
but the $\psi(2S)\pi^0$s are not.
In order to be $J^P = 1^-$, the $\psi(2S)\pi^0$ has to be in a $P$-wave.
In addition, $\psi(2S)\pi^0$ is an isospin one state. Thus, further
suppressions are expected.

In the arguments up to this point, we assume the $\pi^0$s as real.
However, the reconstructed $\pi^0$s can be fake.
Using MC events, we observe that 13\% of $\pi^0$s, found in the
$X(3915)$ candidates, are fake.
This number is considered a lower limit as we found that the abundance of
low-$p_T$ $\pi^0$s is higher in real data than in MC.
Thus, the fraction of fake $\pi^0$s is higher at low $p_T$ than at
high $p_T$.
The observed $\psi(2S)\pi^0/\psi(2S)$ is 6/231, where the $\pi^0$s
are either real or fake.
In summary, it is plausible that the high-$Q^2$ event is
a $\psi(2S)\pi^0$ background.

If we remove the high-$Q^2$ event from the $f_{\rm BW + flat} $ fit,
the result is $a_{\rm BW} = 0.043\pm 0.017$~GeV/$c^2$/(5~MeV/$c^2$) and
$a_{\rm flat} = 0.025\pm 0.036$~/(5~MeV/$c^2$).
From that, we obtain $n_{\rm sig} = 7.9^{+3.1}_{-3.0}$.
As a note, the significance for 8 events is 5.2$\sigma$.

\begin{figure}[t]
  \centering
  \includegraphics*[width=0.47\textwidth]{pub_Q2_mod.pdf}
  \caption{\label{fig:Q2-dist}Measured $Q^2$ distributions:
(a) the number of events per 1~(GeV/$c$)$^2$ and
(b) $Z_{\gamma^*\gamma}{\cal B}(X\rightarrow J/\psi\omega)$.
The solid (red) curve shows the SBG prediction
based on the data of the no-tag two-photon measurement,
$\Gamma_{\gamma\gamma}(0){\cal B}(X\rightarrow J/\psi\omega)=54$~eV/$c^2$,
shown as a small (red) circle.
  }
\end{figure}

In the low-$Q^2$ region, there are 8 events in the
$M(J/\psi\pi^+\pi^-\pi^0)$ range [3.911~GeV/$c^2$; 3.958~GeV/$c^2$].
Figure~\ref{fig:Q2-dist} shows the $Q^2$ distributions both for
the number of events and for
$Z_{\gamma^*\gamma}{\cal B}(X\rightarrow J/\psi\omega)$,
where $Z_{\gamma^*\gamma}$ is a $Q^2$-dependent decay function
defined by Eq.~(\ref{eq:Zgamma}).
The $Z_{\gamma^*\gamma}{\cal B}(X\rightarrow J/\psi\omega)$ distribution
is obtained by multiplying the event distribution by a correction
function further detailed in the Appendix [see Eq.~(\ref{eq:C})].
The integrated yield of the
$Z_{\gamma^*\gamma}{\cal B}(X\rightarrow J/\psi\omega)$
distribution in the $Q^2$ range of 1.5~(GeV/$c$)$^2$ to 10.0~(GeV/$c$)$^2$
is $1.9 \pm 0.9$ times the expectation from the no-tag measurement,
$\Gamma_{\gamma\gamma}(0){\cal B}(X(3915)\rightarrow J/\psi\omega)=54\pm 9$~eV~\cite{PDG,PDG-2},
combined with its extrapolation to the higher $Q^2$ region
using Eq.~(\ref{eq:SBG-0+}).
The averages, $\langle Q^2 \rangle$,
and the root-mean-squared (RMS) values of $Q^2$,
$\sqrt{\langle(Q^2 - \langle Q^2\rangle)^2 \rangle}$, for the
$Z_{\gamma^*\gamma}{\cal B}(X\rightarrow J/\psi\omega)$ distribution
are listed in Table~\ref{tab:Q2-dist}, both for the measurement and for
a theoretical prediction based on a $c\bar{c}$
model [see Eq.~(\ref{eq:SBG-0+}) of the Appendix].
They are obtained from the same range in $Q^2$ as above, i.e.,
$1.5~(\text{GeV}/c)^{2} \le Q^2 \le 10.0~(\text{GeV}/c)^{2}$.
The measured average $Q^2$, $4.5\pm 0.7$~(GeV/$c$)$^2$, agrees with
the theoretical prediction, 4.8~(GeV/$c$)$^2$. Their difference
is approximately 10\% of the RMS widths of their distributions,
which are $1.9\pm 0.8$~(GeV/$c$)$^2$ vs. 2.4~(GeV/$c$)$^2$.
The resolution in $Q^2$  is about 0.03~(GeV/$c$)$^2$ depending
on the tag-electron's scattering angle.
Hence, the measurement is consistent with the prediction both in
the averages and the RMS values of $Q^2$.
In summary, the measured $Q^2$ distribution does not show a
significant shift to lower $Q^2$ as expected from some types
of non-$c\bar{c}$ models, but it agrees with
the SBG $c\bar{c}$ model.

\begin{table}[t]
  \centering
  \caption{\label{tab:Q2-dist}Comparison of the measurement and the SBG model
prediction for the $Q^2$ distribution of
$Z_{\gamma^*\gamma}{\cal B}(X\rightarrow J/\psi\omega)$.
$Q^2$ resolution is estimated to be about $0.03$~(GeV/$c$)$^2$.
Used are the 8 events shown in Fig.~\ref{fig:Q2-dist}.
  }
  \begin{tabular}{lrr}
  \hline
  \hline
  Item & Measurement & SBG model \\
  \hline
  Relative yield & $1.9 \pm 0.9$ & 1.0 \\
  $\langle Q^2 \rangle$ (GeV/$c$)$^2$ & $4.5\pm 0.7$ & 4.8 \\
  $\sqrt{\langle(Q^2 - \langle Q^2\rangle)^2 \rangle}$ (GeV/$c$)$^2$
	& $1.9\pm 0.8$ & 2.4 \\
  \hline
  \hline
  \end{tabular}
\end{table}

\section{\label{uncertainty}Systematic Uncertainties}

The largest uncertainty is associated with
the $\pi^0$ selection efficiency, including the rate of fake $\pi^0$s.
By comparing the number of selected events using the different
$\pi^0$ selection algorithms, we estimate 15\% uncertainty
associated with the $\pi^0$ selection algorithm.
Another uncertainty in $\pi^0$ detection is associated to
fake $\pi^0$s from background photons.
In the data before applying $\pi^0$ selection,
a significant number of low-energy photons, either true or fake,
contaminate.
These photons can produce fake $\pi^0$s.
To estimate the effect of such background photons, we look at variations in
the ratio of events with identified $\pi^0$(s) to all events
observed during the whole data-taking period.
From this, we estimate a 5.6\% uncertainty
after correcting the selection efficiency for the events with fake $\pi^0$.
This effect of background photons is also estimated using MC events
simulated with different background conditions, which gives a 3\% variation.
Conservatively, we use the larger 5.6\% as the systematic
uncertainty in $\pi^0$ identification due to
background photons.

Another large uncertainty is associated with $J/\psi$
identification.
The combined uncertainty in $J/\psi$ ID is 8\%. The largest
contribution, 7\%, to this comes from the difference in the ratio
of the number of $J/\psi$ selected events,
$N(J/\psi\rightarrow e^+e^-)/N(J/\psi\rightarrow \mu^+\mu^-)$,
between the real data and MC.
The other smaller contributions are
the uncertainties in the efficiencies of electron and muon IDs, background
levels and radiative $\gamma$ corrections in the case of
$J/\psi\rightarrow e^+e^-$ and the shapes of the invariant-mass
distributions.
They are estimated by the differences in the efficiencies
between the real data and MC by varying the selection conditions.

The uncertainties in electron tagging, 5\%,
and charged pion ID, 3\%, are estimated by
the difference in the efficiencies between real data and MC
by varying the selection conditions.
To calculate the detection efficiency, we set
the fit region for selecting signal events. Because of
the uncertainty in the $X(3915)$ distribution at or near the lower
boundary of the fit region, 3.888~GeV/$c^2$, detection efficiency
will have an uncertainty, which is estimated to be 3\%.
The uncertainties in the $\omega$ selection, 2\%, and
the $p_T$-and-$E_{\rm obs}^*/E_{\rm exp}^*$ selection specified
by Eq.~(\ref{eq:hyper}), 4\%, are estimated using MC by
varying selection condition.
The uncertainty in the luminosity function,
which is defined by Eq.~(\ref{eq:lum-func}), 3\%, is estimated
from the uncertainties in QED modeling and numerical integration.
The other uncertainties are 2\% for missing $p_T$,
2\% for $||\phi({\rm tag})-\phi({\rm rest})|-\pi|$,
1.8\% for track finding, 1.4\% for luminosity measurement,
1\% for $p_T < 0.2$~GeV/$c$, 1\% for $Q^2$ numerical integration,
1\% for energy dependence in the detection efficiency,
and 0.6\% for MC statistics.

Table~\ref{tab:systematics} lists a summary of systematic uncertainties.
As a total, combined quadratically, uncertainty in the reconstruction
efficiency is 20\%.
\begin{table}[ht]
  \centering
  \caption{\label{tab:systematics}Systematic uncertainties in the
reconstruction efficiency.
  }
  \begin{tabular}{lr}
  \hline
  \hline
  Item & Uncertainty \\
  \hline
  $\pi^0$ selection algorithm		& 15\% \\
  $J/\psi$ ID				&  8\% \\
  Fake $\pi^0$ by background		& 5.6\% \\
  Electron tagging			&  5\% \\
  $p_T$-and-$E^*_{\rm obs}/E^*_{\rm exp}$ selection &  4\% \\
  Charged pion ID			&  3\% \\
  Luminosity function			&  3\% \\
  Efficiency window			&  3\% \\
  $\omega$ selection 			&  2\% \\
  Missing $p_T$				&  2\% \\
  $||\phi({\rm tag})-\phi({\rm rest})|-\pi|$ &  2\% \\
  Luminosity measurement		& 1.4\% \\
  Track finding				& 1.8\% \\
  $p_T < 0.2$~GeV/$c$			&  1\% \\
  $Q^2$ numerical integration		&  1\% \\
  Energy dependence in efficiency	&  1\% \\
  MC statistics				& 0.6\% \\
  \hline
  \textbf{Total} 			& \textbf{20\%} \\
  \hline
  \hline
  \end{tabular}
\end{table}

\section{\label{summary}Summary}

We performed the first measurement of the $Q^2$ distribution of $X(3915)$
production in single-tag two-photon interactions.
For signals,
$7.9\pm 3.1({\rm stat.})\pm 1.5({\rm syst.})$
events are observed,
while the expectation is $4.1\pm 0.7$,
derived from the measured decay width at $Q^2=0$,
$\Gamma_{\gamma\gamma}(0){\cal B}(X\rightarrow J/\psi\omega) = 54\pm 9$~eV,
extrapolated to higher $Q^2$ region using the SBG $c\bar{c}$ model~\cite{SBG}.
The measured shape of $Q^2$ distribution does not show a significant
shift to lower $Q^2$ in contrast to the expectation
from some types of non-$c\bar{c}$ models. But, it agrees with the $X(3915)$
being a $c\bar{c}$ state or a compact non-$c\bar{c}$ state,
albeit within still very limited statistical precision.

\begin{acknowledgements}
\input{ack.txt}
\end{acknowledgements}