Perturbative Corrections for Approximate Inference in Gaussian Latent Variable Models

2013·Arxiv

Abstract

Abstract

Expectation Propagation (EP) provides a framework for approximate inference. When the model under consideration is over a latent Gaussian field, with the approximation being Gaussian, we show how these approximations can systematically be corrected. A perturbative expansion is made of the exact but intractable correction, and can be applied to the model’s partition function and other moments of interest. The correction is expressed over the higher-order cumulants which are neglected by EP’s local matching of moments. Through the expansion, we see that EP is correct to first order. By considering higher orders, corrections of increasing polynomial complexity can be applied to the approximation. The second order provides a correction in quadratic time, which we apply to an array of Gaussian process and Ising models. The corrections generalize to arbitrarily complex approximating families, which we illustrate on tree-structured Ising model approximations. Furthermore, they provide a polynomial-time assessment of the approximation error. We also provide both theoretical and practical insights on the exactness of the EP solution.

Keywords: expectation consistent inference, expectation propagation, perturbation correction, Wick expansions, Ising model, Gaussian process

1. Introduction

Expectation Propagation (EP) (Opper and Winther, 2000, Minka, 2001a,b) is part of a rich family of variational methods, which approximate the sums and integrals required for exact probabilistic inference by an optimization problem. Variational methods are perfectly amenable to probabilistic graphical models, as the nature of the optimization problem often allows it to be distributed across a graph. By relying on local computations on a graph, inference in very large probabilistic models becomes feasible.

Being an approximation, some error may invariably be introduced. This paper is specifically concerned with the error that arises when a Gaussian approximating family is used, and lays a systematic foundation for examining and correcting these errors. It follows on earlier work by the authors (Opper et al., 2009). The error that arises when the free energy (the negative logarithm of the partition function or normalizer of the distribution) is approximated, may for instance be written as a Taylor expansion (Opper et al., 2009, Paquet et al., 2009). A pleasing property of EP is that, at its stationary point, the first order term of such an expansion is zero. Furthermore, the quality of the approximation can then be ascertained in polynomial time by including corrections beyond the first order, or beyond the standard EP solution. In general, the corrections improve the approximation when they are comparatively small, but can also leave a question mark on the quality of approximation when the lower-order terms are large.

The approach outlined here is by no means unique in correcting the approximation, as is evinced by cluster-based expansions (Paquet et al., 2009), marginal corrections for EP (Cseke and Heskes, 2011) and the Laplace approximation (Rue et al., 2009), and corrections to Loopy Belief Propagation (Chertkov and Chernyak, 2006, Sudderth et al., 2008, Welling et al., 2012).

1.1 Overview

EP is introduced in a general way in Section 3, making it clear how various degrees of complexity can be included in its approximating structure. The partition function will be used throughout the paper to explain the necessary machinery for correcting any moments of interest. In the experiments, corrections to the marginal and predictive means and variances are also shown, although the technical details for correcting moments beyond the partition function are relegated to Appendix D. The Ising model, which is cast as a Gaussian latent variable model in Section 2, will furthermore be used as a running example throughout the paper.

The key to obtaining a correction lies in isolating the “intractable quantity” from the “tractable part” (or EP solution) in the true problem. This is done by considering the cumulants of both: as EP locally matches lower-order cumulants like means and variances, the “intractable part” exists as an expression over the higher-order cumulants which are neglected by EP. This process is outlined in Section 4, which concludes with two useful results: a shift of the “intractable part” to be an average over complex Gaussian variables with zero diagonal relation matrix, and Wick’s theorem, which allows us to evaluate the expectations of polynomials under centered Gaussian measures. As a last stage, the “intractable part” is expanded in Sections 5 and 7 to obtain corrections to various orders. In Section 6, we provide a theoretical analysis of the radius of convergence of these expansions.

Experimental evidence is presented in Section 8 on Gaussian process (GP) classification and (non-Gaussian) GP regression models. An insightful counterexample where EP diverges under increasing data, is also presented. Ising models are examined in Section 9.

Numerous additional examples, derivations, and material are provided in the appendices. Details on different EP approximations can be found in Appendix A, while corrections to tree-structured approximations are provided in Appendix B. In Appendix C we analytically show that the correction to a tractable example is zero. The main body of the paper deals with corrections to the partition function, while corrections to marginal moments are left to Appendix D. Finally, useful calculations of certain cumulants appear in Appendix E.

2. Gaussian Latent Variable Models

Let ) be an unobserved random variable with an intractable distribution p(x). In the Gaussian latent variable model (GLVM) considered in this paper, terms are combined over a quadratic exponential

with partition function (normalizer)

This model encapsulates many important methods used in statistical inference. As an example, can encode the covariance matrix of a Gaussian process (GP) prior on latent function observations . In the case of GP classification with a class label on a latent function evaluation , the terms are typically probit link functions, for example

The probit function is the standard cumulative Gaussian density Φ(In this example, the partition function is not analytically tractable but for the one-dimensional case N = 1.

An Ising model can be constructed by letting the terms Dirac delta functions). By introducing the symmetric coupling matrix , an Ising model can be written as

In the Ising model, the partition function Z is intractable, as it sums values of x. In the variational approaches, the intractability is addressed by allowing approximations to Z and other marginal distributions, decreasing the computational complexity from being exponential to polynomial in N, which is typically cubic for EP.

3. Expectation Propagation

An approximation to Z can be made by allowing p(x) in Equation (1) to factorize into a product of . This factorization is not unique, and the structure of the factorization of p(x) defines the complexity of the resulting approximation, resulting in different structures in the approximating distribution. Where GLVMs are concerned, a natural and computationally convenient choice is to use Gaussian factors , and as such, the approximating distribution q(x) in this paper will be Gaussian. Appendix A summarizes a number of factorizations for Gaussian approximations.

The tractability of the resulting inference method imposes a pragmatic constraint on the choice of factorization; in the extreme case p(x) could be chosen as a single factor and inference would be exact. For the model in Equation (1), a three-term product may be factorized as (), which gives the typical GP setup. When a division is introduced and the term product factorizes as (), the resulting free energy will be that of the tree-structured EC approximation (Opper and Winther, 2005). To therefore allow for regrouping, combining, splitting, and dividing terms, a power is associated with each , such that

with intractable normalization (or partition function) Appendix A shows how the introduction of lends itself to a clear definition of tree-structured and more complex approximations.

To define an approximation to , which typically take an exponential family form, are chosen such that

has the same structure as p’s factorization. Although not shown explicitly, a dependence on the same subset of variables . The optimal parameters of the approximations are found through a set of auxiliary tilted distributions, defined by

Here a single approximating term is replaced by an original term . Assuming that this replacement leaves still tractable, the parameters in are determined by the condition that ) should be made as similar as possible. This is usually achieved by requiring that these distributions share a set of generalised moments which usually coincide with the sufficient statistics of the exponential family. For example with sufficient statistics ) we require that

Note that those factors ) which are already in the exponential family, such as the Gaussian terms in examples above, can trivially be solved for by setting partition function associated with this approximation is

Appendix A.2 shows that the moment-matching conditions must hold at a stationary point of log . The EP algorithm iteratively updates the -terms by enforcing q to share moments with each of the tilted distributions ; on reaching a fixed point all moments match according to Equation (7) (Minka, 2001a,b). Although is defined in the terminology of EP, other algorithms may be required to solve for the fixed point, and , as a free energy, can be derived from the saddle point of a set of self-consistent (moment-matching) equations (Opper and Winther, 2005, van Gerven et al., 2010, Seeger and Nickisch, 2010). We next make EP concrete by applying it to the Ising model, which will serve as a running example in the paper. The section is finally concluded with a discussion of the interpretation of EP.

3.1 EP for Ising Models

The Ising model in Equation (3) will be used as a running example throughout this paper. To make the technical developments more concrete, we will consider both the N-variate and bivariate cases. The bivariate case can be solved analytically, and thus allows for a direct comparison to be made between the exact and approximate solutions.

We use the factorized approximation as a running example, dividing p(x) in Equation (3) into N + 1 factors with 1) + (see Appendix A for generalizations). We consider the Gaussian exponential family such that The approximating distribution from Equation (5), ), is thus a full multivariate Gaussian density, which we write as

3.1.1 Moment Matching

The moment matching condition in Equation (7) involves only the mean and variance if q(x) fully factorizes according to We therefore only need to match the mean and variances of marginals of q(x) and the tilted distribution ) in Equation (6). The tilted distribution may be decomposed into a Gaussian and a discrete part as ), where the vector consists of all variables apart from may marginalize out ) in terms of two factors:

where we dropped the dependency of for notational simplicity. Through some manipulation, the tilted distribution is equivalent to

This discrete distribution has mean and variance 1 . By adapting the parameters of ) using for example the EP algorithm, we aim to match the mean and variance of the marginal )) to the mean and variance of ). The reader is referred to Section 9 for benchmarked results for the Ising model.

3.1.2 Analytic Bivariate Case

Here we shall compare the exact result with EP and the correction for the simplest non-trivial model, the N = 2 Ising model with no external field

In order to solve the moment matching conditions we observe that the mean values must be zero because the distribution is symmetric around zero. Likewise the linear term in the approximating factors disappears and we can write

, turns into a second order equation with solution now insert this solution into the expression for the EP partition function in Equation (8). By expanding the result to the second order in , we find that

Comparing with the exact expression

we see that EP gives the correct coefficient, but the coefficient comes out wrong. In Section 4 we investigate how cumulant corrections can correct for this discrepancy.

3.2 Two Explanations Why Gaussian EP is Often Very Accurate

EP, as introduced above, is an algorithm. The justification for the algorithm put forward by Minka and adopted by others (see for example recent textbooks by Bishop 2006, Barber 2012 and Murphy 2012) is useful for explaining the steps in the algorithm but may be misleading in order to explain why EP often provides excellent accuracy in estimation of marginal moments and Z.

The general justification for EP (Minka, 2001a,b) is based upon a minimization of Kullback-Leiber (KL) divergences. Ideally, one would determine the approximating distribution q(x) as the minimizer of KL() in an exponential family of (in our case, Gaussian) densities. Since this is not possible—it would require the computation of exact moments— we instead iteratively minimize “local” KL-divergences KL(), between the tilted distribution , with respect to (appearing in q). This leads to the moment matching conditions in Equation (7). The argument for this procedure is essentially that this will ensure that the approximation q will capture high density regions of the intractable posterior p. Obviously, this argument cannot be applied to Ising models because the exact and approximate distributions are very different, with the former being discrete due to the Dirac -functions that constrain 1 to be binary variables. Even though the optimization still implies moment matching, this discrete-continuous discrepancy makes local KL-divergences KL(

In order to justify the usefulness of EP for Ising models we therefore need an alternative argument. Our argument is entirely restricted to Gaussian EP for our extended definition of GLVMs and do not extend to approximations with other exponential families. In the following, we will discuss these assumptions in inference approximations that preceded the formulation of EP, in order to provide a possibly more relevant justification of the method. Although this justification is not strictly necessary for practically using EP nor corrections to EP, it nevertheless provides a good starting point for understanding both.

The argument goes back to the mathematical analysis of the Sherrington-Kirkpatrick (SK) model for a disordered magnet (a so-called spin glass) (Sherrington and Kirckpatrick, 1975). For this Ising model, the couplings J are drawn at random from a Gaussian distribution. An important contribution in the context of inference for this model (the computations of partition functions and average magnetizations) was the work of Thouless et al. (1977) who derived self-consistency equations which are assumed to be valid with a probability (with respect to the drawing of random couplings) approaching one as the number of variables grows to infinity. These so-called Thouless-Anderson-Palmer (TAP) equations are closely related to the EP moment matching conditions of Equation (7), but they differ by partly relying on the specific assumption of the randomness of the couplings. Selfconsistency equations equivalent to the EP moment matching conditions which avoided such assumptions on the statistics of the random couplings were first derived by Opper and Winther (2000) by using a so-called cavity argument (M´ezard et al., 1987). A new important contribution of Minka (2001a) was to provide an efficient algorithmic recipe for solving these equations.

We will now sketch the main idea of the cavity argument for the GLVM. Let without n”) denote the complement to . Without loss of generality we will take the quadratic exponential term to be written as similar definitions of , the exact marginal distribution of may be written as

It is clear that ) depends entirely on the statistics of the random variable This is the total created by all other ‘magnetic moments’ the ‘cavity’ opened once has been removed from the system. In the context of densely connected models with weak couplings, we can appeal to the central limit theoremapproximate by a Gaussian random variable with mean and variance looking at the influence of the remaining variables , the non-Gaussian details of their distribution have been washed out in the marginalization. Integrating out the Gaussian random variable gives the Gaussian cavity field approximation to the marginal distribution:

This is precisely of the form of the marginal tilted distribution ) of Equation (9) as given by Gaussian EP. In the cavity formulation, q(x) is simply a placeholder for the sufficient statistics of the individual Gaussian cavity fields. So we may observe cases, with the Ising model or bounded support factors being the prime examples, where EP gives essentially correct results for the marginal distributions of the and of the partition function Z, while q(x) gives a poor or even meaningless (in the sense of KL divergences) approximation to the multivariate posterior. Note however, that the entire covariance can be computed simply from a derivative of the free energy (Opper and Winther, 2005) resulting in an approximation of this covariance by that of q(x). This may indicate that a good EP approximation of the free energy may also result in a good approximation to the full covariance. The near exactness of EP (as compared to exhaustive summation) in Section 9 therefore shows the central limit theorem at work. Conversely, mediocre accuracy or even failure of Gaussian EP, as also observed in our simulations in Sections 8.3 and 9, may be attributed to breakdown of the Gaussian cavity field assumption. Exact inference on the strongest couplings as considered for the Ising model in Section 9 is one way to alleviate the shortcoming of the Gaussian cavity field assumption.

4. Corrections to EP

The approximation can be corrected in a principled approach, which traces the following outline:

1. The exact partition function Z is re-written in terms of , scaled by a correction factor . This correction factor R encapsulates the intractability in the model, and contains a “local marginal” contribution by each (see Section 4.1).

2. A “handle” on R is obtained by writing it in terms of the cumulants (to be defined in Section 4.2) of ) from Equations (5) and (6). As share their two first cumulants, the mean and covariance from the moment matching condition in Equation (7), a cumulant expansion of R will be in terms of higher-order cumulants (see Section 4.2).

3. R, defined in terms of cumulant differences, is written as a complex Gaussian average. Each factor contributes a complex random variable in this average (see Section 4.3).

4. Finally, the cumulant differences are used as “small quantities” in a Taylor series expansion of R, and the leading terms are kept (see Sections 5 and 7).

The series expansion is in terms of a complex expectation with a zero “self-relation” matrix, and this has two important consequences. Firstly, it causes all first order terms in the Taylor expansion to disappear, showing that is correct to first order. Secondly, due to Wick’s theorem (introduced in Section 4.4), these zeros will contract the expansion by making many other terms vanish.

The strategy that is presented here can be re-used to correct other quantities of interest, like marginal distributions or the predictive density of new data when p(x) is a Bayesian probabilistic model. These corrections are outlined in Appendix D.

4.1 Exact Expression for Correction

We define the (intractable) correction . We can derive a useful expression for R in a few steps as follows: First we solve for in Equation (6), and substitute this into Equation (4) to obtain

We introduce F(x)

to derive the expression for the correction by integrating Equation (11):

Corrections to the marginal and predictive densities of p(x) can be computed from this formulation. This expression will become especially useful because the terms in F(x) turn out to be “local”, that is, they only depend on the marginals of the variables associated with factor ) depend on the subset the remaining variables. The distributions in Equations (5) and (6) differ only with respect to their marginals on ), and therefore

Now we can rewrite F(x) in terms of marginals:

The key quantity, then, is F, after which the key operation is to compute its expected value. The rest of this section is devoted to the task of obtaining a “handle” on F.

4.2 Characteristic Functions and Cumulants

The distributions present in each of the ratios in F(x) in Equation (14) share their first two cumulants, mean and covariance. Cumulants and cumulant differences are formally defined in the next paragraph. This simple observation has a crucial consequence: As the )’s are Gaussian and do not contain any higher order cumulants (three and above), F can be expressed in terms of the higher cumulants of the ). When the term-product approximation is fully factorized, these are simply cumulants of one-dimensional distributions.

Let be the number of variables in subvector In the examples presented in this work, is one or two. Furthermore, let -dimensional vector (. The characteristic function of

and is obtained through the Fourier transform of the density. Inversely,

The cumulants are the coefficients that appear in the Taylor expansion of log around the zero vector,

By this definition of , the Taylor expansion of log

Some notation was introduced in the above two equations to facilitate manipulating a multivariate series. The vector , denotes a multi-index on the elements of . Other notational conventions that employ instead of are:

For example, when indices 2), and (0, 3).

There are two characteristic functions that come into play in F(x) and R in Equation (13). The first is that of the tilted distribution, log ), and the other is the characteristic function of the EP marginal ), defined as . By virtue of matching the first two moments, and ) being Gaussian with cumulants

contains the remaining higher-order cumulants where the tilted and approximate distributions differ. All our subsequent derivations rest upon moment matching being attained. This especially means that one cannot use the derived corrections if EP has not converged.

4.2.1 Ising Model Example

The cumulant expansion for the discrete distribution in Equation (10) becomes

(we’re compactly writing ), from which the cumulants are obtained as

4.3 The Correction as a Complex Expectation

The expected value of F, which is required for the correction, has a dependence on a product of ratios of distributions In the preceding section it was shown that the contributing distributions share lower-order statistics, allowing a twofold simplification. Firstly, the ratio will be written as a single quantity that depends on , which was introduced above in Equation (17). Secondly, we will show that it is natural to shift integration variables into the complex plane, and rely on complex Gaussian random variables (meaning that both real and imaginary parts are jointly Gaussian). These complex random variables that define the ’s have a peculiar property: they have a zero self-relation matrix! This property has important consequences in the resulting expansion.

4.3.1 Complex Expectations

Assume that ) share the same mean and covariance, and substitute log ) in the definition of in Equation (16) to give

Although the variables have not been introduced as random variables, we find it natural to interpret them as such, because the rules of expectations over Gaussian random variables will be extremely helpful in developing the subsequent expansions. We will therefore write ) as an expectation of exp ) over a density

By substituting log 2 into Equation (18), we see that can be viewed as Gaussian, but not for real random variables! We have to consider Gaussian random variables with a real and an imaginary part with

Figure 1: Equation (20) shifts to the complex plane. In the simplest case the joint density

For the purpose of computing the expectation in Equation (19), is a degenerate complex Gaussian that shifts the coefficients into the complex plane. The expectation of exp ) is therefore taken over Gaussian random variables that have covariance matrix as their (real) covariance! As shorthand, we write

Figure 1 illustrates a simple density ), showing that the imaginary component is a deterministic function of is averaged out of the joint density a circularly symmetric complex Gaussian distribution over It is circularly symmetric as , relation matrix, and covariance matrix(notation k indicates the complex conjugate of k). For the purpose of computing the expected values with Wick’s theorem (following in Section 4.4 below), we only need the relationsfor pairs of factors a and b. All of these will be derived next:

into Equation (12), R is equal to

When x is given, the -variables are independent. However, when they are averaged over -variables become coupled. They are zero-mean complex Gaussians

and are coupled with a zero self-relation matrix! In other words, if expected valuesbetween the variables in the set

Complex Gaussian random variables are additionally characterized bythese expectations are not required for computing and simplifying the expansion of log R in Section 5, and are not needed for the remainder of this paper. Figure 2 illustrates the structure of the resulting relation matrix for two different factorizations of the same distribution. Each factor contributes a variable, such that the tree-structured approximation’s relation matrix will be larger than that of the fully factorized one.

Section 5 shows that when and expanded. In the general case, discussed in Section 7, the inner expectation is first expanded (to treat the powers) before computing an expectation over both cases the expectation will involve polynomials in k-variables. The expected values of Gaussian polynomials can be evaluated with Wick’s theorem.

4.4 Wick’s Theorem

Wick’s theorem provides a useful formula for mixed central moments of Gaussian variables. Let be real or complex centered jointly Gaussian variables, noting that they do not have to be different. Then

where the sum is over all partitions of into disjoint pairs even, then there are (21)!! such partitions.is odd, then there are none, and the expectation in Equation (23) is zero.

Consider the one-dimensional variable Wick’s theorem states that is even, and is odd. In other words, , and so forth.

Figure 2: The relation matrices between for two factorizations of illustration is for , while the bottom illustration is of a tree structure (. The white squares indicate a zero relation matrixwith the diagonal being zero. From the properties of Equation (22) there are additional zeros in the tree structure’s relation matrix, where edge and node factors share variables. The factor is shadowed in grey in the left-hand figures, and can make q(x) densely connected.

5. Factorized Approximations

In the fully factorized approximation, with ), the exact distribution in Equation (13) depends on the ). Following Equa- tion (21), the correction to the free energy

is taken directly over the centered complex-valued Gaussian random variables which have a relations

In the section to follow, all expectations shall be with respect to k, which will be dropped where it is clear from the context.

Thus far, R is re-expressed in terms of site contributions. The expression in Equation (24) is exact, albeit still intractable, and will be treated through a power series expansion. Other quantities of interest, like marginal distributions or moments, can similarly be expressed exactly, and then expanded (see Appendix D).

5.1 Second Order Correction to log R

Assuming that the ’s are small on average with respect to k, Equation (24) is expanded and the lower order terms kept:

log R = log

The simplification in the second line is a result of the variance terms being zero from Equation (25). The single marginal terms also vanish (and hence EP is correct to first order) because both

This result can give us a hint in which situations the corrections are expected to be small:

• Firstly, the could be small for values of where the density of k is not small. For example, under a zero noise Gaussian process classification model, a step function ) times a Gaussian, where the latter often has small variance compared to the mean. Hence, ) should be very close to a Gaussian.

• Secondly, for systems with weakly (posterior) dependent variables we might expect that the log partition function log Z would scale approximately linearly with N, the number of variables. Since terms with m = n vanish in the computation of ln R, there are no corrections that are is sufficiently small as . Hence, the dominant contributions to log Z should already be included in the EP approximation. However, Section 8.3 illustrates an example where this need not be the case.

The expectation , as it appears in Equation (26), is treated by substituting its cumulant expansion ! from Equation (17). Wick’s theorem now plays a pivotal role in evaluating the expectations that appear in the expansion:

The second line above follows from contractions in Wick’s theorem. All the self-pairing terms, when for example one of the ’s is paired with another in Equation (23), are zero because, using Equation (23), has to be paired with some factor , and this is possible only when l = s. Wick’s theorem sums over all pairings, and there are l! ways of pairing a giving the result in Equation (27). Finally, plugging Equation (27) into Equation (26) gives the second order correction

5.1.1 Ising Example Continued

We can now compute the second order log R correction for the N = 2 Ising model example of Section 3.1. The covariance matrix has Σ. The uneven terms in the cumulant expansion derived in Section 4.2.1 disappear because m = 0. The first nontrivial term is therefore . In Section 3.1, we saw that log plus terms of order and higher. To lowest order in J we have Σand thus log which exactly cancels the lowest order error of EP.

5.2 Corrections to Other Quantities

The schema given here is applicable to any other quantity of interest, be it marginal or predictive distributions, or the marginal moments of p(x). The cumulant corrections for the marginal moments are derived in Appendix D; for example, the correction to the marginal mean of an approximation

while the correction to the marginal covariance is

5.3 Edgeworth-Type Expansions

To simplify the expansion of Equation (24), we integrated (combined) degenerate complex Gaussians ) to obtain fully complex Gaussian random variables then relied on

The expectationsals, and this can be employed in an alternative derivation. In particular, one can first make a Taylor expansion of exp ) around zero, giving complex-valued polynomials in When the inner average in Equation (24) is then taken over , a real-valued series of Hermite polynomials in arises. These polynomials are orthogonal under q(x). The series that describes the tilted distribution ) is equal to the product of expansion of polynomials for the higher-cumulant deviation from a Gaussian density. This line of derivation gives an Edgeworth expansion foreach factor’s tilted distribution.

As a second step, Equation (24) couples the product of separate Edgeworth expansions (one for each factor) together by requiring an outer average over q(x). The orthogonality of Hermite polynomials under q(x) now come into play: it allows products of orthogonal polynomials under q(x) to integrate to zero. This is similar to contractions in Wick’s theorem, whereof this work, an example of such a derivation appears in Appendix C.1.

6. Radius of Convergence

We may hope that in practice the low order terms in the cumulant expansions will account already for the dominant contributions. But will such an expansion actually converge when extended to arbitrary orders? While we will leave a more general answer to future research, we can at least give a partial result for the example of the Ising model. Let the diagonal of the covariance matrix of the EP approximation q(x). We prove here that a cumulant expansion for R will converge when the eigenvalues of diagonal values of one—are bounded between zero and two.

In practice we’ve found that even if the largest of these eigenvalues grows with N, the second-order correction gives a remarkable improvement. This, with the results in Figure 6, lead us to believe that the power series expansion is often divergent. It may well be that our expansions are only of an asymptotic type (Boyd, 1999) for which the summation of only a certain number of terms might give an improvement whereas further terms would lead to worse results. It leads to a paradoxical situation, which seems common when interesting functions are computed: On the one hand we may have a series which does not converge, but in many ways is more practical; on the other hand one might obtain an expansion that converges, but only impractically. Quoting George F. Carrier’s rule from Boyd (1999):

Divergent series converge faster than convergent series because they don’t have to converge.

For this, we do not yet have a clear-cut answer.

6.1 A Formal Expression for the Cumulant Expansion to All Orders

To discuss the question when our expansion will converge when extended to arbitrary orders, we introduce a single extra parameter , which controls the strength of the contribution of cumulants. Expanded into a series in powers of , contributions of cumulants of total order l are multiplied by a factor , for example course, at the end of the calculation, we set replacing

in Equation (24). Hence, we define

By working backwards, and expressing everything by the original densities over correction can be written as

where the density ) is a multivariate Gaussian with mean and covariance given by

where . Hence, we see that the expansion in powers of is actually equivalent to an expansion in products of nondiagonal elements of

Noticing that as ) depends on through the density we can see by expressing in terms of eigenvalues and eigenvectors that for any fixed ) is an analytic function of the complex variable z as long as is positive definite. Since

this is equivalent to the condition that the matrix ) is positive definite. Introducing , the eigenvalues of , positive definiteness fails when for the first time 1 + ) is convergent for

Setting z = 1, this is equivalent to the condition

This means that the eigenvalues have to fulfil 0 Unfortunately, we can not conclude from this condition that pointwise convergence of leads to convergence of ) (which is an integral of !). However, in cases where the integral eventually becomes a finite sum, such as the Ising model, pointwise convergence in x leads to convergence of

6.1.1 Ising Model Example

From Section 4.2.1 the tilted distribution for the running example Ising model is

1)], and hence ) is a unit-variance Gaussian,

follows from Equation (31). The arguments of the previous section show that the radius of ) is determined by the condition that the matrix ) is positive definite or the eigenvalues

1 + c, meaning that cumulant expansion for ) is convergent for the N = 2 Ising model. For N > 2, it is easy to show that this is not necessarily true. Consider the ‘isotropic’ Ising model with and zero external field, then Σ1[. The eigenvalues are now 1 + ((the latter with degeneracy 1). For finite c, the largest eigenvalue will scale with N and thus be larger than the upper value of two that would be required for convergence. Scaling with N for the largest eigenvalue of is also observed in the Ising model simulations Section 9.

We conjecture that convergence of the cumulant series for ) also implies convergence of the series for log ) but leave an investigation of this point to future research. We only illustrate this point for the N = 2 Ising model case, where we have the explicit formula

As can be easily seen, an expansion in converges for 1 which gives the same radius of convergence |c| < 1 as for the expansion of R.

7. General Approximations

The general approximations differ from the factorized approximation in that an expansion in terms of expectations under doesn’t immediately arise. Consider R in Equation (21): Its inner expectations are over , and outer expectations are over x. First take the binomial expansion of the inner expectation, and keep it to second order in

Notice that ) can be complex, but , as it appears in the above expansion, is real-valued. Using this result, again expand . The correction to log R, up to second order, is

In the above relation the first-order terms all disappeared as to zero. This is a general case of Equation (26), in which B we show how to use the general result for the case where the factorization is a tree and our factors are edges (pairs) and nodes (single variables).

8. Gaussian Process Results

One of the most important applications of EP is to statistical models with Gaussian process (GP) priors, where x is a latent variable with Gaussian prior distribution with a kernel matrix K as covariance

It is well known that for many models, like GP classification, inference with EP is on par with MCMC ground truth (Kuss and Rasmussen, 2005). Section 8.1 underlines this case, and shows corrections to the partition function on the USPS data set over a range of kernel hyperparameter settings.

A common inference task is to predict the output for previously unseen data. Under a GP regression model, a key quantity is the predictive mean function. The predictive mean is analytically tractable when the latent function is corrupted with Gaussian noise to produce observations This need not be the case; in Section 8.2 we examine the problem of quantized regression, where the noise model is non-Gaussian with sharp discontinuities. We show practically how the corrections transfer to other moments, like the predictive mean. Through it, we arrive at a hypothetical rule of thumb: if the data isn’t “sensible” under the (probabilistic) model of interest, there is no guarantee for EP giving satisfactory inference.

Armed with the rule of thumb, Section 8.3 constructs an insightful counterexample where the EP estimate diverges or is far from ground truth with more data. Divergence in the partition function is manifested in the initial correction terms, giving a test for the approximation accuracy that doesn’t rely on any Monte Carlo ground truth.

8.1 Gaussian Process Classification

The GP classification model arises when we observe N data points with class labels , and model y through a latent function x with a GP prior. The likelihood terms for are assumed to be ), where Φ() denotes the cumulative Normal density.

An extensive MCMC evaluation of EP for GP classification on various data sets was given by Kuss and Rasmussen (2005), showing that the log marginal likelihood of the data can be approximated remarkably well. As shown by Opper et al. (2009), an even more accurate estimation of the approximation error is given by considering the second order correction in Equation (28). For GPC we generally found that the l = 3 term dominates l = 4, and we do not include any higher cumulants here.

Figure 3 illustrates the correction to log R, with l = 3, 4, on the binary subproblem of the USPS 3’s vs. 5’s digits data set, with N = 767. This is the same set-up of Kuss and Rasmussen (2005) and Opper et al. (2009), using the kernel and we refer the reader to both papers for additional and complimentary figures and results. We evaluated Equation (28) on a similar grid of log For the same grid values we obtained Monte Carlo estimates of log Z, and hence log R. The correction, compared to the magnitude of the log Z grids by Kuss and Rasmussen (2005), is remarkably small, and underlines their findings on the accuracy of EP for GPC.

The correction from Equation (28), as computed here, is ), and compares favorably to ) complexity of EP for GPC.

Figure 3: A comparison of log R using a perturbation expansion of Equation (28) against Monte Carlo estimates of log R, using the USPS data set from Kuss and Rasmussen (2005). The second order correction to log R, with l = 3, 4, is used on the left; the right plot uses a Monte Carlo estimate of log R.

8.2 Uniform Noise Regression

We turn our attention to a regression problem, that of learning a latent function x(s) from inputs and matching real-valued observations . A frequent nonparametric treatment assumes that x(s) is a priori drawn from a GP prior with covariance function ), from which a corrupted version y is observed. Analytically tractable inference is no longer possible in this model when the observation noise is non-Gaussian. Some scenarios include that of quantized regression, where is formed by rounding ) to, say, the nearest integer, or where x(s) indicates a robot’s path in a control problem, with conditions to stay within certain “wall” bounds. In these scenarios the latent function reconstructed from by adding sharply discontinuous uniformly random

We now assume an EP approximation ), which can be obtained by using the moment calculations in Appendix E.2. To simplify the exposition of the predictive marginal, we follow the notation of Rasmussen and Williams (2005, Chapter 3) and let (), so that the final EP approximation multiplies into a joint Gaussian N(x ; 0, K).

8.2.1 Making Predictions for New Data

The latent function ) at any new input is obtained by the predictive marginal of ). The marginal )—given below in Equation (34)—is directly obtained from the EP approximation ). However, the correction to its mean, as was given in Equation (29), requires covariances Σ, which are derived here.

Let be a vector containing the covariance function evaluations ). Again following Rasmussen and Williams (2005)’s notation, let ˜be the diagonal matrix containing 1along its diagonal. The EP covariance, on the inclusion of

with . There is no observation associated with in the first line above, and its inclusion has 3. The second line follows by computing matrix partitioned inverses twice on . The joint EP approximation for any new input point is directly obtained as

with the marginal

According to Equation (29), one needs the covariances Σto correct the marginal’s mean; they appear in the last column of in Equation (33). The correction is

The sum over pairs include the added dimension , and thus pairs (The cumulants for this problem, used both for EP and correcting it, are derived in Appendix E.2.

8.2.2 Predictive Corrections

In Figure 4 we investigate the predictive mean correction for two cases, one where the data cannot realistically be expected to appear under the prior, and the other where the prior is reasonable. For , the values of ) are predicted using a GP with squared exponential covariance function

In the first instance, the prior amplitude and lengthscale are deliberately set to values that are too big; in other words, a typical sample from the prior would not match the observed data. We illustrate the posterior marginal ), and using Equations (29) and (30), show visible corrections to its mean and variance.For comparison, Figure 4

Figure 4: Predicting ) with a GP. The “boxed” bars indicate the permissible values; they are linked to observations through the uniform likelihood ]. Due to the ] noise model, ) is ambivalent to where in the “box” ) is placed. A second order correction to the mean of in a dotted line. The lightly shaded function plots the likelihood was also Gaussian with variance matching that of the “box”. In the top figure both the prior amplitude and lengthscale are overestimated. In the bottom figure, were chosen by maximizing log with respect to their values. Notice the smaller EP approximation error.

additionally shows what the predictive mean would have been were observed under Gaussian noise with the same mean and variance as ]: it is substantially different.

Figure 5: Predicting ) with a GP with 41, while the second order correction estimates it as log , the correction to the variance is not as accurate as that on the left. The right correction is log 28, and its discrepancy with log 45 (EP+corr) is much bigger.

In the second instance, log is maximized with respect to the covariance function hyperparameters to get a kernel function that more reasonably describes the data. The correction to the mean of ) is much smaller, and furthermore, generally follows the “Gaussian noise” posterior mean. When the observed data is not typical under the prior, the correction to is substantially bigger than when the prior is representative of the data.

8.2.3 Underestimating the Truth

Under closer inspection, the variance in Figure 4 is slightly underestimated in regions where there are many close box constraints . However, under sparser constraints relative to the kernel width, EP accurately estimates the predictive mean and variance. In Figure 5 this is taken further: for 1], it is clear that q(x) becomes too narrow. The second order correction, on the other hand, provides a much closer estimate to the ground truth.

One might inquire about the behavior of the EP estimate as in Figure 5. In the next section, this will be used as a basis for illustrating a special case where log diverges.

8.3 Gaussian Process in a Box

In the following insightful example—a special case of uniform noise regression—log diverges from the ground truth with more data. Consider the ratio of functions x(s) over [0, 1], drawn from a GP prior with kernel ), such that x(s) lies within the [Figure 6 illustrates three random draws from a GP prior, two of which are not contained in the [] interval. The ratio of functions contained in the interval is equal to the

Figure 6: Samples from a GP prior with kernel two of which are not contained in the [] interval, are shown top left. As N increases in Equation (35), with diverges, while log Z converges to a constant. This is shown top right. The +’s and ’s indicate the inclusion of the fourth (+) and fourth and sixth () cumulants from the 2in Equation (28) (an arrangement by total order would include 3in show the growth for 2correction relative to the exact correction.

normalizing constant of

The fraction of samples from the GP prior that lie inside [] shouldn’t change as the GP is sampled at increasing granularity of inputs s. As Figure 6 illustrates, the MCMC estimate of log Z converges to a constant as . The EP estimate log , on the other hand, diverges to . (The cumulants that are required for the correction in Equation (28), and recipes for deriving them, are given in Appendix E.1.) Of course the correction also depends on the value a chosen. Figure 7 shows that for both correction is zero for large N.

Figure 7: The accurateness of log depends on the size of the [] box relative to with the estimation being exact as . The second order correction for Figure 6’s kernel is illustrated here over varying a’s. The +’s and ’s indicate the inclusion of the 4) cumulants in Equation (28). Of these, the top pair of lines are for N = 100, and the bottom pair for N = 50.

An intuitive explanation, due to Philipp Hennig, takes a one-dimensional model p(x) = 1). A fully-factorized approximation therefore has factors, as removing them doesn’t change p(x). However, each additional I[|x| < a] truncates the estimate, forcing EP to further reduce the variance of q(x). The EP estimate using N factors is correct (see Appendix C for a similar example and analysis), even though the original problem remains unchanged. Even though this immediate solution cannot be applied to Equation (35), the redundancy across factors could be addressed by a principled junction tree-like factorization, where tuples of “neighboring” factors can be co-treated. Although beyond the scope of this paper, Appendix A gives a guideline on how to structure such an approximation.

9. Ising Model Results

This section discusses various aspects of corrections to EP as applied to the Ising model—a Bayesian network with binary variables and pairwise potentials—in Equation (3).

We consider the set-up proposed by Wainwright and Jordan (2006) in which N = 16 nodes are either fully connected or connected to their nearest neighbors in a 4-by-4 grid. The external field (observation) strengths are drawn from a uniform distribution 25. Three types of coupling strength statistics are considered: repulsive (anti-ferromagnetic) tractive (ferromagnetic)

Previously we have shown (Opper and Winther, 2005) that EP/EC gives very competitive results compared to several standard methods. In Section 9.1 we are interested in investigating whether a further improvement is obtained with the cumulant expansion. In

Table 1: Average absolute deviation (AAD) of marginals in a Wainwright-Jordan set-up, comparing loopy belief propagation (LBP), log-determinant relaxation (LD), EC, EC with l = 4 second order correction (EC c), and an EC tree (EC t). Results in bold face highlight best results, while italics indicate where the cumulant expression is less accurate than the original approximation.

Section 9.2, we revisit the correction approach proposed in Paquet et al. (2009) and make and empirical comparison with the cumulant approach.

9.1 Cumulant Expansion

For the factorized approximation we use Equations (26) and (29) for the log Z and marginal corrections, respectively. The expression for the cumulants of the Ising model is given in Section 4.2.1. The derivation of the corresponding tree expressions may be found in Appendices B and E.4.

Table 1 gives the average absolute deviation (AAD) of marginals

while Table 2 gives the absolute deviation of log Z averaged of 100 repetitions. In two cases (Grid, with the EC tree solver. It might be some cases that a solution does not exist but we ascribe numerical instabilities in our implementation as the main cause for these problems. It is currently out of the scope of this work to come up with a better solver. We choose to report the average performance for those runs that could attain a high degree of expectation consistency: . This was 69 out of 100 in the mentioned cases and 100 of 100 in the remaining.

We observe that for the Grid simulations, the corrected marginals in factorized approximation are less accurate than the original approximation. In Figure 8 we vary the coupling strength for a specific set-up (Grid Mixed) and observe a cross-over between the

Table 2: Absolute deviation log partition function in a Wainwright-Jordan set-up, comparing EC, EC with l = 4 second order correction (EC c), EC with a full second order expansion (EC c), EC tree (EC t) and EC tree with l = 4 second order correction (EC tc). Results in bold face highlight best results. The cumulant expression is consistently more accurate than the original approximation.

Figure 8: Error on marginal (left) and log Z (right) for grid and mixed couplings as a function of coupling strength.

correction and original for the error on marginals as the coupling strength increases. We conjecture that when the error of the original solution is high then the number of terms needed in the cumulant correction increases. The estimation of the marginal seems more sensitive to this than the log Z estimate. The tree approximation is very precise for the whole coupling strength interval considered and the fourth order cumulant in the second order expansion is therefore sufficient to get often quite large improvements over the original tree approximation.

9.2 The ε-Expansion

In Paquet et al. (2009) we introduced an alternative expansion for R and applied it to Gaussian processes and mixture models. It is obtained from Equation (12) using a finite series expansion, where the normalized deviation

is treated as the small quantity instead of higher order cumulants. R has an exact representation with 2terms that we may truncate at lowest non-trivial order:

The linear terms are all equal to one because

) is a binary distribution the quadratic term becomes a weighted sum of ratios of Normal distributions:

The final expression for the lowest order approximation to R is then

From Table 2 we observe an improvement over the original factorized approximation and results similar to the cumulant correction to the factorized approximation for all settings. The -expansion may also used to calculate marginals and applied to generalized factorizations. These topics will be studied elsewhere.

10. Future Directions

Corrections to Gaussian EP approximations were examined in this paper. The Gaussian measure allowed for a convenient set of mathematical tools to be employed, mostly because it admits orthogonality of a set of polynomials, the Hermite polynomials, which allowed a clean simplification of many expressions. So far we have restricted ourselves to expansions to low orders in cumulants. Our results indicate that these first corrections to EP can already provide useful information about the quality of the EP solution. Small corrections typically show that EP is fairly accurate and the corrections improve on that. On the other hand, large corrections indicate that the EP approximation performs poorly. The low order corrections can yield a step in the right direction but in general their result may not be trusted and alternatives to the Gaussian EP approximation should be considered. It will be interesting to develop similar expansions to EP approximations with other exponential families besides the Gaussian one.

Can we expect that higher order terms in the cumulant expansion will give more reliable approximations? Before such a question could be attacked one first would need to decide in which order the terms of the expansion should be evaluated in order to obtain the most dominant contributions. For example, we might think of trying to first compute all terms in the second order expansion of the exponential in Equation (26), and then move on to higher orders. An alternative is to sort the expansion by the total sum of the orders of cumulants involved. This is in fact possible by introducing a suitable expansion parameter (which is later set equal to one) such that the formal Taylor series with respect to this parameter yields the desired expansion. However, it is not clear yet if and when such a power series expansion would actually converge. It may well be that our expansions are only of an asymptotic type (Boyd, 1999) for which the summation of only a certain number of terms might give an improvement whereas further terms would lead to worse results.

We expect that such questions could at least be answered for toy models such as the Gaussian process in a box model of Section 8.3. Our results for the latter example (together with the related uniform noise regression case) indicates that EP may not be understood as an off the shelf method for approximately calculating arbitrary high dimensional sums or integrals. One may conjecture that its quality strongly depends on