On Energy-Based Models with Overparametrized Shallow Neural Networks

Energy-based models (EBMs) are a simple yet powerful framework for generative modeling. They are based on a trainable energy function which defines an associated Gibbs measure, and they can be trained and sampled from via well-established statistical tools, such as MCMC. Neural networks may be used as energy function approximators, providing both a rich class of expressive models as well as a flexible device to incorporate data structure. In this work we focus on shallow neural networks. Building from the incipient theory of overparametrized neural networks, we show that models trained in the so-called"active"regime provide a statistical advantage over their associated"lazy"or kernel regime, leading to improved adaptivity to hidden low-dimensional structure in the data distribution, as already observed in supervised learning. Our study covers both maximum likelihood and Stein Discrepancy estimators, and we validate our theoretical results with numerical experiments on synthetic data.


Introduction
A central problem in machine learning is to learn generative models of a distribution through its samples. Such models may be needed simply as a modeling tool in order to discover properties of the data, or as a way to generate new samples that are similar to the training samples. Generative models come in various flavors. In some cases very few assumptions are made on the distribution and one simply tries to learn generator models in a black-box fashion [Goodfellow et al., 2014, Kingma andWelling, 2013], while other approaches make more precise assumptions on the form of the data distribution. In this paper, we focus on the latter approach, by considering Gibbs measures defined through an energy function f , with a density proportional to exp{−f (x)}. Such energy-based models (EBMs) originate in statistical physics [Ruelle, 1969], and have become a fundamental modeling tool in statistics and machine learning [Wainwright and Jordan, 2008, Ranzato et al., 2007, LeCun et al., 2006, Du and Mordatch, 2019, Song and Kingma, 2021. If data is assumed to come from such a model, the learning algorithms then attempt to estimate the energy function f . The resulting learned model can then be used to obtain new samples, typically through Markov Chain Monte Carlo (MCMC) techniques.
In this paper, we study the statistical problem of learning such EBMs from data, in a non-parametric setting defined by a function class F, and with possibly arbitrary target energy functions. If we only assume a simple Lipschitz property on the energy, learning such models will generally suffer from the curse of dimensionality [von Luxburg and Bousquet, 2004], in the sense that an exponential number of samples in the dimension is needed to find a good model. However, one may hope to achieve better guarantees when additional structure is present in the energy function.
The standard method to train EBMs is maximum likelihood estimation. One generic approach for this is to use gradient descent, where gradients may be approximated using MCMC samples from the current trained model. Such sampling procedures may be difficult in general, particularly for complex energy landscapes, thus we also consider different estimators based on un-normalized measures which avoid the need of sampling. We focus here on approaches based on minimizing Stein discrepancies Mackey, 2015, Liu andWang, 2016], which have recently been found to be useful in deep generative models [Grathwohl et al., 2020], though we note that alternative approaches may be used, such as score matching [Hyvärinen, 2005, Song and Kingma, 2021, Song and Ermon, 2019, Block et al., 2020.
Our main focus is to study the resulting estimators when using gradient-based optimization over infinitelywide neural networks in different regimes, showing the statistical benefits of the 'feature learning' regime when the target models have low-dimensional structure, thus extending the analogous results for supervised least-squares [Bach, 2017a] and logistic [Chizat and Bach, 2020] regression. More precisely, we make the following contributions: • We derive generalization bounds for the learned measures in terms of the same metrics used for training (KL divergence or Stein discrepancies). Using and extending results from the theory of overparametrized neural networks, we show that when using energies in the class F 1 we can learn target measures with certain low-dimensional structure at a rate controlled by the intrinsic dimension rather than the ambient dimension (Corollary 1 and Corollary 2).
• We show in experiments that while F 1 energies succeed in learning simple synthetic distributions with low-dimensional structure, F 2 energies fail (Sec. 6).

Related work
A recent line of research has studied the question of how neural networks compare to kernel methods, with a focus on supervised learning problems. Bach [2017a] studies two function classes that arise from infinite-width neural networks with different norms penalties on its weights, leading to the two different spaces F 1 and F 2 , and shows the approximation benefits of the F 1 space for adapting to low-dimensional structures compared to the (kernel) space F 2 , an analysis that we leverage in our work. The function space F 1 was also studied by Ongie et al. [2019], Savarese et al. [2019], Williams et al. [2019] by focusing on the ReLU activation function. More recently, this question has gained interest after several works have shown that wide neural networks trained with gradient methods may behave like kernel methods in certain regimes [see, e.g., Jacot et al., 2018]. Examples of works that compare 'active/feature learning' and 'kernel/lazy' regimes include [Chizat and Bach, 2020, Ghorbani et al., 2019, Wei et al., 2020, Woodworth et al., 2020. We are not aware of any works that study questions related to this in the context of generative models in general and EBMs in particular.
Other related work includes the Stein discrepancy literature. Although Stein's method [Stein, 1972] dates to the 1970s, it has been popular in machine learning in recent years. Gorham and Mackey [2015] introduced a computational approach to compute the Stein discrepancy in order to assess sample quality. Later, Chwialkowski et al. [2016] and  introduced the more practical kernelized Stein discrepancy (KSD) for goodness-of-fit tests, which was also studied by Gorham and Mackey [2017]. Liu and Wang [2016] introduced SVGD, which was the first method to use the KSD to obtain samples from a distribution, and Barp et al. [2019] where the first to employ KSD to train parametric generative models. More recently, Grathwohl et al. [2020] used neural networks as test functions for Stein discrepancies, which arguably yields a stronger metric, and have shown how to leverage such metrics for training EBMs. The empirical success of their method provides an additional motivation for our theoretical study of the F 1 Stein Discrepancy (Subsec. 4.2).
Finally, another notable paper close in spirit to our goal is Block et al. [2020], which provides a detailed theoretical analysis of a score-matching generative model using Denoising Autoencoders followed by Langevin diffusion. While their work makes generally weaker assumptions and also includes a non-asymptotic analysis of the sampling algorithm, the resulting rates are unsuprisingly cursed by dimension. Our focus is on the statistical aspects which allow faster rates, leaving the quantitative computational aspects aside.
The estimated distribution is simply νf , and samples can be obtained by the MCMC algorithm of choice.
An alternative estimator is the one that arises from minimizing the Stein discrepancy (SD) corresponding to a function class H. If H is a class of functions from K to R d+1 , the Stein discrepancy Mackey, 2015, Liu et al., 2016] for H is a non-symmetric functional defined on pairs of probability measures over K as where A ν : K → R (d+1)×(d+1) is the Stein operator. In order to leverage approximation properties on the sphere, we will consider functions h defined on K = S d . In this case, the Stein operator is defined by A ν h(x) := (s ν (x) − d · x)h(x) + ∇h(x) (see Lemma 5), where s ν (x) = ∇ log( dν dτ (x)) is named the score function. The term d · x is important for the spherical case in order to have SD H (ν, ν) = 0, while it does not appear when considering K = R d . The Stein discrepancy estimator iŝ H0 ≤ 1} for some reproducing kernel Hilbert space (RKHS) H 0 with kernel k with continuous second order partial derivatives, there exists a closed form for the problem (2) and the corresponding object is known as kernelized Stein discrepancy (KSD) [Liu et al., 2016, Gorham andMackey, 2017]. For K = S d , the KSD takes the following form (Lemma 6): x )), and we useũ ν (x, x ) to denote the sum of the first three terms (remark that the fourth term does not depend on ν). One KSD estimator that can be used iŝ The optimization problem for this estimator is convex (Sec. 5), but it is biased. On the other hand, the estimatorf is unbiased, but the optimization problem is not convex.

Neural network energy classes
We are interested in the cases in which F is one of two classes of functions related to shallow neural networks, as studied by Bach [2017a].
Feature learning regime. F is the ball B F1 (β) of radius β > 0 of F 1 , which is the Banach space of Kernel regime. F is the ball B F2 (β) of radius β > 0 of F 2 , which is the (reproducing kernel) Hilbert space of functions f : K → R such that for some absolutely continuous ρ ∈ M(S d ) with dρ dτ ∈ L 2 (S d ) (wherẽ τ the uniform probability measure over S d ), we have that for all Remark that since |h(θ)|dτ (θ) ≤ ( |h(θ)| 2 dτ (θ)) 1/2 by the Cauchy-Schwarz inequality, we have F 2 ⊂ F 1 and B F2 ⊂ B F1 . The TV norm in F 1 acts as a sparsity-promoting penalty, which encourages the selection of few well-chosen neurons and may lead to favorable adaptivity properties when the target has a low-dimensional structure. In particular, Bach [2017a] shows that single ReLU units belong to F 1 but not to F 2 , and their L 2 approximations in F 2 have exponentially high norm in the dimension. Ever since, several works have further studied the gaps arising between such nonlinear and linear regimes [Wei et al., 2019, Ghorbani et al., 2020, Malach et al., 2021. In App. D, we present dual characterizations of the maximum likelihood F 1 and F 2 EBMs as entropy maximizers under L ∞ and L 2 moment constraints (an infinite-dimensional analogue of Della Pietra et al. [1997]; see also Mohri et al. [2012], Theorem 12.2).
The ball radius β acts as an inverse temperature. The low temperature regime β 1 corresponds to expressive models with lower approximation error but higher statistical error: the theorems in Sec. 4 provide bounds on the two errors and the results of optimizing such bounds w.r.t. β. In the following, we will assume that the set K ⊂ R d+1 is compact. We note that there are two interesting choices for K: (i) for K = S d , we obtain neural networks without bias term; and (ii) for K = K 0 × {R}, where K 0 ⊂ R d with norm bounded by R, we obtain neural networks on K 0 with a bias term.

Statistical guarantees for shallow neural network EBMs
In this section, we present our statistical generalization bounds for various EBM estimators based on maximum likelihood and Stein discrepancies, highligting the adaptivity to low-dimensional structures that can be achieved when learning with energies in F 1 . All the proofs are in App. A.

Guarantees for maximum likelihood EBMs
The following theorem provides a bound of the KL divergence between the target probability measure and the maximum likelihood estimator in terms of a statistical error and an approximation error.
Theorem 1. Assume that the class F has a (distribution-free) Rademacher complexity bound R n (F) ≤ βC √ n and L ∞ norm uniformly bounded by β. Given n samples {x i } n i=1 from the target measure ν, consider the maximum likelihood estimator (MLE)ν := νf , wheref is the estimator defined in (1). With probability at least 1 − δ, we have If dν dτ (x) = e −g(x) / K e −g(y) dτ (y) for some g : K → R, i.e. −g is the log-density of ν up to a constant term, then with probability at least 1 − δ, Equation (7) follows from using a classical argument in statistical learning theory. To obtain equation (8) we bound the last term of (7) by 2 inf f ∈F g − f ∞ using Lemma 1 in App. A. We note that other metrics than L ∞ may be used for the approximation error, such as the Fisher divergence, but these will likely lead to similar guarantees under our assumptions. Making use of the bounds developed by Bach [2017a], Corollary 1 below applies (8) to the case in which F is the F 1 ball B F1 (β) for some β > 0 and the energy of the target distribution is a sum of Lipschitz functions of orthogonal projection to low-dimensional subspaces.
Assumption 1. The target probability measure ν is absolutely continuous w.r.t. the uniform probability measure τ over K and it satisfies ∀x , where ϕ j are (ηR −1 )-Lipschitz continuous functions on the R-ball of R k such that ϕ j ∞ ≤ η, and U j ∈ R k×d with orthonormal rows.
Assume that Assumption 1 holds. Then, we can choose β > 0 such that with probability at least 1 − δ we have where the notationÕ indicates that we overlook logarithmic factors and constants depending only on the dimension k.
Remarkably, Corollary 1 shows that for our class of target measures with low-dimensional structure, the KL divergence between ν andν decreases as n − 1 k+3 . That is, the rate "breaks" the curse of dimensionality since the exponent only depends on the dimension k of the low-dimensional spaces, not to the ambient dimension d. This can be seen as an alternative, more structural approach to alleviate dimension-dependence compared to other standard assumptions such smoothness classes for density estimation [e.g., Singh et al., 2018, Tsybakov, 2008. As discussed earlier, a motivation for Assumption 1 comes from Markov Random Fields, where each ϕ j corresponds to a local potential defined on a neighborhood determined by U j . Note that the bound scales linearly with respect to the number of local potentials J. As our experiments illustrate (see Sec. 6), it is easy to construct target energies that are much better approximated in F 1 than in F 2 . Indeed, we find that the test error tends to decrease more quickly as a function of the sample size when training both layers of shallow networks rather than just the second layer, which corresponds to controlling the F 1 norm.

Guarantees for Stein Discrepancy EBMs
We now consider EBM estimators obtained by minimizing Stein discrepancies, and establish bounds on the Stein discrepancies between the target measure and the estimated one. As in Subsec. 4.1, we begin by providing error decompositions in terms of estimation and approximation error. The following theorem applies to the Stein discrepancy estimator when the set of test functions H is the unit ball of the space of F d+1 in a mixed F/ 2 norm, with F = F 1 or F 2 . For F 1 , we will denote this particular setting as F 1 -Stein discrepancy, or F 1 -SD. Although F 1 -SD has not been studied before to our knowledge, the empirical work of Grathwohl et al. [2020] does use Stein discrepancies with neural network test functions, which provides practical motivation for considering such a metric.
≤ 1}, we have that for the estimatorν defined in (3), with probability at least 1 − δ, where C 2 is a universal constant and ∇f denotes the Riemannian gradient of f .
Notice that unlike in Theorem 1, the statistical error terms in Theorem 2 depend on the ambient dimension d. While we do not show that this dependence is necessary, studying this question would be an interesting future direction. Remark as well the similarity of the approximation term with the term 2 inf f ∈F g − f ∞ from equation (8), albeit in this case it involves the L ∞ norm of the gradients. Furthermore, note that the only assumption on the set F is a uniform L ∞ bound on F 1 , while Theorem 1 also requires a more restrictive Rademacher complexity bound on F. This illustrates the fact that the Stein discrepancy is a weaker metric than the KL divergence.
In Theorem 3 we give an analogous result for the unbiased KSD estimator (6), under the following reasonable assumptions on the kernel k, which follow .
Assumption 2. The kernel k has continuous second order partial derivatives, and it satisfies that for any non-zero function g ∈ L 2 (S d ), Theorem 3. Let K = S d . Assume that the class F is such that sup f ∈F { ∇f ∞ } ≤ βC 1 . Let KSD be the kernelized Stein discrepancy for a kernel that satisfies Assumption 2. If we take n samples {x i } n i=1 of a target measure ν with almost everywhere differentiable log-density, and consider the unbiased KSD estimator (6), we have with probability at least 1 − δ, The statistical error term in Theorem 3 is obtained using the expression of the variance of the estimator (6) . Note that Assumption 2 is fulfilled, for example, for the radial basis function (RBF) kernel where the notationÕ indicates that we overlook logarithmic factors and constants depending only on the dimension. (ii) Whenν is the unbiased KSD estimator (6) and the assumptions of Theorem 3 hold, β > 0 can be chosen so that with probability at least 1 − δ we have that KSD(ν,ν) is upper-bounded bỹ .
Noticeably, the rates in Corollary 2 are also of the form O(n − 1 k+3 ), which means that just as in Corollary 1, the low-dimensional structure in the target measure helps in breaking the curse of dimensionality.
Proof sketch. The main challenge in the proof of Corollary 2 is to bound the approximation terms in Theorem 2 and Theorem 3. To do so, we rely on Lemma 7 in App. A, which shows the existence ofĝ in a ball of F 2 such that sup x∈S d ∇ĝ(x)−∇g(x) 2 has a certain bound when g is bounded and has bounded and Lipschitz gradient. Lemma 7 might be of independent interest: in particular, it can be used to obtain a similar adaptivity result for score-matching EBMs, which optimize the Fisher divergence E x∼ν [ ∇ log( dν dp (x)) − ∇f (x) 2 ].

Algorithms
This section provides a description of the optimization algorithms used for learning F 1/2 -EBMs using the estimators studied in Sec. 4, namely maximum likelihood, KSD, and F 1 -SD.

Algorithms for F 1 EBMs
We provide the algorithms for the three models using a common framework. We define the function Φ : R × R d+1 → F 1 as Φ(w, θ)(x) = wσ( θ, x ). Given a convex loss R : F 1 → R, we consider the problem for some λ > 0. It is known [e.g., Neyshabur et al., 2015] that, since |w| 2 + θ 2 2 ≥ 2|w| θ 2 with equality when moduli are equal, this problem is equivalent to And by the definition of the F 1 norm, this is equivalent to inf f ∈F1 R (f ) + λ f F1 , which is the penalized form of inf f ∈B F 1 (β) R (f ) for some β > 0. Our F 1 EBM algorithms solve problems of the form (9) for different choices of R, or equivalently, minimize the functional R over an F 1 ball. The functional R takes the following forms for the three models considered: (i) Cross-entropy: We have that R(f ) = 1 , which is convex (and differentiable) because the free energy obeys such properties [e.g., by adapting Wainwright and Jordan, 2008, Prop 3.1 to the infinite-dimensional case].
(ii) Stein discrepancy: the estimator (5) corresponds to , which is convex as the supremum of convex (linear) functions.
(iii) Kernelized Stein discrepancy: we have R(f ) = 1 n 2 n i,j=1ũ ν f (x i , x j ), which is convex (in fact, it is quadratic in ∇f ).
In order to optimize (9), we discretize measures in P(R d+2 ) as averages of point masses 1 m m i=1 δ (w (i) ,θ (i) ) , each point mass corresponding to one neuron. Furthermore, we define the function G : Then, as outlined in Algorithm 1, we use gradient descent on G to optimize the parameters of the neurons, albeit possibly with noisy estimates of the gradients.
Algorithm 1 Generic algorithm to train F 1 EBMs Computing an estimate the gradient of G involves computing the gradient of R 1 m m i=1 Φ(w (i) , θ (i) ) . Denoting by z i = (w (i) , θ (i) ), z = (z i ) m i=1 and by ν z the Gibbs measure corresponding to the energy The expectation under ν z is estimated using MCMC samples of the EBM. Thus, the quality of gradient estimation depends on the performance of the MCMC method of choice, which can suffer for non-convex energies and low temperatures.
The gradient estimation involves d + 1 optimization procedures over balls of F 1 to compute h j , which we solve using Algorithm 1. Thus, the algorithm operates on two timescales.
(iii) Kernelized Stein discrepancy: Using (4), the gradient of R(f z ) with respect to z i takes the expression E x,x ∼νn [∇ zi u νz (x, x )], which can be developed into closed form. The only issue is the quadratic dependence on the number of samples.

Algorithms for F 2 EBMs
Considering convex losses R : F 1 → R as in Subsec. 5.1, the penalized form of the problem To optimize this, we discretize the problem: we take m samples (θ (i) ) m i=1 of the uniform measure τ that we keep fixed, and then solve the random features problem Remark that this objective function is equivalent to the objective function G((w (i) , θ (i) ) m i=1 ) in equation (10) when (θ (i) ) m i=1 are kept fixed. Thus, we can solve (11) by running Algorithm 1 without performing gradient descent updates on (θ (i) ) m i=1 . That is, while for the F 1 EBM training both the features and the weights are learned via gradient descent, for F 2 only the weights are learned.

Qualitative convergence results
The overparametrized regime corresponds to taking a large number of neurons m. In the limit m → ∞, under appropriate assumptions the empirical measure dynamics corresponding to the gradient flow of G((w (i) , θ (i) ) m i=1 ) converge weakly to the mean-field dynamics Mei et al. [2018], Chizat and Bach [2018], Rotskoff and Vanden-Eijnden [2018]. Leveraging a result from Chizat and Bach [2018] we argue informally that in the limit m → ∞, t → ∞, with continuous time and exact gradients, the gradient flow of G converges to the global optimum of F over P(R d+2 ) (see more details in App. B).
In contrast with this positive qualitative result, we should mention a computational aspect that distinguishes these algorithms from their supervised learning counterparts: the Gibbs sampling required to estimate the gradient at each timestep. A notorious challenge is that for generic energies (even generic energies in F 1 ), either the mixing time of MCMC algorithms is cursed by dimension Bakry et al. [2014] or the acceptance rate is exponentially small. The analysis of the extra assumptions on the target energy and initial conditions that would avoid such curse are beyond the scope of this work, but a framework based on thermodynamic integration and replica exchange [Swendsen and Wang, 1986] would be a possible route forward.
Experimental setup. We generate data on the sphere S d from teacher models by using a simple rejection sampling strategy, given an estimate of the minimum of f * (which provides an estimated upper bound on the unnormalized density e −f * for rejection sampling). This minimum is estimated using gradient descent with many random restarts from uniform points on the sphere. For different numbers of training samples, we run our gradient-based algorithms in F 1 and F 2 with different choices of step-sizes and regularization parameters λ, using m = 500 neurons. We report test metrics after selecting hyperparameters on a validation set of 2000 samples. For computing gradients in maximum likelihood training, we use a simple Metropolis-Hastings   Figure 3: Test metrics obtained for MLE, KSD and F 1 -SD training on a four-neuron teacher with weights w * 1 , w * 2 = 7.5 and w * 3 , w * 4 = −7.5. (top) Test performance measured with cross-entropy estimates. (bottom) Test KSD and F 1 -SD for models trained with the same metric. For reference, the black discontinuous lines show the teacher KSD and F 1 -SD of the teacher model w.r.t. 5000 and 2000 test samples, respectively. Confidence estimates are over 10 different data samplings.
algorithm with uniform proposals on the sphere. To obtain non-negative test KL divergence estimates, which are needed for the log-log plots, we sample large numbers of points uniformly on the hypersphere, and compute the KL divergence of the restriction of the EBMs to these points. The sampling techniques that we use are effective for the toy problems considered, but more refined techniques might be needed for more complex problems in higher dimension or lower temperatures.
Learning planted neuron distributions in hyperspheres. We consider the task of learning planted neuron distributions in d = 15 and d = 10. Remark that in this setting, when F = B F1 (β) with β large enough there is no approximation error. We compare the behavior of F 1 and F 2 models with different estimators in Figures 1, 2 and 3, corresponding to models with J = 1, 2, 4 teacher neurons, respectively. The error bars show the average and standard deviation for 10 runs. In the three figures, the top plot in the first column represents the test KL divergence of the F 1 and F 2 EBMs trained with maximum likelihood for an increasing number of samples, showcasing the adaptivity of F 1 to distributions with low-dimensional structure versus the struggle of the F 2 model. In Figures 1 and 2 the bottom plot in the first column shows the same information for a teacher with the same structure but different values for the output weights. We observe that the separation between the F 1 and the F 2 models increases when the teacher models have higher weights.
In the three figures, the plots in the second column show the test KL divergence and test KSD, respectively, for EBMs trained with KSD (with RBF kernel with σ 2 = 1). We observe that we are able to train EBMs successfully by optimizing the KSD; even though maximum likelihood training is directly optimizing the KL divergence, the test KL divergence values we obtain for the KSD-trained models are on par, or even slightly better, comparing at equal values of n. It is also worth noticing that in Figure 1, we observe a separation between F 1 and F 2 in the KL divergence plot, but not in the KSD plot. It seems that in this particular instance, although the training is successful, the KSD is too weak of a metric to tell that the F 1 EBMs are better than F 2 EBMs.
In the three figures, the plots in the third column show the test KL divergences and test F 1 -SD for EBMs  trained with F 1 -SD. Remark that the error bars are wider due to the two timescale algorithm used for F 1 -SD, which seems to introduce more variability. While the plots only go up to n = 3000, the test cross-entropy curves show a separation between F 1 and F 2 very similar to maximum likelihood training when comparing at equal values of n.
App. C contains additional experiments for the cases J = 1, w * 1 = 10 and J = 2, w * i = −2.5, training with KSD and F 1 -SD. Figure 4 shows a 3D visualization of the teacher and trained models, energies and densities corresponding to two teacher neurons with negative weights in d = 3. Since the dimension is small and the temperature is not too small, we used train and test sizes for which the statistical error due to train and test samples is negligible. Interestingly, while the F 1 model achieves a KL divergence close to zero at the end of training ( Figure 5), in Figure 4 we see that the positions of the neurons of the trained model do not match the teacher neurons. In fact, there are some neurons with positive weights in the high energy region. This effect might be linked with the fact that there is a constant offset of around 0.3 between the teacher energy and the trained energy. The offset is not reflected in the Gibbs measures of the models, which are invariant to constant terms.

3D visualizations and time evolution in
Figure 5 also shows that for this particular instance, the convergence is polynomial in the iteration number. We attach a video of the training dynamics: https://github.com/CDEnrich/ebms_shallow_nn/blob/main/ KLfeatzunnorm1.mp4.

Conclusions and discussion
We provide statistical error bounds for EBMs trained with KL divergence or Stein discrepancies, and show benefits of using energy models with infinite-width shallow networks in in "active" regimes in terms of adaptivity to distributions with low-dimensional structure in the energy. We empirically verify that networks in "kernel" regimes perform significantly worse in the presence of such structures, on simple teacher-student experiments.
A theoretical separation result in KL divergence or SD between F 1 and F 2 EBMs remains an important open question: one major difficulty for providing a lower bound on the performance for F 2 is that L 2 (or L ∞ ) approximation may be not be appropriate for capturing the hardness the problem, since two log-densities differing substantially in low energy regions can have arbitrarily small KL divergence. Another direction for future work is to apply the theory of shallow overparametrized neural networks to other generative models such as GANs or normalizing flows.
On the computational side, in App. B we leverage existing work to state qualitative convergence results in an idealized setting of infinite width and exact gradients, but it would be interesting to develop convergence results for maximum likelihood that take the MCMC sampling into account, as done for instance by Bortoli et al. [2020] for certain exponential family models. In our setting, this would entail identifying a computationally tractable subset of F 1 energies. A more ambitious and long-term goal is to instead move beyond the MCMC paradigm, and devise efficient sampling strategies that can operate outside the class of log-concave densities, as for instance Gabrié   A Proofs of Sec. 4 Theorem 1. Assume that the class F has a (distribution-free) Rademacher complexity bound R n (F) ≤ βC √ n and L ∞ norm uniformly bounded by β. Given n samples {x i } n i=1 from the target measure ν, consider the maximum likelihood estimator (MLE)ν := νf , wheref is the estimator defined in (1). With probability at least 1 − δ, we have If dν dτ (x) = e −g(x) / K e −g(y) dτ (y) for some g : K → R, i.e. −g is the log-density of ν up to a constant term, then with probability at least 1 − δ, Proof. In the first place, remark that for all ν 1 , ν 2 ∈ P(K) that are absolutely continuous w.r.t. p, we have is the cross-entropy and H(ν 1 ) is the differential entropy. Hence, for all ν 1 , ν 2 , ν 3 ∈ P(K), Secondly, notice that for any ν ∈ P(K) and measurable f : Thus, if we apply (13) on ν and its empirical version ν n = 1 n n i=1 δ xi , we obtain that with probability at least 1 − δ, for all f ∈ F: where we have used the Rademacher generalization bound (Mohri et al. [2012], Theorem 3.3) and the Rademacher complexity bound from the assumption of the theorem.
For the proof of (8) we apply Lemma 1 into (7).
Lemma 1. Let g : K → R be such that dν dτ (x) = e −g(x) / K e −g(y) dτ (y), i.e. −g is the log-density of ν up to a constant term. Then, Proof. Notice that ν = ν g . Thus, for any f ∈ F, Here, we bound and applying Lemma 2 to f and g, we obtain log Plugging these two bounds into (15), we obtain D KL (ν||ν f ) ≤ 2 f − g ∞ , which yields the result.
We do not claim that the upper-bound in Lemma 1 is tight; it might be possible to provide a bound involving a weaker metric. Regardless, it suffices for our purposes.
Lemma 2. Let f : K → R, g : K → R be measurable functions. For some α ∈ [0, 1], Proof. We define the function where p α (y) is the density of the Gibbs probability measure corresponding to the energy αf + (1 − α)g. We make use of the mean value theorem: Lemma 3 (Approximation of Lipschitz functions by F 2 balls, Proposition 6 of Bach [2017a]). For δ greater than a constant depending only on d, for any function f :

Rδ η
Proof. From Bach [2017a]. Notice that the factor in the bound is Rδ log δ η . The R factor stems from the fact that we consider the neural network features to lie in S d , while Bach [2017a] considers them in the hypersphere of radius R −1 .
Lemma 4 (Rademacher complexity bound for B F1 , Section 5.1 of Bach [2017a]; Kakade et al. [2009]). Suppose that K ⊆ {x ∈ R d+1 | x 2 ≤ R}. The Rademacher complexity of the function class B F1 is bounded by Assume that Assumption 1 holds. Then, we can choose β > 0 such that with probability at least 1 − δ we have where the notationÕ indicates that we overlook logarithmic factors and constants depending only on the dimension k.

Rβ ηJ
The first inequality holds because for all x ∈ K, U x 2 ≤ x 2 ≤ R by the fact that U has orthonormal rows, and the second inequality holds by (16). Thus, Notice that the assumptions of Theorem 1 are fulfilled: the Rademacher complexity bound for B F1 (Lemma 4) implies that R n (B F1 (β)) ≤ βR √ n and it is also easy to check that If we minimize the right-hand side w.r.t. β (disregarding the log factor), we obtain that the optimal value is Lemma 5 (Stein operator for functions on S d ). For a probability measure ν on the sphere S d with a continuous and almost everywhere differentiable density dν dτ , the Stein operator A ν is defined as for any h : S d → R d+1 that is continuous and almost everywhere differentiable, where ∇ denotes the Riemannian gradient. That is, for any h : S d → R d+1 that is continuous and almost everywhere differentiable, the Stein identity holds: Now, if we take the inner product of the right-hand side with the canonical basis vector e k ∈ R d+1 , we obtain where in the second equality we used that (I − xx ) the projection matrix to the tangent space of S d at x, in the third equality we used that it is symmetric, and in the last equality we used integration by parts on S d (∇· denotes the Riemannian divergence).
To compute ∇ · (e k − x k x), remark that by the invariance to change of basis it is equal to the divergence of the function g : R d+1 → R d+1 defined as x → e k − x k x x 2 , when restricted to S d . And we have For x ∈ S d , the right-hand side simplifies to −(d + 1)x k + 2x k − x k = −dx k , which means that the right-hand side of (18) becomes Lemma 6 (Kernelized Stein discrepancy for probability measures on S d ). For K = S d , and ν 1 , ν 2 ∈ P(K) with continuous, almost everywhere differentiable log-densities, the kernelized Stein discrepancy KSD(ν 1 , ν 2 ) is equal to where Proof. The argument for the first equality is from Theorem 3.8 of , but we rewrite it with our notation. Using the Stein identity, which holds by Lemma 5, we have Given the form of the Stein operator for functions on S d (Lemma 5), the proof of the second equality of (19) is a straightforward analogy of the proof of Theorem 3.6 of , which is for the Stein operator for functions on R d .
where H i are normed spaces of functions from S d to R. Assume that the following Rademacher complexity type bounds hold for 1 ≤ i ≤ d + 1: If we take n samples {x i } n i=1 of a target measure ν with almost everywhere differentiable log-density, and consider the Stein Discrepancy estimator (SDE)ν := νf , wheref is the estimator defined in (3), we have that with probability at least 1 − δ, SD H (ν,ν) is upper-bounded by Proof. Notice that by the definition of the Stein operator, Thus, using a union bound, we obtain that and through a change of variables, that means that with probability at least 1 − δ, All that is left is to upper-bound E Sn Φ i (S n ) for any i using Rademacher complexity bounds: and this is upper-bounded by By Talagrand's Lemma (Mohri et al. [2012], Theorem 5.7) and the uniform L ∞ bound on {∇ i f |f ∈ F} (notice that y → (β∇ i f (x j ) + dx j,i )y has Lipschitz constant uniformly upper-bounded by β∇ i f (x j ) + dx j,i ∞ , which means that the assumptions of Talagrand's Lemma are fulfilled), we have where we used the Rademacher complexity bound of B Hi . Using the Rademacher complexity bound of ∇ i h i as well, we conclude that the right-hand side of (20) can be upper-bounded by 2(βC1+Rd)C2+2C3 √ n . Thus, with probability at least 1 − δ, for all h ∈ H, We conclude the proof with an argument similar to the one of Theorem 1: In the second equality we use the definition of the Stein discrepancy (equation (2)). The inequality follows from (21) applied on ν f and onν = νf . The last equality holds because of the definition off and the definition of C 4 .
≤ 1}, we have that for the estimatorν defined in (3), with probability at least 1 − δ, where C 2 is a universal constant and ∇f denotes the Riemannian gradient of f .
which follows from the definition of the Stein operator and the Stein identity: E ν1 [A ν1 h(x)] = 0. Thus, for any ν ∈ P(K), Moreover, by Lemma 4: And where the last inequality follows from the Rademacher complexity bound on the hyperplane hypothesis, which is obtained through a VC dimension argument (Bach [2017a], Section 5.1; Bartlett and Mendelson [2002], Theorem 6). Moreover, h ∞ ≤ 1 and ∇ i h ∞ ≤ 1 for all h ∈ F 1 . The proof concludes by plugging (22), , all the upper-bounds of the proof hold for H = B F d+1 2 as well.
Theorem 5. Let K = S d . Let KSD be the kernelized Stein discrepancy for a positive definite kernel k with continuous second order partial derivatives, such that for any non-zero function g ∈ L 2 (S d ), of a target measure ν with almost everywhere differentiable log-density, and consider the unbiased KSD estimator (6), we have with probability at least 1 − δ, Proof. For the kernelized Stein discrepancy estimator we can write The third equality holds because of the definition ofũ ν in terms of u ν . In the first inequality we have used that for anyν (different from ν) with almost-everywhere differentiable log-density, Moreover, using the argument of Theorem 5.1 of , by Lemma 6, Theorem 3. Let K = S d . Assume that the class F is such that sup f ∈F { ∇f ∞ } ≤ βC 1 . Let KSD be the kernelized Stein discrepancy for a kernel that satisfies Assumption 2. If we take n samples {x i } n i=1 of a target measure ν with almost everywhere differentiable log-density, and consider the unbiased KSD estimator (6), we have with probability at least 1 − δ, Proof. We apply Theorem 5. We can bound , Lemma 7. For a function g : S d → R, we define the partial derivative ∂ i g : S d → R as the restriction to S d of the partial derivative of the polynomial power series extension of g to R n (i.e. the extension of a spherical harmonic to R n is the polynomial whose restriction to S d is equal to the spherical harmonic (Atkinson and Han [2012], Definition 2.7)). We denote by ∂g = (∂ i g) d+1 i=1 the vector of partial derivatives of g. The Riemannian gradient ∇g : S d → R d+1 , which is intrinsic (does not depend on the extension chosen), fulfills That is, ∇g(x) is the projection of ∂g(x) to the tangent space of S d at x.
For δ greater than a constant depending only on d, for any function g : S d → R such that for all x, y ∈ S d we have |g(x)| ≤ η and |g(x) − g(y)| ≤ η x − y 2 , and ∇g(x) 2 ≤ η and ∇g(x) − ∇g(y) 2 ≤ L x − y 2 , and g is even, there existsĝ ∈ F 2 such that ĝ F2 ≤ δ and where C(d) are constants depending only on the dimension d.
Proof. We will use some ideas and notation of the proof of Prop. 3 of Bach [2017a]. We can decompose . g k is the k-th spherical harmonic of g and P k is the k-th Legendre polynomial in dimension d + 1. Analogously, for any i between 1 and d + 1 for some r ∈ (0, 1). Equivalently,ĝ(x) = k,λ k =0 r k g k (x). Since g k is a homogeneous polynomial of degree k (Atkinson and Han [2012], Definition 2.7), we have thatĝ(x) = k,λ k =0 g k (rx) = g(rx).
With this choice ofĝ, the first equation of (25) holds by Prop. 3 of Bach [2017a].
Using this characterization ofĝ, by the chain rule we compute the Riemannian gradient The polynomial power series extension ∇g of ∇g is by definition equal to ∇g( Since the terms of k≥0 (∂g) k (x) − (∂g) k (x), x x are polynomials on x, this expression is equal to the polynomial power series of ∇g by uniqueness of the polynomial power series. Thus, for all x ∈ R d+1 , The second equality follows from Lemma 8, which states that ∂(g k ) = (∂g) k . Hence, by (26), we have r∂g(rx) − r ∂g(rx), rx rx = r ∇g(rx) = r k≥0 (∇g) k (rx) = r k≥0 r k (∇g) k (x). Thus, in analogy with Bach [2017a], we have Hence, keeping the analogy with Bach [2017a] (and Bourgain and Lindenstrauss [1988], Equation 2.13), we obtain that In the last equality we have used that ∂ i g is L-Lipschitz by assumption. And In the second equality we have used that ∇g(x) − ∇ĝ(x) is orthogonal to x (because it belongs to the tangent space at x), and the Pythagorean theorem. As in Bach [2017a], for δ > 0 large enough the argument is concluded by taking 1 − r = (C 1 (d)η/δ) 2/(d+1) ∈ (0, 1), which means that the (square root of the) error in the right-hand side of (27) Using that g is η-Lipschitz, by the argument of Bach [2017a] we have that ĥ L 2 (S d ) ≤ C 1 (d)η(1 − r) (−d−1)/2 , where C 1 (d) is a constant that depends only on d and consequently ĝ F2 ≤ C 1 (d)η(1 − r) (−d−1)/2 . And for our choice of r, this bound becomes ĝ F2 ≤ C 1 (d)η((C 1 (d)η/δ) 2/(d+1) ) (−d−1)/2 = δ.
Lemma 8. For g : S d → R with spherical harmonic decomposition g(x) = k≥0 g k (x) and with partial derivative with spherical harmonic decomposition Proof. Remark that the spherical harmonics on S d can be characterized as the restrictions of the homogeneous harmonic polynomials on R d+1 (Atkinson and Han [2012], Definition 2.7). k-th degree homogeneous polynomials are of sums of monomials of the form α i1,...,ir x i1 1 · · · · · x ir r , where r l=1 i l = k, and harmonic polynomials are those such that ∆p = d+1 i=1 ∂ 2 p ∂x 2 i = 0. Thus, for all k ≥ 0, g k can be seen as the restrictions to S d of homogeneous harmonic polynomials of degree k.
Notice that the i-th partial derivative of a homogeneous harmonic polynomial p of degree k is a homogeneous harmonic polynomial of degree k − 1. That is because by commutation of partial derivatives, we have Thus, ∂ i (g k ) are homogeneous harmonic polynomials of degree k − 1, which means that their restrictions to S d are spherical harmonics. Since ∂ i g(x) = ∂ i ( k≥0 g k (x)) = k≥0 ∂ i (g k )(x) and the spherical harmonic decomposition is unique, ∂ i (g k ) must be precisely the spherical harmonic components of ∂ i g.
Corollary 2. Let F = B F1 (β). Let Assumption 3 hold. (i) Whenν is the F 1 -SD estimator (2) and the assumptions of Theorem 2 hold, we can choose the inverse temperature β > 0 such that with probability at least 1 − δ we have that SD B d+1 where the notationÕ indicates that we overlook logarithmic factors and constants depending only on the dimension. (ii) Whenν is the unbiased KSD estimator (6) and the assumptions of Theorem 3 hold, β > 0 can be chosen so that with probability at least 1 − δ we have that KSD(ν,ν) is upper-bounded bỹ Proof. We will use Theorem 2. Let g : Letφ j : S k → R be the restriction of ϕ j to S k . By Lemma 7, there existsψ j : S k → R such thatψ j ∈ F 2 and ψ j F2 ≤ β/J, and Moreover, if we denote by ψ j : {x ∈ R k+1 | x 2 ≤ 1} → R the 1-homogeneous extension ofψ j , we can write the (Euclidean) gradient of ψ j at the point rx (with r ∈ [0, 1], x ∈ S d ) in terms of the (Riemannian) gradient ofψ j at x: Thus, by Equation 28, and renaming C(k), This also shows thatg j has F 1 norm g j F1 ≤ ψ j F2 ≤ β/J, which would mean thatg = J j=1g j ∈ F 1 and g F1 ≤ β. Moreover, The first inequality holds because the Riemannian gradient is the orthogonal projection of the Euclidean gradient of the extension, and orthogonal projections are 1-Lipschitz. The following equality holds because U j has orthonormal rows. The second inequality holds because for all x ∈ S d , U j x 2 ≤ x 2 = 1 by the fact that U j has orthonormal rows, and the third inequality holds by (28).
Thus, for part (i), we have Plugging this into Theorem 2 and using that If we optimize this bound with respect to β as in the proof of Corollary 1, we obtain and the optimal β is 2B √ n/(A(k + 1)) k+1 k+3 , where For part (ii), we plug into Theorem 3, and we obtain (using the notation of Theorem 3) that with probability at least 1 − δ, where A, B, D are defined appropriately. If we set β to minimize Aβ+B n 1/4 + Dβ −2/(k+1) , we obtain β = 2n 1/4 D A(k+1) k+1 k+3 , and the right-hand side becomes
For any m > 0, we define the m-particle gradient flow t → u m (t) = ((w Theorem 6. [Chizat and Bach [2018], Thm. 3.3; informal] Let R be a convex differentiable loss defined on a Hilbert space with differential dR Lipschitz on bounded sets and bounded on sublevel sets which satisfies a technical Sard-type regularity assumption. Let (µ t ) t≥0 be a Wasserstein gradient flow corresponding to F in (9), such that the support of µ 0 is contained in B(0, r b ) and separates the spheres r a S d+1 and r b S d+1 for some 0 < r a < r b . If (µ t ) t converges to µ ∞ in W 2 , then µ ∞ is a global minimizer of F . Moreover, if (µ m,t ) t≥0 is the empirical measure of (u m (t)) t≥0 and µ m,0 → µ 0 weakly, we have lim t,m→∞ F (µ m,t ) = lim m,t→∞ F (µ m,t ) = F (µ ∞ ).
Theorem 6 states that when the number m of particles (read neurons) goes to infinity, the function value of the gradient flow of the function G((w (i) , θ (i) ) m i=1 ) converges to a global optimum of F over P(R d+2 ). Remark that Algorithm 1 corresponds to the gradient descent algorithm on G((w (i) , θ (i) ) m i=1 ) with noisy gradient estimates. Thus, in the small stepsize and exact gradient limits, the iterates of Algorithm 1 approximate the gradient flow of G((w (i) , θ (i) ) m i=1 ). This reasoning provides an informal justification that Algorithm 1 should have a sensible behavior in the appropriate limits.
Observation 1. While Theorem 6 assumes that R is defined on a Hilbert space, this assumption is not convenient in our case because R(f ) = 1 n n i=1 f (x i ) + log K e −f (x) dx is not well defined on L 2 (R d ), as it involves pointwise evaluations. However, following the argument of Chizat and Bach [2018], up to the technical Sard-type regularity assumption, it suffices to show that F (µ) is a convex differentiable loss with a first variation δ δµ F (µ) such that |y| 2 ϕ(y/|y|) dµ(y) and · BL is the bounded Lipschitz norm.
To apply Theorem 6, we must check that the two statements in Observation 1 hold. Since for the maximum likelihood loss we have: we obtain that for all (w, θ) ∈ R d+2 and µ, µ ∈ P(R d+2 ), This shows the first point in Observation 1. For the second point, we have the following bound: F 2 is an RKHS with kernel k(x, y) = S d σ( x, θ )σ( y, θ )dp(θ) and for the ReLU unit this kernel has a closed-form expression [Cho and Saul, 2009]. Thus, one approach to optimize EBMs with energies over F 2 -balls is to apply the representer theorem and to write an optimizer f ∈ F 2 as f (·) = n i=1 α i k(x i , ·) for some α ∈ R n , as well as . Then, f becomes a finite-dimensional linear function of α, and thus any loss F that is convex in f is also convex on α. However, this approach scales quadratically with the number of samples and in practical terms, it is quite far from the way neural networks are typically trained.
The approach that we use to optimize EBMs over F 2 -balls is to sample random features (θ i ) m i=1 on S d from a probability measure with density q(·) and consider an approximate kernel k m (x, and Recht, 2008, Bach, 2017b. The functions in the finite dimensional RKHS H m with kernel k m are of the form h( Thus, learning a distribution with log-densities restricted in a ball of H m reduces to learning the outer layer weights (w i ) n i=1 . Namely, for R as in Subsec. 5.1, we optimize the loss which is convex. The gradient flow for G (with scaled gradient m∇ i G((w j ) m j=1 )) is and we can approximate it by gradient descent, which converges to the optimum (w i ) m i=1 if the gradients are exact and the stepsize is well chosen.
The connection between learning in H m balls and learning in F 2 balls is not straightforward. Applying Proposition 2 of Bach [2017b] and making use of the eigenvalue decay of the F 2 kernel on S d [Bach, 2017a], for an appropriate choice of q we have that for all f ∈ B F2 , there existsf ∈ H m with f Hm ≤ 2 such that f −f L 2 (p) ≤ O (m/ log(m)) −(d+3)/2 . This L 2 error bound is sufficient to produce a quantitative result for least squares regression. However, for the three losses considered in this paper we would need bounds for f −f ∞ and ∇f − ∇f ∞ , which do not seem to be available (Bach [2017b] does provide a bound on f −f ∞ , but under the assumption that kernel eigenfunctions have a common L ∞ norm bound, which does not hold for spherical harmonics in S d ).
Nonetheless, a mean-field qualitative approach analogous to the F 1 case is still possible (see Proposition 2.6 of Chizat and Bach [2018]). The learning objective in F 2 can be written as and the mean-field dynamics is  If we choose q(·) = 1, we have that (29) is the m-particle approximation of (30). Let h be the global minimizer of F , which is reached at a linear rate by (30) because F is strongly convex. Skipping through the details, the argument of Lemma C.15 of Chizat and Bach [2018] could be adapted to yield:

C Additional experiments
In this section, we show plots corresponding to additional experiments. Figure 6 shows results for KSD and F 1 -SD training in the case J = 1, w * 1 = 10. Compared to the plots for J = 1, w * 1 = 2 shown in Figure 1, the separation between the F 1 and F 2 EBMs becomes much more apparent. Figure 7 shows results for KSD and F 1 -SD training in the case J = 2, w * 1 = −2.5. The separation between the F 1 and F 2 EBMs is smaller than in the case J = 2, w * 1 = −5 shown in Figure 2.

D Duality theory for F 1 and F 2 MLE EBMs
In this section we present the dual problems of min f ∈F H(ν n , ν f ) (i.e. problem (1)) for the cases F = F 1 , F 2 (Subsec. D.1), F = B F1 (β) (Subsec. D.2) and F = B F2 (β) (Subsec. D.3). The dual problems take the form of entropy maximization under hard constraints, L ∞ and L 2 moment penalizations, respectively. The tools used involve a generalized minimax theorem and Fenchel duality, which was also used for results of the same flavor in finite dimensions (c.f. Mohri et al. [2012]). The proofs are in App. E.  The teacher model has two neurons of negative weights w * 1 , w * 2 = −2.5 and random positions in the hypersphere. (Right column) Test cross entropy and test F 1 -SD for model trained with F 1 -SD at different training sample sizes n, in d = 10. The teacher model has two neurons of negative weights w * 1 , w * 2 = −2.5 and random positions in the hypersphere.

D.1 Duality for the unconstrained problem
Consider the following entropy maximization problem under generalized moment constraints: recalling that τ is the uniform probability measure over K and letting β > 0 be arbitrary. The constraints in this problem can be interpreted either (i) as an equality constraint in C(S d ), i.e., the set of continuous functions on S d , or (ii) as an equality constraint in L 2 (S d ), i.e., the set of square-integrable functions on S d . Each interpretation yields a different dual problem.
By the Riesz-Markov-Kakutani representation theorem, the set of signed Radon measures M(S d ) can be seen as the continuous dual of C(S d ). Hence, in the case (i), the Lagrangian for problem (31) is L 1 : M(K)×M(S d )× which is equivalent to the MLE problem (1) when F = F 1 .
Letτ to denote the uniform probability measure over S d . In the case (ii), the Lagrangian for problem (31) is which is equivalent to the MLE problem (1) when F = F 2 .

D.2 Duality for the F 1 -ball constrained problem
Using ν n = 1 n n i=1 δ xi to denote the empirical measure, consider the following problem, which can be seen as an L ∞ -penalized version of (31): As shown in Theorem 8, the dual of this problem is a modified version of (32) in which µ is constrained to have TV norm bounded by 1: Remark that by the definition of F 1 , the problem (35) is equivalent to MLE problem (1) in the case F = B F1 (β).
Theorem 8. The problem (35) is the dual of the problem (34), and strong duality holds. Moreover, the solution ν of (35) is unique and its density satisfies where µ is a solution of (34) and Z β is a normalization constant.

D.3 Duality for the F 2 -ball constrained problem
The following problem can be seen as an L 2 -penalized version of (31): And as shown in Theorem 9, the dual of this problem is a modified version of (33) in which h is constrained to have L 2 norm bounded by 1: Remark that by the definition of F 2 , the problem (35) is equivalent to MLE problem (1) in the case F = B F2 (β).
Theorem 9. The problem (37) is the dual of the problem (36), and strong duality holds. Moreover, the solution ν of (37) is unique and its density satisfies where h is a solution of (36) and Z β is a normalization constant.

E Proofs of App. D
Theorem 10. [Kneser [1952]] Let X be a non-empty compact convex subset of a locally convex topological vector space space E and Y a non-empty convex subset of a locally convex topological vector space space F . Let the function f : X × Y → R be such that: (i) For each y ∈ Y , the function x → f (x, y) is upper semicontinuous and concave, (ii) For each x ∈ X, the function y → f (x, y) is convex.
Observation 2. Notice that for any functional f : M(K) → R, we have Theorem 7. Strong duality holds between the entropy maximization problem (31) and each of the two dual problems (32)-(33).
Proof. We start with (32). First, we prove that it is indeed the dual problem of (31). As stated in the main text, the problem (31) admits a Lagrangian L 1 : The Lagrange dual function is where we have used that at the optimal ν, the first variation of L 1 w.r.t. ν must be zero: The Lagrange dual problem is sup µ,λ,g≥0 and the right-hand side is precisely (32). In the second equality we used that the optimal choice for g is g = 0.
In the third equality we used that the optimal λ must satisfy the first-order optimality condition: To prove strong duality, we need to show that L 1 (ν, µ, g, λ).
If we defineL 1 : P(K) × M(S d ) → R asL 1 (ν, µ) = L 1 (ν, µ, 0, 0), we have that the assumptions of Theorem 10 hold for −L 1 . Indeed, by Lemma 9 we have that −L 1 (·, µ) is a concave and upper semicontinuous function of ν. And by Prokhorov's theorem, P(K) is a compact subset of the locally convex topological vector space of signed Radon measures with the topology of weak convergence (tightness follows from the fact that K is compact). Thus, On the one hand, notice that by Observation 2, where we have used Lemma 10 in the first equality, Observation 2 in the second equality and the definition of L 1 in the third equality. Thus, the strong duality (40) follows from plugging (42) and (43) into (41).
To show that (33) is also a dual problem of (31), we consider the Lagrangian L 2 : M(K)×L 2 (S d )×C(K)×R → R defined as The reasoning to obtain the dual problem (33) is analogous. Strong duality in this case can be stated as L 1 (ν, h, g, λ).
Notice that C(K) ⊆ M(K) * by the fact that a vector space has a natural embedding in its continuous bidual (but the continuous bidual is in general larger). Through this identification, (46) implies that we can write A * µ(x) = S d σ( θ, x )dµ(θ). Our goal now is to compute the convex conjugates f * and g * . By the argument of Lemma B.37 of Mohri et al. [2012], which works in the infinite-dimensional case as well, the convex conjugate f * : C(K) → R ∪ {+∞} is shown to be: Remark that f * has domain M(K) * , which is larger than C(K). However, knowing the restriction of f * to C(K) will suffice for our purposes.
Moreover, g * : M(S d ) → R ∪ {+∞} fulfills: In the first equality we have used the definition of g, in the fourth equality we have used Theorem 10 (remark that {µ ∈ M(S d ) : |µ | T V ≤ 1} is compact in the weak convergence topology), and in the fifth equality we have used that sup ϕ∈C(S d ) S d ϕ(θ) d(µ − µ )(θ) = +∞ unless µ = µ .
With these definitions, notice that problem (34) can be rewritten as inf ν∈M(K) {f (ν) + g(Aν)} and problem (35) can be rewritten as sup µ∈M(S d ) {−f * (−A * µ) − g * (µ)}. Thus, strong duality between (34) and (35) follows from Fenchel strong duality, which holds by checking condition 2 of Theorem 11. We have to see that A dom f ∩ cont g = ∅. Consider ϕ(·) = K σ( ·, x ) dν(x) ∈ C(S d ) for some ν ∈ P(K) absolutely continuous w.r.t. τ . Then, we have that ϕ ∈ A dom f . Moreover, since g is a continuous function (in the supremum norm topology), cont g = C(S d ) and hence ϕ ∈ cont g as well, which means that the intersection is not empty.
Notice that in our case d * = sup µ∈M(S d ) {−f * (−A * µ) − g * (µ)} ∈ R, which by Theorem 11 implies that the supremum is attained: let µ be one maximizer. We show that p * = inf ν∈M(K) {f (ν) + g(Aν)} = inf ν∈P(K) {f (ν) + g(Aν)} admits a minimizer by the direct method of the calculus of variations: notice that f and g • A are lower semicontinuous in the topology of weak convergence (f by Lemma 9 and g • A because it is a maximum of continuous functions, and thus its sublevel sets are closed because they are the intersection of closed sublevel sets), and P(K) is compact.
Theorem 9. The problem (37) is the dual of the problem (36), and strong duality holds. Moreover, the solution ν of (37) is unique and its density satisfies where h is a solution of (36) and Z β is a normalization constant.
Proof. The proof is largely analogous to the proof of Theorem 8. We use Theorem 11 with X = M(K) as before, and Y = L 2 (S d ), the Hilbert space of square-integrable functions on S d under the base measureτ , which is of course self-dual. We define f as before, and g : L 2 (S d ) → R ∪ {+∞} as , and consequently, g * : L 2 (S d ) → R ∪ {+∞} fulfills: and using Theorem 10 once more, this is equal to: With these definitions, notice that problem (36) can be rewritten as inf ν∈M(K) {f (ν) + g(Aν)} and problem (37) can be rewritten as sup ψ∈L 2 (S d ) {−f * (−A * ψ) − g * (ψ)}. The rest of the proof is analogous.