Figure 1d shows loops corresponding to the higher order Pfaffian terms on our illustrative example. The first and second subsets of triplets (d.1 and d.2) include summation over two loops whereas the remaining Pfaffian terms include uniquely one loop.

Exhaustive enumeration of all the subsets of triplets leads to a 2|T | **time algorithm, which**
is prohibitive. However, many triplet combinations may lead to forbidden configurations. Experimentally, we found that a principled way to look for higher order Pfaffian terms with large contribution is to search first for pairs of triplets, then groups of four, and so on. For large graphs, this also becomes intractable. Actually, the **problem is very similar to the**
problem of selecting loop terms rC **with largest contribution. The advantage of the Pfaffian**
representation, however, is that Z∅**is always the Pfaffian term that accounts for the largest**
number of loop terms and is the most contributing term in the series. In this work we do not derive any heuristic for searching Pfaffian terms with larger contributions. Instead, in Section 3.1 we study the full Pfaffian series and subsequently **we restrict ourselves on the**
accuracy of Z∅.

Algorithm 2 describes the full procedure to compute all terms using the representation of expression (6). The main loop of the algorithm can be interrupted at any time, thus leading to a sequence of algorithms producing corrections incrementally.

Algorithm 2 **Pfaffian series**
Arguments: **Forney graph** G
1: z **:= 0.**
2: for all (Ψ ∈ T ) do 3: Build extended graph GextΨ
applying rules of Figure 2.

4: Set Pfaffian orientation in GextΨ
according to Algorithm 1 5: Build matrices Aˆ and Bˆ.

6: Compute Pfaffian with sign correction zΨ **according to Equation (3).**
7: z := z + zΨQa∈Ψ µa;¯a.

8: **end for**
9: **RETURN** Z
BP · z

## 3. Experiments

In this Section we study numerically the proposed algorithm. To facilitate the evaluation and the comparison with other algorithms we focus on the binary Ising model, a particular case of the model (1) where factors only depend on the disagreement between two variables and take the form fa (σab, σac) = exp Ja;{ab,ac}σabσac**. We consider also nonzero local**
potentials parametrized by fa (σab) = exp Ja;{ab}σab**in all variables so that the model**
becomes planar-intractable.

We create different inference problems by choosing different interactions {Ja;{ab,ac}} and local field parameters {Ja;{ab}}**. To generate them we draw independent samples from a**
Normal distribution {Ja;{ab,ac}} ∼ N (0, β/2) and {Ja;{ab}} ∼ N (0, β**Θ), where Θ and** β determine how difficult the inference problem is. Generally, **for Θ = 0 the planar problem is** tractable. For Θ > 0, small values of β **result in weakly coupled variables (easy problems)**