Solving the Rubik's cube with stepwise deep learning

This paper explores a novel technique for learning the fitness function for search algorithms such as evolutionary strategies and hillclimbing. The aim of the new technique is to learn a fitness function (called a Learned Guidance Function) from a set of sample solutions to the problem. These functions are learned using a supervised learning approach based on deep neural network learning, that is, neural networks with a number of hidden layers. This is applied to a test problem: unscrambling the Rubik's Cube using evolutionary and hillclimbing algorithms. Comparisons are made with a previous LGF approach based on random forests, with a baseline approach based on traditional error‐based fitness, and with other approaches in the literature. This demonstrates how a fitness function can be learned from existing solutions, rather than being provided by the user, increasing the autonomy of AI search processes.

Another approach focuses on how the fitness function itself can be transformed to avoid these problems. This has a long history in the evolutionary computation, typified by work on fitness scaling (Grefenstette, 1986;Hopgood & Mierzejewska, 2009;Kreinovich et al., 1993;Ware et al., 2003). The principle aim of fitness scaling is to prevent premature convergence of the search algorithm, by composing a scaling function with the fitness function that does not change the ranking of points in the search space but ensures a more even distribution of the fitness values allocated to those points. A more recent version of this transformation approach is exemplified by geometric semantic genetic programming (GSGP) (Moraglio et al., 2012), which attempts to reconfigure the problem so that a much simpler search process such as hillclimbing can be used.
The aim of these approaches is to create a transformed problem that can be solved more effectively. However, these transformations sometimes involve some kind of tradeoff. For example, in basic GSGP, this tradeoff is against the size of the solution, though more recent implementations have used a caching strategy to make implementation more efficient (Vanneschi et al., 2013). This idea of reconfiguring the fitness function prior to the main evolutionary algorithms being run is one source of inspiration for the work in this paper; this has been explored elsewhere in evolutionary computation in work showing how a good choice of genotype-phenotype mapping can be used to create a smoother landscape (Asselmeyer et al., 1996).
A domain-specific smoothing the search landscape is in the form of pattern databases (Culberson & Schaeffer, 1998). These consist of patterns in the search space such that any point that matches the pattern has the same cost of solution-typically, these represent symmetries of the underlying problem. If a solution of a particular cost is known for one problem that matches the pattern, then any other solution matching the pattern will have at most that cost to solve because all of the moves to the solution can be similarly transformed. As with the approaches discussed earlier, the idea here is to transform the search space. The key difference is that pattern databases are based on domain knowledge. Some work has used learning methods to generalize from pattern databases-for example, by using neural networks to learn how pattern databases can be combined (Mehdi Samadi & Felner, 2008).
Another important source of inspiration is the view that traditional fitness functions take a very narrow view of the problem; while a traditional fitness function is a good guide as to which elements of the population to choose for the next generation, it is a very simple representation of the complexity of a problem. Instead, it is argued, rather than a fitness function that returns a single number or a ranking, we should be using more complex fitness drivers that give us more information about the population member, allowing a more directed application of operators (Krawiec, 2016;Krawiec et al., 2016). However, such fitness drivers can require more domain-specific knowledge than a traditional fitness function. One of the aims of this paper is to give a generic method by which more information about problems can be incorporated into the evolutionary search, in this case by pre-training.
A more fundamental problem for evolutionary algorithms is that for some problems, defining the fitness function is difficult, because each problem has a different goal state. Call these non-oracular problems. As an example, consider the protein-folding problem in bioinformatics (Dobson, 2003). Biological proteins consist of a sequence of amino acids, which then fold into a three-dimensional shape, which is (with a few exceptions such as prions) entirely dependent on the sequence. To define this as a traditional evolutionary search is problematic, because we do not have access to a measure of how far a particular configuration is from the solution-indeed, if we did know what configuration we were searching for, we would have solved the problem! Therefore, evolutionary computing approaches to these types of problems have focused on learning parameters in, or functional forms of, a domain-specific model (Widera et al., 2010).
Another potential advantage to pre-training for simplifying the fitness landscape is that more extensive computational effort can be expended during an early training phase, and then when evolution is applied to a specific problem, the evolutionary algorithm can run in fewer generations because more domain-specific information has been encoded into the fitness function. This may be of importance in some application where running a traditional evolutionary algorithm might be infeasible because of the need for a large population and many generations to escape local minima, whereas a smaller population and fewer generations might be needed for the simpler function. This is particularly important where the evolutionary algorithm needs to run in a time-constrained situation.

| DEEP LEARNED GUIDANCE FUNCTIONS
Fitness functions in evolutionary learning are provided as part of the problem definition. Typically, these are then used directly-individuals are evaluated using the fitness function, and operators in the search are used to avoid problems in the fitness landscape such as local minima. However, an alternative approach has been applied in both evolutionary learning (Szubert et al., 2013) and reinforcement learning (Erez & Smart, 2008), where the fitness function is shaped so that it more directly represents routes through the fitness landscape from an arbitrary point to the desired target.
A form of this called Learned Guidance Functions (LGFs) was introduced by (Johnson, 2018). The input to this is a search space and set of existing solution trajectories for the problem. For example, in the protein folding problem these would be sequences of points in the space of threedimensional structures, going from a sequence to a completely folded structure. For an image denoising problem, this would be a sequence of images from a clean image to a very noisy one. These are an example of True Distance Heuristics (Sturtevant et al., 2009), but with a particular layered structure and the use of a predictive function to give the heuristic value rather than a look-up table.
These solution trajectories can be obtained in a number of ways. For some problems, we will have access to a set of already-solved examples.
For others, we can construct artificial examples by starting from a solved state and carrying out a number of moves from that solved state to generate trajectories in reverse (a similar approach has been called Autodidactic Iteration in McAleer et al., 2018).
These trajectories can then be used to create a training set for a supervised learning problem. Each trajectory will consist of a number of states in the search space of the problem, each of which is paired with a number that is the number of steps away from the target that it took to get to that state.
These pairs then become the training set: so, the task for the supervised learning problem is to build a model that takes an arbitrary state of the system and assigns a number predicting how many steps it will take to get to the target state. The LGF is the model learned from this supervised learning process.
This LGF then be used as a ranking function in an evolutionary algorithm. Take each member of the population, and apply the LGF to it. Then select the individuals that will form the parents of the next population from the lowest scoring ones on the LGF-these are the ones that are being predicted as being closest to the solution.
A similar approach has been taken in the recent papers by McAleer et al. (2018) and Agostinelli et al. (2019), though these are grounded in a reinforcement learning approach rather than a supervised learning approach. Compared with the work in this paper, their algorithm learns a mapping from points in the state space of the cube to a pair consisting of a value and a policy. The approach that we take in this paper can be seen as the equivalent of learning just the value function. Indeed, the greedy experiment in McAleer et al. (2018, fig. 5a) demonstrates a similar percentage-solved behaviour to the experiments in this paper. However, their full approach, which combines value function approximation and policy learning, and then uses an Monte Carlo Tree Search (Browne et al., 2012) approach, continues to work for problems of much higher complexity that the experiments in this paper, but at the cost of much longer search times (around 10 min as shown by Figure 5b and the description on Page 6 of McAleer et al., 2018) compared to a few seconds for our approach. In summary, their approach is scalable to larger problems (more scrambles), but at the cost of increased computation time, whereas the approach in this paper offers a quicker solution approach to simpler problem instances but fails for larger problems.

| Formalization
Now we formalize this idea. Consider a search space S consisting of a set of points, which is the node set of a directed graph M S , which represent the possible moves (mutations) from each point in the search space. Identify one or more of these as goal states; these might be the only goal states, or they might represent a sample of the class of states that the eventual problem is trying to solve. X can now be used as a training set for a supervised learning algorithm. The trained model from that supervised learning algorithm, L : S ! ℤ ≥0 , is a function that takes a set in the search space and predicts how many moves are needed to get to the goal state. This function will be used as an alternative kind of fitness function in the experiments below.

| EXAMPLE: APPLYING DEEP LGFS TO THE RUBIK'S CUBE
In Johnson (2018) we applied the LGF to the problem of unscrambling the Rubik's Cube. We used a number of classifiers from the scikit-learn library (scikit-learn, n.d.) to implement LGFs, and demonstrated that (1) the LGF can learn to recognize the number of turns that have been made to a cube to a decent level of accuracy; and, (2) that this LGF can then be used to unscramble particular states of the cube in a sensible number of moves. Unscrambling is not one of the non-oracular problems, because the goal state is known, but it has a complex fitness landscape with many local minima, and so is a good test for these kind of algorithms.
The search space C consists of all possible configurations of coloured facelets on the six faces of the cube, each of which has a 3 × 3 set of facelets. The move set M is notated by a list of twelve 90 moves, Singmaster, 1981), which are functions from C to C. We use the notation m(c) to denote the application of move m M to the cube c, returning the new state of the cube.
An earlier paper by the author (Johnson, 2018) applied a number of learning algorithms to the problem of learning LGFs for the Rubik's cube, with random forests demonstrating itself to be the best approach. That paper did not use any deep learning (Goodfellow et al., 2017) approaches-in this paper, we extend the work by using deep learning.

| Constructing the LGF
The LGF for this problem is constructed as follows (pseudocode in Algorithm 1). For n s iterations, start with a solved cube and make n ℓ − 1 moves.
Each time a move is made (and in the initial state), the pair consisting of the current state and the number of moves made to get to that state is added to the training set. This is illustrated in Figure 1.
The LGF is then constructed from this training set by applying a supervised learning algorithm, specifically a deep neural network implemented in the Keras framework (Keras, n.d.) on TensorFlow (TensorFlow, n.d.). The specific network used is illustrated in Figure 2. This is a fairly standard deep learning network, with dropout (Srivastava et al., 2014) used to encourage generalization and prevent overfitting. The categorical crossentropy function was used for the loss function, and the Adam optimizer was applied. Future work will apply meta-learning of parameters and network shape to optimize the model produced (Hutter et al., 2019).
Once an LGF is learned, it can be applied to the task at hand, which is to take a scrambled state of the cube and move through the search space with the aim of finding the solved state. This is done using a variant on evolution strategies (ES). The initial state of the cube is duplicated to fill the population. Then, in each generation a number of mutants are generated by making a random move for each member of the population.
Any solutions that are predicted by the LGF to be closer to the solution than the current one are placed in an intermediate population pool, and a new generation created by uniform random sampling with replacement from this pool to bring the population up to full size. This is repeated until one of three states occurs: (1) the solution is found; (2) none of the mutants produce any improvement, in which case the algorithm is restarted;

Algorithm 1
Training set construction for the Rubik's cube 1: procedure ConstructTrainingSetRubik (n s , n m ) This is summarized in pseudocode in Algorithm 2, where the inputs are: n ℓ , the problem size (number of scrambling twists given); n p , the population size; θ the maximum number of generations; and, L the LGF function used. In the results tables, this is referred to as Deep LGF + ES

(ES refers to Evolution Strategies).
When Algorithm 2 discovers a solution, it is then easy to reconstruct a path through the various population iterations to create a solution trajectory. Each population consists of solutions that are mutations (moves) from the previous population, so by starting from the solution, identifying which cube it was a mutation of, and so on back to the original scrambled state, we can construct the trajectory (Figure 3).
A second approach is to use the same fitness function, but to remove the population-based aspect of it. So, instead of having a population of cube states, we start with a single cube state, and then at each iteration generate all possible 12 moves from that state, and use the deep network to predict how many scrambles it is from solved. If one of these predictions is closer to the solved state, then that move is made. Otherwise, a random (neutral) move is made. This is summarized in pseudocode in Algorithm 3, with the input parameters being a subset of those in Algorithm 2.
In the results tables, this is referred to as Deep LGF + Hillclimbing.

| Sources of error
Note that if a perfect LGF existed for a problem, we could solve the problem in a minimal number of steps. Starting from an arbitrary scrambled state, we can examine all possible moves from that state. At least one of these will be closer in terms of number of moves to the target state, and so we can move the state of the system to the state which is closest, and repeat until we reach the target state.
In practice, there are two forms of error. The first is in the formation of the training set for the problem. A particular sequence of scrambling moves of length n might, nonetheless, end up with the cube in a state which could have been reached using fewer moves. A simple example of this is where one move is followed by a move which is the inverse of that move (this is explored in more detail by Johnson (2018)).
The second is where the model makes the wrong prediction. For these reasons, the fitness landscape created by a real LGF will still have local minima.
F I G U R E 2 Keras deep learning network used for training

| EXPERIMENTS AND RESULTS
The experiments were carried out as follows. LGF function used L is the result of running Algorithm 1 with trajectory length n m and 100,000 trajectories, then using that as the training set for the Keras network in Figure 2 with 50 epochs.
The total time to run all of these experiments was under 3 h, not including time to train the models (training time was between 11-59 s per epoch depending on the size of the model).
Results for the unscrambling experiments are presented in two tables. Table 1 shows for each (n ℓ , n m ) pair the percentage of times that the problem was solved. Table 2 shows the number of generations taken by successful algorithms to unscramble the cube.
There are a number of observations. Firstly, for the smaller problem sizes, a solution to the problem is frequently found; for problems below size 9, at least half of the attempts are successful, and it is very reliable for small problem sizes. Secondly, the size of the model makes little difference-using a larger model than the problem size is of little value. Thirdly, the number of generations needed is small for the lower problem sizes, but increases for large problem sizes; this may be an effect of more re-initialisations needing to be carried out.
Fourthly, note that some of the average lengths in Table 2 are shorter than the problem size. This is because the problems were constructed by scrambling randomly n ℓ times, but no check was made to ensure that the resulting state could not be solved in less than n ℓ moves; indeed, doing such a check is rather complex. Therefore, the starting state for some runs may contain a problem that can be solved in fewer than n ℓ moves.
Tables 3 and 4 present the results for the Deep LGF combined with hillclimbing. The total time to run all of these experiments was considerably than the previous experiments, at around 52 h, not including time to train the models (the models were reused from the previous experiment).
The results in terms of percentages solved are overall worse for the hillclimbing approach rather than the evolutionary approach. The lengths to solution are larger in the hillclimbing process for simpler problems, but are in many cases better for more complex problems.
Tables 5 and 6 compare the results of the experiments in this paper (a preliminary version of which as presented by Johnson (2019)) to two experiments in a previous paper (Johnson, 2018). The main experiments in the current paper (Deep LGF + ES) varied from the experiments in this earlier paper (Random Forest LGF + Hillclimbing) in two main ways. Firstly, the models were trained using a random forest classifier (the implementation in the scikit-learn package (scikit-learn, n.d.)). The tables give the results for models trained on examples with up to 13 moves. Secondly, the unscrambling in the earlier paper was based on only the hill-climbing approach, whereas in this paper the evolution strategy approach is also used.

F I G U R E 3 Reconstructing a single trajectory from the evolutionary search once the solution has been found
These tables also contain a comparison with a baseline experiment also described in detail in the earlier paper (Johnson, 2018) (Error + Hillclimbing). This also uses a simple hill-climbing method, but the choice of moves is made by choosing the move that maximizes the number of correct facelets. This is a traditional error-based fitness function.
It is notable that the percentage of successes in the Deep LGF + ES approach is considerably higher than the Random Forest LGF + Hillclimbing approach. However, the length of the solutions found by the new approach is much larger for larger problem sizes. This may well reflect the use of reinitialisation in the latter approach; in the earlier paper, a search that did not terminate was considered a failure. The evolutionary approach using the Deep LGF is largely a better performer than the hillclimbing approach. All of these new methods clearly outperform the traditional error-based fitness measure, demonstrating the value of this pre-training step.

| WHAT IS BEING LEARNED BY THE DEEP NETWORK?
What is being learned by the network? Is it learning the structure of solutions, or relationships between fixed values in the network? This next set of experiments will investigate this.
Recent work by Geirhos et al. (2020) that demonstrates that deep learning learns so-called "short-cuts" to classifying a particular class of items in the training set. For example, a deep learning network will misclassify an image of a cow when presented on a beach rather than a field, whereas a human would not make the same mistake (Beery et al., 2018), and an object recognition task will prioritize the position of an object in the visual field over the content of the objects in that field (Geirhos et al., 2020). Deep learning demonstrates cognitive biases in learning (Geirhos et al., 2019), that might not always align with human cognitive biases.
T A B L E 4 Average number of generations needed to solve problem of each size using trained model of each size   Size of problem   2  3  4  5  6  7  8  9  10  11 12 13 Size of model We would expect an intelligent solution to a problem such as learning an LGF for the Rubik's Cube to learn the structure of the problem, rather than learning relationships between colours at specific positions. It seems reasonable that a system that has learned tout court that a cube with four sides each with a single row of miscoloured faces is likely to represent a single scramble has found a more intelligent and generalisable solution to the problem than one that has learned that relationship based on lines of specific colours, and is unable to generalize when those colours are permuted.
In this next set of experiments, we will explore these ideas in the kinds of models used in this paper. This will be done by creating a restricted training set, T a based solely on anticlockwise scrambles. This will then be tested on two test sets: one generated using anticlockwise scrambles, and generated using clockwise scrambles.
The hypothesis being tested is that the deep learning system will not be able to develop a structural understanding from the restricted training set. If it did, it should be able to transfer that understanding from the anticlockwise training to the clockwise test.
In detail, a solved cube c was generated, and the pair (c, 0) added to the (initially empty) training set T a . This was then scrambled 10 times (n = 1, …, 10) using just anticlockwise scrambles, and after each of these scrambles the pair (c, n) added to T a . This process was repeated 100,000 times to produce the training set. A model M was then generated using T a and the network from earlier (Figure 2), again trained over 50 epochs.
The testing was as following. For each of n = 2, …, 9, a testing set was created of 100 cubes generated using random anticlockwise turns only, as was a second testing set in the same way except that clockwise turns were created. Predictions were then made using M, and a count made of how often the model made a correct prediction. The results are presented in Table 7.
Notably, the performance is overall considerably worse for the clockwise test set, though this evens out for more complex problems (7,8,9 scrambles) where the overall performance is lower. This demonstrates, certainly for the regime seen at smaller problem size, that the learning does not generalize from the anticlockwise to the clockwise.
This points towards alternative approaches based on representations that have a more explicit representation of these structural relationships. For example, it would be interesting to explore whether a representation that contains explicit representations of relationships between facelets on the cube, such as program code manipulated with an approach such as genetic programming (GP) (Poli et al., 2008), would make this kind of generalization more possible.
More generally, an interpretable representation (Murdoch et al., 2019) such as those generated by GP (at least for small program sizes) would not only be useful for human interpretability, but would also expose the key meaningful ideas used in classification. This would then provide a representation that could be manipulated by another AI system, facilitating transfer learning (Pan & Yang, 2010) to other related problems.

| OTHER RELATED WORK
There are similarities between this approach and the idea of a learned value function in reinforcement learning (Sutton & Barto, 2018 In the metaheuristics literature the idea of learning from a set of already-solved problems is explored in the idea of target analysis (Glover & Greenberg, 1989). This takes a set of solved problems from a problem class, and uses those known solutions to set the parameters of a metaheuristic. This is different to our approach in that it still relies on the metaheuristic operators to avoid local optima in the landscape, whereas the approach in this paper uses those already-solved problems to construct a new landscape based on a metric which is designed to have fewer such local optima. The idea of learning a metric from a large set of examples is explored in the literature on metric learning (Kulis, 2013), and it would be interesting to explore further connections between metric learning and the idea of constructing new fitness functions.
It should be noted that there are algorithms specifically for solving the Rubik's Cube, as summarized in the book by Slocum et al. (2011). However, comparisons with these are less relevant to this paper, which was using the Rubik's Cube as an example to see whether a learning algorithm could learn naïvely from it. The importance of these methods that can learn without explicit human knowledge has been emphasized as an important route towards artificial general intelligence (Silver et al., 2017). More widely, this approach points towards a human-like approach where the machine learns concepts such as similarity from existing solutions. There is an interesting continuum of learning methods from much contemporary machine learning that attempts to generalize from examples (whether labelled examples as in supervised learning or reward as in reinforcement learning), through the approaches in this paper that learn from solution trajectories, through to expert systems where the inputs data explicitly attempts to articulate what aspects of successful solutions should be used by the machine.

| SUMMARY AND FUTURE WORK
We have developed the idea of deep learning for learned guidance functions, and shown how these can then be used as fitness drivers in evolutionary computation. This has been applied to a case study of solving a Rubik's Cube, and shown to have a advantages in terms of frequency of finding a solution and the size of the models needed when compared to a random forest-based LGF; however, the number of generations needed is, for more complex problems, larger. Compared to the work by Agostinelli et al. (2019), the results for smaller problems are comparable but quicker to compute, but the combination of policy and value learning approach in that paper allows reliable solution of more complex problems compared to the approach in this paper that relies solely on value function approximation.
There are a number of areas for future work. Firstly, there is much of scope for optimizing the deep learning system using automated machine learning approaches both to optimize the parameters and the structure of the system (Hutter et al., 2019). Secondly, there are a number of further experiments that would investigate the behaviour further: investigating the frequency of and impact of the reinitialisation in this method, using measures of landscape smoothness to understand the effect of the LGF on the landscape, and experimenting with different population sizes.
Finally, there are a large number of other problems to which this approach could be applied, for example, protein folding, and de-noising of audio and video files, audio transcription (Souto-Rico et al., 2020), etc.

CONFLICT OF INTEREST
The author declares that there is no conflict of interest that could be perceived as prejudicing the impartiality of the research reported.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available from the corresponding author upon reasonable request. The code that supports the findings of this study is available in the supplementary material of this article.


extfieldsofstudy:
Computer Science
Mathematics

text:
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.