Does Loopy BP give the same solutions as a Gibbs sampler? The literature in MCMC and LBP never refer to the fact that the two methods look (on expectation) exactly the same.  To illustrate, first consider a simple Ising model, that is, a graphical model where all the variables are boolean.  We can denote these variables as $\{ X_i | X \in \{0,1\}, i \in (1,...,n)\}$ and we can write the joint distribution in some factored form as: 
$$
P(X) = \frac{1}{Z}\prod_{\alpha} \Psi_\alpha(X_\alpha)
$$
Where $\alpha$ is some subset of the variables and $\Psi_\alpha$ is a potential function over these variables.  Assuming we are interested in computing the marginal distribution of a given set of variables we could run Loopy BP to compute these or we could run a Gibbs sampler to simulate from $P$ to provide samples $\{X^{(1)},...,X^{(m)}\}$ and compute the marginal as: $X_i = \frac{1}{m}\sum_{j=1}^m X^{(j)}_i$.
For now, consider running an infinite number of Gibbs samplers at the same time, this would normally be computationally intractable since we would need to maintain an infinite number of sample sets. However, for discrete distributions we can represent these samples efficiently with a single vector. In our Ising model, for example, we can represent the sets of samples with a set of 2D vectors $\{V_i | V \in \mathbb{R}^2, i \in (1,...,n)\}$.  The first component in the vector $V_i$ gives the proportion of the chains where $X_i = 1$ and the second component gives the proportion where $X_i = 0$ (of course we only actually need 1 dimension to represent this, but never mind). So, assuming we initialize the variables our Gibbs samplers uniformly then all of these vectors would equal $[0.5, 0.5]$. We can then continue with the regular Gibbs updates: loop through each of the variables (vectors) and re-sample them according to the neighbouring potentials (conditionals).  After a burn in period these vectors will represent the set of samples which are drawn from the stationary distribution $P$, we can then compute the marginal distribution for each of the variables by averaging these vectors.  
Now, other than the averaging step, this procedure is identical to loopy BP.  And, in fact, the averaging is not likely to matter since these values tend to converge in mosts cases anyway. 
If in fact, loopy BP is exactly the same as an "infinite" Gibbs sampler why is the literature on these two methods so different.  All the Loopy BP analysis seems to be concerned with how it is minimizing the Bethe free energy while the Gibbs literature is focussed on mixing rates and the ergodicity of the Markov chain. Also, would simply averaging the messages in a Loopy BP inference procedure provide a correct estimate in cases where the updates are oscillating? Could the advancements in the Gibbs sampling literature, such as block Gibbs, be used in Loopy BP schemes to speed convergence? Lastly, if Loopy BP converges to an incorrect solution, does that imply that the Gibbs chain is non-ergodic? 
 A: They are different
Consider BP update equations
$$M_{i}(x)\propto \prod_j m_{ij}(x)$$ 
$$m_{ij}(x) \propto \sum_y \psi(x,y) \prod_{k \ne i} m_{jk}(y)$$
First note is that quantities that are updated are on edges, rather than nodes. Now, note the effect of $k \ne i$ part. This makes sure that updates are "non-backtracking". So, when quantities going into node $i$ change, this flows into node $j$, but then it doesn't flow back into node $i$ on the third step. In contrast, with Gibbs sampling you can update node $i$ at time 1, that affects neighboring value of node $j$ at time 2, and this change comes back to affect node $i$ at time step 3. Since this immediate feedback loop is not present in BP, on a graph without cycles, the values stop changing after a finite number of updates.
Naive Mean field keeps states on graph nodes unlike BP, so NMF updates are more similar to Gibbs sampling. 
Another thing to note is that unless graph is a tree, belief propagation converges to the wrong answer. Consider ferromagnetic Ising model on a loop graph, 0 magnetic field. If you you make coupling strength very large, belief propagation will converge either to 0 or to 1, whereas correct answer is 0.5
As far as averaging BP values when BP is oscillating -- yes, people have done that and that helps with convergence, however that's not too useful because convergence and quality of results of BP are correlated. In other words, you could force BP to converge, but then the results will be inaccurate.
Block Gibbs is similar in spirit to Cluster Belief Propagation and Generalized Belief Propagation, which allow one to update several variables simultaneously.
For analysis of BP convergence you can consider analysis by Joris Mooij (
Sufficient conditions for convergence of Loopy Belief Propagation). Basically belief propagation gets into problems when branching factor of non-backtracking random walk on the graph gets too large. On other hand, polynomial mixing time of Gibbs sampler is more closely related to branching factor of a self-avoiding walk on the graph which is always smaller.
