Pseudo-Spherical Contrastive Divergence
Year: 2021

Abstract:
Energy-based models (EBMs) offer flexible distribution parametrization. However, due to the intractable partition function, they are typically trained via contrastive divergence for maximum likelihood estimation. In this paper, we propose pseudo-spherical contrastive divergence (PS-CD) to generalize maximum likelihood learning of EBMs. PS-CD is derived from the maximization of a family of strictly proper homogeneous scoring rules, which avoids the computation of the intractable partition function and provides a generalized family of learning objectives that include contrastive divergence as a special case. Moreover, PS-CD allows us to flexibly choose various learning objectives to train EBMs without additional computational cost or variational minimax optimization. Theoretical analysis on the proposed method and extensive experiments on both synthetic data and commonly used image datasets demonstrate the effectiveness and modeling flexibility of PS-CD, as well as its robustness to data contamination, thus showing its superiority over maximum likelihood and f -EBMs.

heading: 1 Introduction
text: Energy-based models (EBMs) provide a unified framework for generative and discriminative learning by capturing dependencies between random variables with an energy function. Due to the absence of the normalization constraint, EBMs offer much more flexibility in distribution parametrization and architecture design compared to properly normalized probabilistic models such as autoregressive models [53, 23], flow-based models [14, 15, 46] and sum-product networks [72]. Recently, deep EBMs have achieved considerable success in realistic image generation [17, 66, 12, 30], molecular modeling [90] and model-based planning [16], thanks to modern deep neural networks [54, 49, 35] for parametrizing expressive energy functions and improved Markov Chain Monte Carlo (MCMC) techniques [62, 75, 40, 17, 66] for efficiently sampling from EBMs.
Training EBMs consists of finding an energy function that assigns low energies to correct configurations of variables and high energies to incorrect ones [55], where a central concept is the loss functional that is used to measure the quality of the energy function and is minimized during training. The flexibility of EBMs does not come for free: it makes the design of loss functionals particularly challenging, as it usually involves the partition function that is generally intractable to compute. As a result, EBMs are typically trained via CD [37], which belongs to the “analysis by synthesis” scheme [31] and performs a sampling-based estimation of the gradient of KL between data distribution and energy-based distribution. Since different loss functionals will induce different solutions in practice (when the model is mis-specified and data is finite) and KL may not provide the right inductive
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
bias [25, 91], inspired by the great success of implicit generative models [28, 67, 3], [89] proposed a variational framework to train EBMs by minimizing general f -divergences [10]. Although this framework enables us to specify various modeling preferences such as diversity/quality tradeoff, they rely on learning additional components (variational functions) within a minimax framework, where the optimization is complicated by the notion of Nash equilibrium and local optimality [42] and suffers from instability and non-convergence issues [59]. Along this line, noise contrastive estimation (NCE) [34] can train EBMs with a family of loss functionals induced by different Bregman divergences. However, in practice, NCE usually relies on carefully-designed noise distribution such as context-dependent noise distribution [41] or joint learning of a flow-based noise distribution [21].
In this paper, we draw inspiration from statistical decision theory [13] and propose a novel perspective for designing loss functionals for training EBMs without involving auxiliary models or variational optimization. Specifically, to generalize maximum likelihood training of EBMs, we focus on maximizing pseudo-spherical scoring rules [76, 27], which are strictly proper such that the data distribution is the unique optimum and homogeneous such that they can be evaluated without the knowledge of the normalization constant. Under the “analysis by synthesis” scheme used in CD and f -EBM [89], we then derive a practical algorithm termed Pseudo-Spherical Contrastive Divergence (PS-CD). Different from f -EBM, PS-CD enables us to specify flexible modeling preferences without requiring additional computational cost or unstable minimax optimization. We provide a theoretical analysis on the sample complexity and convergence property of PS-CD, as well as its connections to maximum likelihood training. With experiments on both synthetic data and commonly used image datasets, we show that PS-CD achieves significant sample quality improvement over conventional maximum likelihood training and competitive performance to f -EBM without expensive variational optimization. Based on a set of recently proposed generative model evaluation metrics [61], we further demonstrate the various modeling tradeoffs enabled by PS-CD, justifying its modeling flexibility. Moreover, PS-CD is also much more robust than CD in face of data contamination.

heading: 2 Preliminaries

heading: 2.1 Energy-Based Distribution Representation and Sampling
text: Given a set of i.i.d. samples {xi}Ni=1 from some unknown data distribution p(x) defined over the sample space X ⊂ Rm, the goal of generative modeling is to learn a θ-parametrized probability distribution qθ(x) to approximate the data distribution p(x). In the context of energy-based modeling, instead of directly parametrizing a properly normalized distribution, we first parametrize an unnormalized energy function Eθ : X → R, which further defines a normalized probability density via the Boltzmann distribution:
qθ(x) = qθ(x)
Zθ = exp(−Eθ(x)) Zθ , (1)
where Zθ := ∫ X exp(−Eθ(x))dx is the partition function (normalization constant). In this paper, unless otherwise stated, we will use q to denote an unnormalized density and q to denote the corresponding normalized distribution. We also assume that the exponential of the negative energy belongs to the L1 space, E := { Eθ : X → R : ∫ X exp(−Eθ(x))dx <∞ } , i.e., Zθ is finite.
Since energy-based models (EBMs) represent a probability distribution by assigning unnormalized scalar values (energies) to the data points, we can use any model architecture that outputs a bounded real number given an input to implement the energy function, which allows extreme flexibility in distribution parametrization. However, it is non-trivial to sample from an EBM, usually requiring MCMC [75] techniques. Specifically, in this work we consider using Langevin dynamics [62, 88], a gradient-based MCMC method that performs noisy gradient descent to traverse the energy landscape and arrive at the low-energy configurations:
x̃t = x̃t−1 − 2 ∇xEθ(x̃t−1) +
√ zt, (2)
where zt ∼ N (0, I). The distribution of x̃T converges to the model distribution qθ(x) ∝ exp(−Eθ(x)) when → 0 and T → ∞ under some regularity conditions [88]. In order to sample from an energy-based distribution efficiently, many scalable techniques have been proposed such as learning non-convergent, non-persistent, short-run MCMC [66] and using a sample replay buffer to improve mixing time and sample diversity [17]. In this work, we leverage these recent advances when we need to obtain samples from an EBM.

heading: 2.2 Maximum Likelihood Training of EBMs via Contrastive Divergence
text: The predominant approach to training explicit density generative models is to approximately minimize the KL divergence between the (empirical) data distribution and model distribution. Minimizing KL divergence is equivalent to the following maximum likelihood estimation (MLE) objective:
min θ LMLE(θ; p) = min θ −Ep(x) [log qθ(x)] = min θ Ep(x)[Eθ(x)] + logZθ. (3)
Because of the intractable partition function (an integral over the sample space), we cannot directly optimize the above MLE objective. To tackle this issue, [37] proposed contrastive divergence (CD) algorithm as a convenient way to estimate the gradient of LMLE(θ; p) using samples from qθ: ∇θLMLE(θ; p) = Ep(x)[∇θEθ(x)] +∇θ logZθ = Ep(x)[∇θEθ(x)]− Eqθ(x)[∇θEθ(x)], (4)
which can be interpreted as decreasing the energies of real data from p and increasing the energies of fake data generated by qθ. As discussed above, evaluating Equation (4) typically relies on MCMC methods such as the Langevin dynamics sampling procedure defined in Equation (2) to produce samples from the model distribution qθ, which induces a surrogate gradient estimation:
∇θLCD−K(θ; p) = Ep(x)[∇θEθ(x)]− EqKθ (x)[∇θEθ(x)], (5)
where qKθ denotes the distribution after K steps of MCMC transitions from an initial distribution (typically data distribution or uniform distribution), and Equation (4) corresponds to LCD−∞.

heading: 2.3 Strictly Proper Scoring Rules
text: Stemming from statistical decision theory [13], scoring rules evaluate the quality of probabilistic forecasts by assigning numerical scores based on the predictive distributions and the events that materialize. Formally, consider a compact sample space X . Let M be a space of all locally 1- integrable non-negative finite measures and P be a subspace consisting of all probability measures on the sample space X . A scoring rule S(x, q) specifies the utility of forecasting using a probability forecast q ∈ P for a given sample x ∈ X . With slightly abused notation, we write the expected score of S(x, q) under a reference distribution p as:
S(p, q) := Ep(x)[S(x, q)]. (6)
Definition 1 (Proper Scoring Rules [26]). A scoring rule S : X × P → R is called proper relative to P if the corresponding expected score satisfies:
∀p, q ∈ P, S(p, q) ≤ S(p, p). (7) It is strictly proper if the equality holds if and only if q = p.
In prediction and elicitation problems, strictly proper scoring rules encourage the forecaster to make honest predictions based on their true beliefs [22]. In estimation problems, where we want to approximate a distribution p with another parametric distribution qθ, strictly proper scoring rules provide attractive learning objectives:
arg max qθ∈Pθ S(p, qθ) = arg max qθ∈Pθ
Ep(x)[S(x, qθ)] = p (when p ∈ Pθ). (8)
When a scoring rule S is strictly proper relative to P , the associated generalized entropy function and divergence function are defined as:
G(p) := sup q∈P
S(p, q) = S(p, p), D(p, q) := S(p, p)− S(p, q) ≥ 0. (9)
G(p) is convex and represents the maximally achievable utility, while D(p, q) is the Bregman divergence [6] associated with the convex function G and the equality holds only when p = q.
Next, we introduce a specific kind of scoring rules that are particularly suitable for learning unnormalized statistical models. Definition 2 (Homogeneous Scoring Rules [69]). A scoring rule is homogeneous if it satisfies (here the domain of the score function is extended to X ×M):
∀λ > 0,x ∈ X , S(x, q) = S(x, λ · q). (10)
Since scaling the model distribution q by a positive constant λ does not change the value of a homogeneous scoring rule, such homogeneity allows us to evaluate it without computing the intractable partition function of an energy-based distribution. Thus, strictly proper and homogeneous scoring rules are natural candidates for new training objectives of EBMs. Example 1. A notable example of scoring rules is the widely used logarithm score: S(x, q) = log q(x). The associated generalized entropy is the negative Shannon entropy: G(p) = Ep(x)[log p(x)], and the associated Bregman divergence is the KL divergence: D(p, q) = Ep(x) [log(p(x)/q(x))]. From Definitions 1 and 2, we know that the logarithm score is strictly proper but not homogeneous. Specifically, for a θ-parametrized energy-based distribution qθ = qθ/Zθ, since S(x, qθ) = S(x, Zθ · qθ) = S(x, qθ) + logZθ 6= S(x, qθ) and logZθ cannot be ignored during the optimization of θ, we need to use tailored methods such as contrastive divergence [37] or doubly dual embedding [11] to tackle the intractable partition function.

heading: 3 Training EBMs by Maximizing Homogeneous Scoring Rules
text: In this section, we derive a new principle for training EBMs from the perspective of optimizing strictly proper homogeneous scoring rules. All proofs for this section can be found in Appendix B.

heading: 3.1 Pseudo-Spherical Scoring Rule
text: In this section, we introduce the pseudo-spherical scoring rule, which is a representative family of strictly proper homogeneous scoring rules that have great potentials for training deep energy-based models and allow flexible and convenient specification of modeling preferences, yet have not been explored before in the context of energy-based generative modeling. Definition 3 (Pseudo-Spherical Scoring Rule [76, 27]). For γ > 0, the pseudo-spherical scoring rule is defined as:
S(x, q) := q(x)γ ( ∫ X q(y) γ+1dy) γ γ+1 = q(x)γ ( ∫ X q(y) γ+1dy) γ γ+1 =
( q(x)
‖q‖γ+1
)γ (11)
where ‖q‖γ+1 := (∫ X q(y) γ+1dy ) 1 γ+1 .
The expected pseudo-spherical score under a reference distribution p is defined as:
Sps(p, q) := Ep(x)[S(x, q)] = Ep(x)[q(x)γ ] ( ∫ X q(y) γ+1dy) γ γ+1
(12)
Example 2. The classic spherical scoring rule [19] is a special case in the pseudo-spherical family, which corresponds to γ = 1:
S(x, q) = q(x) ( ∫ X q(y) 2dy) 1 2 = q(x) ‖q‖2 (13)
The family of pseudo-spherical scoring rules is appealing because it introduces a different and principled way for assessing a probability forecast. For example, the spherical scoring rule has an interesting geometric interpretation. Suppose the sample space X contains n mutually exclusive and exhaustive outcomes. Then a probability forecast can be represented as a vector q = (q1, . . . , qn). Let vector p = (p1, . . . , pn) represent the oracle probability forecast. The expected spherical score can be written as:
S(p, q) = Ep(x)[S(x, q)] = ∑ i piqi√∑ i q 2 i = ‖p‖2 〈p, q〉 ‖p‖2‖q‖2 = ‖p‖2 cos(∠(p, q)) (14)
where 〈p, q〉 and ∠(p, q) denote the inner product and the angle between vectors p and q respectively. In other words, when we want to evaluate the expected spherical score of a probability forecast q under real data distribution p using samples, the angle between p and q is the sufficient statistics. Since we know that both p and q belong to the probability simplex P = {v| ∑ x∈X v(x) = 1 and ∀x ∈ X ,v(x) ≥ 0.}, the expected score will be minimized if and only if the angle is zero, which implies p = q. More importantly, since all we need to do is to minimize the angle of
the deviation, we are allowed to scale q by a constant. Specifically, when q is an energy-based distribution q = ( exp(−E1)∑ i exp(−Ei) , . . . , exp(−En)∑ i exp(−Ei) ) , we can instead evaluate and minimize the angle between data distribution p and the unnormalized distribution q = (exp(−E1), . . . , exp(−En)), since ∠(p, q) = ∠(p, q). More generally, we have the following theorem to justify the use of pseudo-spherical scoring rules for training energy-based models: Theorem 1 ([26, 70]). Pseudo-spherical scoring rule is strictly proper and homogeneous.
As the original definition of pseudo-spherical scoring rule (Equation (11)) takes the form of a fraction, for computational considerations, in this paper we instead focus on optimizing its composite scoring rule (Definition 2.1 in [44]): Definition 4 (γ-score [20]). For the expected pseudo-spherical score Sps(p, q) defined in Equation (12) with γ > 0, the expected γ-score is defined as:
Sγ(p, q) := 1
γ log(Sps(p, q)) =
1 γ log ( Ep(x)[q(x)γ ] ) − log(‖q‖γ+1) (15)
Since 1γ log(u) is strictly increasing in u, Sγ(p, q) is a strictly proper homogeneous composite score:
arg max q∈P Sγ(p, q) = arg max q∈P
1 γ log(Sps(p, q)) = arg max q∈P Sps(p, q) = p. (16)

heading: 3.2 Pseudo-Spherical Contrastive Divergence
text: Suppose we parametrize the energy-based model distribution as qθ ∝ qθ = exp(−Eθ) and we want to minimize the negative γ-score in Equation (15):
min θ Lγ(θ; p) = min θ − 1 γ
log ( Ep(x)[qθ(x)γ ] ) + log(‖qθ‖γ+1) (17)
In the following theorem, we derive the gradient of Lγ(θ; p) with respect to θ: Theorem 2. The gradient of Lγ(θ; p) with respect to θ can be written as:
∇θLγ(θ; p) = − 1
γ ∇θ log
( Ep(x)[exp(−γEθ(x))] ) − Erθ(x)[∇θEθ(x)] (18)
where the auxiliary distribution rθ is also an energy-based distribution defined as:
rθ(x) := qθ(x) γ+1∫ X qθ(x) γ+1dx = exp(−(γ + 1)Eθ(x))∫ X exp(−(γ + 1)Eθ(x))dx .
In App. B.1, we provide two different ways to prove the above theorem. The first one is more straightforward and directly differentiates through the term log(‖qθ‖γ+1). The second one leverages a variational representation of log(‖qθ‖γ+1), where the optimal variational distribution happens to take an analytical form of r∗θ(x) ∝ qθ(x)γ+1, thus avoiding the minimax optimization in other variational frameworks such as [89, 11, 12]. The main challenge in maximizing γ-score is that it is generally intractable to exactly compute the gradient of the second term in Equation (15).
During training, estimating the second term of Equation (18) requires us to obtain samples from the auxiliary distribution rθ ∝ exp(−(γ + 1)Eθ), while at test time, we want to sample from the model distribution qθ ∝ exp(−Eθ) that approximates the data distribution. Due to the restrict regularity conditions on the convergence of Langevin dynamics, in practice, we found it challenging to use the iterative sampling process in Equation (2) with a fixed number of transition steps and step size to produce samples from rθ and qθ simultaneously, as the temperature γ + 1 in rθ simply amounts to a linear rescaling to the energy function during training. Thus for generality, as in contrastive divergence [37, 17, 66] and f -EBM [89], we make the minimal assumption that we only have a sampling procedure to produce samples from qθ for both learning and inference procedures.
In this case, we can leverage the analytical form of rθ and self-normalized importance sampling [68] (which has been used to derive gradient estimators in other contexts such as importance weighted autoencoder [7, 18]) to obtain a consistent estimation of Equation (18):
Algorithm 1 Pseudo-Spherical Contrastive Divergence. 1: Input: Empirical data distribution pdata. Pseudo-spherical scoring rule hyperparameter γ. 2: Initialize energy function Eθ. 3: repeat 4: Draw a minibatch of samples {x+1 , . . . ,x + N} from pdata.
5: Draw a minibatch of samples {x−1 , . . . ,x − N} from qθ ∝ exp(−Eθ) (e.g., using Langevin dynamics with a sample replay buffer). 6: Update the energy function by stochastic gradient descent:
̂∇θLNγ (θ; p) = −∇θ 1
γ log
( 1
N N∑ i=1 exp(−γEθ(x+i ))
) − ∑N i=1 exp(−γEθ(x − i ))∇θEθ(x − i )∑N
i=1 exp(−γEθ(x − i ))
7: until Convergence
Theorem 3. Let x+1 , . . . ,x + N be i.i.d. samples from p(x) and x − 1 , . . . ,x − N be i.i.d. samples from qθ(x) ∝ exp(−Eθ(x)). Define the gradient estimator as:
̂∇θLNγ (θ; p) := −∇θ 1
γ log
( 1
N N∑ i=1 exp(−γEθ(x+i ))
) − ∑N i=1 ωθ(x − i )∇θEθ(x − i )∑N
i=1 ωθ(x − i )
(19)
where the self-normalized importance weight ωθ(x−i ) := rθ(x − i )/qθ(x − i ) = exp(−γEθ(x − i )). Then the gradient estimator converges to the true gradient (Equation (18)) in probability, i.e., ∀ > 0:
lim N→∞
P (∥∥∥ ̂∇θLNγ (θ; p)−∇θLγ(θ; p)∥∥∥ ≥ ) = 0.
We summarize the pseudo-spherical contrastive divergence (PS-CD) training procedure in Algorithm 1. In Appendix A, we also provide a simple PyTorch implementation for stochastic gradient descent (SGD) with the gradient estimator in Equation (19).

heading: 3.3 Connections to Maximum Likelihood Estimation and Extension to γ < 0
text: From Equation (9) in Section 2.3, we know that γ-score induces the following Bregman divergence (the divergence function associated with proper composite scoring rule is analogously defined in Def. 2.1 in [44]):
Dγ(p, qθ) = Sγ(p, p)− Sγ(p, qθ) and maximizing γ-score is equivalent to minimizing Dγ(p, qθ). In the following lemma, we show that when γ → 0, Dγ(p, qθ) will recover the KL divergence between p and qθ, and the gradient of PS-CD will recover the gradient of CD. Lemma 1. When γ → 0, we have:
lim γ→0 Dγ(p, qθ) = DKL(p‖qθ); lim γ→0 ∇θLγ(θ; p) = ∇θLMLE(θ; p).
Inspired by [86, 56] that generalize Rényi divergence beyond its definition to negative orders, we now consider the extension of γ-score with γ < 0 (although it may not be strictly proper for these γ values). The following lemma shows that maximizing such scoring rule amounts to maximizing a lower bound of logarithm score (MLE) with an additional Rényi entropy regularization. Lemma 2. When −1 ≤ γ < 0, we have:
Sγ(p, q) ≤ Ep(x)[log q(x)] + γ
γ + 1 Hγ+1(q)
whereHγ+1(q) is the Rényi entropy of order γ + 1.

heading: 4 Theoretical Analysis
text: In this section, to gain a deeper understanding of our PS-CD algorithm and how the proposed estimator behaves, we analyze its sample complexity and convergence property under certain conditions. All the proofs for this section can be found in Appendix C.

heading: 4.1 Sample Complexity
text: We start with analyzing the sample complexity of the consistent gradient estimator in Equation (19), that is how fast it approaches the true gradient value. We first make the following assumption: Assumption 1. The energy function is bounded by K and the gradient is bounded by L:
∀x ∈ X , θ ∈ Θ, |Eθ(x)| ≤ K, ‖∇θEθ(x)‖ ≤ L.
This assumption is typically easy to satisfy because in practice we always use a bounded sample space (e.g. normalizing images to [0,1] or truncated Gaussian) to ensure stability. For example, in image modeling experiments, we use L2 regularization on the outputs of the energy function (hence bounded energy values), as well as normalized inputs and spectral normalization [60] for the neural network that realizes the energy function (hence Lipschitz continuous with bounded gradient).
With vector Bernstein inequality [47, 32], we have the following theorem showing a sample complexity of O
( log(1/δ)
2
) such that the estimation error is less than with probability at least 1− δ:
Theorem 4. For any constants > 0 and δ ∈ (0, 1), when the number of samples N satisfies:
N ≥ 32L 2e8γK (1 + 4 log(2/δ))
2
we have:
P (∥∥∥ ̂∇θLNγ (θ; p)−∇θLγ(θ; p)∥∥∥ ≤ ) ≥ 1− δ.

heading: 4.2 Convergence
text: Typically, convergence of SGD are analyzed for unbiased gradient estimators, while the gradient estimator in PS-CD is asymptotically consistent but biased. Building on the sample complexity bound above and prior theoretical works for analyzing SGD [8, 24], we analyze the convergence of PS-CD. For notational convenience, we use L(θ) to denote the loss function Lγ(θ; p) = −Sγ(p, qθ). Besides Assumption 1, we further make the following assumption on the smoothness of L(θ): Assumption 2. The loss function L(θ) is M -smooth (with M > 0):
∀θ1,θ2 ∈ Θ, ‖∇L(θ1)−∇L(θ2)‖ ≤M‖θ1 − θ2‖.
This is a common assumption for analyzing first-order optimization methods, which is also used in [24, 8]. Also this is a relatively mild assumption since we do not require the loss function to be convex in θ. Since in non-convex optimization, the convergence criterion is typically measured by gradient norm, following [64, 24], we use ‖∇L(θ)‖ ≤ ξ to judge whether a solution θ is approximately a stationary point. Theorem 5. For any constants α ∈ (0, 1) and δ ∈ (0, 1), suppose that the step sizes satisfy ηt < 2(1 − α)/M and the sample size Nt used for estimating ĝt is sufficiently large (satisfying Equation (36)). Let L∗ denote the minimum value of L(θ). Then with probability at least 1− δ, the output of Algorithm 2 (in Appendix C.2), θ̂, satisfies (constant C := αML2e4γK):
E[‖∇L(θ̂)‖2] < 2(L(θ1)− L∗) + 12C
∑T t=1 η
2 t∑T
t=1(2(1− α)ηt −Mη2t )
The above theorem implies the following corollary, which shows a typical convergence rate of O(1/ √ T ) for non-convex optimization problems:
Corollary 1. Under the conditions in Theorem 5 except that we use constant step sizes: ηt = min{(1− α)/M, 1/ √ T} for t = 1, . . . , T . Then with probability at least 1− δ, we have (constant C := αML2e4γK):
E[‖∇L(θ̂)‖2] < 2M(L(θ1)− L ∗) (1− α)2T + 2(L(θ1)− L∗) + 12C (1− α) √ T
In Appendix C.2, we discuss more on the strongly convex (Theorem 7) and convex cases (Theorem 8).

heading: 5 Related Work
text: Direct KL Minimization. Under the “analysis by synthesis” scheme, [37] proposed Contrastive Divergence (CD), which estimates the gradient of the log-partition function (arising from KL) using samples from some MCMC procedure. To improve the mixing time of MCMC, [17] proposed to employ Persistent CD and a replay buffer to store intermediate samples from Markov chains throughout training, and [66] proposed to learn non-convergent short-run MCMC. Both approaches (long-run and short-run MCMC) work well with PS-CD in our experiments. PS-CD may also benefit from recent advances on unbiased MCMC [40, 73], which we leave as interesting future work.
Fenchel Duality. By exploiting the primal-dual view of KL, recent works [11, 12, 2] proposed to cast maximum likelihood training of EBMs as minimax problems, which introduce a dual sampler and are approximately solved by alternating gradient descent ascent updates. Along this line, to allow flexible modeling preferences, [89] proposed f -EBM to enable the use of any f -divergence to train EBMs, which also relies on minimax optimization. By contrast, in this work, we leverage the analytical form of the optimal variational distribution and self-normalized importance sampling to reach a framework that requires no adversarial training and has no additional computational cost compared to CD while allowing flexible modeling preferences. Besides convenient optimization, PS-CD and f -EBM trains EBMs with two different families of divergences (hence complementary) with KL being the only shared one, since any pseudo-spherical scoring rule corresponds to a Bregman divergence (Section 2.3) and the only member in f -divergence that is also Bregman divergence is α-divergence (with KL as special case) (Theorem 4 in [1]).
Homogeneous Scoring Rules. [84] proposed to learn unnormalized statistical models on discrete sample space by maximizing γ-score, which uses empirical data distribution (p̂(x) = nx/n, where nx is the number of appearance of x in the dataset and n is the total number of data) as a surrogate to the real data density p and relies on a localization trick to bypass the computation of ‖qθ‖γ+1. Consequently, it is only amenable to finite discrete sample space such as natural language [51], whereas PS-CD is applicable to any unnormalized probabilistic model in continuous domains. Another popular homogeneous scoring rule is the Hyvärinen score, which gives rise to the score matching objective [39] for EBM training. However, score matching and its variants [87, 82] have difficulties in low data density regions and do not perform well in practice when training EBMs on high-dimensional datasets [80]. Moreover, since the score matching objective involves the Hessian of log-density functions that is generally expensive to compute [58], methods such as approximate propagation [45], curvature propagation [58] and sliced score matching [83] are needed to approximately compute the trace of the Hessian.
Noise Contrastive Estimation. Another principle for learning EBMs is Noise Contrastive Estimation (NCE) [34], where an EBM is learned by contrasting a prescribed noise distribution with tractable density against the unknown data distribution. Using various Bregman divergences, NCE can be generalized to a family of different loss functionals [33, 85]. However, finding an appropriate noise distribution for NCE is highly non-trivial. In practice, NCE typically works well in conjunction with a carefully-designed noise distribution such as context-dependent noise distribution [41] or joint learning of a flow-based noise distribution [21].
In this work, we focus on generalizing maximum likelihood by deriving novel training objectives for EBMs without involving auxiliary models (e.g., the variational function in [89], the flow-based noise distribution in [21] and the amortized sampler in [50, 12, 29]).

heading: 6 Experiments
text: In this section, we demonstrate the effectiveness of PS-CD on several 1-D and 2-D synthetic datasets as well as commonly used image datasets.
Setup. The 2-D synthetic datasets include Cosine, Swiss Roll, Moon, Mixture of Gaussian, Funnel and Rings, which cover different modalities and geometries (see Figure 2 in App. D.1 for illustration). To test the practical usefulness, we use MNIST [54], CIFAR-10 [48] and CelebA [57] in our experiments for modeling natural images. Following [80], for CelebA, we first center-crop the images to 140 × 140, then resize them to 64 × 64. More experimental details about the data processing, model architectures, sampling strategies and additional experimental results can be found in App. D.
Effects of Different γ Values. To illustrate the modeling flexibility brought by PS-CD and provide insights on the effects of different γ values, we first conduct a 1-D synthetic experiment similar to the one in [89]. As shown in Figure 1, when fitting a mixture of Gaussian with a single Gaussian (i.e., model mis-specification case), the family of PS-CD offers flexible tradeoffs between quality and diversity (i.e., mode collapse vs. mode coverage). Although in the well-specified case these objectives induce the same optimal solution, in this example, we can see that a larger γ leads to higher entropy. More importantly, compared to f -EBM [89]
that also provides similar modeling flexibility and includes CD as a special case, PS-CD does not require expensive and unstable minimax optimization (no additional computational cost compared to CD). In App. D.2, we further visualize the objective landscapes for different γ values. As shown in
Figure 3 and 4, when the model is well-specified, different objectives will induce the same optimal solution since they are strictly proper; when the model is mis-specified (corresponding to practical scenarios), different objectives will exhibit different modeling preferences.
Furthermore, to better demonstrate the modeling flexibility brought by PS-CD in high-dimensional case, we conduct experiments on CIFAR-10 using a set of indicative and reliable metrics (Density, Coverage, Precision, Recall) proposed by [61] to evaluate the effects of γ from various perspectives. Please refer to App. D.3 for experimental results (Table 4) and detailed discussions.
Table 1: FID scores for CIFAR-10 conditional and CelebA unconditional image generation. We list comparisons with results reported by CD [17], Noise-Conditional Score Network (NCSN) [81] and f -EBMs [89]. γ = 1.0 corresponds to maximizing spherical scoring rule (Example 2).
Method FID
CIFAR-10 (32× 32) Conditional Contrastive Divergence (KL) 37.90 f -EBM (KL) 37.36 f -EBM (Reverse KL) 33.25 f -EBM (Squared Hellinger) 32.19 f -EBM (Jensen Shannon) 30.86
Pseudo-Spherical CD (γ = 2.0) 33.19 Pseudo-Spherical CD (γ = 1.0) 29.78 Pseudo-Spherical CD (γ = 0.5) 35.02 Pseudo-Spherical CD (γ = −0.5) 27.95
CelebA (64× 64) Contrastive Divergence (KL) 26.10 NCSN (w/o denoising) 26.89 NCSN (w/ denoising) 25.30 NCSNv2 (w/o denoising) 28.86 NCSNv2 (w/ denoising) 10.23 Pseudo-Spherical CD (γ = 1.0) 24.76 Pseudo-Spherical CD (γ = −0.5) 20.35
2-D Synthetic Data. For quantitative evaluation of the 2-D synthetic data experiments, we follow [79] and report the maximum mean discrepancy (MMD, [5]) between the generated samples and validation samples in
Table 3 in App. D.1, which demonstrates that PS-CD outperforms its CD counterpart on all but the Funnel dataset. From the histograms of samples shown in
Figure 2 in App. D.1, we can also have similar observations. For example, CD fails to place high densities in the center of the right mode in MoG, while PS-CD places the modes correctly.
Image Generation. In Figure 5 in App.D.4, we show MNIST, CIFAR-10 and CelebA samples produced by PS-CD (with γ = 1.0), which demonstrate that our approach can produce highly realistic images with simple model architectures. As suggested in [17], we use Fréchet Inception Distance (FID) [36] as the quantitative evaluation metric for CIFAR-10 and CelebA, as Langevin dynamics converge to local minima that artificially inflate Inception Score [77]. From Table 11, we can see that various members (different γ values) in the family of PS-CD can outperform CD significantly, and more surprisingly, PS-CD also shows competitive performance to the recently proposed f -EBMs, without requiring expensive minimax optimization. While our method currently does not outperform the stateof-the-art image generation methods such as improved denoising score matching [81], which relies on carefully selected noise schedule and specially designed noise-dependent score network (modified U-Net architecture, hence not directly comparable to our results), we think that our work opens up a new research direction by bridging statistical decision theory (homogeneous proper scoring rules) and deep energy-based generative modeling. Moreover,
1For CelebA dataset, we reproduced the short-run MCMC method [66] using our code base. Moreover, the f -EBM paper only reported results on CelebA 32x32 and we empirically found it is not comparable to our method in CelebA 64x64 (higher resolution), indicating better scalability of PS-CD to high-dimensional case.
under the setting of simple model architectures and the same hyperparameter configuration (e.g., batch size, learning rate, network structure, etc.), our empirical results suggest clear superiority of PS-CD over traditional CD and recent f -EBMs.
OOD Detection & Robustness to Data Contamination. We further test our methods on out-ofdistribution (OOD) detection tasks. For the conditional CIFAR-10 model, we follow the evaluation protocol in [17] and use s(x) = maxy∈Y −E(x, y) as the score for detecting outliers. We use SVHN [65], Textures [9], Uniform/Gaussian Noise, CIFAR-10 Linear Interpolation and CelebA as the OOD datasets. We summarize the results in Table 5 in App. D.6, from which we can see that PS-CD consistently outperforms CD and other likelihood-based models.
Inspired by the OOD detection performance and previous work on robust parameter estimation under data contamination [43], we further test the robustness of CD and PS-CD on both synthetic and natural image datasets. Specifically, suppose p(x) is the underlying data distribution and there is another contamination distribution ω(x), e.g. uniform noise. In generative modeling under data contamination, our model observes i.i.d. samples from the contaminated distribution p̃(x) = cp(x) + (1 − c)ω(x), where 1 − c ∈ [0, 1/2) is the contamination ratio. A theoretical advantage of pseudo-spherical score is its robustness to data contamination: the optimal solution of minθ Sps(p̃, qθ) is close to that of minθ Sps(p, qθ) under some conditions (e.g. the density of ω(x) mostly lies in the region for which the target density p(x) is small) [20, 43]. From Table 2, we can see that CD suffers from data contamination severely: as the contamination ratio increases, the performance degrades drastically. By contrast, PS-CD shows good robustness against data contamination and a larger γ leads to better robustness. For example, PS-CD with γ = 1.0 can properly approximate the target distribution when the contamination ratio is 0.2, while PS-CD with γ = 2.0 can do so when the contamination ratio is up to 0.3.
We conduct similar experiments on MNIST and CIFAR-10 datasets, where we use uniform noise as the contamination distribution and the contamination ratio is 0.1 (i.e. 10% images in the training set are replaced with random noise). After a warm-up pretraining2 (when the model has some OOD detection ability), we train the model with the contaminated data and measure the training progress. We observe that CD gradually generates more random noise and diverge after a few training steps, while PS-CD is very robust. As shown in Table 6 in App. D.6, for a slightly pre-trained unconditional CIFAR-10 model (a simple 5-layer CNN with FID of 68.77), we observe that the performance of CD degrades drastically in terms of FID, while PS-CD can continuously improve the model even using the contaminated data. We provide visualizations and theoretical explanations in App. D.6.

heading: 7 Conclusion
text: From the perspective of maximizing strictly proper homogeneous scoring rules, we propose pseudospherical contrastive divergence (PS-CD) to generalize maximum likelihood estimation of energybased models. Different from prior works that involve joint training of auxiliary models [89, 21, 50, 12, 29], PS-CD allows us to specify flexible modeling preferences without additional computational cost compared to contrastive divergence. We provide a theoretical analysis on the sample complexity and convergence property of the proposed method, as well as its connection to maximum likelihood. Finally, we demonstrate the effectiveness of PS-CD with extensive experiments on both synthetic data and commonly used image datasets.
2Note that it is impossible for a randomly initialized model to be robust to data contamination since without additional inductive bias, it will simply treat the contaminated distribution p̃(x) as the target.

heading: Acknowledgements
text: This research was supported by NSF(#1651565, #1522054, #1733686), ONR (N000141912145), AFOSR (FA95501910024), ARO (W911NF-21-1-0125) and Sloan Fellowship.

Review summaries:
The authors propose a spectrum of objectives for learning energy-based models that generalise maximum likelihood. These objectives--known as pseudo-spherical scoring rules--are similar to (persistent) contrastive divergence in terms of their implementation and computational cost, but have an extra degree of flexibility that can boost downstream performance for e.g. image generation.

This paper provides an efficient framework to estimate energy based model (EBM).  While conventional KL or general f-divergence-based contrastive divergence learning suffers from instability and non-convergence issues, the proposed framework namely pseudo-spherical contrastive divergence (PS-CD) can be implemented as SGD type algorithm whose sample complexity and convergence property are provided.  PS-CD includes the conventional CD as a special case and thus provides gradient-based optimization framework for CD.  Experiments to evaluate the effectiveness of the proposed approach are conducted with both synthetic and image benchmark datasets. 

This paper presents a generalization of the Contrastive Divergence that uses pseudo-spherical scoring rules to define a new family of loss functions for generative modeling. This family of objective functions has an auxiliary loss parameter that can balance the trade-off between diversity and quality. The method is applied to learn EBMs on image benchmark datasets.

The paper introduces a class of divergence to learn EBMs. The class of divergence includes KL as a special case. The paper approximates the gradient of the class of divergence by self-normalized importance sampling and analyzes the convergence when optimizing with the approximate gradient. Experiments are conducted to evaluate the influence of $\gamma$.