Computer Science
Computer Science-2021.gz:143491
Computer Science
90bf69b7238d274c4de9459185c6a6d53302c595
Pseudo-Spherical Contrastive Divergence

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.


Introduction
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 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.

Energy-Based Distribution Representation and Sampling
Given a set of i.i.d. samples {x i } N i=1 from some unknown data distribution p(x) defined over the sample space X ⊂ R m , 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: 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 L 1 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: where z t ∼ N (0, I). The distribution ofx 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.

Maximum Likelihood Training of EBMs via Contrastive Divergence
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: 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 L MLE (θ; p) using samples from q θ : 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: where q K θ denotes the distribution after K steps of MCMC transitions from an initial distribution (typically data distribution or uniform distribution), and Equation (4) corresponds to L CD−∞ .

Strictly Proper Scoring Rules
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 1integrable 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: 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).
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: When a scoring rule S is strictly proper relative to P, the associated generalized entropy function and divergence function are defined as: 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): 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: , and the associated Bregman divergence is the KL divergence: . 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 θ ) + log Z θ = S(x, q θ ) and log Z θ 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.

Training EBMs by Maximizing Homogeneous Scoring Rules
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.

Pseudo-Spherical Scoring Rule
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: where q γ+1 := X q(y) γ+1 dy 1 γ+1 .
The expected pseudo-spherical score under a reference distribution p is defined as: Example 2. The classic spherical scoring rule [19] is a special case in the pseudo-spherical family, which corresponds to γ = 1: 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 = (q 1 , . . . , q n ). Let vector p = (p 1 , . . . , p n ) represent the oracle probability forecast. The expected spherical score can be written as: 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(−E 1 ), . . . , exp(−E n )), 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 S ps (p, q) defined in Equation (12) with γ > 0, the expected γ-score is defined as: Since 1 γ log(u) is strictly increasing in u, S γ (p, q) is a strictly proper homogeneous composite score:

Pseudo-Spherical Contrastive Divergence
Suppose we parametrize the energy-based model distribution as q θ ∝ q θ = exp(−E θ ) and we want to minimize the negative γ-score in Equation (15): 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: where the auxiliary distribution r θ is also an energy-based distribution defined as: 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) 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: . Define the gradient estimator as: where the self-normalized importance weight . Then the gradient estimator converges to the true gradient (Equation (18)) in probability, i.e., ∀ > 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).

Connections to Maximum Likelihood Estimation and Extension to γ < 0
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.

in [44]):
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: 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: where H γ+1 (q) is the Rényi entropy of order γ + 1.

Theoretical Analysis
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.

Sample Complexity
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: 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 L 2 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 complex- 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: we have:

Convergence
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): 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 N t 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 := αM L 2 e 4γK ): 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: Then with probability at least 1 − δ, we have (constant C := αM L 2 e 4γK ): In Appendix C.2, we discuss more on the strongly convex (Theorem 7) and convex cases (Theorem 8).

Related Work
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) = n x /n, where n x 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]).

Experiments
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.   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.   [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 1 1 , 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, Table 2: Robustness to data contamination on Gaussian datasets. The data distribution is N (−1, 0.5) and the contamination distribution is N (2, 0.05). We measure the KL divergence between clean target distribution p and converged model distribution q θ , D KL (p q θ ). 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) = max y∈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.

Contamination
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 distributioñ is the contamination ratio. A theoretical advantage of pseudo-spherical score is its robustness to data contamination: the optimal solution of min θ S ps (p, q θ ) is close to that of min θ S ps (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 pretraining 2 (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.

Conclusion
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.

A Implementation of PS-CD Algorithm
In the following, we provide a simple PyTorch [71] implementation of Step 6 in Algorithm 1, where we directly differentiate through log(mean(exp(·))) for the first term and perform stop gradient operation on the self-normalized importance weight (exp(−γE θ (x − i ))) for the second term. Mathematically, these two implementations are equivalent and we use different implementations here to illustrate the difference between their derivations. We further apply L 2 regularization on the outputs of the energy function to stabilize training.

B Proofs for Section 3 B.1 Proof for Theorem 2
In this section, we provide two different ways to prove Theorem 2. 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 for KL and f -divergences [89,11,12] and revealing the elegance of PS-CD. Theorem 2. For an energy-based distribution q θ ∝ q θ = exp(−E θ ), the gradient of the negative γ-score L γ (θ; p) = −S γ (p, q θ ) with respect to θ can be written as: where the auxiliary distribution r θ is also an energy-based distribution defined as: Proof. First proof: direct differentiation. From Equation (17), we have: Second proof: a variational representation with optimal variational distribution taking analytical form. The main challenge is that the log q θ γ+1 term in L γ (θ; p) is generally intractable to compute. To solve this issue, we introduce the following variational representation: Let ∆ X denote the set of all normalized probability density functions on sample space X . With Fenchel duality, we have: where the maximum is attained at Proof. With Jensen's inequality, we have: The equality holds if and only if: As r * ∈ ∆ X is a normalized distribution, we have: Now, suppose we parametrize the energy-based model distribution as q θ = exp(−E θ ) and the variational distribution as r ψ . By plugging the variational representation in Lemma 1 into γ-score (Equation (15)), we obtain the following minimax formulation to minimize the negative γ-score: where the game value function L γ (θ, ψ; p) is defined as (s.t. L γ (θ, ψ * (θ); p) = −S γ (p, q θ )): By Lemma 1, we know that r ψ * (θ) ∝ q γ+1 θ .
Thus we have:

B.2 Proof for Theorem 3
Theorem 3 (Consistent Gradient Estimation). Let x + 1 , . . . , x + N be i.i.d. samples from p(x) and . Define the gradient estimator as: where the self-normalized importance weight . Then the gradient estimator converges to the true gradient in probability: Proof. First, let us write ∇ θ L N γ (θ; p) and ∇ θ L γ (θ; p) as: For the first term in Equation (24), since {x + i } N i=1 are i.i.d. samples from p(x), by weak law of large numbers, the numerator and denominator of the first term in Equation (24) converges to the numerator and denominator of the first term in Equation (23) in probability. By Slutsky's theorem (i.e., for random variables X N , X, Y N , Y , if X N p → X, Y N p → Y and X, Y are constants, then X N /Y N p → X/Y ), the first term of Equation (24) converges to the first term of Equation (23) in probability.
Let us use Z r and Z q to denote the partition function for r θ and q θ . The second term of Equation (24) can be written as: , the numerator of Equation (25) converges to E r θ (x) [∇ θ E θ (x)] in probability, while the denominator of Equation (25) converges to 1 in probability (E q θ (x) [r θ (x)/q θ (x)] = 1). By Slutsky's theorem, the second term of Equation (24) converges to the second term of Equation (23) in probability. Furthermore, since convergence in probability is also preserved under addition transformation, the gradient estimator in Equation (24) converges to the true gradient in Equation (23) in probability.

B.3 Connections to Maximum Likelihood Estimation and Extension to γ < 0
Lemma 2. Let D γ (p, q) be the divergence corresponding to γ-scoring rule. Then, we have: lim γ→0 D γ (p, q) = D KL (p q) Proof. As introduced in Equation (9) in Section 2.3, the divergence corresponding to the γ-scoring rule is: When γ → 0, with Taylor series, we know that: Therefore, we have: The above lemma implies that the KL divergence minimization (maximum likelihood estimation) is a special case of γ-divergence minimization when γ → 0, which also implies the following corollary: Corollary 1. When γ → 0, the gradient of pseudo-spherical contrastive divergence is equal to the gradient of contrastive divergence: lim Proof. This is a direct consequence of Lemma 2. It can also be verified by checking the PS-CD gradient in Equation (18) (when γ → 0, r θ = q θ ∝ exp(−E θ )): Inspired by [86,56] that generalize Rényi divergence beyond its definition to negative orders, we consider the extension of γ-scoring rule with γ < 0 (although it is no longer strictly proper for these γ values) and show that maximizing such scoring rule is equivalent to maximizing a lower bound of logarithm scoring rule (MLE) with an additional Rényi entropy regularization. Lemma 3. When −1 ≤ γ < 0, we have: where H γ+1 (q) is the Rényi entropy of order γ + 1.
Proof. As a generalization to Shannon entropy, the Rényi entropy of order α is defined as: With Jensen's inequality, for −1 ≤ γ < 0, we have:

C Theoretical Analysis
In this section, we provide a theoretical analysis on the sample complexity of the gradient estimator, as well as the convergence property of stochastic gradient descent with consistent (but biased given finite samples) gradient estimators as presented in Algorithm 1.

C.1 Sample Complexity
We start with analyzing the sample complexity of the consistent gradient estimator, that is how fast it approaches the true gradient value or how many samples we need in order to empirically estimate the gradient at a given accuracy with a high probability.
We first make the following assumption, which is similar to the one used in [4,47]: Assumption 1. The energy function is bounded by K and the gradient is bounded by L (with K > 0 and L > 0): The assumption is typically easy to enforce in practice. For example, in the experiments we use L 2 regularization on the outputs of the energy function, as well as normalized inputs and spectral normalization [60] for the neural network that realizes the energy function. Theorem 4. Under Assumption 1, given any constants > 0 and δ ∈ (0, 1), when the number of samples N satisfies: N ≥ 32L 2 e 8γK (1 + 4 log(2/δ)) 2 we have: Proof. For notation simplicity, we use p N to denote the empirical distribution of First, we observe that: where the partition functions of r θ and q θ cancel out. Based on Equation (23) and (24), with triangle inequality, the estimation error can be upper bounded as: Now we introduce the following lemma that will provide us a probability upper bound that an empirical mean of independent random variables deviates from its expected value more than a certain amount.
Lemma 4 (Vector Bernstein Inequality [47,32]). Let X 1 , . . . , X N be independent vector-valued random variables. Assume that each one is centered, uniformly bounded and the variance is also bounded: ∀i, E[X i ] = 0 and X i ≤ µ and E[ X i 2 ] ≤ σ 2 Define X := 1 N (X 1 + . . . + X N ). Then we have for 0 < t < σ 2 /µ: We then define the following vector-valued random variable: From Equation (27), we know that: With σ 2 := 4L 2 e 4γK , from Lemma 4, we know that: To obtain a sample complexity bound such that the probability bound in Equation (29) is less than 1 − √ 1 − δ, we need to solve for N : Solving Equation (30) gives us: for δ ∈ (0, 1], we use the following slightly weaker bound such that it looks cleaner: Since Equation (28) is an upper bound of ∆ p , we know that when the sample size satisfies Equation (32), we have: Similarly, we can obtain the same sample complexity bound for ∆ q θ such that: From Equation (26), we know that ∆ p + ∆ q θ is an upper bound of the gradient estimation error. Also note that the event ∆ p ≤ 2 and the event ∆ q θ ≤ 2 are independent from each other (the samples for p N and the samples for q N θ are independent samples from p and q θ respectively). Thus when the sample size satisfies Equation (32), we have:

C.2 Convergence of Pseudo-Spherical Contrastive Divergence Algorithm
In this section, we analyze the convergence property of the PS-CD algorithm presented in Algorithm 1. For notation simplicity, we define g as the true gradient in Equation (18) and g as the gradient estimator in Equation (19). We further use L(θ) to denote the loss function L γ (θ, ψ * (θ); p) = −S γ (p, q θ ).
Let us consider the following stochastic gradient descent (SGD) update rule: where η t is the step size at step t, T is the total number of steps and g t := ∇ θ L N γ (θ; p) θ=θt is the consistent (but biased) gradient estimation of g t at step t. Note that θ t and g t are random variables that depend on the previous history g 1 , . . . , g t−1 . For brevity, in the following we will omit such dependency in the notations.
Most works for analyzing the convergence behavior of SGD relies on the assumption that the gradient estimator g t is asymptotically unbiased, e.g., [63,52,78,24,74], while in our case the gradient estimator is not unbiased but consistent (see Section 1.2 in [8] for a detailed discussion on the distinctions between unbiasedness and consistency). Therefore, in this work we generalize the theory developed in [24] and [8] to analyze the convergence rate for PS-CD.
Besides Assumption 1 used for analyzing the sample complexity of the gradient estimator, we further make the following assumption: Assumption 2. The loss function L(θ) is M -smooth (with M > 0): This is a common assumption used for analyzing first-order optimization methods, which is also used in [24,8]. Also note that 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. Now, let us consider Algorithm 2, which is a variant of SGD that allows early stopping before reaching the iteration limit T according to some probability distribution p Z over iteration indexes Obtain the gradient estimator g t with a sample batch size of N t .
Note that this is equivalent (more efficient in terms of computation) to running the algorithm to the iteration limit T and then selecting the final solution from {θ 1 , . . . , θ T } according to distribution p Z .
We have the following theorem that characterizes the convergence property of Algorithm 2: Theorem 5. Under Assumptions 1 and 2, for arbitrary constants α ∈ (0, 1) and δ ∈ (0, 1), suppose that the step sizes satisfy η t < 2(1 − α)/M and the probability distribution over iteration indexes is chosen to be: , t = 1, . . . , T and the sample size N t used for estimating g t satisfies: Denote by L * the minimum value of L(θ). Then with probability at least 1 − δ, we have: Proof. First, with Assumption 1 (|E θ (x)| ≤ K, ∇ θ E θ (x) ≤ L), we can bound the norm of the true gradient as: From Theorem 4, we know that when sample size at each step satisfies Equation (36), we have: Thus, with probability at least 1 − δ, we have: where the last inequality is due to the fact that L * ≤ L(θ T +1 ).

Dividing both sides by
By definition of p Z in Equation (35), which is used to select a final solution among {θ 1 , . . . , θ T }, we know that: Now, let us consider a simple case where we use a constant step size, which gives us the following corollary: Corollary 2. Under the conditions in Theorem 5 except that we use constant step sizes: then with probability at least 1 − δ, we have: Proof. Since we are using a constant step size, by Theorem 5, we know that: By Equation (40), we have: Note that an alternative result for constant step sizes can also be obtained from Theorem 6 in [8]: Theorem 6. Under Assumptions 1 and 2, for arbitrary constants α ∈ (0, 1) and δ ∈ (0, 1), suppose we use constant step sizes: and the sample size N t used for estimating g t satisfies Equation (36), then with probability at least 1 − δ, we have: Although both Corollary 2 and Theorem 6 give a convergence rate of O(1/ √ T ), the strategy in Theorem 6 requires extra computational effort to compute ∇L(θ t ) for t = 1, . . . , T in order to select the solution with the minimum gradient norm. Since ∇L(θ) cannot be computed exactly, Monte Carlo estimation will incur additional approximation error. By contrast, the strategy in our analysis does not have such issues and Theorem 5 provides a general analysis on the convergence rate of the randomized SGD algorithm with consistent but biased gradient estimators, which allows using different step sizes.
For example, starting from Equation (37) and with the fact that: Assumption 4. The loss function L(θ) is convex and the parameter space has finite diameter D: sup θ1,θ2∈Θ θ 1 = D. Let θ * ∈ arg min θ∈Θ L(θ). Theorem 8. Under Assumptions 1 and 4, for arbitrary constant δ ∈ (0, 1), suppose we use decreasing step sizes η t = 1/ √ t for t = 1, . . . , T and the sample size N t used for estimating g t satisfies: then with probability at least 1 − δ, we have: The above theorems follow from the sample complexity bound in Theorem 4 and the results in [8] (Theorem 2 with constant ρ = J/2 and Theorem 5 with constants ρ = 1, c = 1), which we refer to for a detailed proof. Note that the condition on J in Theorem 7 is optional and without the condition, we can obtain the same convergence rate at the cost of a cumbersome form in the R.

D.2 Understanding the Effects of Different γ Values in 1-D Examples
In this section, we aim to provide insights on the effects of different γ values with 1-D toy experiments. Specifically, we use an EBM with a quadratic energy function (corresponding to a Gaussian distribution): where µ and σ are two trainable parameters.
First, we show that when the real data distribution is also a Gaussian distribution such that the model is well-specified, then different γ values will induce the same optimal distribution since they are strictly proper. To verify this property, we visualize the objective landscape in Figure 3.
Second, we use the same quadratic energy function to fit a mixture of Gaussians. We visualize the objective landscapes in Figure 4, which shows that when the model is mis-specified, different objectives will exhibit different modeling preferences (inducing different solutions). This corresponds to the practical scenarios, where such property enables us to flexibly specify different inductive biases to make tradeoff among various modeling factors such as diversity/quality.

D.3 Understanding the Effects of Different γ Values in Image Generation
Although FID has been the most popular evaluation metric for image generative models, it is problematic since it summarizes the difference between two distributions into a single number and fails to separate important aspects such as fidelity and diversity [61]. To better demonstrate the modeling flexibility brought by the proposed PS-CD framework, we conduct experiments on CIFAR-10 dataset using a set of more indicative and reliable metrics proposed by [61] to evaluate the effects of γ from various perspectives. Table 4: Effects of γ on CIFAR-10 image generation. We use the same image embeddings (activations of a pre-trained inception network) to compute these metrics and FID to ensure consistency. We briefly introduce these metrics here and refer to [61] for accurate descriptions and mathematical definitions. Denote data distribution as P (X) and model distribution as Q(X). Based on manifold estimation, Precision is defined as the portion of Q(X) that can be generated by P (X) and Recall is symmetrically defined as the portion of P (X) that can be generated by Q(X); Density improves upon Precision to count how many real-sample neighbourhood spheres contain a certain fake sample; Coverage improves upon Recall to measure the fraction of real samples whose neighbourhoods contain at least one fake sample. From Table 4, we have some interesting observations: (1) PS-CD with γ = −0.5 and γ = 1.0 get best FID scores because they can simultaneously achieve good balance among these metrics (e.g., high Density and Coverage score); (2) By contrast, PS-CD with γ = 2.0 achieves the highest Density score but a relatively low Coverage score, which potentially leads to a slightly worse FID; (3) Many members in the PS-CD family showed superior performance over traditional contrastive divergence in most metrics, demonstrating the potential of our method. Just like various f -divergences used in generative modeling, different members in the PS-CD family can represent complicated inductive bias in practice (although being strictly proper in model well-specified case). Since these single-valued evaluation metrics measure the generative performance in a complicated way, we think it is normal that the change of γ is not monotone to the change of each metric. For specific application scenarios, we may mainly care about a certain metric and we should choose γ accordingly.

Density Coverage Precision
We would like to emphasize that, a major contribution of our paper is opening the door to a new family of EBM training objectives and enabling us to flexibly specify modeling preferences, without introducing additional computational cost compared to CD (unlike adversarial training in f -EBM).   To keep a fair comparison, all the compared methods use the same architecture to implement the energy function, except that f -EBMs require an additional variational function that uses the same architecture as the energy function. The architectures used for CIFAR-10 (32 × 32) and CelebA (64 × 64) datasets are shown in Figure 6. We use leaky-ReLU non-linearity with default leaky factor 0.2 throughout the architectures (between all the convolution layers). Following [17,89], we apply spectral normalization/L 2 regularization (on the outputs of the models) with coefficient 1.0 to improve the stability.

D.4 Image Generation Samples for PS-CD
For CIFAR-10, to keep a fair comparison, we use the same sampling strategy as [17,89], where a sample replay buffer is employed to improve the mixing of Langevin dynamics. Specifically, we use 60 steps Langevin dynamics together with a sample replay buffer of size 10000 to produce samples in the training phase. In each Langevin step, we use a step size of 10.0 and a random noise with standard deviation of 0.005.
For CelebA, which has a higher data dimension, we use the sampling strategy in [66] to improve the efficiency of sampling, where we always start the Markov chains from a fixed uniform distribution and run a fixed number of Langevin steps (100) with a constant step size.
For all the experiments, we use Adam optimizer to optimize the parameters of the energy function. In each training iteration, we use a batch size of 128 for CIFAR-10 and 64 for CelebA. We run the PS-CD algorithms for about 50K iterations of parameter updates for CIFAR-10 and about 100K iterations for CelebA.
For computational cost, the CIFAR-10 experiments take about 48 hours on 4 Titan Xp GPUs, while the CelebA experiments take about 16 hours since we learn non-convergent short-run MCMC.  To show the practical advantage of PS-CD in face of data contamination, we further conduct experiments on MNIST and CIFAR-10 datasets, where we use random 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 pretraining (when the model has some OOD detection ability), we train the model with the contaminated data and measure the training progress of CD and PS-CD.

D.6 OOD Detection & Robustness to Data Contamination
As shown in Figure 7, CD gradually generates more and more random noise and diverge after a few training steps, while PS-CD is very robust. In particular, as shown in Table 6, for a slightly pretrained 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 believe that robustness to data contamination is a valuable property for modern deep generative models and there is actually a natural interpretation for the robustness of PS-CD. Compared to CD, there is an extra weight term before the gradient of the energy: exp(−γE θ (xi)) j exp(−γE θ (xj ) ∇ θ E θ (x i ) (the first term in Eq. (19)). Suppose x i ∼ ω is a noise data from the contaminated distributionp in a batch of samples, for a model with OOD detection ability, it will assign a much higher energy to x i than normal data and the weight before ∇ θ E θ (x i ) will be close to zero. In short, PS-CD naturally integrates the OOD detection ability of EBMs into the training process, which then leads to robustness to data contamination.