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.

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.

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.

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.

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.

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.

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).

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.

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.

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 the fact that such sums or integrals may or may not have an interpretation in terms of a proper statistical inference model which contain data that are highly probable with respect to the model. It would be interesting to see if one can develop a theory for the average case performance of EP under such statistical assumptions of the data.

As p(x) is a latent Gaussian model, the g-terms in Equation (5) are chosen in this paper to give a Gaussian approximation

The sufficient statistics ) and natural parameters of the Gaussian are defined as

where . There exists a bijection between the canonical parameters and natural parameters, such that the mean and covariance can be determined with

In Equation (1) we can define ), such that it is essentially a rescaling of factor . In the Ising model in Equation (3), this means that . In the Gaussian process classification model in Equation (2), this implies that

It remains to define a suitable factorization for the term-product ). This factoriza- tion can be fully factorized, factorized over disjoint sets of variables, factorized as a tree, or follow more arbitrary factorizations (see the simple example in Appendix C). A few such factorizations are given below in increasing orders of complexity. In each case we do not Furthermore, even though the term factorization may be chosen to fully factorize, q(x) may be fully connected through the inclusion of

A common factorization of ) is to set ). The natural parameters of are chosen to be ), corresponding to (). For clarity the other and Λ parameters in are not shown, as they are clamped at zero. This gives an approximation q(x) that is defined by

As a second step the N variables can be subdivided into factorization over terms couples pairs of variables through

In this case each factor will have a contribution to the overall ap- proximation, and, as is a function of two variables, it is parameterized by the “correlated Gaussian form” ). By symmetry Λ. The result- ing q(x) is defined in terms of these disjoint sets with

A tree structure factorization can be defined by extending the above “disjoint pairs” case to allow for overlaps between terms. Let G define a spanning tree structure over all x, and let define the edges in the tree. Let be the number of edges emanating from node in the graph. Through a clever regrouping of terms into a “junction tree” form with

the term-approximation will be tree-structured. In this example the powers are 1 for edge factors ) for node factors ) be parameterized by , as was done in the two examples above. Using

the resulting q(x) has parameter vector

It is useful to note that the form of the tree-structured approximation given here is that used by Opper and Winther (2005); it approximates the “junction tree” form using a Power EP factorization (Minka, 2004). The factorization and stationary condition is different from that of Tree EP (Minka and Qi, 2004).

The EP moment matching conditions from Equation (7) are uniquely met at the stationary point of log in Equation (8), and are shown here. Consider the logarithm of the normalizer,

Using the sufficient statistics and natural parameters defined above, the two normalizers that constitute Equation (36) are

Using these definitions, the derivatives of the terms in Equation (36) with respect to some EP factor c’s parameters

When , the following therefore holds:

Let D be a square matrix where the values in column ; all the rows in D are equal and it is singular. Furthermore, let . By stacking all the a column vector , the above set of equalities lead to a system of equations

(The Kronecker product is only required as the sufficient statistics’ differences dimensionality “dim”, usually larger than one.) As is nonsingular, it is solved by , and hence

The choice of parameterization of might give an overcomplete representation, and the exact moment-matching conditions might have more than one unique solution. However, this does not invalidate that at the stationary point of Equation (36), all moment-matching conditions must hold.

Let the factorization of the term-product ) take the form of a tree G with edges , as is described in Appendix A.1.3. The number connections to a node or

vertex n shall be denoted by . From Equation (32) the second order expansion is

where the inner expectations are over , while the outer expectations are over The edge-edge, edge-node, and node-node expectations that are needed in Equation (37) are given in the following three sections.

The edge-edge expectation provides a beautiful illustration of the combinatorics that may be involved in Wick’s theorem. For , the following expectation needs to be evaluated:

The vectors that are summed over to get . From the independence of

and therefore

Wick’s theorem is again instrumental in computing , as all possible pairings of the random variables ) need to be included. As expansion of Equation (39) occur when all the variables in are paired. This immediately means that , as there will be some remaining variables in ) that can’t be paired and have to be self-paired with zero expectation.

Given l = s, evaluate the expectation in Equation (39). We introduce the “pairing count” vector with elements and constraint number of pairings of count the number of pairings of there are terms, the sum of its outgoing pairings should equal

A furthermore requirement is that

where be as in the Wick expansion below. Define B to be the set of all such ’s, and let ) count the number of permuted configurations for a given pairing . From Wick’s theorem the expected value is equal to the sum over all possible pairings

A simple scheme to enumerate all

so that The remaining components of are uniquely determined from

How many permuted pairings ) are there?

1. There areways of choosing ’s, and then ways of choosing pair with.

2. This leaves a remaining (’s, that need to be paired with (There are such pairings.

3. There are also ’s, that need to be paired with variables. There areways of picking a , and a further ()! ways of arranging the remaining

4. Finally, the () remaining need to be coupled with the remaining ’s, and there are (())! such arrangements.

Multiplying the possible pairings from the four steps above gives

which adds up to the total number of possible pairings !. A further useful simplification is , and is used below.

The absence of any self-interacting loops from Wick’s theorem lets the drop away in Equation (38), as all terms are zero except for when l = s. Substituting into Equation (38) gives the final result,

The derivation for the edge-node expectations is similar to that of the edge-edge case,

where the expectations in the last line are again over is evaluated with Wick’s theorem, there are zero relation of ensures that the only non-zero terms in the Wick sum are those where all the ’s are paired with ’s; in other words, when l = s. There are l! possible pairings, which cancels l! in the denominator.

The above edge-node expectation is for any edge and node in the tree, but notice that it simplifies greatly when the edge is a connection to node is the edge variable corresponding to . In this case the covariance with respect to the opposite pair is zero, with ’s will have a non-zero contribution to the sum, namely when

The node-node expectation is given in Equation (27), and is also used for

The following example illustrates a tractable one-dimensional model with two factors. It is shown analytically that the correction to log must be zero, and that the result is reflected in the higher-order terms in Equation (32), which are also zero.

Consider the factorization of a probit term with a Gaussian prior into

where Φ(x) is the cumulative Gaussian density function, and be computed exactly, but for the sake of example p(x) will be approximated with

Choose approximation has parameter vector . The EP fixed point is defined by . (For example, subtracting at the fixed point will leave equal to a scaled version of the prior ). The factor ) is hence incorporated into the prior, giving . By a symmetric argument, .) Although it is trivial to show that will be equal to the true partition function Z, we shall prove it by showing that the correction term is log R = 0.

In this section a transformation of variables from will be used to make the derivation slightly simpler, and therefore

Below we analytically show that the correction log R is zero, and hence that

where ) is a shorthand for

Because , the cumulants will be the same for all . Furthermore, are both distributed according to the same density. Now define, using e1 +

In the second line above a transformation of variables was made in the integral, with u = , such that . The Jacobian 1ensures proper normalization so that the average is over 1). In the last line degree l,

which can be obtained for any real y and integer l = 0, 1, 2, . . . from the average

The remarkable property , ensures that (41). Furthermore, ) follows from the equivalence in cumulants roots in Equation (40) disappear to give , proving that R = 1 in Equation (40).

The second order expansion in Equation (32) in Section 7 evaluates to zero, as the matching cumulants and equal distributions of ensure that

Corrections to the marginal distributions follow from a similar derivation to that of the normalizing constant. As a simplification, let the Gaussian approximation be centred with ), and assume that q(x) is arises from the fully factorized approximation in Section 5. In this appendix corrections will be computed for the mean , and variance

A further simplification that will be employed in the following section is a change of variables

which is zero-mean complex Gaussian random variable with a relation . Following Equation (24), the correction reads

The lowest order correction to the EP marginal’s mean follows from the result in Equation (13):

pears as

second order term also disappears asThese equivalences can be seen by taking ) (and also its derivative) as a expansion over powers of should be zero. Hence

The derivative of the characteristic function, as required in Equation (42), is

The expectations for in Equation (42) evaluate to

with the second term disappearing as is always self-paired in Wick’s theorem. Finally, by substituting Equation (43) into (42), the correction to the mean is

The correction to the second moments follow the same recipe as that of the marginal mean in Appendix D.1. We proceed by first treating

Reapplying the recipe gives the correction to the covariance:

Much of this paper hinges on cumulants beyond the second order. These are frequently more cumbersome to obtain than the initial moments that are required by EP. This appendix provides details of the cumulants used in this paper.

The cumulants of a distribution ) can be obtained from its moments through

they are derived for doubly-truncated Gaussian distributions in Appendices E.1 and E.2. One might also directly take derivatives of the cumulant generating function, and the cumulants of a Probit-times-Gaussian distribution, common to GP classification models, are derived this way in Appendix E.3.

The tree-structured approximation in Sections 7 and 9.1, and Appendices A.1.3 and B, require cumulants over two variables. They are presented in Appendix E.4 for the Ising model.

Consider the centered distribution ). The odd moments of this tilted distributions are, by symmetry,

with the Probit function being Φ(. Subscripts n are dropped where they are clearly implied by their context. To get the even moments, consider

Using the partition function, we get

The same calculation from Appendix E.1 can be repeated to get the moments of the non-centered truncated Gaussian ). The subscripts n are dropped where evident. The partition function is

Figure 9: The moments of ), as a function of Gaussian variance , the moments converge to that of a uniform distribution.

By again taking increasing derivatives of ) with respect to , the moments solved for are

Finally,

As Figure 9 illustrates, these moments will converge to that of a uniform distribution as the Gaussian’s variance grows large.

Figure 10: The third and fourth cumulants of the density in Appendix E.3. The step function Θ(x), with m = v = 0, is taken as an example here. The third cumulant is always positive, while the fourth cumulant is positive only when

EP approximations to Probit regression models, and Gaussian process classification models in general (see Section 8.1), depend on the moments of We introduce 0 so that the likelihood can become a step function at v = 0, for example. We shall obtain the cumulants by taking derivatives of the characteristic function. The characteristic function of ), as described by Equation (15), is

with

The cumulants are determined from the derivatives of log ) at zero; a lengthy calculation shows that they are

where

We need some third and fourth order two-variable cumulants and thus generalize the results of Section 4.2 to the bivariate case. To do this we can exploit the cumulant generating property of log denote the joint order cumulant of variable one and

two, respectively. We can generate this cumulant from derivatives of log

We can also express this as a recursion in terms of cumulants:

By explicit calculation for a bivariate binary distribution we get the first two orders’ cumulants: is equal to the covariance between the two variables (to be matched with q(x)). The fact that we can write in terms of the first order cumulant shows that we can express all order cumulants in terms of the first and second order cumulant for example:

Using the same recursion it is easy to show: 24

D. Barber. Bayesian Reasoning and Machine Learning. Cambridge University Press, 2012.

C. M. Bishop. Pattern Recognition and Machine Learning. Springer, 2006.

S. Blinnikov and R. Moessner. Expansions for nearly Gaussian distributions. Astronomy and Astrophysics Supplement Series, 130:193–205, 1998.

J. P. Boyd. The devil’s invention: Asymptotic, superasymptotic and hyperasymptotic series. Acta Applicandae Mathematicae, 56:1–98, 1999.

M. Chertkov and V. Y. Chernyak. Loop series for discrete statistical models on graphs. Journal of Statistical Mechanics: Theory and Experiment, 2006:P06009, 2006.

B. Cseke and T. Heskes. Approximate marginals in latent Gaussian models. Journal of Machine Learning Research, 12:417–457, 2011.

M. Kuss and C. E. Rasmussen. Assessing approximate inference for binary Gaussian process classification. Journal of Machine Learning Research, 6:1679–1704, 2005.

M. M´ezard, G. Parisi, and M. A. Virasoro. Spin Glass Theory and Beyond, volume 9 of Lecture Notes in Physics. World Scientific, 1987.

T. P. Minka. Expectation propagation for approximate Bayesian inference. In UAI 2001, pages 362–369, 2001a.

T. P. Minka. A family of algorithms for approximate Bayesian inference. PhD thesis, MIT Media Lab, 2001b.

T. P. Minka. Power EP. Technical Report MSR-TR-2004-149, Microsoft Research Ltd, 2004.

T. P. Minka and Y. Qi. Tree-structured approximations by expectation propagation. In Advances in Neural Information Processing Systems 16. 2004.

K. P. Murphy. Machine Learning: A Probabilistic Perspective. The MIT Press, 2012.

M. Opper, U. Paquet, and O. Winther. Improving on expectation propagation. In Advances in Neural Information Processing Systems 21, pages 1241–1248. 2009.

M. Opper and O. Winther. Gaussian processes for classification: Mean field algorithms. Neural Computation, 12:2655–2684, 2000.

M. Opper and O. Winther. Expectation consistent approximate inference. Journal of Machine Learning Research, 6:2177–2204, 2005.

U. Paquet, M. Opper, and O. Winther. Perturbation corrections in approximate inference: Mixture modelling applications. Journal of Machine Learning Research, 10:935–976, 2009.

C. E. Rasmussen and C. K. I. Williams. Gaussian Processes for Machine Learning. The MIT Press, 2005.

H. Rue, S. Martino, and N. Chopin. Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 71(2):319–392, 2009.

M. W. Seeger and H. Nickisch. Fast convergent algorithms for expectation propagation approximate Bayesian inference. Arxiv preprint arXiv:1012.3584, 2010.

D. Sherrington and S. Kirckpatrick. Solvable model of a spin-glass. Phys. Rev. Lett., 35 (26):1792–1796, December 1975.

E. Sudderth, M. Wainwright, and A. Willsky. Loop series and Bethe variational bounds in attractive graphical models. In Advances in Neural Information Processing Systems 20, pages 1425–1432. 2008.

D. J. Thouless, P. W. Anderson, and R. G. Palmer. Solution of a ‘solvable model of a spin glass’. Phil. Mag., 35:593, 1977.

M. A. J. van Gerven, B. Cseke, F. P. de Lange, and T. Heskes. Efficient Bayesian multivariate fMRI analysis using a sparsifying spatio-temporal prior. NeuroImage, 50(1):150–161, 2010.

M. J. Wainwright and M. I. Jordan. Log-determinant relaxation for approximate inference in discrete Markov random fields. IEEE Transactions on Signal Processing, 54(6):2099– 2109, 2006.

M. Welling, A. Gelfand, and A. Ihler. A cluster-cumulant expansion at the fixed points of belief propagation. In Uncertainty in Artificial Intelligence (UAI). 2012.