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 generatespseudo-samples in a deterministic fashion with the property of asymptotically matching the empirical moments of the data (akin to maximum entropy). The herding algorithm generates pseudo-samples with the following simple recursion:
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, andis 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 fixed-step-size subgradient ascent on the zero-temperature 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 pseudo-samples 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 auto-correlations amongst the pseudo-samples and was used by Chen et al. (2010) to sub-select a small collection of representative “super-samples” 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. Frank-Wolfe 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 non-uniform weights, one based on line-search, one based on an active-set 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 line-search variant in finite-dimensional settings and show how the conditions assumed by Chen et al. (2010) in fact never hold in the infinite-dimensional 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 low-dimensional 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 well-defined.
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 infinite-dimensional kernels on in Section 5.1); it is also true if contains both an injective feature map
and its self-tensor-product, which is the case in the graphical model examples in Section 5.2.
Mean element and expectation.
We consider a fixed probability distributionover . 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:
For all functions , we then have
and similarly . We thus get, using Cauchy-Schwarz 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 worst-case 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 well-established 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 Bayes-Hermite 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, quasi-Monte Carlo methods consider sequences of so-called “quasi-random” quadrature points so that the empirical average approaches the integral. These quasi-random 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 Franke-Wolfe algorithms
Given a smooth (twice continuously differentiable) convex function on a compact convex set with gradient , Frank-Wolfe 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 upper-bound 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 minimum-norm-point algorithm will find that minimizes on the convex hull of all previously visited points, using a specific active-set 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 Frank-Wolfe algorithm
To relate herding with conditional gradient algorithms, we consider the following optimization problem:
with the trivial unique solution . A conditional gradient algorithm to solve this optimization problem with stepsize use the iterates
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 step-sizes , if we assume that we start the algorithm with (which implies ), then we get that which thus leads to non-uniform 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 step-size case.
A second choice of step-size 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 non-uniform 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 log-likelihood of an exponential model that he called the “tipi function”. Our interpretation of the procedure as a Frank-Wolfe algorithm might therefore appear as a conflicting interpretation at first sight. To establish a natural connection between these two interpretations, we can compute the Fenchel-dual 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 non-differentiable, 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
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 finite-dimensional settings, this assumption is always satisfied under reasonable conditions, while it is never satisfied in a large class of infinite-dimensional 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., ).111It 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:
Assume that is finite-dimensional, 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.
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 non-zero constant functions in ; otherwise for some and thus the orthogonality condition implies that . Both cases imply , hence is a norm.
Assume is a compact subspace of , and that the kernel is a continuous function on . If is infinite-dimensional, 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 Lipschitz-continuous 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 infinite-dimensional 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 finite-dimensional settings, one needs to compute ; while in an infinite-dimensional 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 NP-hard 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.
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 set-up. In this section, we consider and the norm on the infinite-dimensional 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 non-uniform density (left). We compare the following:
cg-1/(t+1): conditional gradient procedure with , which is the original herding procedure of Welling (2009a), leading to uniform weights.
cg-l.search: conditional gradient procedure with line search (with non-uniform weights).
min-norm-point: Minimum-norm point algorithm. This leads to non-uniform weights.
random: Random selection of points (from ), averaged over 10 replications.
: a classical quasi-random sequence, with uniform weights. For non-uniform densities, we first apply the inverse cumulative distribution function.
For all of these (except for min-norm-point), we also consider an extra a posteriori projection step (denoted by the -proj suffix), which computes the optimal set of non-uniform weights by finding the best approximation of in the convex hull of the points selected by the algorithm. We can draw the following conclusions:
Min-norm-point algorithms always perform best.
Conditional gradient with line search is performing slightly worse than regular herding. (Note that we are in the infinite-dimensional 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 min-norm-point algorithm.
5.2 Estimation on graphical models
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 set-up, we can compute the expectation
in 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 finite-dimensional 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 min-norm-point procedure over the two other procedures (which include regular herding). Note that the line-search 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 pseudo-sample 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 well-known 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 non-uniform weights)—the min-norm-point 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 line-search 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 non-independent 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.
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 pseudo-sample 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).
- 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. Super-samples 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. Quasi-random sequences and their discrepancies. SIAM Journal on Scientific Computing, 15(6):1251–1279, 1994.
- O’Hagan (1991) O’Hagan, A. Bayes-Hermite 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(1-2):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.