1 Introduction
The herding algorithm has recently been presented by Welling (2009b)
as a computationally attractive alternative method for learning in intractable Markov random fields models (MRF). Instead of first estimating the parameters of the MRF by maximum likelihood / maximum entropy (which requires approximate inference to estimate the gradient of the partition function), and then sampling from the learned MRF to answer queries, herding directly generates
pseudosamples in a deterministic fashion with the property of asymptotically matching the empirical moments of the data (akin to maximum entropy). The herding algorithm generates pseudosamples with the following simple recursion:(1)  
where is the observation space; is a feature map from to
, which could be viewed as the vector of sufficient statistics for some exponential family, and
is a mean vector to match (the empirical moment vector of the same family). Unlike in frequentist learning of MRFs, the parameter never converges to a point in herding and actually follows a “weakly chaotic” walk (Welling & Chen, 2010).The herding updates can be motivated from two different perspectives. From the learning perspective, the herding updates can be derived by performing fixedstepsize subgradient ascent on the zerotemperature limit of the annealed likelihood function of the MRF—called the “tipi function” by Welling (2009b). From this perspective, herding was later generalized to MRFs with latent variables (Welling, 2009a) as well as discriminative MRFs (Gelfand et al., 2010).
From the moment matching perspective, which has been explored more in details by Chen et al. (2010), the herding updates can be derived as an effective way to choose greedily pseudosamples in order to quickly decrease the moment discrepancy (Chen et al., 2010). Under suitable regularity conditions, decreases as for the herding updates—this is faster than i.i.d. sampling from the distribution generating (e.g., the training data) which would yield the slower rate. This faster rate has been explained by negative autocorrelations amongst the pseudosamples and was used by Chen et al. (2010) to subselect a small collection of representative “supersamples” from a much larger set of i.i.d. samples. We make the following contributions:

We show that herding as described by Eq. (1) is equivalent to a specific type of conditional gradient algorithm (a.k.a. FrankWolfe algorithm) for the problem of estimating the mean . This provides a novel understanding and another explicit cost function that herding is minimizing.

This interpretation yields improvements, for the task of estimating means, with other faster variants of the conditional gradient algorithm, which lead to nonuniform weights, one based on linesearch, one based on an activeset algorithm.

Based on existing results from convex optimization, we extend and improve the convergence results of herding. In particular, we provide a linear convergence rate for the linesearch variant in finitedimensional settings and show how the conditions assumed by Chen et al. (2010) in fact never hold in the infinitedimensional setting.

We run experiments that show that algorithms which estimate faster the mean than herding generate samples that are not better (and typically worse) than the ones obtained with herding when evaluated in terms of the ability to approximate a sample with large entropy, a property which (if or when satisfied by herding) could be the basis for an interpretation of herding as a learning algorithm (Welling, 2009b).
2 Mean estimation
We start with a similar setup as Chen et al. (2010), where herding can be interpreted as a way to approximate integrals of functions in a reproducing kernel Hilbert space (RKHS). We consider a set and a mapping from to a RKHS . Through this mapping, all elements of may be identified with real functions on defined by , for . We denote by the associated positive definite kernel. Note that the mapping may be explicit (classically in lowdimensional settings) or implicit—where the kernel trick can be used, see Section 4.3 and Chen et al. (2010).
Throughout the paper, we assume that the data is uniformly bounded in feature space, that is, for all , , for some ; this condition is needed for the updates of Eq. (1) to be welldefined.
We denote by the marginal polytope (Wainwright & Jordan, 2008; Chen et al., 2010), i.e., the convexhull of all vectors for . Note that for any , we have
and that for all and (i.e., all functions with finite norm are bounded).
Extreme points of the marginal polytope. In all the cases we consider in Section 5, it turns out that all points of the form , are extreme points of (see an illustration in Figure 1). This is indeed always true when is constant for all (for example for our infinitedimensional kernels on in Section 5.1); it is also true if contains both an injective feature map
and its selftensorproduct
, which is the case in the graphical model examples in Section 5.2.Mean element and expectation.
We consider a fixed probability distribution
over . Following Chen et al. (2010), we denote by the mean element (see, e.g., Smola et al., 2007):Note that in the learning perspective, is the empirical distribution on the data and so is the corresponding empirical moment vector to match. To approximate this mean, we consider points combined linearly with positive weights that sum to one. These define , the associated weighted empirical distribution, and the approximating mean:
(2) 
For all functions , we then have
and similarly . We thus get, using CauchySchwarz inequality,
and controlling is enough to control the error in computing the expectation for all with finite norm. Note that a random i.i.d. sample from would lead to an expected worstcase error which is less than —a classical result based on Rademacher averages (see, e.g. Boucheron et al., 2005).
In this paper, we will try to find a good estimate of based on a weighted set of points from , generalizing Chen et al. (2010), and show how this relates to herding.
3 Related work
This paper brings together three lines of work, namely the approximation of integrals, herding and convex optimization. The links between the first two were clearly outlined by Chen et al. (2010), while the present paper provides the novel interpretation of herding as a wellestablished convex optimization algorithm.
3.1 Quadrature/cubature formulas
The evaluation of expectation, or equivalently of integrals, is a classical problem in numerical analysis. When the input space is a compact subset of and is the density of the distribution with respect to the Lebesgue measure, then this is equivalent to evaluating the integral . Quadrature formulas are aimed at computing such integrals as a weighted combinations of values of at certain points, which is exactly the problem we consider in Section 2.
Although a thorough review of quadrature formulas is outside of the scope of this paper, we mention two methods which are related to our approach. First, given a positive definite kernel and a given set of points (typically sampled i.i.d. from a given distribution), the BayesHermite quadrature of O’Hagan (1991) essentially computes an orthogonal projection of onto the affine hull of this set of points. This does not lead to positive quadrature weights, and one may thus replace the affine hull by the convex hull to obtain such nonnegative weights, which we do in our experiments in Section 5.
Moreover, quasiMonte Carlo methods consider sequences of socalled “quasirandom” quadrature points so that the empirical average approaches the integral. These quasirandom sequences are such that the approximation error goes down as (up to logarithmic terms) for functions of bounded variation, as opposed to for a random sequence. In simulations, we used a Sobol sequence (see, e.g., Morokoff & Caflisch, 1994).
3.2 FrankeWolfe algorithms
Given a smooth (twice continuously differentiable) convex function on a compact convex set with gradient , FrankWolfe algorithms are a class of iterative optimization algorithms that only require (in addition to the computation of the gradient ) to be able to optimize linear functions on . The first class of such algorithms is often referred to as conditional gradient algorithms: given an iterate , the minimum of over is computed, and the next iterate is taken on the segment between and , i.e., , where . There are two natural choices for , (a) simply taking and (b) performing a line search to find the point in the segment with optimal value (either for or for a quadratic upperbound of ). These are illustrated in Figure 2, and convergence rates are detailed in Section 4.2. Moreover, for quadratic functions, the conditional gradient algorithm with step sizes has a dual interpretation as subgradient ascent (see, e.g., Bach, 2011), which we outline in Section 4.1.
Finally, in order to minimize the number of iterations, a variant known as the minimumnormpoint algorithm will find that minimizes on the convex hull of all previously visited points, using a specific activeset algorithm (see Bach, 2011, Sec. 6, for details). For convex sets with finitely many extreme points, it converges in a finite number of iterations with higher (though still polynomial) iteration computational cost (Wolfe, 1976).
4 Herding as a FrankWolfe algorithm
To relate herding with conditional gradient algorithms, we consider the following optimization problem:
(3) 
with the trivial unique solution . A conditional gradient algorithm to solve this optimization problem with stepsize use the iterates
(4) 
But these updates are exactly the same as herding via the change of variable . Indeed, the minimizer of a linear function of a convex set can be restricted to be an extreme point of ; this implies that for a certain . The herding updates in Eq. (1) are thus equivalent to the conditional gradient minimization of with step size given by .
With this choice of step size, we get , that is and we thus get uniform weights in Eq. (2).
For general stepsizes , if we assume that we start the algorithm with (which implies ), then we get that which thus leads to nonuniform weights in Eq. (2), though they still sum to one. The conditional gradient updates in Eq. (4) can thus be seen as a generic algorithm to obtain a weighted set of points to approximate , and traditional herding is the stepsize case.
A second choice of stepsize for is to use a line search, which leads in this setting (where the global unconstrained minimum happens to belong to ) to This leads to a variant of herding with nonuniform weights.
We finally comment on the initialization for the updates in Eq. (4). In the kernel herding algorithm of Chen et al. (2010), the authors use as this is required to interpret the herding updates as greedily minimizing (with the additional assumptions that is constant). In our setting, this corresponds to choosing (which might be outside of , though this is not problematic in practice). Another standard choice (for MRFs in particular) is to use (), which means that the first point is chosen randomly from the extreme points of —this is the scheme we used. As is common in convex optimization, we didn’t see any qualitative difference in our experiments between the two types of initialization.
4.1 Dual problem and subgradient descent
Welling (2009b) proposed originally an algorithmic interpretation of herding as performing subgradient ascent with constant step size on a function obtained as the zero temperature limit of the loglikelihood of an exponential model that he called the “tipi function”. Our interpretation of the procedure as a FrankWolfe algorithm might therefore appear as a conflicting interpretation at first sight. To establish a natural connection between these two interpretations, we can compute the Fencheldual optimization problem to Eq. (3). Indeed, we have (with standard arguments for swapping the min and max operations):
The dual function is strongly concave and nondifferentiable, and a natural algorithm to maximize it is thus subgradient ascent with a step size equal to , which is known to be equivalent to the primal conditional gradient algorithm with step sizes (Bach, 2011, App. A). It is therefore not surprising that herding is equivalent to subgradient ascent with a decreasing stepsize on this function (with the identification ). The presence of the squared norm which is added to the “tipi function” merely reflects the change of scaling between and . It is worthwhile mentioning that other functions, like Bregman divergences, would have led to a different dual function hence adding a different term than a squared norm to the “tipi function”.
4.2 Convergence analysis
Without further assumptions on the problem, then the two choices of step sizes lead to a convergence rate of the form (Dunn, 1980; Bach, 2011):
where is diameter of the marginal polytope. Note that the convergence in does not lead to improved estimation of integrals over random sampling. Moreover, without further assumptions, current theoretical analysis does not allow to distinguish between the two forms of conditional gradient algorithms (although they differ a bit in practice, see Section 5).
However, if we assume that within the affine hull of , there exists a ball of center and radius that is included in (i.e., is in the relative interior of ), then both choices of step sizes yield faster rates than random sampling. For the version with line search, we actually obtain a linear convergence rate (Beck & Teboulle, 2004):
For the version without line search (), Chen et al. (2010) shows the slower convergence rate in :
Concerning the assumption that is in the relative interior of , we now show that in finitedimensional settings, this assumption is always satisfied under reasonable conditions, while it is never satisfied in a large class of infinitedimensional settings (namely for Mercer kernels).
We first provide an equivalent definition of this condition which is used in the proofs. Let be the affine hull of , the orthogonal projection of on , and the space of directions (or difference space) of (i.e., ).^{1}^{1}1It turns out that is a function taking a constant value since the orthogonality condition yields , i.e., for all , by the reproducing property of . Now there exists so that any element of which is at distance less than of is in if and only if , the maximum of over and is less than the maximum of over . Given the properties of and , this is equivalent to:
(5)  
Proposition 1
Assume that is finitedimensional, that is a compact topological measurable space with a continuous kernel function, and that the distribution has full support on . Then so that Eq. (5) is true.
Proof sketch.
It is sufficient to show that is a norm on : as all norms are equivalent in finite dimensions, we get for some , yielding Eq. (5). is convex and positively homogeneous by construction. Now implies that , and thus for in the support of (assumed to be ) using the fact that is continuous (since the kernel is continuous), and so is a constant function. We then have two possibilities: either , in which case one can show that there is no nonzero constant functions in ; otherwise for some and thus the orthogonality condition implies that . Both cases imply , hence is a norm.
Proposition 2
Assume is a compact subspace of , and that the kernel is a continuous function on . If is infinitedimensional, then there exists no so that Eq. (5) is true.
Proof sketch. We can apply Mercer theorem to the kernel of the projection onto the orthogonal of . This kernel is also a Mercer kernel, and we get an orthonormal basis of with countably manyeigenvalues that are summable. Moreover, is known to be an orthonormal basis of the associated feature space (Cucker & Smale, 2002), and for all , , with uniform convergence. This implies that for , we have , and .
If we assume that there exists so that Eq. (5) is true, then we have for all ,
. Since is compact, there exists a cover of with finitely many balls of radius . Let be the finite set of centers.
Since all functions are Lipschitzcontinuous with constant , then for all ,
. Since is finite, there exists so that for infinitely many values of ; this contradicts the summability of . Hence the result.
The last proposition shows that the current theory only supports the slower rates of for the two conditional gradient algorithms in infinitedimensional settings. On the other hand, we note that, in some cases, traditional herding performs empirically better without known theoretical justification (see Section 5).
4.3 Computational issues
In order to run a herding algorithm, there are two potential computational bottlenecks:
Computing expectations : in a learning context (empirical moment matching), these are done through an empirical average. In an integral evaluation context, in finitedimensional settings, one needs to compute ; while in an infinitedimensional setting, following Chen et al. (2010), expectations of the form need to be computed. This can be done for some pairs of kernels/distributions, like the ones we choose in Section 5, but not in general.
Minimizing with respect to : in general, this computation may be relatively hard (it is for example NPhard in the context of the graphical models we consider in Section 5). In practice, Chen et al. (2010) and Welling (2009a) perform local search, while another possibility is to perform the minimization through exhaustive search in a finite sample. Note that a convex relaxation through variational methods (Wainwright & Jordan, 2008) could provide an interesting alternative.
5 Experiments
The goals of these simulations are (a) to compare the different algorithms aimed at estimating integrals, i.e., assess herding for the task of mean estimation (Section 5.1 and Section 5.2), and (b) to briefly study the empirical relationship with maximum entropy estimation in a learning context (Section 5.3).
5.1 Kernel herding on
Problem setup. In this section, we consider and the norm on the infinitedimensional space of functions with zero mean and whose th derivative exists and is in . As shown by Wahba (1990), the associated kernel is equal to , where is the th Bernoulli polynomial, with and .
We consider either the uniform density on , or a randomly selected density of the form , for which all required expectations may be computed in closed form. In particular, the mean element is computed as which may be computed in closed form by expanding all terms in the Fourier basis. As for the optimization step, it consists in minimizing over the interval which can be done efficiently with exhaustive search.
Comparison of mean estimation procedures. We compare in Figure 3 two kernels, i.e., with (left and middle plots) and (right plot), the following mean estimation procedures, and plot , for two densities, the uniform density (middle and right) and a randomly selected nonuniform density (left). We compare the following:

cg1/(t+1): conditional gradient procedure with , which is the original herding procedure of Welling (2009a), leading to uniform weights.

cgl.search: conditional gradient procedure with line search (with nonuniform weights).

minnormpoint: Minimumnorm point algorithm. This leads to nonuniform weights.

random: Random selection of points (from ), averaged over 10 replications.

sobol
: a classical quasirandom sequence, with uniform weights. For nonuniform densities, we first apply the inverse cumulative distribution function.
For all of these (except for minnormpoint), we also consider an extra a posteriori projection step (denoted by the proj suffix), which computes the optimal set of nonuniform weights by finding the best approximation of in the convex hull of the points selected by the algorithm. We can draw the following conclusions:

Minnormpoint algorithms always perform best.

Conditional gradient with line search is performing slightly worse than regular herding. (Note that we are in the infinitedimensional setting and so they both have as theoretical rate.)

The extra projection step always significantly improves performance, and sometimes enough that random selection of point combined with a reprojection outperforms regular herding (at least for , i.e., with a smaller feature space).

On the right plot, it turns out that the Sobol sequence is known to achieve the optimal rate of for for the associated Sobolev space (Wahba, 1990). In this situation, regular herding empirically achieves the same rate; however, the theoretical analysis provided in the present paper or by Chen et al. (2010) does not allow to explain or support this statement theoretically.
Estimation from a finite sample. In Figure 4, we compare three of the previously mentioned herding procedures when all quantities are computed from a random sample of size . In plain, testing errors are computed (using exact expectations) while in dashed, the training errors are computed. All methods eventually fit the empirical mean, with no further progress on the testing error, this behavior happening faster with the minnormpoint algorithm.
5.2 Estimation on graphical models
We consider
and a random variable computed as the sign (in
) of a Gaussian random vector in , together with composed of and of all of its pairwise products . In this setup, we can compute the expectationin closed form, as the mass assigned to the positive orthant by a bivariate Gaussian distribution with correlation
, which is equal to (Abramowitz & Stegun, 1964). We are here in the finitedimensional setting and the faster rates derived in Section 4.2 apply.We generated samples from such a distribution and performed herding with exact expectations but with minima with respect to computed over this sample (by exhaustive search over the sample). We plot results in Figure 5, where we see the superiority of the minnormpoint procedure over the two other procedures (which include regular herding). Note that the linesearch algorithm is slower than the rule, which seems to contradict the bounds. The bounds depend on the distance between the mean and the boundary of the marginal polytope. If this is too small (much like if the strong convexity constant is too small for gradient descent), the linear convergence rate can only be seen for larger numbers of iterations.
5.3 Herding and maximum/minimum entropy
Given a moment vector obtained from the empirical mean of on data, the goal of herding is to produce a pseudosample whose moments match without having to estimate the canonical parameters of the corresponding model. A natural candidate for such a distribution is the maximum entropy distribution and we will compare it to the results of herding in cases where it can be easily computed, namely for (with ) and either or . In this setup, following Welling (2009a), the distribution on is estimated by the empirical distribution .
Learning independent bits. We first consider and some specific feasible moment . It is wellknown that the maximum entropy distribution is the one with independent bits. In the top panels of Figure 6, we compare the norm between the maximum entropy probability vector and the one estimated by two versions of herding, namely conditional gradient with stepsize (regular herding with uniform weights) and with line search (with nonuniform weights)—the minnormpoint algorithm leads to quantitatively similar results. We show results in Figure 6 for a mean vector which is a random uniform vector in (left plots), and for a mean which is random with uniform values in (middle plots), and for mean values which is are random with uniform values in (right plots). Note that the difference between rational and irrational means was already brought up by Welling & Chen (2010) through the link between herding and Sturmian sequences.
For each of the mean vector, we plot in the top plots, the error in estimating the full maximum entropy distribution (a vector of size ), and in the bottom plots, the error in estimating the feature means (a vector of size ). We can draw the following conclusions:

For a random vector (left plots), then regular herding (with no line search) empirically converges to the maximum entropy distribution.

For rational ratios between the means (but irrational means, middle plots), then there is no convergence to the maximum entropy distribution.

For rational means (right), there is no convergence either, but the behavior is more erratic.

The linesearch procedure never converges to the maximum entropy procedure. On the opposite, it happens to be close to a minimum entropy solution, where many states have probability zero.
Experiments considered in Figure 6 considered a single draw of the mean vector , but similar empirical conclusions may be drawn from other random samples, and we conjecture that for almost surely all random vectors (which would avoid rational ratios between mean values), then regular herding converges to the maximum entropy distribution. The next experiment shows that this is not the case in general.
Learning nonindependent bits. We now consider , and a certain random feasible moment . As before, we compare the norm between the maximum entropy probability vector and the one estimated by the two versions of herding. We present results in Figure 7 for a mean vector obtained by the corresponding exponential family distribution with zero weights for the features and constant weights on the features . We see that the herding procedures, while leading to a consistent estimation of the mean vector, does not converge to the maximum entropy distribution and other unreported experiments have led to similar results.
6 Conclusion
We showed that herding generates a sequence of points which give in sequence the descent directions of a conditional gradient algorithm minimizing the quadratic error on the moment vector. Therefore, if herding is only assessed in terms of its ability to approximate the moment vector, it is outperformed by other more efficient algorithms. Clearly, herding was originally defined with another goal, which was to generate a pseudosample whose distribution could approach the maximum entropy distribution with a given moment vector. Our experiments suggest empirically, that while this is the case in certain cases, herding fails in other case, which are not chosen to be particularly pathological. This probably prompts for a further study of herding.
Our experiments also show that algorithms that are more efficient than herding at approximating the moment vector fail more blatantly to approach a maximum entropy distribution and they present characteristics which would rather suggest a minimization of the entropy. This suggests the question of whether there is a tradeoff between approximating most efficiently the mean vector and approximating well the maximum entropy distribution, or if the two goals are in fact rather aligned? In any case, we hope that formulating herding as an optimization problem can help form a better understanding of its goals and its properties.
Acknowledgements. We thank Ferenc Huszár and Zoubin Ghahramani for helpful discussions. This work was supported by the European Research Council (SIERRA Project) and the city of Paris (“Research in Paris” program).
References
 Abramowitz & Stegun (1964) Abramowitz, M. and Stegun, I.A. Handbook of mathematical functions. Dover publications, 1964.
 Bach (2011) Bach, F. Learning with submodular functions: A convex optimization perspective. Technical Report 1111.6453, Arxiv, 2011.
 Beck & Teboulle (2004) Beck, A. and Teboulle, M. A conditional gradient method with linear rate of convergence for solving convex linear systems. Math. Meth. Op. Res., 59(2):235–247, 2004.
 Boucheron et al. (2005) Boucheron, S., Bousquet, O., and Lugosi, G. Theory of classification: A survey of some recent advances. ESAIM Probability and statistics, 9:323–375, 2005.
 Chen et al. (2010) Chen, Y., Welling, M., and Smola, A. Supersamples from kernel herding. In Proc. UAI, 2010.
 Cucker & Smale (2002) Cucker, F. and Smale, S. On the mathematical foundations of learning. Bull. AMS, 39(1), 2002.
 Dunn (1980) Dunn, J. C. Convergence rates for conditional gradient sequences generated by implicit step length rules. SIAM J. Control & Opt., 18:473–487, 1980.

Gelfand et al. (2010)
Gelfand, A., van der Maaten, L., Chen, Y., and Welling, M.
On herding and the perceptron cycling theorem.
In Adv. NIPS, 2010.  Morokoff & Caflisch (1994) Morokoff, W.J. and Caflisch, R.E. Quasirandom sequences and their discrepancies. SIAM Journal on Scientific Computing, 15(6):1251–1279, 1994.
 O’Hagan (1991) O’Hagan, A. BayesHermite quadrature. J. Stat. Planning & Inference, 29(3):245–260, 1991.
 Smola et al. (2007) Smola, A., Gretton, A., Song, L., and Schölkopf, B. A Hilbert space embedding for distributions. In Algorithmic Learning Theory, pp. 13–31. Springer, 2007.
 Wahba (1990) Wahba, Grace. Spline Models for Observational Data. SIAM, 1990.
 Wainwright & Jordan (2008) Wainwright, M.J. and Jordan, M.I. Graphical models, exponential families, and variational inference. Found. Trends Mach. Learn., 1(12):1–305, 2008.
 Welling (2009a) Welling, M. Herding dynamic weights for partially observed random field models. In Proc. UAI, 2009a.
 Welling (2009b) Welling, M. Herding dynamical weights to learn. In Proc. ICML, 2009b.
 Welling & Chen (2010) Welling, M. and Chen, Y. Statistical inference using weak chaos and infinite memory. J. Physics: Conf. Series, 233, 2010.
 Wolfe (1976) Wolfe, P. Finding the nearest point in a polytope. Math. Progr., 11(1):128–149, 1976.