Variational Potential Flow: A Novel Probabilistic Framework for Energy-Based Generative Modelling

## Abstract

Abstract Energy based models (EBMs) are appealing for their generality and simplicity in data likelihood modeling, but have conventionally been difficult to train due to the unstable and time-consuming implicit MCMC sampling during contrastive divergence training.
In this paper, we present a novel energy-based generative framework, Variational Potential Flow ( VAPO ), that entirely dispenses with implicit MCMC sampling and does not rely on complementary latent models or cooperative training. The VAPO framework aims to learn a potential energy function whose gradient (flow) guides the prior samples, so that their density evolution closely follows an approximate data likelihood homotopy. An energy loss function is then formulated to minimize the Kullback-Leibler divergence between density evolution of the flow-driven prior and the data likelihood homotopy.
Images can be generated after training the potential energy, by initializing the samples from Gaussian prior and solving the ODE governing the potential flow on a fixed time interval using generic ODE solvers. Experiment results show that the proposed VAPO framework is capable of generating realistic images on various image datasets. In particular, our proposed framework achieves competitive FID scores for unconditional image generation on the CIFAR-10 and CelebA datasets.

## 1 Introduction

In recent years, deep generative modeling has garnered significant attention for unsupervised learning of complex, high-dimensional data distributions [].
In particular, probabilistic generative models such as variational autoencoders [], normalizing flows [], score-matching or diffusion models [, , ]. Poisson flow [, ], and energy-based models (EBMs) [, ] aim to maximize the likelihood (probability density) underlying the data. By design, these probabilistic frameworks enhance training stability, accelerate model convergence, and reduce mode collapse compared to generative adversarial networks [], albeit at the cost of a slow sampling procedure and poor model scalability [].
Among these frameworks, EBMs have emerged
as a flexible and expressive class of probabilistic generative models [, , , , , , , , ].
EBMs model high-dimensional data space with a network-parameterized energy potential function that assigns data regions with energy that is directly (or inversely) proportional to the unnormalized data likelihood []. This provides a natural interpretation of the network model in the form of an energy landscape, thereby endowing EBMs with inherent interpretability.

Deep EBMs are particularly appealing since they impose no restrictions on the network architecture, potentially resulting in high expressiveness []. Moreover, they are more robust and generalize well to out-of-distribution samples [, ] as regions with high probability under the model but low probability under the data distribution are explicitly penalized during training. Additionally, EBMs, which trace back to Boltzmann machines [], have strong ties to physics models and can thus borrow insights and techniques from statistical physics for their development and analysis [].
On these grounds, EBMs have been applied across a diverse array of applications apart from image modelling, including text generation [, ], point cloud synthesis [], scene graph generation [], anomaly detection [, ], earth observation [], robot learning [, ], trajectory prediction [, ], and molecular design [, ].

Despite a number of desirable properties, deep EBMs require implicit Langevin Markov Chain Monte Carlo (MCMC) sampling during the contrastive divergence training. MCMC sampling in a high-dimensional setting, however, has shown to be challenging due to poor mode mixing and excessively long mixing time [, , , , , ]. As result, energy potential functions learned with non-convergent MCMC do not have valid steady-states, in the sense that MCMC samples can differ greatly from data samples [].
Current deep EBMs are thus plagued by high variance training and high computational complexity due to MCMC sampling.
In view of this, recent works have explored learning complementary latent model to amortize away the challenging MCMC sampling [, , , , ], or cooperative learning where model-generated samples serve as initial points for subsequent MCMC revision in the latent space [, ].
While such approaches alleviate the burden of MCMC sampling, it comes at the expense of the inherent flexibility and composability of EBMs [].
Moreover, co-optimizing
multiple models
adds complexity [, ] to the implementation of these approaches.

In this paper, we introduce Variational Potential Flow (VAPO), a novel energy-based generative framework that eliminates the need for implicit MCMC sampling and complementary models. At the core of VAPO lies the construction of a homotopy (smooth path) that bridges the prior distribution with the data likelihood. Subsequently, a potential flow with model-parameterized potential energy function is designed to guide the evolution of prior sample densities along this approximate data likelihood homotopy. Applying a variational approach to this path-matching strategy ultimately yields a probabilistic Poisson’s equation, where the weak solution corresponds to minimizing the energy loss function of our proposed VAPO.

Our contributions are summarized as follows:

- •
We introduce VAPO, a novel energy-based generative framework that entirely dispenses with the unstable and inefficient implicit MCMC sampling. Our proposed framework learns a potential energy function whose gradient (flow) guides the prior samples, ensuring that their density evolution path closely follows the approximate data likelihood homotopy.
- •
We derive an energy loss function for VAPO by constructing a variational formulation of the intractable homotopy path-matching problem. Solving this energy loss objective is equivalent to minimizing the Kullback-Leibler divergence between density evolution of the flow-driven prior and the approximate data likelihood homotopy.
- •
To assess the effectiveness of our proposed VAPO for image generation, we conduct experiments on the CIFAR-10 and CelebA datasets and benchmark the performances against state-of-the-art generative models. Our proposed framework achieves competitive FID scores of 0.0 and 0.0 for unconditional image generation on CIFAR-10 and CelebA, respectively.

## 2 Background and Related Works

In this section, we provide an overview of EBMs, particle flow, and the deep Ritz method, collectively forming the cornerstone of our proposed VAPO framework.

### 2.1 Energy-Based Models (EBMs)

Denote $\bar{x}\in\Omega\subseteq\mathbb{R}^{n}$ as the training data, EBMs approximate the data likelihood ${p}_{\text{data}}(\bar{x})$ via defining a Boltzmann distribution, as follows:

$$ $\displaystyle\begin{split}{p}_{\theta}(x)=\frac{e^{\Phi_{\theta}(x)}}{\int_{ \Omega}e^{\Phi_{\theta}(x)}\;dx}\end{split}$ (1) $$

where $\Phi_{\theta}$ is an energy function modelled by deep neural networks.
Given that the denominator of (1), i.e., the partition function, is analytically intractable for high-dimensional data, EBMs perform the maximum likelihood estimation (MLE) by minimizing the negative log likelihood loss $\mathcal{L}_{\mathrm{MLE}}(\theta)=\operatorname{\mathbb{E}}_{{p}_{\text{data}
}(\bar{x})}[\log{p}_{\theta}(\bar{x})]$ and approximate its gradient via the contrastive divergence []:

$$ $\displaystyle\begin{split}\nabla_{\theta}\mathcal{L}_{\mathrm{MLE}}= \operatorname{\mathbb{E}}_{{p}_{\text{data}}(\bar{x})}\!\big{[}\nabla_{\theta} \Phi_{\theta}(\bar{x})\big{]}-\operatorname{\mathbb{E}}_{{p}_{\theta}(x)}\! \big{[}\nabla_{\theta}\Phi_{\theta}(x)\big{]}\end{split}$ (2) $$

However, EBMs are computationally intensive due to the implicit MCMC generating procedure, required for generating negative samples $x\in\Omega\sim{p}_{\theta}(x)$ for gradient computation (2) during training.

### 2.2 Particle Flow

Particle flow, initially introduced by the series of papers [], is a class of nonlinear Bayesian filtering (sequential inference) methods that aim to approximate the posterior distribution $p(x_{t}|\bar{x}_{0:t})$ of the state of system given the observations. While particle flow methods are closely related to normalizing flows [] and neural ordinary differential equations [], these latter frameworks do not explicitly accommodate a Bayes update.
In particular, particle flow performs the Bayes update $p(x_{t}|\bar{x}_{0:t})\propto p(x_{t}|\bar{x}_{0:t-1})\,p(\bar{x}_{t}|x_{t},
\bar{x}_{0:t-1})$ by subjecting prior samples $x_{t}\sim p(x_{t}|\bar{x}_{0:t-1})$ to a series of infinitesimal transformations through the ordinary differential equation (ODE) $\frac{dx}{d\tau}=v(x,\tau)$ parameterized by a flow velocity (field) function $v(x,\tau)$, in a pseudo-time interval $\tau\in[0,1]$ in between sampling time steps. The flow velocity is designed such that the driven Kolmogorov forward path evolution (Fokker–Planck dynamics, see (13)) of the sample particles, coincides with a data log-homotopy (smooth path) that inherently perform the Bayes update.
Despite its efficacy in time-series inference [, , ] and resilience to the curse of dimensionality [], particle flow has yet to be explored in generative modelling for high-dimensional data.

### 2.3 Deep Ritz Method

The deep Ritz method is a deep learning-based variational numerical approach, originally proposed in [], for solving scalar elliptic partial differential equations (PDEs) in high dimensions.
Consider the following Poisson’s equation, fundamental to many physical models:

$$ $\displaystyle\begin{split}\Delta u(x)=f(x),\quad x\in\Omega\end{split}$ (3) $$

subject to boundary condition

$$ $\displaystyle\begin{split}u(x)=0,\quad x\in\partial\Omega\end{split}$ (4) $$

where $\Delta$ is the Laplace operator, and $\partial\Omega$ denotes the boundary of $\Omega$. For a Sobolev function $u\in\mathcal{H}^{1}_{0}(\Omega)$ (see Proposition 2 for definition) and square-integrable $f\in L^{2}(\Omega)$, the variational principle ensures that a weak solution $u^{*}$ of the Euler-Lagrange boundary value equation (3)-(4) is equivalent to the variational problem of minimizing the Dirichlet energy [], as follows:

$$ $\displaystyle\begin{split}u^{*}=\arg\underset{v}{\min}\int_{\Omega}\bigg{(}\, \frac{1}{2}\|\nabla v(x)\|^{2}-f(x)v(x)\bigg{)}\;dx\end{split}$ (5) $$

where $\nabla$ denotes the Del operator (gradient). In particular, the deep Ritz method parameterizes the trial energy function $v$ using neural networks, and performs the optimization (5) via stochastic gradient descent. Due to its versatility and effectiveness in handling high-dimensional PDE systems, the deep Ritz method is predominantly applied for finite element analysis []. In [], the deep Ritz method is used to solve the probabilistic Poisson’s equation resulting from the feedback particle filter [].
Nonetheless, the method has not been explored for generative modelling.

## 3 Variational Energy-Based Potential Flow

In this section, we introduce a novel generative modelling framework, Variational Energy-Based Potential Flow (VAPO), drawing inspiration from both particle flow and the calculus of variations. First, we establish a homotopy that transforms a prior to the data likelihood and derive the evolution of the prior in time. Then, we design an energy-generated potential flow and a weighted Poisson’s equation that aligns the evolving density distribution of transported particles with the homotopy-driven prior. Subsequently, we formulate a variational loss function where its optimization with respect to the flow-generating potential energy is equivalent to solving the Poisson’s equation. Finally, we describe the model architecture that is used to parameterize the potential energy function and the backward ODE integration for generative sampling.

### 3.1 Bridging Prior and Data Likelihood: Log-Homotopy Transformation

Let $\bar{x}\in\Omega$ denote the training data, ${p}_{\text{data}}(\bar{x})$ be the data likelihood, $x\in\Omega$ denote the approximate data samples.
To achieve generative modelling, our objective is to closely approximate the training data $\bar{x}$ with the data samples $x$. On this account, we define a conditional data likelihood $p(\bar{x}|x)=\mathcal{N}(\bar{x};x,\Pi)$ with isotropic Gaussian noise with covariance $\Pi=\text{diag}(\sigma^{2})$ and standard deviation $\sigma\in\Omega$. This is equivalent to considering a state space model $x=\bar{x}+\nu$, where $\nu\in\Omega\sim\mathcal{N}(\nu;0,\Pi)$. Here, we set a small $\sigma$ so that $x$ closely resembles the training data $\bar{x}$.

Subsequently, consider a conditonal (data-conditioned) density function $\rho:\Omega^{2}\times[0,1]\rightarrow\mathbb{R}$, as follows:

$$ $\displaystyle\begin{split}\rho(x;\bar{x},t)=\frac{e^{f(x;\bar{x},t)}}{\int_{ \Omega}e^{f(x;\bar{x},t)}\;dx}\end{split}$ (6) $$

where $f:\Omega^{2}\times[0,1]\rightarrow\mathbb{R}$ is a log-linear function:

$$ $\displaystyle\begin{split}f(x;\bar{x},t)=\log q(x)+t\,\log p(\bar{x}|x)\end{split}$ (7) $$

parameterized by the auxiliary time variable $t\in[0,1]$, and we let $q(x)=\mathcal{N}(x;0,\Lambda)$ be a isotropic Gaussian prior density with covariance $\Lambda=\text{diag}(\omega^{2})$ and standard deviation $\omega\in\Omega$. Here, $\text{diag}(\cdot)$ denotes the diagonal function.
By construction, we have
$\rho(x;\bar{x},0)=q(x)$ at $t=0$, and $\rho(x;\bar{x},1)=p(x|\bar{x})$ at $t=1$ since we have

$$ $\displaystyle\begin{split}\rho(x;\bar{x},1)&=\frac{e^{f(x;\bar{x},1)}}{\int_{ \Omega}e^{f(x;\bar{x},1)}\;dx}=\frac{p(\bar{x}|x)\,q(x)}{{p}_{\text{data}}( \bar{x})}=\frac{p(\bar{x},x)}{{p}_{\text{data}}(\bar{x})}=p(x|\bar{x})\end{split}$ (8) $$

where we have used the fact that ${p}_{\text{data}}(\bar{x})=\int_{\Omega}e^{f(x;\bar{x},1)}\,dx=\int_{\Omega}p(
x,\bar{x})\,dx$.
Therefore, the conditional density function $\rho(x;\bar{x},t)$ here (6) essentially represents a density homotopy between the prior $q(x)$ and the posterior $p(x|\bar{x})$.

In particular, the density function $\rho(x;\bar{x},t)$ also defines a conditional (data-conditioned) homotopy between the prior $q(x)$ and the exact posterior $p(x|\bar{x})$, the latter of which gives a maximum a posteriori (Bayesian) estimate of the approximate data samples after observing true training data.

To obtain an estimate of the intractable data likelihood for generative sampling, we then consider a (approximate) data likelihood homotopy
$\bar{\rho}:\Omega\times[0,1]\rightarrow\mathbb{R}$ as follows:

$$ $\displaystyle\begin{split}\bar{\rho}(x;t)&=\int_{\Omega}{p}_{\text{data}}(\bar {x})\,\rho(x;\bar{x},t)\;d\bar{x}\end{split}$ (9) $$

Considering this, it remains that $\bar{\rho}(x;0)=q(x)$ at $t=0$. Furthermore, given that we have
$\bar{\rho}(x;1)=\int_{\Omega}{p}_{\text{data}}(\bar{x})\,p(x|\bar{x})\;d\bar{x
}=\bar{p}(x)$ at $t=1$,
the data likelihood homotopy $\bar{\rho}(x;t)$ here inherently performs a kernel density approximation of the true data likelihood, using the normalized kernel $p(x|\bar{x})$ obtained from the conditional homotopy $\rho(x;\bar{x},1)$ at $t=1$.
Therefore, the approximate data likelihood $\bar{p}(x)$ acts as a continuous interpolation of the data likelihood ${p}_{\text{data}}(x)$, represented by Dirac delta function $\delta(x-\bar{x})$ centered on the discrete training data $\bar{x}$.

Nevertheless, the conditional homotopy (8) is intractable due to the normalizing constant in the denominator. This intractability rules out a close-form solution of the data likelihood homotopy (9), thus it is not possible to sample directly from the data likelihood estimate. Taking this into account, we introduce the potential flow method
in the following section, where we model the evolution of the prior samples (particles) instead, such that their distribution adheres to the data likelihood homotopy.

### 3.2 Modelling Potential Flow in a Homotopy Landscape

Our aim is to model the flow of the prior particles in order for their distribution to follow the data likelihood homotopy and converge to the data likelihood. To accomplish this, we first derive the evolution of the latent prior density with respect to time in the following proposition.

###### Proposition 1 .

Consider the data likelihood homotopy $\bar{\rho}(x;t)$ in (9) with Gaussian conditional data likelihood $p(\bar{x}|x)=\mathcal{N}(\bar{x};x,\Pi)$. Then, its evolution in time $t\in[0,1]$ is given by the following PDE:

$$ $\displaystyle\begin{split}\frac{\partial\bar{\rho}(x;t)}{\partial t}=-\,\frac{ 1}{2}\,\operatorname{\mathbb{E}}_{{p}_{\text{data}}(\bar{x})}\Big{[}\rho(x; \bar{x},t)\,\big{(}\gamma(x,\bar{x})-\bar{\gamma}(x,\bar{x})\big{)}\Big{]}\end {split}$ (10) $$

where

$$ $\displaystyle\begin{split}\gamma(x,\bar{x})=(x-\bar{x})^{T}\,\Pi^{-1}\,(x-\bar {x})\end{split}$ (11) $$

is the innovation term in the conditional data likelihood, and $\bar{\gamma}(x,\bar{x})=\operatorname{\mathbb{E}}_{\rho(x;\bar{x},t)}[\gamma(x
,\bar{x})]$ denotes the expectation of the innovation with respect to the conditional homotopy and on the latent variables.

###### Proof.

Refer to Appendix A.1.
∎

Figure: Figure 1: A planar visualization of the potential-generated field (represented by coloured arrows) that transports the prior particles towards the approximate data likelihood (represented by the blue contour).
Refer to caption: extracted/5745897/VectorField_2.png

Our proposed potential flow method involves subjecting the latent prior samples to a potential-generated velocity field, such that the flow trajectories of these sample particles $x(t)$ within the interval $t\in[0,1]$ are governed by the following ordinary differential equation (ODE):

$$ $\displaystyle\frac{d(x(t))}{dt}=\nabla\Phi(x(t))$ (12) $$

where $\Phi:\Omega\rightarrow\mathbb{R}$ is a scalar potential energy function. Therefore, $\nabla\Phi\in\Omega$ is the velocity vector field generated by the potential energy, and $\nabla$ denotes the Del operator (gradient) with respect to the data samples $x$. The scalar potential (12) is a result of the Helmholtz decomposition of a vector field $V(x(t))$ and disregarding the solenoidal (rotational) component.
Henceforth, the time variable of $x(t)$ is implicitly assumed, thus we omit it for simplicity.

Considering a potential flow of the form (12), a direct consequence is that the (approximate likelihood) density ${\rho}^{\Phi}$ of the flow-driven prior samples evolves according to a Fokker–Planck (Kolmogorov forward) equation as follows:

$$ $\displaystyle\begin{split}\frac{\partial{\rho}^{\Phi}(x;t)}{\partial t}=-\, \nabla\cdot\Big{(}{\rho}^{\Phi}(x;t)\,\nabla\Phi(x)\Big{)}\end{split}$ (13) $$

where $\nabla\cdot$ denotes the divergence operator.
In particular, the Fokker–Planck equation (13) exemplifies a form of continuity (transport) equation commonly used for modelling fluid advection. In this analogy, ${\rho}^{\Phi}$ corresponds to fluid density, $\Phi$ represents a field-driving potential energy, and its gradient $\nabla\Phi$ acts as the resulting conservative (irrotational) velocity field.

The goal of our proposed framework is to model the potential energy function in the potential flow (12), such that the progression of the prior density subject to potential flow emulates the evolution of the data likelihood homotopy.
In particular, we seek to solve the problem of minimizing the Kullback-Leibler divergence (KLD) in the following proposition.

###### Proposition 2 .

Consider a potential flow of the form (12) and given that $\Phi\in\mathcal{H}^{1}_{0}(\Omega,\bar{\rho})$, where $\mathcal{H}^{n}_{0}$ denotes the (Sobolev) space of $n$-times differentiable functions that are compactly supported, and square-integrable with respect to data likelihood homotopy $\bar{\rho}(x;t)$.

Then, the problem of solving for the optimal potential energy function $\Phi(x)$ that satisfies the following probabilistic (density-weighted) Poisson’s equation:

$$ $\displaystyle\begin{split}&\nabla\cdot\Big{(}\bar{\rho}(x;t)\,\nabla\Phi(x) \Big{)}=\frac{1}{2}\,\operatorname{\mathbb{E}}_{{p}_{\text{data}}(\bar{x})} \Big{[}\rho(x;\bar{x},t)\,\big{(}\gamma(x,\bar{x})-\bar{\gamma}(x,\bar{x})\big {)}\Big{]}\end{split}$ (14) $$

is equivalent to minimizing the KLD $\operatorname{\mathcal{D}_{\mathrm{KL}}}\big{[}{\rho}^{\Phi}(x;t)\,\|\,\bar{
\rho}(x;t)\big{]}$ between the flow-driven prior ${\rho}^{\Phi}(x;t)$ and the data likelihood homotopy $\bar{\rho}(x;t)$ at time $t$.

###### Proof.

Refer to Appendix A.2.
∎

In hindsight, the left-hand side of the probabilistic Poisson’s equation (LABEL:eq:unmarginalized_PDE_equation) resembles the evolution of the flow-driven prior given by the Fokker-Plank equation (13). In addition, the right-hand side resembles the evolution of data likelihood homotopy given by PDE (10), with the conditional homotopy $\rho(x;\bar{x},t)$ replaced by flow-driven prior ${\rho}^{\Phi}(x;t)$. Therefore, the probabilistic Poisson’s equation is an attempt to solve the approximation
$\frac{\partial{\rho}^{\Phi}(x;t)}{\partial t}\equiv\frac{\partial\bar{\rho}(x;
t)}{\partial t}$.

Nevertheless, explicitly solving the probabilistic Poisson’s equation (LABEL:eq:unmarginalized_PDE_equation) is challenging in a high-dimensional setting. Numerical methods that approximate the solution often do not scale well with the data dimension. For example, the Galerkin approximation requires a selection of the basis functions, which becomes non-trivial when the dimensionality is high []. The diffusion map-based algorithm, on the other hand,
requires a large number of particles, which grows exponentially with respect to the dimensionality, in order to achieve error convergence []. Taking this into consideration, we propose an energy loss function in the following section, where we cast the Poisson’s equation as a variational problem compatible with stochastic gradient descent.

###### Proposition 1 .

###### Proof.

###### Proposition 2 .

###### Proof.

### 3.3 Variational Energy Loss Function Formulation: Deep Ritz Approach

In this section, we introduce an energy method which presents a variational formulation of the probabilistic Poisson’s equation.
Given that the aim is to minimize the divergence between the data likelihood homotopy and the flow-driven prior and directly solving the probabilistic Poisson’s equation is difficult, we first consider a weak formulation of (LABEL:eq:unmarginalized_PDE_equation) as follows:

$$ $\displaystyle\begin{split}\int_{\Omega}\bigg{(}&\,\frac{1}{2}\,\operatorname{ \mathbb{E}}_{{p}_{\text{data}}(\bar{x})}\Big{[}\rho(x;\bar{x},t)\,\big{(} \gamma(x,\bar{x})-\bar{\gamma}(x,\bar{x})\big{)}\Big{]}-\nabla\cdot\Big{(}\bar {\rho}(x;t)\,\nabla\Phi(x)\Big{)}\bigg{)}\,\Psi(x)\;dx\,=\,0\end{split}$ (15) $$

where the equation must hold for all differentiable trial functions $\Psi$. In the following proposition, we introduce an energy loss objective that is equivalent to solving this weak formulation of the probabilistic Poisson’s equation.

###### Proposition 3 .

The variational problem of minimizing the following loss function:

$$ $\displaystyle\begin{split}\mathcal{L}(\Phi;t)=&\;\frac{1}{2}\,\operatorname{ \mathrm{Cov}}_{\rho(x;\bar{x},t)\,{p}_{\text{data}}(\bar{x})}\big{[}\Phi(x), \gamma(x,\bar{x})\big{]}+\frac{1}{2}\,\operatorname{\mathbb{E}}_{\bar{\rho}(x; t)}\Big{[}\big{\|}\nabla\Phi(x)\big{\|}^{2}\Big{]}\end{split}$ (16) $$

with respect to the potential energy $\Phi$,
is equivalent to solving the weak formulation (15) of the probabilistic Poisson’s equation (LABEL:eq:unmarginalized_PDE_equation).
Here, $\|\cdot\|$ denotes the Euclidean norm, and $\operatorname{\mathrm{Cov}}$ denotes the covariance.

Furthermore, the variational problem (16) has a unique solution if for all energy functions $\Phi\in\mathcal{H}^{1}_{0}(\Omega;\bar{\rho})$, the data likelihood homotopy $\bar{\rho}$ satisfy the Poincaré inequality:

$$ $\displaystyle\begin{split}\operatorname{\mathbb{E}}_{\bar{\rho}(x;t)}\Big{[} \big{\|}\nabla\Phi(x)\big{\|}^{2}\Big{]}\geq\lambda\,\operatorname{\mathbb{E}} _{\bar{\rho}(x;t)}\Big{[}\big{\|}\Phi(x)\big{\|}^{2}\Big{]}\end{split}$ (17) $$

for some positive scalar constant $\lambda>0$ (spectral gap).

###### Proof.

Refer to Appendix A.3.
∎

In sum, leveraging Propositions 2 and 3, we reformulate the intractable task of minimizing the KLD between flow-driven prior and data likelihood homotopy equivalently as a variational problem with energy loss function (16). By optimizing the potential energy function $\Phi$ with respect to the energy loss and transporting the prior samples through the potential flow ODE (12), the prior particles follow a trajectory that accurately approximates the data likelihood homotopy. In doing so, the potential flow $\nabla\Phi$ drives the prior samples to posterior regions densely populated with data, thus enabling us to perform generative modelling.

The minimum covariance objective in (16) plays an important role by ensuring that the normalized innovation is inversely proportional to the potential energy. As a result, the potential-generated velocity field $\nabla\Phi$ consistently points in the direction of greatest potential ascent, thereby driving the flow of prior particles towards high likelihood regions of the true posterior, as illustrated in Figure 1. In other words, the potential energy is conjugate to the approximate data likelihood $\bar{p}(\bar{x})$, analogous to Hamiltonian fluid mechanics []. It is worth noticing that instead of being an ad hoc addition, the L2 regularization term $\operatorname{\mathbb{E}}_{\bar{\rho}(x;t)}[\|\nabla\Phi(x)\|^{2}]$ on the velocity field in (16) arises, from first-principle derivation, as a direct consequence of considering the data likelihood homotopy.

Given that the aim is to solve the probabilistic Poisson’s equation (LABEL:eq:unmarginalized_PDE_equation) for all $t$, we include an auxiliary time integral to the energy loss function (16) as follows:

$$ $\displaystyle\begin{split}\mathcal{L}^{\mathrm{VAPO}}(\theta)=\int_{\mathbb{R} }\,\mathcal{L}(\theta;t)\;dt=\operatorname{\mathbb{E}}_{\,\mathcal{U}(t;0,1)} \big{[}\mathcal{L}(\theta;t)\big{]}\end{split}$ (18) $$

where we have applied Monte Carlo integration, and $\mathcal{U}(a,b)$ denotes the uniform distribution over interval $[a,b]$.
In addition, the data likelihood homotopy may not satisfy the Poincaré inequality (17). Hence, we include the right-hand side of the inequality to the loss function (16) to enforce uniqueness of its minimizer. This addtional L2 loss also regularize the energy function, preventing its values from exploding. The spectral gap constant $\lambda$ is left as a training hyperparameter.

In addition, the energy loss (16) requires us to sample from the conditional and data likelihood density homotopies. By design, both the prior $q(x)=\mathcal{N}(x;0,\Lambda)$ and the conditional data likelihood $p(\bar{x}|x)=\mathcal{N}(\bar{x};x,\Pi)$ are assumed to be Gaussian. As a consequence, the Bayes update (6) results in a Gaussian density $\rho(x;\bar{x},t)=\mathcal{N}\big{(}x;\mu(\bar{x},t),\Sigma(\bar{x},t)\big{)}$, from which the time-varying mean and covariance can be derived using the Bayes’ theorem [], as follows:

$$ $\displaystyle\mu(\bar{x},t)=t\,\Sigma(t)\,\Pi^{-1}\,\bar{x},\qquad\Sigma(t)= \big{(}\Lambda^{-1}\,+\,t\,\Pi^{-1}\big{)}^{-1}$ (19) $$

Therefore, to sample from $\rho(x;\bar{x},t)$ or $\rho(x;t)$, we first sample data $x$ from ${p}_{\text{data}}(\bar{x})$ and compute the mean and covariance according to (19). Then, we can generate samples of the approximate data $x$ using the reparameterization trick $x=\mu(\bar{x},t)+\sqrt{\Sigma}(t)\,\epsilon$, where $\epsilon\sim\mathcal{N}(\epsilon;0,I)$ and $\sqrt{\Sigma}$ is the square root decomposition of $\Sigma$, i.e., $\Sigma=\sqrt{\Sigma}\sqrt{\Sigma}^{T}$. A detailed derivation of (19) is provided in Appendix A.4.

Nevertheless, parameterizing the conditional homotopy using mean and covariance (19) causes it to converge too quickly to the posterior $\rho(x;\bar{x},1)=p(x|\bar{x})$. As a consequence, most samples are closely clustered around the observed data. To mitigate this issue, a strategy is to slow down its convergence by reparameterizing it with $t+\varepsilon=e^{\tau}$, where $\tau\in[\ln\varepsilon,\ln(1+\varepsilon)]$. This time reparameterization compels $t+\varepsilon$ to follow a log-uniform (reciprocal) distribution $\mathcal{R}(t+\varepsilon;\varepsilon,1+\varepsilon)$ defined over the interval $[\varepsilon,1+\varepsilon]$. Here, the hyperparameter $\varepsilon$ is a small positive constant that determines the sharpness of the log-uniform density, and the rate at which its tail decays to zero.

Incorporating all of the above considerations, the final energy loss function becomes:

$$ $\displaystyle\begin{split}\mathcal{L}_{\mathrm{VAPO}}(\theta)&=\frac{1}{2}\, \operatorname{\mathbb{E}}_{\mathcal{R}(t+\varepsilon;\,\varepsilon,1+ \varepsilon)}\big{[}\mathcal{L}(\theta;t)\big{]}\end{split}$ (20) $$

where

$$ $\displaystyle\begin{split}\mathcal{L}(\theta;t)=&\;\operatorname{\mathrm{Cov}} _{\rho(x;\bar{x},t)\,{p}_{\text{data}}(\bar{x})}\big{[}\Phi_{\theta}(x),\gamma (x,\bar{x})\big{]}\\ &+\operatorname{\mathbb{E}}_{\bar{\rho}(x;t)}\Big{[}\|\nabla\Phi_{\theta}(x) \big{\|}^{2}\Big{]}+\lambda\operatorname{\mathbb{E}}_{\bar{\rho}(x;t)}\Big{[} \big{\|}\Phi_{\theta}(x)\big{\|}^{2}\Big{]}\end{split}$ (21) $$

The algorithm for training VAPO is shown in Algorithm 1.

Figure: Algorithm 1 VAPO Training

###### Proposition 3 .

###### Proof.

### 3.4 Energy Parameterization and ODE Sampling

To implement stochastic gradient descent on top of the energy loss function (18), we adopt the deep Ritz approach and in particular, we model the potential energy function $\Phi_{\theta}$ as deep neural networks with parameters $\theta$. Here, we restrict our model architecture to convolutional and fully-connected layers, which are shown to satisfy the universal approximation property within weighted Sobolev spaces, i.e., the neural network model densely approximates functions in $\mathcal{H}^{1}_{0}(\Omega,\bar{\rho})$ and enables model convergence.

**Table 1: Comparison of FID scores on unconditional CIFAR-10 image generation. FID baselines are obtained from [].**
| Models | FID $\downarrow$ | Models | FID $\downarrow$ |
| --- | --- | --- | --- |
| EBM-based methods | Other likelihood-based methods |  |  |
| EBM-SR [] | 44.5 | VAE [] | 78.4 |
| JEM [] | 38.4 | PixelCNN [] | 65.9 |
| EBM-IG [] | 38.2 | PixelIQN [] | 49.5 |
| EBM-FCE [] | 37.3 | ResidualFlow [] | 47.4 |
| CoopVAEBM [] | 36.2 | Glow [] | 46.0 |
| CoopNets [] | 33.6 | DC-VAE [] | 17.9 |
| Divergence Triangle [] | 30.1 | GAN-based methods |  |
| VERA [] | 27.5 | WGAN-GP [] | 36.4 |
| EBM-CD [] | 25.1 | SN-GAN [] | 21.7 |
| GEBM [] | 19.3 | SNGAN-DDLS [] | 15.4 |
| HAT-EBM [] | 19.3 | BigGAN [] | 14.8 |
| CF-EBM [] | 16.7 | Score-based and Diffusion methods |  |
| CoopFlow [] | 15.8 | NCSN [] | 25.3 |
| CLEL-base [] | 15.3 | NCSN-v2 [] | 10.9 |
| VAEBM [] | 12.2 | DDPM Distil. [] | 9.36 |
| DRL [] | 9.58 | DDPM [] | 3.17 |
| VAPO (Ours) | 16.6 | NCSN++[] | 2.20 |

**Table 2: Comparison of FID scores on unconditional CelebA $64\times 64$. FID baselines obtained from [].**
| Models | FID $\downarrow$ |
| --- | --- |
| NCSN [] | 25.3 |
| NCSN-v2 [] | 10.2 |
| EBM-Triangle [] | 24.7 |
| EBM-SR [] | 23.0 |
| Divergence Triangle [] | 18.2 |
| CoopNets [] | 16.7 |
| VAPO (Ours) | 14.5 |

After modelling the potential energy $\Phi_{\theta}$, the gradient $\nabla\Phi_{\theta}$ can be used to generate approximate data samples from the potential flow ODE (12). Given that the ODE is defined on the interval $[0,1]$ by construction, its starting and terminal time is predetermined and therefore known, in contrast to most flow-based generative frameworks. On top of that, the potential flow ODE is compatible with general-purpose ODE solvers, such as the explicit and implicit Runge-Kutta methods of different orders and the forward and backward Euler methods, which can readily be employed for sampling.

## 4 Experiments

In this section, we show that VAPO is an effective generative model for images. In particular, Section 4.1 demonstrates that VAPO is capable of generating realistic unconditional images on the well-known CIFAR-10 and CelebA datasets. Section 4.2 demonstrates that VAPO is capable of performing smooth interpolation between two generated samples.
Implementation details, including model architecture and training, numerical ODE solver, datasets and FID evaluation are provided in Appendix B.
Apart from that, we also show that VAPO exhibits extensive mode coverage and robustness to anomalous data, as well as generalizing well to unseen data without over-fitting. Specifically, Appendix C.1 evaluates model over-fitting and generalization based on the energy histogram of CIFAR-10 train and test sets and the nearest neighbors of generated samples. Appendix C.2 examines robustness to anomalous data by assessing its performance on out-of-distribution (OOD) detection
on various image datasets.

Figure: Figure 2: Generated samples on unconditional CIFAR-10 $32\times 32$ (left) and CelebA $64\times 64$ (right).
Refer to caption: extracted/5745897/uncurated.png

### 4.1 Unconditional Image Generation

Figure 2 shows the uncurated and unconditional image samples generated from the learned energy model on the CIFAR-10 and CelebA $64\times 64$ datasets. More generated samples are provided in Appendix D. The samples are of decent quality and resemble the original datasets despite not having the highest fidelity as achieved by state-of-the-art models. Tables 1 and 2 summarize the quantitative evaluations of our proposed VAPO model in terms of FID [] scores on the CIFAR-10 and CelebA datasets. On CIFAR-10, VAPO achieves a competitive FID that is better than the majority of existing EBM-based generative models. Having dispensed with the implicit MCMC sampling, VAPO still outperforms most of the EBM approaches without relying on complementary latent models or cooperative training.
On CelebA, VAPO obtains an FID that outperforms some existing EBMs but falls short compared to [] and state-of-the-art models.

### 4.2 Image Interpolation

Figure 3 shows the interpolation results between pairs of generated CelebA samples, where it demonstrates that VAPO is capable of smooth and semantically meaningful image interpolation. To perform interpolation for two samples $x_{1}(1)$ and $x_{2}(1)$, we construct a spherical interpolation between the initial Gaussian noise $x_{1}(0)$ and $x_{2}(0)$, and subject them to sampling over the potential flow ODE. More interpolation results on CIFAR-10 and CelebA are provided in Appendix D.

Figure: Figure 3: Interpolation results between the leftmost and rightmost generated CelebA $64\times 64$ samples.
Refer to caption: extracted/5745897/celeba_interp_small.png

## 5 Limitations and Future Work

We propose VAPO, a novel energy-based generative modelling framework without the need for expensive and unstable MCMC runs amidst training. Despite the improvement over the majority of existing EBMs, there is still a large performance gap between VAPO and the state-of-the-art score-based (or diffusion) and Poisson flow models [, , ]. To close this gap, diffusion recovery likelihood [, ], which is shown to be more tractable than marginal likelihood, can be incorporated into the VAPO framework for a more controlled diffusion-guided energy optimization. The dimensionality augmentation technique of [, ] can also be integrated given that fundamentally, both Poisson flow and VAPO aim to model potential field governed by a Poisson’s equation. On top of that, the scalability of VAPO to higher resolution images and its generalizability to other data modalities have yet to be validated. In addition, the current VAPO framework does not allow for class-conditional generation. Moreover, the training of VAPO requires a large number of iterations to converge and thus warrants improvement. These important aspects are earmarked for future extensions of our work.

## Appendix A Proofs and Derivations

### A.1 Proof of Proposition 1

###### Proof.

Differentiating the conditional homotopy $\rho(x;\bar{x},t)$ in (6) with respect to $t$, we have

$$ $\displaystyle\begin{split}&\frac{\partial\rho(x;\bar{x},t)}{\partial t}=\frac{ 1}{\int_{\Omega}e^{f(x;\bar{x},t)}\;dx}\,\frac{\partial[e^{f(x;\bar{x},t)}]}{ \partial t}-\frac{e^{f(x;\bar{x},t)}}{[\int_{\Omega}e^{f(x;\bar{x},t)}\;dx]^{2 }}\,\frac{\partial[\int_{\Omega}e^{f(x;\bar{x},t)}\;dx]}{\partial t}\\ &=\frac{1}{\int_{\Omega}e^{f(x;\bar{x},t)}\;dx}\,\frac{\partial[e^{f(x;\bar{x} ,t)}]}{\partial f}\,\frac{\partial f(x;\bar{x},t)}{\partial t}-\frac{e^{f(x; \bar{x},t)}}{[\int_{\Omega}e^{f(x;\bar{x},t)}\;dx]^{2}}\,\int_{\Omega}\frac{ \partial[e^{f(x;\bar{x},t)}]}{\partial f}\,\frac{\partial f(x;\bar{x},t)}{ \partial t}\;dx\\ &=\frac{e^{f(x;\bar{x},t)}}{\int_{\Omega}e^{f(x;\bar{x},t)}\;dx}\,\frac{ \partial f(x;\bar{x},t)}{\partial t}-\frac{e^{f(x;\bar{x},t)}}{\int_{\Omega}e^ {f(x;\bar{x},t)}\;dx}\,\int_{\Omega}\frac{e^{f(x;\bar{x},t)}}{\int_{\Omega}e^{ f(x;\bar{x},t)}\;dx}\,\frac{\partial f(x;\bar{x},t)}{\partial t}\;dx\\ &=\rho(x;\bar{x},t)\left(\frac{\partial f(x;\bar{x},t)}{\partial t}-\int_{ \Omega}\rho(x;\bar{x},t)\,\frac{\partial f(x;\bar{x},t)}{\partial t}\;dx\right )\\ &=-\,\frac{1}{2}\,\rho(x;\bar{x},t)\Bigg{(}(x-\bar{x})^{T}\,\Pi^{-1}\,(x-\bar{ x})-\int_{\Omega}\rho(x;\bar{x},t)\,(x-\bar{x})^{T}\,\Pi^{-1}\,(x-\bar{x})\;dx \Bigg{)}\end{split}$ (22) $$

where we used the quotient rule in the first equation, and chain rule in the second equation.

Let $\gamma(x,\bar{x})=(x-\bar{x})^{T}\,\Pi^{-1}\,(x-\bar{x})$ and using the fact that

$$ $\displaystyle\begin{split}\frac{\partial\bar{\rho}(x;t)}{\partial t}=\frac{ \partial\int_{\Omega}\rho(x;\bar{x},t)\,{p}_{\text{data}}(\bar{x})\;d\bar{x}}{ \partial t}=\int_{\Omega}\frac{\partial\rho(x;\bar{x},t)}{\partial t}\,{p}_{ \text{data}}(\bar{x})\;d\bar{x}\end{split}$ (23) $$

we can substitute (22) into (23) to get

$$ $\displaystyle\begin{split}\frac{\partial\bar{\rho}(x;t)}{\partial t}=-\,\frac{ 1}{2}\int_{\Omega}{p}_{\text{data}}(\bar{x})\,\rho(x;\bar{x},t)\,\bigg{(} \gamma(x,\bar{x})-\int_{\Omega}\rho(x;\bar{x},t)\,\gamma(x,\bar{x})\;dx\bigg{) }\;d\bar{x}\end{split}$ (24) $$

Given that both $\rho(x;\bar{x},t)$ and ${p}_{\text{data}}(x)$ are normalized (proper) density functions, writing (24) in terms of expectations
yields the PDE in (10).
∎

###### Proof.

### A.2 Proof of Proposition 2

Here, we used the Einstein tensor notation interchangeably with the conventional notation for vector dot product and matrix-vector multiplication in PDE.

Given that the context is clear, we write $g_{t}(\cdot)$ in place of time-varying functions $g(\cdot,t)$. For brevity, we will also omit the time index $t$, and write $g$ in place of $g_{t}(x)$.

###### Proof.

Applying the forward Euler method to the particle flow ODE (12) using step size $\Delta_{t}$, we obtain:

$$ $\displaystyle\begin{split}x_{t+\Delta_{t}}&=\alpha_{t}(x_{t})=x_{t}+\Delta_{t} \,u(x_{t})\end{split}$ (25) $$

where

$$ $\displaystyle\begin{split}u(x_{t})&=\nabla\Phi(x_{t})\end{split}$ (26) $$

where we denote $x_{t}$ as the discretizations random variables $x(t)$.

Assuming that the $\alpha_{t}:\Omega\rightarrow\Omega$ is a diffeomorphism (bijective function with differentiable inverse), the push-forward operator $\alpha_{t\#}:\mathbb{R}\rightarrow\mathbb{R}$ on density function ${\rho}^{\Phi}_{t}\mapsto{\rho}^{\Phi}_{t+\Delta_{t}}:=\alpha_{t\#}{\rho}^{\Phi
}_{t}$ is defined by:

$$ $\displaystyle\begin{split}\int_{\Omega}{\rho}^{\Phi}_{t+\Delta_{t}}(x)\,g(x)\; dx=\int_{\Omega}\alpha_{t\#}{\rho}^{\Phi}_{t}(x)\,g(x)\;dx=\int_{\Omega}{\rho} ^{\Phi}_{t}(x)\,g\big{(}\alpha_{t}(x)\big{)}\;dx\end{split}$ (27) $$

for any measurable function $g$.

Associated with the change-of-variables formula (27) is the following density transformation:

$$ $\displaystyle\begin{split}{\rho}^{\Phi}_{t+\Delta_{t}}\big{(}\alpha_{t}(x)\big {)}&=\frac{1}{|D\alpha_{t}|}\,{\rho}^{\Phi}_{t}(x)\end{split}$ (28) $$

where $|D\alpha_{t}|$ denotes the Jacobian determinant of $\alpha_{t}$.

From (10) and (24), we have

$$ $\displaystyle\begin{split}\frac{\partial\ln\bar{\rho}_{t}(x)}{\partial t}&= \frac{1}{\bar{\rho}_{t}(x)}\,\frac{\partial\bar{\rho}_{t}(x)}{\partial t}=-\, \frac{1}{\bar{\rho}_{t}(x)}\,\frac{1}{2}\,\operatorname{\mathbb{E}}_{{p}_{ \text{data}}(\bar{x})}\Big{[}\rho_{t}(x,\bar{x})\,\big{(}\gamma(x,\bar{x})- \bar{\gamma}(x,\bar{x})\big{)}\Big{]}\end{split}$ (29) $$

Applying the forward Euler method to (29), we obtain

$$ $\displaystyle\begin{split}\ln\bar{\rho}_{t+\Delta_{t}}(x)\geq\ln\bar{\rho}_{t} (x)-\frac{\Delta_{t}}{2}\,\frac{1}{\bar{\rho}_{t}(x)}\,\operatorname{\mathbb{E }}_{{p}_{\text{data}}(\bar{x})}\Big{[}\rho_{t}(x,\bar{x})\,\big{(}\gamma(x, \bar{x})-\bar{\gamma}(x,\bar{x})\big{)}\Big{]}\end{split}$ (30) $$

Applying the change-of-variables formula (27) and density transformation (28), then substituting (30) into the KLD $\operatorname{\mathcal{D}_{\mathrm{KL}}}\big{[}{\rho}^{\Phi}_{t+\Delta_{t}}\,
\|\,\bar{\rho}_{t+\Delta_{t}}\big{]}$ at time $t+\Delta_{t}$, we have

$$ $\displaystyle\begin{split}&\operatorname{\mathcal{D}_{\mathrm{KL}}}\big{[}{ \rho}^{\Phi}_{t+\Delta_{t}}(x)\,\|\,\bar{\rho}_{t+\Delta_{t}}(x)\big{]}=\int_{ \Omega}{\rho}^{\Phi}_{t}(x)\,\ln\Bigg{(}\frac{{\rho}^{\Phi}_{t+\Delta_{t}}\big {(}\alpha_{t}(x)\big{)}}{\bar{\rho}_{t+\Delta_{t}}\big{(}\alpha_{t}(x)\big{)}} \Bigg{)}\;dx\\ &=\int_{\Omega}{\rho}^{\Phi}_{t}(x)\,\bigg{(}\ln{\rho}^{\Phi}_{t}(x)-\ln|D \alpha_{t}|-\ln\bar{\rho}_{t}\big{(}\alpha_{t}(x)\big{)}\\ &\qquad\qquad\qquad+\frac{\Delta_{t}}{2}\,\frac{1}{\bar{\rho}_{t}\big{(}\alpha _{t}(x)\big{)}}\,\operatorname{\mathbb{E}}_{{p}_{\text{data}}(\bar{x})}\bigg{[ }\rho_{t}\big{(}\alpha_{t}(x),\bar{x}\big{)}\,\Big{(}\gamma\big{(}\alpha_{t}(x ),\bar{x}\big{)}-\bar{\gamma}\big{(}\alpha_{t}(x),\bar{x}\big{)}\Big{)}\bigg{] }+C\bigg{)}\;dx\\ \end{split}$ (31) $$

Consider minimizing the KLD (LABEL:eq:min_KLD_problem_expand) with respect to ${\alpha_{t}}$ as follows:

$$ $\displaystyle\begin{split}&\min_{\alpha_{t}}\;\;\operatorname{\mathcal{D}_{ \mathrm{KL}}}(\alpha_{t})\\ &=\min_{\alpha_{t}}\;\;\underbrace{\frac{\Delta_{t}}{2}\int_{\Omega}{\rho}^{ \Phi}_{t}(x)\,\frac{1}{\bar{\rho}_{t}\big{(}\alpha_{t}(x)\big{)}}\, \operatorname{\mathbb{E}}_{{p}_{\text{data}}(\bar{x})}\bigg{[}\rho_{t}\big{(} \alpha_{t}(x),\bar{x}\big{)}\,\Big{(}\gamma\big{(}\alpha_{t}(x),\bar{x}\big{)} -\bar{\gamma}\big{(}\alpha_{t}(x),\bar{x}\big{)}\Big{)}\bigg{]}\;dx}_{\mathcal {D}^{\mathrm{KL}}_{1}(\alpha_{t})}\\ &\qquad\qquad\!\!\!-\underbrace{\int_{\Omega}{\rho}^{\Phi}_{t}(x)\,\ln\bar{ \rho}_{t}(\alpha_{t}(x))\;dx}_{\mathcal{D}^{\mathrm{KL}}_{2}(\alpha_{t})}\,-\, \underbrace{\int_{\Omega}{\rho}^{\Phi}_{t}(x)\,\ln|D\alpha_{t}|\;dx}_{\mathcal {D}^{\mathrm{KL}}_{3}(\alpha_{t})}\end{split}$ (32) $$

where we have neglected the constant terms that do not depend on ${\alpha_{t}}$.

To solve the optimization (32), we consider the following optimality condition in the first variation of $\operatorname{\mathcal{D}_{\mathrm{KL}}}$:

$$ $\displaystyle\mathcal{I}(\alpha,\nu)=\frac{d}{d\epsilon}\,\operatorname{ \mathcal{D}_{\mathrm{KL}}}\big{(}\alpha(x)+\epsilon\,\nu(x)\big{)}\,\bigg{|}_{ \epsilon=0}=0$ (33) $$

which must hold for all trial function $\nu(x)$.

Taking the variational derivative of the first functional $\mathcal{D}^{\mathrm{KL}}_{1}$ in (32), we have

$$ $\displaystyle\begin{split}&\mathcal{I}^{1}(\alpha,\nu)=\frac{d}{d\epsilon}\, \mathcal{D}^{\mathrm{KL}}_{1}(\alpha+\epsilon\nu)\,\bigg{|}_{\epsilon=0}\\ &=\frac{\Delta}{2}\int_{\Omega}{\rho}^{\Phi}(x)\,\frac{d}{d\epsilon}\bigg{\{} \frac{1}{\bar{\rho}(\alpha+\epsilon\nu)}\,\operatorname{\mathbb{E}}_{{p}_{ \text{data}}(\bar{x})}\Big{[}\rho(\alpha+\epsilon\nu,\bar{x})\,\big{(}\gamma( \alpha+\epsilon\nu,\bar{x})-\bar{\gamma}(\alpha+\epsilon\nu,\bar{x})\big{)} \Big{]}\bigg{\}}\,\Bigg{|}_{\epsilon=0}\;dx\\ &=\frac{\Delta}{2}\int_{\Omega}{\rho}^{\Phi}(x)\,D\,\bigg{\{}\frac{1}{\bar{ \rho}(x)}\,\operatorname{\mathbb{E}}_{{p}_{\text{data}}(\bar{x})}\Big{[}\rho(x ,\bar{x})\,\big{(}\gamma(x,\bar{x})-\bar{\gamma}(x,\bar{x})\big{)}\Big{]}\bigg {\}}\;\nu\;dx\end{split}$ (34) $$

where $Dg:=\nabla^{T}g$ denotes the Jacobian of function $g(x)$ with respect to $x$.

A Taylor series expansion of the derivative $\frac{\partial g}{\partial x_{i}}(\alpha)$ with respect to $x_{i}$ yields

$$ $\displaystyle\begin{split}\frac{\partial g(\alpha)}{\partial x_{i}}&=\frac{ \partial g(x+\Delta u)}{\partial x_{i}}=\frac{\partial g(x)}{\partial x_{i}}+ \Delta\,\sum_{j}\frac{\partial^{2}g(x)}{\partial x_{i}\,\partial x_{j}}\,u_{j} +O(\Delta^{2})\\ \end{split}$ (35) $$

Using the Taylor series expansion (35), (34) can be written in tensor notation as follows:

$$ $\displaystyle\begin{split}\mathcal{I}^{1}(\alpha,\nu)&=\frac{\Delta}{2}\int_{ \Omega}{\rho}^{\Phi}(x)\,\sum_{i}\,\frac{\partial}{\partial x_{i}}\bigg{\{} \frac{1}{\bar{\rho}(x)}\,\operatorname{\mathbb{E}}_{{p}_{\text{data}}(\bar{x}) }\Big{[}\rho(x;\bar{x})\,\big{(}\gamma(x,\bar{x})-\bar{\gamma}(x,\bar{x})\big{ )}\Big{]}\bigg{\}}\,\nu_{i}\;dx\;+\;O(\Delta^{2})\end{split}$ (36) $$

Taking the variational derivative of the second functional $\mathcal{D}^{\mathrm{KL}}_{2}$ in (32) yields

$$ $\displaystyle\begin{split}\mathcal{I}^{2}(\alpha,\nu)&=\frac{d}{d\epsilon}\, \mathcal{D}^{\mathrm{KL}}_{2}(\alpha+\epsilon\nu)\,\bigg{|}_{\epsilon=0}\\ &=\int_{\Omega}{\rho}^{\Phi}(x)\,\frac{d}{d\epsilon}\ln\bar{\rho}(\alpha+ \epsilon\nu)\,\bigg{|}_{\epsilon=0}\;dx\\ &=\int_{\Omega}{\rho}^{\Phi}(x)\,\frac{1}{\bar{\rho}(\alpha)}\,\nabla\bar{\rho }(\alpha)\cdot\nu\;dx\\ &=\int_{\Omega}{\rho}^{\Phi}(x)\,\nabla\ln\bar{\rho}(\alpha)\cdot\nu\;dx\\ \end{split}$ (37) $$

where we have used the derivative identity $d\ln g=\frac{1}{g}\,dg$ to obtain the second equation.

Using the Taylor series expansion (35), (37) can be written in tensor notation as follows:

$$ $\displaystyle\begin{split}\mathcal{I}^{2}(\alpha,\nu)&=-\int_{\Omega}{\rho}^{ \Phi}(x)\,\sum_{i}\Bigg{(}\frac{\partial\ln\bar{\rho}(x)}{\partial x_{i}}- \Delta\,\sum_{j}\,\frac{\partial^{2}\ln\bar{\rho}(x)}{\partial x_{i}\,\partial x _{j}}\,u_{j}\Bigg{)}\,\nu_{i}\;dx\;+\;O(\Delta^{2})\\ &=-\int_{\Omega}{\rho}^{\Phi}(x)\,\sum_{i}\Bigg{(}\frac{\partial\ln\bar{\rho}( x)}{\partial x_{i}}-\Delta\,\sum_{j}\,\frac{\partial^{2}\ln\bar{\rho}(x)}{ \partial x_{i}\,\partial x_{j}}\,u_{j}\Bigg{)}\,\nu_{i}\;dx\;+\;O(\Delta^{2}) \end{split}$ (38) $$

Similarly, taking the variational derivative of the $\mathcal{D}^{\mathrm{KL}}_{3}$ term in (32), we have

$$ $\displaystyle\begin{split}\mathcal{I}^{3}(\alpha,\nu)&=\frac{d}{d\epsilon}\, \mathcal{D}^{\mathrm{KL}}_{3}(\alpha+\epsilon\nu)\,\bigg{|}_{\epsilon=0}\\ &=\int_{\Omega}{\rho}^{\Phi}(x)\,\frac{d}{d\epsilon}\ln\big{|}D(\alpha+ \epsilon\nu)\big{|}\,\bigg{|}_{\epsilon=0}\;dx\\ &=\int_{\Omega}{\rho}^{\Phi}(x)\,\frac{1}{\big{|}D\alpha\big{|}}\,\frac{d}{d \epsilon}\big{|}D(\alpha+\epsilon\nu)\big{|}\,\bigg{|}_{\epsilon=0}\;dx\\ &=\int_{\Omega}{\rho}^{\Phi}(x)\,\operatorname{\mathrm{tr}}\big{(}{D\alpha}^{- 1}D\nu\big{)}\;dx\end{split}$ (39) $$

where we have used the following Jacobi’s formula:

$$ $\displaystyle\frac{d}{d\epsilon}\big{|}D(\alpha+\epsilon\nu)\big{|}\,\bigg{|}_ {\epsilon=0}=\left|D\alpha\right|\operatorname{\mathrm{tr}}\left(D\alpha^{-1}D \nu\right)$ (40) $$

to obtain the last equation in (39).

The inverse of Jacobian $D\alpha^{-1}$ can be expanded via Neuman series to obtain

$$ $\displaystyle\begin{split}D\alpha^{-1}&=\big{(}\,\mathrm{I}+\Delta\,Du\,\big{) }^{-1}=\mathrm{I}-\Delta\,Du\;+\;O(\Delta^{2})\end{split}$ (41) $$

Substituting in (41) and using the Taylor series expansion (35), (37) can be written in tensor notation as follows:

$$ $\displaystyle\begin{split}\mathcal{I}^{3}(\alpha,\nu)&=\int_{\Omega}\sum_{i} \Bigg{(}{\rho}^{\Phi}(x)\,\frac{\partial\nu_{i}}{\partial x_{i}}-\Delta\sum_{j }\,{\rho}^{\Phi}(x)\,\frac{\partial u_{j}}{\partial x_{i}}\,\frac{\partial\nu_ {i}}{\partial x_{j}}\Bigg{)}\;dx\,+\,O(\Delta^{2})\\ &=\int_{\Omega}\sum_{i}\Bigg{(}\frac{\partial{\rho}^{\Phi}(x)}{\partial x_{i}} \,\nu_{i}-\Delta\sum_{j}\,\frac{\partial}{\partial x_{j}}\bigg{\{}{\rho}^{\Phi }(x)\,\frac{\partial u_{j}}{\partial x_{i}}\bigg{\}}\,\nu_{i}\Bigg{)}\;dx\,+\, O(\Delta^{2})\\ &=\int_{\Omega}\sum_{i}\Bigg{(}\frac{\partial{\rho}^{\Phi}(x)}{\partial x_{i}} -\Delta\sum_{j}\,\frac{\partial}{\partial x_{j}}\bigg{\{}{\rho}^{\Phi}(x)\, \frac{\partial u_{j}}{\partial x_{i}}\bigg{\}}\Bigg{)}\,\nu_{i}\;dx\,+\,O( \Delta^{2})\end{split}$ (42) $$

where we have used integration by parts to obtain the second equation.

Taking the limit $\lim\Delta\rightarrow 0$, the terms $O(\Delta^{2})$ that approach zero exponentially vanish.
Subtracting (36) by (38) and (42) then equating to zero, we obtain the first-order optimality condition (33) as follows:

$$ $\displaystyle\begin{split}\int_{\Omega}\bar{\rho}(x)\,\sum_{i}\Bigg{(}\sum_{j} &-\,\frac{\partial}{\partial x_{i}}\bigg{\{}\frac{1}{\bar{\rho}(x)}\,\frac{ \partial}{\partial x_{j}}\Big{\{}\bar{\rho}(x)\,u_{j}\Big{\}}\bigg{\}}\\ &+\frac{1}{2}\,\frac{\partial}{\partial x_{i}}\bigg{\{}\frac{1}{\bar{\rho}(x)} \,\operatorname{\mathbb{E}}_{{p}_{\text{data}}(\bar{x})}\Big{[}\rho(x;\bar{x}) \,\big{(}\gamma(x,\bar{x})-\bar{\gamma}(x,\bar{x})\big{)}\Big{]}\bigg{\}}\Bigg {)}\;\nu_{i}\;dx=0\end{split}$ (43) $$

where we have assumed that ${\rho}^{\Phi}(x)\equiv\bar{\rho}(x)$ holds, and used the following identities:

$$ $\displaystyle\begin{split}\frac{\partial\ln\bar{\rho}(x)}{\partial x_{i}}&= \frac{1}{\bar{\rho}(x)}\,\frac{\partial\bar{\rho}(x)}{\partial x_{i}}\\ \frac{\partial^{2}\ln\bar{\rho}(x)}{\partial x_{i}\,\partial x_{j}}&=\frac{ \partial}{\partial x_{i}}\bigg{(}\frac{1}{\bar{\rho}(x)}\frac{\partial\bar{ \rho}(x)}{\partial x_{j}}\bigg{)}\end{split}$ (44) $$

Given that $\nu_{i}$ can take any value, the equation (43) holds (in the weak sense) only if the terms within the round bracket vanish.
Integrating this term with respect to the $x_{i}$, we are left with

$$ $\displaystyle\begin{split}\sum_{j}\,\frac{\partial}{\partial x_{j}}\Big{\{} \bar{\rho}(x)\,u_{j}\Big{\}}=\frac{1}{2}\,\operatorname{\mathbb{E}}_{{p}_{ \text{data}}(\bar{x})}\Big{[}\rho(x;\bar{x})\,\big{(}\gamma(x,\bar{x})-\bar{ \gamma}(x,\bar{x})\big{)}\Big{]}\,+\,\bar{\rho}(x)\,C\end{split}$ (45) $$

which can also be written in vector notation as follows:

$$ $\displaystyle\begin{split}&\nabla\cdot\big{(}\bar{\rho}(x)\,u\big{)}=\frac{1}{ 2}\,\operatorname{\mathbb{E}}_{{p}_{\text{data}}(\bar{x})}\Big{[}\rho(x;\bar{x })\,\big{(}\gamma(x,\bar{x})-\bar{\gamma}(x,\bar{x})\big{)}\Big{]}\,+\,\bar{ \rho}(x)\,C\end{split}$ (46) $$

To find the scalar constant $C$, we integrate both sides of (46) to get

$$ $\displaystyle\begin{split}\int_{\Omega}\nabla\cdot\big{(}\bar{\rho}(x)\,u\big{ )}\;dx&=\frac{1}{2}\,\int_{\Omega}\operatorname{\mathbb{E}}_{{p}_{\text{data}} (\bar{x})}\Big{[}\rho(x;\bar{x})\,\big{(}\gamma(x,\bar{x})-\bar{\gamma}(x,\bar {x})\big{)}\Big{]}\;dx\,+\,\int_{\Omega}\bar{\rho}(x)\,C\;dx\\ &=\frac{1}{2}\,\int_{\Omega}\operatorname{\mathbb{E}}_{{p}_{\text{data}}(\bar{ x})}\Big{[}\rho(x;\bar{x})\,\big{(}\gamma(x,\bar{x})-\bar{\gamma}(x,\bar{x}) \big{)}\Big{]}\;dx\,+\,C\end{split}$ (47) $$

Applying the divergence theorem to the left-hand side of (47), we have

$$ $\displaystyle\begin{split}\int_{\Omega}\nabla\cdot\big{(}\bar{\rho}(x)\,u\big{ )}\;dx=\int_{\partial\Omega}\bar{\rho}(x)\,u\cdot\hat{n}\;dx\end{split}$ (48) $$

where $\hat{n}$ is the outward unit normal vector to the boundary $\partial\Omega$ of $\Omega$.

Given that $\bar{\rho}(x)$ is a normalized (proper) density with compact support (vanishes on the boundary), the term (48) becomes zero and we obtain $C=0$.
Substituting this and $u(x)=\nabla\Phi_{\theta}(x)$ into (46), we arrive at the PDE

$$ $\displaystyle\begin{split}\nabla\cdot\big{(}\bar{\rho}_{t}(x)\,\nabla\Phi(x) \big{)}=\frac{1}{2}\,\operatorname{\mathbb{E}}_{{p}_{\text{data}}(\bar{x})} \Big{[}\rho_{t}(x,\bar{x})\,\big{(}\gamma(x,\bar{x})-\bar{\gamma}(x,\bar{x}) \big{)}\Big{]}\end{split}$ (49) $$

Assume that the base case ${\rho}^{\Phi}_{0}(x)\equiv\bar{\rho}_{0}(x)$ holds, and that there exists a solution to (49) for every $t$. The proposition follows by the principle of induction.

∎

###### Proof.

### A.3 Proof of Proposition 3

###### Proof.

The energy loss function in (16) can generally be written as follows:

$$ $\displaystyle\begin{split}\mathcal{L}(\Phi,t)=&\ \frac{1}{2}\,\operatorname{ \mathbb{E}}_{\rho(x;\bar{x},t)\,{p}_{\text{data}}(x)}\Big{[}\Phi(x)\,\big{(} \gamma(x,\bar{x})-\bar{\gamma}(x,\bar{x})\big{)}\Big{]}+\frac{1}{2}\, \operatorname{\mathbb{E}}_{\bar{\rho}(x;t)}\,\Big{[}\big{\|}\nabla\Phi(x)\big{ \|}^{2}\Big{]}\end{split}$ (50) $$

where we have assumed, without loss of generality, that a normalized potential energy $\bar{E}_{\theta}(x;t)=0$. For an unnormalized solution $\Phi(x)$, we can always obtain the desired normalization by subtracting its mean.

The optimal solution $\Phi$ of the functional (50) is given by the first-order optimality condition:

$$ $\displaystyle\mathcal{I}(\Phi,\Psi)$ $\displaystyle=\frac{d}{d\epsilon}\,\mathcal{L}(\Phi(x)+\epsilon\Psi(x),t)\, \bigg{|}_{\epsilon=0}=0$ (51) $$

which must hold for all trial function $\Psi$.

Taking the variational derivative of the particle flow objective (51) with respect to $\epsilon$, we have

$$ $\displaystyle\begin{split}&\mathcal{I}(\Phi,\Psi)=\frac{d}{d\epsilon}\, \mathcal{L}(\Phi+\epsilon\Psi)\,\bigg{|}_{\epsilon=0}\\ &=\frac{1}{2}\,\int_{\Omega\times\Omega}{p}_{\text{data}}(x)\,\rho(x;\bar{x}) \,\big{(}\gamma(x,\bar{x})-\bar{\gamma}(x,\bar{x})\big{)}\,\frac{d}{d\epsilon} (\Phi+\epsilon\Psi)\;d\bar{x}\,dx\\ &\quad\,+\frac{1}{2}\,\int_{\Omega}\bar{\rho}(x)\,\frac{d}{d\epsilon}\big{\|} \nabla(\Phi+\epsilon\Psi)\big{\|}^{2}\;dx\\ &=\frac{1}{2}\,\int_{\Omega\times\Omega}{p}_{\text{data}}(x)\,\rho(x;\bar{x}) \,\big{(}\gamma(x,\bar{x})-\bar{\gamma}(x,\bar{x})\big{)}\,\Psi\;d\bar{x}\,dx \;+\;\int_{\Omega}\bar{\rho}(x)\,\nabla\Phi\cdot\nabla\Psi\;dx\end{split}$ (52) $$

Given that $\Phi\in\mathcal{H}^{1}_{0}(\Omega;\bar{\rho})$, the energy values vanish on the boundary $\partial\Omega$. Therefore, the second summand of the last expression in (52) can be written, via multivariate integration by parts, as

$$ $\displaystyle\begin{split}&\int_{\Omega}\bar{\rho}(x)\,\nabla\Phi\cdot\nabla \Psi=-\int_{\Omega}\nabla\cdot\big{(}\bar{\rho}(x)\,\nabla\Phi\big{)}\,\Psi\; dx\end{split}$ (53) $$

By substituting (53) into (52), we get

$$ $\displaystyle\begin{split}\mathcal{I}(\Phi,\Psi)&=\frac{1}{2}\,\int_{\Omega}\, \int_{\Omega}{p}_{\text{data}}(x)\rho(x;\bar{x})\,\big{(}\gamma(x,\bar{x})- \bar{\gamma}(x,\bar{x})\big{)}\,\Psi\;d\bar{x}\,dx\;-\;\int_{\Omega}\nabla \cdot\big{(}\bar{\rho}(x)\,\nabla\Phi\big{)}\,\Psi\;dx\\ &=\int_{\Omega}\bigg{(}\,\frac{1}{2}\,\int_{\Omega}{p}_{\text{data}}(x)\rho(x; \bar{x})\,\big{(}\gamma(x,\bar{x})-\bar{\gamma}(x,\bar{x})\big{)}\;d\bar{x}\,- \,\int_{\Omega}\nabla\cdot\big{(}\bar{\rho}(x)\,\nabla\Phi\big{)}\bigg{)}\Psi \;dx\end{split}$ (54) $$

and equating it to zero, we obtain the weak formulation (15) of the probabilistic Poisson’s equation.

Given that the Poincaré inequality (17) holds, [, Theorem 2.2] presents a rigorous proof of existence and uniqueness for the solution $\Phi$ to the weak formulation (15), based on the Hilbert-space form of the Riesz representation theorem.
∎

###### Proof.

### A.4 Derivation of Time-varying Mean and Variance in ( 19 )

Given the following marginal Gaussian distribution for $z$ and a conditional Gaussian distribution for $x$ given $x$, as defined in Section 3.1:

$$ $\displaystyle q(x)=\mathcal{N}(x;0,\Lambda)$ (55a) $\displaystyle p(\bar{x}|x)=\mathcal{N}(\bar{x};x,\Pi)$ (55b) $$

The posterior distribution of $x$ given $\bar{x}$ is obtained via Bayes’ theorem as

$$ $\displaystyle\begin{split}p(x|\bar{x})=\frac{p(\bar{x}|x)\,q(x)}{\int_{\Omega} p(\bar{x}|x)\,q(x)\;dx}=\mathcal{N}(x;\mu,\Sigma)\end{split}$ (56) $$

and remains a Gaussian, whose mean and variance are given by:

$$ $\displaystyle\mu(\bar{x})=\Sigma\,\Pi^{-1}\,\bar{x}$ (57a) $\displaystyle\Sigma=\big{(}\Lambda^{-1}\,+\,\Pi^{-1}\big{)}^{-1}$ (57b) $$

In fact, the conditional homotopy (6) can be written as

$$ $\displaystyle\begin{split}\rho(x;\bar{x},t)=\frac{p(\bar{x};t|x)\,q(x)}{\int_{ \Omega}p(\bar{x};t|x)\,q(x)\;dx}\end{split}$ (58) $$

where

$$ $\displaystyle\begin{split}p(\bar{x};t|x)=\mathcal{N}(x;\mu,\frac{1}{t}\Pi)\end {split}$ (59) $$

Notice that the terms involving $t$ in the numerator and denominator of (58) cancel each other out.
Substituting the variance of (59) into (58) and using (56)-(57), we obtain (19).

## Appendix B Experimental Details

### B.1 Model architecture

Our network structure is based on the Wide ResNet []. We adopt the same model hyperparameters used in [] for different datasets. In particular, the number of downsampled resolutions increases with the image size of the dataset and the number of ResBlocks in each resolution varies. Nevertheless, there are a few major differences between our network model and the ones used in []:

- 1.
We replace LeakyReLU activations with Gaussian Error Linear Unit (GELU) activations [], which we found improves training stability and convergence.
- 2.
We do not use spectral normalization []; instead, we use weight normalization with data-dependent initialization [].
- 3.
Following [, ], our training includes an additional spectral regularization loss which penalizes the spectral norm of each convolutional layer in the Wide ResNet to regularize the sharpness of the energy model.
- 4.
We drop the time embedding to render the time variable implicit in our energy model.

### B.2 Training

We use the Lamb optimizer [] and a learning rate of $0.001$ for all the experiments. We find that Lamb performs better than Adam over large learning rates. Following [], we set a smaller $\beta_{1}$ of $0.5$ in Lamb for the more high-resolution CelebA and LSUN images to help stabilize training.
For CIFAR-10, CelebA and LSUN, we use a batch size of 256, 128 and 64, respectively.
For all experiments, we set a spectral gap constant $\lambda$ of $0.001$, and a sharpness constant $\varepsilon$ of $0.0001$ in our training. Here, we set the standard deviation $\omega$ of the prior density to be $1$ so that the data likelihood homotopy is variance-preserving. Also, we set the standard deviation $\sigma$ of conditional data likelihood to be $0.01$ so that the difference between samples $x$ and data $\bar{x}$ is indistinguishable to human eyes [].
All models are trained for 1.3M iterations on a single NVIDIA A100 (80GB) GPU.

### B.3 Numerical Solver

In our experiments, the default solver of ODEs used is the black box solver in the Scipy library with the RK45 method [] following []. Since time variable $t$ is implicit in our energy model $\Phi(x(t))$, we can set a longer ODE interval, allowing the additional ODE iterations to further refine the samples within regions of high likelihood and improve the quality of generated images. We observe that setting a terminal time $t_{\text{end}}$ of $1.625$ for the RK45 ODE solver gives the best results.

### B.4 Datasets

We use the following datasets in our experiments: CIFAR-10 [], CelebA [] and LSUN []. CIFAR-10 is of resolution $32\times 32$, and contains $50,000$ training images and $10,000$ test images. CelebA contains $202,599$ face images, of which $162,770$ are training images and $19,962$ are test images. For processing, we first clip each image to $178\times 178$ and then resize it to $64\times 64$. For LSUN, we use the church outdoor and bedroom categories, which contain $126,227$ and $3,033,042$ training images respectively. Both categories contain $300$ test images. For processing, we first crop each image to a square image whose side is of length which is the minimum of the height and weight, and then we resize it to $64\times 64$ or $128\times 128$. For resizing, we set the anti-alias to True. We apply horizontal random flip as data augmentation for all datasets during training.

### B.5 Quantitative Evaluation

We employ the FID and inception scores as quantitative evaluation metrics for assessing the quality of generated samples. For CIFAR-10, we compute the Frechet distance between $50,000$ samples and the pre-computed statistics on the training set [13]. For CelebA $64\times 64$, we follow the setting in [] where the distance is computed between $5,000$ samples and the pre-computed statistics on the test set. For model selection, we follow [] and pick the checkpoint with the smallest FID scores, computed on 2,500 samples every 10,000 iterations.

### B.6 Potential Societal Impact

Generative models is a rapidly growing field of study with overarching implications in science and society. Our work proposes a new generative model VAPO that allows image generation via efficient and adaptive sampling.
The usage of the proposed model could have both positive and negative outcomes depending on the downstream application. For example, VAPO can be used to efficiently produce high-quality image/audio samples via the fast backward ODE. At the same time, it could promote deepfake technology and undermine social security. Generative models are also vulnerable to backdoor adversarial attacks on publicly available training data. Addressing the above concerns requires further collaborative research efforts aimed at mitigating misuse of AI technology.

## Appendix C Mode Evaluation

In this section, we evaluate the mode coverage and over-fitting of the proposed VAPO.

### C.1 Model Over-fitting and Generalization

To assess over-fitting, Figure 4 plots the histogram of the energy outputs on the CIFAR-10 train and test dataset. The energy histogram shows that the learned energy model assigns similar energy values to both train and test set images. This indicates that VAPO generalizes well to unseen test data and extensively covers all the modes in the training data.

Figure: Figure 4: Histogram of energy output for CIFAR-10 train and test set.
Refer to caption: extracted/5745897/train_eval_energy_histogram_cropped.png

In addition, Figure 5 presents the nearest neighbors of the generated samples in the train set of CIFAR-10. It shows that nearest neighbors are significantly different from the generated samples, thus suggesting that our models do not over-fit the training data and generalize well across the underlying data distribution.

Figure: Figure 5: Generated samples and their five nearest neighbours in the CIFAR-10 train set based on pixel distance.
Refer to caption: extracted/5745897/knn4.png

### C.2 Out-of-Distribution Detection

We evaluate robustness of our proposed VAPO model to anomalous data by assessing its performance on unsupervised out-of-distribution (OOD) detection. Given that potential energy is conjugate to the approximate data likelihood, the energy model can be used to distinguish between the in-distribution and out-distribution samples based on the energy values it assigns.
In particular, the energy model trained on CIFAR-10 train set is used for assigning normalized energy values to in-distribution samples (CIFAR-10 test set) and out-distribution samples from various other image datasets. The area under the receiver operating characteristic curve (AUROC) is used as a quantitative metric to determine the efficacy of the VAPO model in OOD detection, where a high AUROC score indicates that the model correctly assigns low energy to out-distribution samples.

Table 3 compares the AUROC scores of VAPO with various likelihood-based and EBM-based models. The result shows that VAPO performs exceptionally well on the CIFAR-10 interpolated dataset. However, its performance is average on CIFAR-100 and SVHN. This suggests that the perturbation of training data using the data likelihood homotopy may not sufficiently explore the data space in comparison to MCMC methods. The investigation into the underlying cause is left for future work.

**Table 3: Comparison of AUROC scores $\uparrow$ for OOD detection on several datasets.**
| Models | CIFAR-10<br>interpolation | CIFAR-100 | SVHN |
| --- | --- | --- | --- |
| PixelCNN | 0.71 | 0.63 | 0.32 |
| GLOW | 0.51 | 0.55 | 0.24 |
| NVAE | 0.64 | 0.56 | 0.42 |
| EBM-IG | 0.70 | 0.50 | 0.63 |
| VAEBM | 0.70 | 0.62 | 0.83 |
| CLEL | 0.72 | 0.72 | 0.98 |
| DRL | - | 0.44 | 0.88 |
| VAPO (Ours) | 0.78 | 0.50 | 0.61 |

## Appendix D Additional Results

Figures 6 and 7 show additional examples of image interpolation on CIFAR-10 and CelebA $64\times 64$, respectively.
Figures 8 and 9 show additional uncurated examples of unconditional image generation on CIFAR-10 and CelebA $64\times 64$, respectively.

Figure: Figure 6: Additional interpolation results on unconditional CelebA $64\times 64$.
Refer to caption: extracted/5745897/cifar10_interp_big.png

Figure: Figure 7: Additional interpolation results on unconditional CelebA $64\times 64$.
Refer to caption: extracted/5745897/celeba_interp_big.png

Figure: Figure 8: Additional uncurated samples on unconditional CIFAR-10 $32\times 32$.
Refer to caption: extracted/5745897/cifar10_big.png

Figure: Figure 9: Additional uncurated samples on unconditional CelebA $64\times 64$.
Refer to caption: extracted/5745897/celeba_big_2.png

## Peer Reviews

confidence: 3
contribution: 2
initial_rating: 3
presentation: 1
questions: -	The reviewer observed from the appendix that the authors used a lot of tricks for the network architecture and optimizer. Is this fair compared to the baselines presented in Table 1.
soundness: 1
strengths: The paper aims to tackle an important problem that is reduce the complexity of energy-based generative models due to the reliance on MCMC sampling or auxiliary latent models.
summary: The paper introduces a new framework for energy-based generative models called Variational Potential Flow (VAPO). This approach eliminates the need for implicit MCMC sampling and does not depend on auxiliary latent models or cooperative training methods. The VAPO framework achieves this by learning a potential energy function path, where the gradient flow guides prior samples along an approximate data likelihood homotopy. Additionally, the authors develop an energy loss function through a variational formulation that leverages the KL divergence between the density evolution of the flow-driven prior and the data likelihood homotopy. The framework is tested on CIFAR-10 and CelebA datasets for unconditional image generation.
weaknesses: -	The presentation could be improved, as the authors included extensive mathematical content in the main paper. The reviewer recommends simplifying this mathematical material within the main text and focusing on presenting high-level concepts and core results to enhance accessibility for readers.
-	The paper allocates considerable space to unnecessary mathematical content, which detracts from the clarity and quality of the experimental section. Additionally, the proposed method underperforms compared to many other energy-based model (EBM) approaches. Given that the paper's objective is to simplify and stabilize EBM components, such as MCMC sampling and auxiliary latent models, the authors should include a comprehensive analysis. This should involve demonstrating aspects like training convergence and complexity to better support the method's effectiveness and contributions.
-	In its current form, the paper presents extensive mathematical content along with some experimental results, but lacks thorough analysis. Consequently, it is unclear how beneficial the proposed tools truly are. The authors should consider demonstrating the proposed framework on simple, synthetic datasets to provide clearer insights into its behavior and effectiveness.

---

confidence: 4
contribution: 3
initial_rating: 6
presentation: 3
questions: * Could the questions about hyperparameters in the first weakness be addressed?
* Could the questions about the relative performance of the proposed method and diffusion models parameterized by potential energy functions be addressed?

Technical:
* About the equation (6): if I am understanding correctly, joint probability of $x$ and $\bar{x}$ satisfies: $p( \bar{x} | x) q(x) = p(\bar{x}, x) = p (x | \bar{x}) p_\text{data} (\bar{x})$. But why does $p( \bar{x} | x) q(x) = p (x | \bar{x}) p_\text{data} (\bar{x})$?
* In (41) second and third equality, the signs be reversed after integration by parts, right?
* I am not sure how (46) and (47) imply $C=0$. It looks like it should instead be: $C = - (1/2) \int E [ \rho(x; \bar{x} (\gamma (x, \bar{x}) - \bar{\gamma} (x, \bar{x}))] dx$.
soundness: 3
strengths: * This work brings an interesting and novel perspective to the study of generative modeling with potential energy. The density homotopy introduced in Section 3.1 is a fresh and potentially fruitful direction for using smooth interpolation between unnormalized densities in a manner similar to interpolation between data and prior samples in diffusion models and flow matching. Using this homotopy as the training target for a SDE with a learnable drift potential seems like a natural and effective choice.
* The theoretical background has a high degree of technical merit and thorough presentation. It is notable (and somewhat surprising) that the time evolution of the density homotopy has the simple and intuitive form in (8) and that the objective in (15) can be used to train a potential that can be used to match the homotopy over time. This could inspire future research in a similar direction.
* The experimental results show strong performance among EBM methods.
summary: This work presents a novel framework for learning energy-based SDEs that can transport samples from a latent prior distribution to a distribution that closely approximates data samples (up to a small additive noise). First, a conditional density homotopy between the prior distribution of noise and the posterior distribution of a slightly perturbed sample given a data sample is defined, along with an unconditional homotopy defined by taking the product of the data density and conditional homotopy and integrating w.r.t. the data variable. This defines a continuous sequence of densities interpolating between the prior and approximate data distribution which will be the target of the generative model. The generative model is an energy function whose gradients are trained to act as the drift term of an SDE with predetermined diffusion coefficients. The evolution of this SDE is intended to match the density homotopy, so that the gradients of the energy function can produce samples along the homotopy trajectory that start from the prior and end at the approximate data distribution. To achieve this, several propositions are used to present analytical forms for the evolution of the homotopy over time, establish a Poisson equation which is equivalent to minimizing the KL-divergence between the homotopy densities and the optimal trained model, and establish a practical loss function that can enforce the Poisson equation and be used to train a neural network. Experiments are conducted for unconditional image generation on CIFAR-10 and Celeb-A which show strong performance among EBMs.
weaknesses: * There are a few hyperparameters which are not fully explored, such as the noise schedule $\beta (t)$ (left constant), perturbation level $\sigma$, spectral gap $\lambda$, and $\varepsilon$. Some discussion of why these were chosen and the sensitivity to these hyperparameters would be helpful. In particular, the $\beta$ schedule is crucial for diffusion models. Why is it not essential here? Could the results be improved with a better schedule?
* Although experimental results are strong among EBMs, they lag behind the diffusion and other SOTA methods. Furthermore, this work also performs slightly worse than the related work [a] which learns a potential energy and uses the SDE (10) to draw samples, which is based on a straightforward diffusion objective (which is actually a component of the proposed loss). In general, the paper would be strengthened by comparing the results from the proposed method with a more straightforward diffusion model using a potential energy network instead of a score network (such as [a]), since both cases are essentially training a network so that (10) transports samples from a latent prior to data.

[a] https://openreview.net/forum?id=9AS-TF2jRNb

---

confidence: 4
contribution: 3
initial_rating: 6
presentation: 3
questions: Please check the weakness part.
soundness: 3
strengths: 1. The paper proposes a novel framework for learning EBM. VAPO has the advantage of more efficient training without the needs of expensive MCMC sampling.
2. The paper provides a sound and detailed derivation and formulation. I appreciate its theoretical contribution.
3. The generative performance of VAPO is competitive among the EBM baselines.
summary: This paper introduces VAPO, a new framework to train Energy-Based Model(EBM). Different from traditional EBM training algorithms. VAPO does not require expensive MCMC sampling at each training step. Instead, VAPO defines an interpolation path $p(x, t)$ between a noise distribution $q(x)$ at $t=0$ and distribution $\bar{p}(x)$ at $t=1$. $\bar{p}(x)$ is the kernel density approximation of the true data likelihood. The authors then introduce a potential function $\Phi(x_t)$, whose gradient guides the transport equation of the defined density path. The potential function is learned by constructing a variational formulation of the homotopy path-matching problem. After the model is learned, $\nabla\Phi(x_t)$ is used to guide the data transformation in ODE during the test time to generate new samples. Experimental results show that VAPO achieves competitive generative performance among EBM baselines.
weaknesses: My primary concern lies in the validity of the learned energy (or potential) function $\Phi(x_t)$. While the generative capability highlights one aspect of Energy-Based Models (EBMs), the ability to perform accurate density estimation is equally crucial. Although I acknowledge the competitive generative performance of VAPO, the relationship between the proposed potential function $\Phi(x_t)$ and the log data density remains unclear to me. Specifically:

1. Theoretically, what is the interpretation of the learned $\Phi(x_t)$? Are there any proofs ensuring that $\Phi(x_t)$ converges to the log data density? The paper demonstrates that the learned potential function $\Phi(x_t)$ can guide an ODE to generate valid samples via its gradient. However, a correct gradient field does not necessarily imply a valid energy function.

2. Experimentally, most of the results in the paper focus on demonstrating generative performance. I recommend that the authors include additional evaluations to demonstrate the validity of the learned energy function..

---

confidence: 4
contribution: 2
initial_rating: 3
presentation: 2
questions: 1. What assumptions are made about the potential flow in the theoretical analysis? For example, are there any assumptions about the target distribution? Are these assumptions too restrictive to be realistically met in experiments, or how are they satisfied in practice?

2. In Line 199, what is $\bar{p}(x)$ and why is it stated that $\bar{p}(x)$ acts as a continuous interpolation of the data likelihood $p_{\text{data}}(x)$? Could you clarify what is meant by "continuous interpolation" in this context and how $\bar{p}(x)$ relates to $p_{\text{data}}(x)$?
soundness: 3
strengths: 1. The VAPO method uses a potential energy function to guide the flow of prior samples toward the data likelihood, which is a fresh approach combining elements from particle flow and the Deep Ritz method. Experiments on CIFAR-10 and CelebA show that VAPO achieves competitive performance on FID scores relative to existing EBM-based approaches.
summary: This paper introduces a new generative model framework called VAPO. This framework learns a potential energy function to guide prior samples through a density evolution that approximates the data likelihood, bypassing the need for implicit and unstable MCMC sampling. VAPO applies Poisson's equation and the deep Ritz method to define and solve the flow of prior samples, ensuring they align with the data likelihood homotopy. VAPO is tested on image datasets like CIFAR-10 and CelebA, achieving competitive FID scores for unconditional image generation compared to state-of-the-art models.
weaknesses: 1. VAPO-A adopts the architecture from VAEBM [1], yet its FID scores fall short of baseline performance. This raises concerns regarding the overall advantage of VAPO. If the authors aim to address the high variance training, computational complexity, and low flexibility of MCMC and strengthen the claim of stability and efficiency of VAPO, the authors should conduct targeted experiments that compare training variance or computation time with MCMC-based methods. For example, record the loss value for each training epoch and calculate the standard deviation, measure the average time taken to generate a single sample, and record the total time and memory consumption required for the model to converge and compare it with MCMC methods. This would highlight the practical benefits of VAPO beyond standard performance metrics like FID.

[1] Zhisheng Xiao, Karsten Kreis, Jan Kautz, and Arash Vahdat. {VAEBM}: A symbiosis between variational autoencoders and energy-based models. In International Conference on Learning Representations, 2021.

2. VAPO incorporates computationally costly methods, such as numerical SDE solver. The authors also acknowledge that VAPO requires a large number of training iterations to converge, which might be a significant drawback for scaling this model to more complex or higher-dimensional datasets. While the elimination of MCMC improves training stability, the optimization process remains computationally expensive.

3. Previous approaches using MCMC methods and those based on flow models should be discussed together in the Related Work section and clearly referenced in the main text. This would provide readers with a comprehensive context for understanding the advancements presented in VAPO. Additionally, Algorithm 1, which outlines the algorithmic workflow, should be included in the main body of the paper rather than placed in the appendix. This is a critical component for grasping the methodology and overall contributions of the work. I recommend placing the Related Work section as Section 2 and positioning Algorithm 1 before Section 4 to enhance the flow and accessibility of the content.

4. Some sections, particularly those involving complex mathematical derivations, such as Sections 3.2 and 3.3, are difficult to follow due to the density of the technical content. These parts could be made more accessible by including additional explanatory text or visual aids to help guide readers through the more intricate aspects of the methodology. Moreover, the sketch proof of the main theorem, which is essential for understanding the core contributions, should be included in the main text.

5. The citation formatting throughout the paper is often incorrect. For instance, many references are not properly enclosed in brackets, which disrupts the flow of reading and makes it difficult to distinguish between the main text and references.