Inference and Sampling of $K_{33}$-free Ising Models

TL;DR

This paper introduces polynomial-time algorithms for inference and sampling in a broad class of Ising models, including those with $K_{33}$-free topologies, extending beyond planar graphs.

Contribution

It extends tractable inference and sampling algorithms to $K_{33}$-free Ising models, generalizing planar cases to models with complex topologies.

Findings

Polynomial-time algorithms for $K_{33}$-free Ising models.

Extension of tractability from planar to $K_{33}$-free topologies.

Efficient sampling and inference in models with unbounded genus.

Abstract

We call an Ising model tractable when it is possible to compute its partition function value (statistical inference) in polynomial time. The tractability also implies an ability to sample configurations of this model in polynomial time. The notion of tractability extends the basic case of planar zero-field Ising models. Our starting point is to describe algorithms for the basic case computing partition function and sampling efficiently. To derive the algorithms, we use an equivalent linear transition to perfect matching counting and sampling on an expanded dual graph. Then, we extend our tractable inference and sampling algorithms to models, whose triconnected components are either planar or graphs of $O(1)$ size. In particular, it results in a polynomial-time inference and sampling algorithms for $K_{33}$ (minor) free topologies of zero-field Ising models - a generalization of planar graphs with a potentially unbounded genus.

Equations

$$\mathbb{P}(S=X)=\frac{1}{Z}\exp \left(\sum _{e=\{v,w\}\in E}{J}_e{x}_v{x}_w\right),$$

$$Z=\sum _{X\in \{-1,+1{\}}^N}\exp \left(\sum _{e=\{v,w\}\in E}{J}_e{x}_v{x}_w\right).$$

$$\forall e^{\ast }\in E^{\ast }:c_{e^{\ast }}=\{\begin{array}{ll}\exp (2J_{g(e^{\ast })}),& e^{\ast }\in E_I^{\ast }\\ 1,& e^{\ast }\in E_C^{\ast }\end{array}$$

$$\mathbb{P}(M(S)=E^{\prime })=\frac{1}{Z^{\ast }}\prod _{e^{\ast }\in E^{\prime }}{c}_{e^{\ast }},$$

$$Z^{\ast }=\sum _{E^{\prime }\in \text{PM}(G^{\ast })}\prod _{e^{\ast }\in E^{\prime }}{c}_{e^{\ast }}=\frac{1}{2}Z\exp \left(\sum _{e\in E}{J}_e\right)$$

$$Z=2^{-h}{Z}_1{Z}_2...Z_h.$$

$$\pi_{a}(x^{\prime},x^{\prime\prime})=\sum_{\begin{subarray}{c}x_{p}=x^{\prime},x_{t}=x^{\prime\prime}\\ \forall u\in V^{T}_{a}\setminus e_{\mathcal{V}}:\,x_{u}=\pm 1\end{subarray}}\exp\biggl{(}\sum_{\begin{subarray}{c}e=\{v,w\}\\ e\in E^{T}_{a}\setminus\{e_{\mathcal{V}}\}\end{subarray}}J_{e}x_{v}x_{w}\biggr{)},$$

$${\pi }_a(x^{\prime },x^{\prime \prime })$$

$$=Z_a{\mathbb{P}}_a(x_1=x^{\prime },x_2=x^{\prime \prime }).$$

$${\mathbb{P}}_a(x_1=x_2)$$

$$=\frac{1}{Z^{\ast }}\sum _{\begin{array}{c}{E}^{\prime }\in \text{PM}(G^{\ast }),\\ e_{\mathcal{V}}^{\ast }\in E^{\prime }\end{array}}\prod _{e^{\ast }\in E^{\prime }}{c}_{e^{\ast }}.$$

$$\{E^{\prime }\in \text{PM}(G^{\ast })|e_{\mathcal{V}}^{\ast }\in E^{\prime }\}=\{E^{\prime \prime }\cup \{e_{\mathcal{V}}^{\ast }\}|E^{\prime \prime }\in \text{PM}(G_{\mathcal{V}}^{\ast })\}.$$

$${\mathbb{P}}_a(x_1=x_2)$$

$${\pi }_a(+1,+1)$$

$${\pi }_a(+1,-1)$$

$$\displaystyle\pi_{a}(x^{\prime},x^{\prime\prime})=\sum_{\begin{subarray}{c}x_{p}=x^{\prime},x_{t}=x^{\prime\prime},\\ \forall u\in V_{a}\setminus e_{\mathcal{V}}:x_{u}=\pm 1\end{subarray}}\biggl{[}\exp\biggl{(}\sum_{\begin{subarray}{c}e=\{v,w\}\\ e\in E_{a}\setminus E_{\mathcal{V}}\end{subarray}}J_{e}x_{v}x_{w}\biggr{)}$$

$$\displaystyle\cdot\prod_{i=1}^{q}\pi_{c_{i}}(x_{p^{i}},x_{t^{i}})\biggr{]}.$$

$$\displaystyle\pi_{a}(x^{\prime},x^{\prime\prime})=\sum_{\begin{subarray}{c}x_{p}=x^{\prime},x_{t}=x^{\prime\prime},\\ \forall u\in V_{a}\setminus e_{\mathcal{V}}:x_{u}=\pm 1\end{subarray}}\exp\biggl{(}\sum_{\begin{subarray}{c}e=\{v,w\}\\ e\in E_{a}\setminus E_{\mathcal{V}}\end{subarray}}J_{e}x_{v}x_{w}$$

$$\displaystyle+\sum_{i=1}^{q}B_{i}x_{p^{i}}x_{t^{i}}\biggr{)}\cdot\exp\biggl{(}\sum_{i=1}^{q}A_{i}\biggr{)}.$$

$${\pi }_a$$

$$\displaystyle\cdot\sum_{\begin{subarray}{c}x_{p}=x^{\prime},x_{t}=x^{\prime\prime},\\ \forall u\in V_{a}\setminus e_{\mathcal{V}}:x_{u}=\pm 1\end{subarray}}\exp\biggl{(}\sum_{e=\{v,w\}\in E_{a}}J_{e}x_{v}x_{w}\biggr{)}.$$

$$\pi_{a}(x^{\prime},x^{\prime\prime})=\exp\biggl{(}\sum_{i=1}^{q}A_{i}\biggr{)}\cdot Z_{a}\mathbb{P}_{a}(x_{p}=x^{\prime},x_{t}=x^{\prime\prime})$$

$$\displaystyle Z=\sum_{X\in\{-1,+1\}^{N}}\exp\biggl{(}\sum_{e=\{v,w\}\in E}J_{e}x_{v}x_{w}\biggr{)}=$$

$$\displaystyle\sum_{\forall u\in V_{a}:x_{u}=\pm 1}\biggl{[}\exp\biggl{(}\sum_{\begin{subarray}{c}e=\{v,w\}\\ e\in E_{a}\setminus E_{\mathcal{V}}\end{subarray}}J_{e}x_{v}x_{w}\biggr{)}\cdot\prod_{i=1}^{q}\pi_{c_{i}}(x_{p^{i}},x_{t^{i}})\biggr{]}.$$

$$\{E^{\prime }\in \text{PM}(G^{\ast })|e_{\mathcal{V}}^{\ast }\notin E^{\prime }\}=\text{PM}({\overline{G}}_{\mathcal{V}}^{\ast }).$$

$$\mathbb{P}(M(S)=E^{\prime })$$

$$=\frac{2}{Z}\exp \left(\sum _{e=\{v,w\}\in E}{J}_e{x}_v^{\prime }{x}_w^{\prime }\right)$$

$$=\frac{2}{Z}\exp \left(\sum _{e^{\ast }\in E^{\prime }\cap E_I^{\ast }}2J_{g(e^{\ast })}-\sum _{e\in E}{J}_e\right)$$

$$=\frac{2}{Z}\exp \left(-\sum _{e\in E}{J}_e\right)\prod _{e^{\ast }\in E^{\prime }\cap E_I^{\ast }}{c}_{e^{\ast }}$$

$$=\frac{2}{Z}\exp \left(-\sum _{e\in E}{J}_e\right)\prod _{e^{\ast }\in E^{\prime }}{c}_{e^{\ast }}$$

$$=\frac{1}{Z^{\ast }}\prod _{e^{\ast }\in E^{\prime }}{c}_{e^{\ast }}$$

$$Z$$

$$\displaystyle=2\sum_{X\in\mathcal{C}_{+}}\biggl{[}\exp\left(\sum_{e=\{v,w\}\in E_{1}}J_{e}x_{v}x_{w}\right)\cdot\exp\left(\sum_{e=\{v,w\}\in E_{2}}J_{e}x_{v}x_{w}\right)\biggr{]}$$

$$=2\sum _{X_1\in {\mathcal{C}}_{+}^1}\exp \left(\sum _{e=\{v,w\}\in E_1}{J}_e{x}_v{x}_w\right)\cdot \sum _{X_2\in {\mathcal{C}}_{+}^2}\exp \left(\sum _{e=\{v,w\}\in E_2}{J}_e{x}_v{x}_w\right)$$

$$=\frac{1}{2}{Z}_1{Z}_2$$

$$\mathbb{P}(S=X)$$

$$=2\frac{1}{Z_1}\exp \left(\sum _{e=\{v,w\}\in E_1}{J}_e{x}_v{x}_w\right)\cdot \frac{1}{Z_2}\exp \left(\sum _{e=\{v,w\}\in E_2}{J}_e{x}_v{x}_w\right)$$

$$=2{\mathbb{P}}_1(S_1=X_1){\mathbb{P}}_2(S_2=X_2)$$

$$={\mathbb{P}}_1(S_1=X_1)\frac{{\mathbb{P}}_2(S_2=X_2)}{{\mathbb{P}}_2(s_1=x_1)}$$

Full text

1 Introduction

Computing the partition function of the Ising model is generally intractable, even an approximate solution in the special anti-ferromagnetic case of arbitrary topology would have colossal consequences in the complexity theory [\citeauthoryearJerrum and SinclairJerrum and Sinclair1993]. Therefore, a question of interest—rather than addressing the general case—is to look after tractable families of Ising models. In the following, we briefly review tractability related to planar graphs and graphs embedded in surfaces of small genus.

Related work. Onsager [\citeauthoryearOnsagerOnsager1944] gave a closed-form solution for the partition function in the case of a homogeneous interaction Ising model over an infinite two-dimensional square grid without a magnetic field. This result has opened an exciting era of phase transition discoveries, which is arguably one of the most significant contributions in theoretical and mathematical physics of the 20th century. Then, Kac and Ward [\citeauthoryearKac and WardKac and Ward1952] showed in the case of a finite square lattice that the problem of the partition function computation is reducible to a determinant. Kasteleyn [\citeauthoryearKasteleynKasteleyn1963] generalized the results to the case of an arbitrary inhomogeneous interaction Ising model over an arbitrary planar graph. Kasteleyn’s construction was based on mapping of the Ising model to a perfect matching (PM) model with specially defined weights over a modified graph. Kasteleyn’s construction was also based on the so-called Pfaffian orientation, which allows counting of PMs by finding a single Pfaffian (or determinant) of a matrix. Fisher [\citeauthoryearFisherFisher1966] simplified Kasteleyn’s construction such that the modified graph remained planar. Transition to PM is fruitful because it extends planar zero-field Ising model inference to models embedded on a torus [\citeauthoryearKasteleynKasteleyn1963] and, in fact, on any surface of small (orientable) genus $g$, but with a price of the additional, multiplicative, and exponential in genus, $4^{g}$, factor in the algorithm’s run time [\citeauthoryearGallucio and LoeblGallucio and Loebl1999].

A parallel way of reducing the planar zero-field Ising model to a PM problem consists of constructing a so-called expanded dual graph [\citeauthoryearBieche, Uhry, Maynard, and RammalBieche et al.1980, \citeauthoryearBarahonaBarahona1982, \citeauthoryearSchraudolph and KamenetskySchraudolph and Kamenetsky2009]. This approach is more natural and interpretable because there is a one-to-one correspondence between spin configurations and PMs on the expanded dual graph. An extra advantage of this approach is that the reduction allows one to develop an exact efficient sampling. Based on linear algebra and planar separator theory [\citeauthoryearLipton and TarjanLipton and Tarjan1979], Wilson introduced an algorithm [\citeauthoryearWilsonWilson1997] that allows one to sample PMs over planar graphs in $O(N^{\frac{3}{2}})$ time. The algorithms were implemented in [\citeauthoryearThomas and MiddletonThomas and Middleton2009, \citeauthoryearThomas and MiddletonThomas and Middleton2013] for the Ising model sampling, however, the implementation was limited to only the special case of a square lattice. In [\citeauthoryearThomas and MiddletonThomas and Middleton2009] a simple extension of the Wilson’s algorithm to the case of bounded genus graphs was also suggested, again with the $4^{g}$ factor in complexity. Notice that imposing zero field condition is critical, as otherwise, the Ising model over a planar graph is NP-hard [\citeauthoryearBarahonaBarahona1982]. On the other hand, even in the case of zero magnetic field Ising models over general graphs are difficult [\citeauthoryearBarahonaBarahona1982].

Contribution. In this manuscript, we discuss tractability related to the Ising model with zero magnetic fields over graphs more general than planar. Our construction is related to graphs characterized in terms of their excluded minor property. Planar graphs are characterized by excluded $K_{5}$ minor and $K_{33}$ minor (Wagner’s theorem [\citeauthoryearDiestelDiestel2006], Chapter 4.4). Therefore, instead of attempting to generalize from planar to graphs embedded into surfaces of higher genus, it is natural to consider generalizations associated with a family of graphs excluding $K_{5}$ minor or $K_{33}$ minor.

In this manuscript, we show that $K_{33}$-free zero-field Ising models are tractable in terms of inference and sampling and give a tight asymptotic bound, $O(N^{\frac{3}{2}})$, for both operations. For that purpose, we use graph decomposition into triconnected components—the result of recursive splitting by pairs of vertices, disconnecting the graph. Indeed, the $K_{33}$-free graphs are simple to work with because their triconnected components are either planar or $K_{5}$ graphs [\citeauthoryearHallHall1943]. Therefore, the essence of our construction is to decompose the inference task in Ising over a $K_{33}$-free graph into a sequential dynamic programming evaluation over planar or $K_{5}$ graphs in the spirit of [\citeauthoryearStraub, Thierauf, and WagnerStraub et al.2014]. Notice that the triconnected classification of the tractable zero-field Ising models is complementary to the aforementioned small genus classification. We illustrate the difference between the two classifications with an explicit example of a tractable problem over a graph with genus growing linearly with graph size.

Structure. The manuscript is organized as follows. Sections 2 and 3, respectively, establish notations and pose problems of inference and sampling. Section 4 presents transition from the zero-field Ising model to an equivalent tractable perfect matching (PM) model. This provides a description of a $O(N^{\frac{3}{2}})$ inference and sampling method in planar models, which is new (to the best of our knowledge), and it sets the stage for what follows. Section 5 discusses a scheme for polynomial inference and sampling in zero-field models over graphs with triconnected components that are either planar or of $O(1)$ size. Section 6 applies this scheme to $K_{33}$-free zero-field Ising models, resulting in tight asymptotic bounds, which appear to be equivalent to those in the planar case. Section 7 describes benchmarks justifying correctness and efficiency of our algorithm. Technical proofs of statements given throughout the manuscript can be found in the supplementary material.

2 Definitions and Notations

Let $V=\{v_{1},...,v_{|V|}\}$ be a finite set of vertices, a multiset $E$ consisting of $e\subseteq V$, $|e|=2$ be edges, then we call $G=(V,E)$ a graph. We call $G$ normal, if $E$ is a set (i.e., there are no multiple edges in $G$).

A tree is a connected graph without cycles. For $V^{\prime}\subseteq V$, let $G(V^{\prime})$ denote a graph $(V^{\prime},\{\{v,w\}\in E\,|\,v\in V^{\prime},w\in V^{\prime}\})$. Let $H=(V_{H},E_{H})$ be a graph. Then $H$ is a subgraph of $G$, if $V_{H}\subseteq V,E_{H}\subseteq E$. Vertex $v\in V$ is an articulation point of $G$, if $G(V\setminus\{v\})$ is disconnected. $G$ is biconnected if there are no articulation points in $G$. Biconnected component is a maximal subgraph of $G$ without an articulation point.

The graph $G$ is planar if it can be drawn on a plane without edge intersections. The corresponding drawing is referred to as planar embedding of $G$. When no ambiguity arises, we do not distinguish planar graph $G$ from its embedding.

A set $E^{\prime}\subseteq E$ is called a perfect matching (PM) of $G$, if edges of $E^{\prime}$ are disjoint and their union equals $V$. $\text{PM}(G)$ denotes the set of all PMs of $G$. $K_{p}$ denotes a complete (normal) graph on $p$ vertices, and $K_{33}$ denotes a utility graph. Triple bond is a graph of two vertices and three edges between them. Multiple bond is a graph of two vertices and at least three edges between them.

3 Problem Setup

Let $G=(V,E)$ be a normal graph, $|V|=N$. For each $v\in V$, define a random binary variable (a spin) $s_{v}\in\{-1,+1\}$, $S=(s_{v_{1}},...,s_{v_{N}})$. Subscript $i$ will be used as shorthand for $v_{i}$, for brevity, thus $S=(s_{1},...,s_{N})$. For each $e\in E$, define a pairwise interaction $J_{e}\in\mathbb{R}$. We associate assignment $X=(x_{1},...,x_{N})\in\{-1,+1\}^{N}$ to vector $S$ with probability as follows:

[TABLE]

where

[TABLE]

The probability distribution (1) defines the so-called zero-field (or pairwise) Ising model, and $Z$ is called the partition function (PF) of the zero-field Ising (ZFI) model. Notice that $\mathbb{P}(S=X)=\mathbb{P}(S=-X)$.

Given a ZFI model, our goal is to find $Z$ (inference) and draw samples from the model efficiently.

4 Reducing Planar ZFI Model to PM Model

In this section, we consider a special case of planar graph $G$ and introduce a transition from the ZFI model to the perfect matching (PM) model on a different planar graph.

We assume that the planar embedding of $G$ is given (and if not, it can be found in $O(N)$ time [\citeauthoryearBoyer and MyrvoldBoyer and Myrvold2004]). We follow [\citeauthoryearSchraudolph and KamenetskySchraudolph and Kamenetsky2009] in constructions discussed in this section.

4.1 Expanded Dual Graph

First, triangulate $G$ by adding new edges $e$ to $E$ such that $J_{e}=0$. (The triangulation does not change probabilities of the spin assignments.) Graph $G$ is generated (use the same notation as for the original graph for convenience) and is biconnected with every face, including lying on the boundary, forming a triangle. Complexity of the triangulation procedure is $O(N)$, see [\citeauthoryearSchraudolph and KamenetskySchraudolph and Kamenetsky2009] for an example.

Second, construct a new graph, $G_{F}=(V_{F},E_{F})$, where each vertex $f$ of $V_{F}$ is a face of $G$, and there is an edge $e=\{f_{1},f_{2}\}$ in $E_{F}$ if and only if $f_{1}$ and $f_{2}$ share an edge in $G$. By construction, $G_{F}$ is planar, and it is embedded in the same plane as $G$, so that each new edge $e=\{f_{1},f_{2}\}\in E_{F}$ intersects the respective old edge. Call $G_{F}$ a dual graph of $G$. Since $G$ is triangulated, each $f\in V_{F}$ has degree 3 in $G_{F}$.

Third, obtain a planar graph $G^{*}=(V^{*},E^{*})$ and its embedding from $G_{F}$ by substituting each $f\in V_{F}$ by a $K_{3}$ triangle so that each vertex of the triangle is incident to one edge, going outside the triangle (see Figure 1 for an illustration). Call $G^{*}$ an expanded dual graph of $G$.

Newly introduced triangles of $G^{*}$, substituting $G_{F}$’s vertices, are called Fisher cities [\citeauthoryearFisherFisher1966]. We refer to edges outside triangles as intercity edges and denote their set as $E^{*}_{I}$. The set $E^{*}\setminus E^{*}_{I}$ of Fisher city edges is denoted as $E^{*}_{C}$. Notice that $e^{*}\in E^{*}_{I}$ intersects exactly one $e\in E$ and vice versa, which defines a bijection between $E^{*}_{I}$ and $E$; denote it by $g:E^{*}_{I}\to E$. Observe also that $|E^{*}_{I}|=|E|\leq 3N-6$, where $N$ is the size of $G$. Moreover, $E^{*}_{I}$ is a PM of $(V^{*},E^{*})$, and thus $|V^{*}|=2|E^{*}_{I}|=O(N)$. Since $G^{*}$ is planar, one also finds that $|E^{*}|=O(N)$. Constructing $G^{*}$ takes efforts of $O(N)$ complexity.

4.2 Perfect Matching (PM) Model

For $X\in\{-1,+1\}^{N}$, let $I(X)$ be a set $\{e\in E^{*}_{I}\,|\,g(e)=\{v,w\},x_{v}=x_{w}\}$. Each Fisher city is incident to an odd number of edges in $I(X)$. Thus, $I(X)$ can be uniquely completed to a PM by edges from $E^{*}_{C}$. Denote the resulting PM by $M(X)\in\text{PM}(G^{*})$ (see Figure 1 for an illustration). Let $\mathcal{C}_{+}=\{+1\}\times\{-1,+1\}^{N-1}$.

Lemma 1**.**

$M$* is a bijection between $\mathcal{C}_{+}$ and $\text{PM}(G^{*})$.*

Define weights on $G^{*}$ according to

[TABLE]

Lemma 2**.**

For $E^{\prime}\in\text{PM}(G^{*})$ holds

[TABLE]

where

[TABLE]

is the PF of the PM distribution (PM model) defined by (2).

Second transition of (3) reduces the $Z$ computation to solve for $Z^{*}$. Furthermore, only two equiprobable spin configurations $X^{\prime}$ and $-X^{\prime}$ (one of which is in $\mathcal{C}_{+}$) correspond to $E^{\prime}$, and they can be recovered from $E^{\prime}$ in $O(N)$ steps, thus resulting in the statement that one samples from (1) if sampling from (2) is known.

The PM model can be defined for an arbitrary graph $\hat{G}=(\hat{V},\hat{E}),\hat{N}=|\hat{V}|$ with positive weights $c_{e},e\in E^{\prime}$, as a probability distribution over $\hat{M}\in\text{PM}(\hat{G})$: $\mathbb{P}(\hat{M})\propto\prod_{e\in\hat{M}}c_{e}$.

Our subsequent derivations are based on the following:

Theorem 1**.**

Given the PM model defined on planar graph $\hat{G}$ of size $\hat{N}$ with positive edge weights $\{c_{e}\}$, one can find its partition function and sample from it in $O(\hat{N}^{\frac{3}{2}})$ time.

Algorithms, constructively proving the theorem, are directly inferred from [\citeauthoryearWilsonWilson1997, \citeauthoryearThomas and MiddletonThomas and Middleton2009], with minor changes/generalizations. Hence, we outline them in the supplementary material.

Corollary**.**

Inference and sampling of the PM model on $G^{*}$ (and, hence, the ZFI model on $G$) take $O(N^{\frac{3}{2}})$ time.

5 Dynamic Programming within Triconnected Components

Starting with this section, we present new results. We describe a general algorithm that allows us to perform inference and sampling from the ZFI model in the case where the triconnected components of the underlying graph are either planar or of $O(1)$ size.

5.1 Decomposition into Biconnected Components

Consider a ZFI model (1) over a normal graph $G=(V,E)$, $|V|=N$. If $G$ is disconnected, then distribution (1) is decomposed into a product of terms associated with independent ZFI models over the connected components of $G$. Hence, we assume below, without loss of generality, that $G$ is connected.

Let $G_{1},...,G_{h}$ be biconnected components of $G$. They form a tree if an edge is drawn between $G_{i}$ and $G_{j}$ whenever $G_{i}$ and $G_{j}$ share an articulation point. A simple reduction (see supplementary material) shows that inference and sampling on $G$ are reduced to a series of inference and sampling on ZFI models induced by subgraphs $G_{1},...,G_{h}$.

Lemma 3**.**

Let $Z_{1},...,Z_{h}$ be partition functions of ZFI models induced by $G_{1},...,G_{h}$. Then,

[TABLE]

Sampling from $\mathbb{P}(S=X)$ is reduced to a series of sampling on $G_{1},...,G_{h}$ and $O(N)$ post-processing.

Observe also that all the articulation points and the biconnected components of $G$ can be found in $O(N+|E|)$ steps [\citeauthoryearHopcroft and TarjanHopcroft and Tarjan1973a]. Therefore later on, we assume without loss of generality that $G$ is biconnected.

5.2 Biconnected Graph as a Tree of Triconnected Components

In this subsection we follow [\citeauthoryearHopcroft and TarjanHopcroft and Tarjan1973b, \citeauthoryearGutwenger and MutzelGutwenger and Mutzel2001], see also [\citeauthoryearMaderMader2008] to define the tree of triconnected components. Following discussions of the previous subsection, one considers here a biconnected $G$.

Let $v,w\in G$. Divide $E$ into equivalence classes $E_{1},...,E_{k}$ so that $e_{1},e_{2}$ are in the same class if they lie on a common simple path that has $v,w$ as endpoints. $E_{1},...,E_{k}$ are referred to as separation classes. If $k\geq 2$, then $\{v,w\}$ is a separation pair of $G$, unless (a) $k=2$ and one of the classes is a single edge or (b) $k=3$ and each class is a single edge. Graph $G$ is called triconnected if it has no separation pairs.

Let $\{v,w\}$ be a separation pair in $G$ with equivalence classes $E_{1},...,E_{k}$. Let $E^{\prime}=\cup_{i=1}^{l}E_{l},E^{\prime\prime}=\cup_{i=l+1}^{k}E_{l}$ be such that $|E^{\prime}|\geq 2$, $|E^{\prime\prime}|\geq 2$. Then, graphs $G_{1}=(\cup_{e\in E^{\prime}}e,E^{\prime}\cup\{e_{\mathcal{V}}\}),G_{2}=(\cup_{e\in E^{\prime\prime}}e,E^{\prime\prime}\cup\{e_{\mathcal{V}}\})$ are called split graphs of $G$ with respect to $\{v,w\}$, and $e_{\mathcal{V}}$ is a virtual edge, which is a new edge between $v$ and $w$, identifying the split operation. Due to the addition of $e_{\mathcal{V}}$, $G_{1}$ and $G_{2}$ are not normal in general.

Split $G$ into $G_{1}$ and $G_{2}$. Continue splitting $G_{1},G_{2}$, and so on, recursively, until no further split operation is possible. The resulting graphs are split components of $G$. They can either be $K_{3}$ (triangles), triple bonds, or triconnected normal graphs.

Let $e_{\mathcal{V}}$ be a virtual edge. There are exactly two split components containing $e_{\mathcal{V}}$: $G_{1}=(V_{1},E_{1})$ and $G_{2}=(V_{2},E_{2})$. Replacing $G_{1}$ and $G_{2}$ with $G^{\prime}=(V_{1}\cup V_{2},(E_{1}\cup E_{2})\setminus\{e_{\mathcal{V}}\})$ is called merging $G_{1}$ and $G_{2}$. Do all possible mergings of the cycle graphs (starting from triangles), and then do all possible mergings of multiple bonds starting from triple bonds. Components of the resulting set are referred to as the triconnected components of $G$. We emphasize again that some graphs (i.e., cycles and bonds) in the set of triconnected components are not necessarily triconnected.

Lemma 4**.**

[\citeauthoryearHopcroft and TarjanHopcroft and Tarjan1973b]** Triconnected components are unique for $G$. Total number of edges within the triconnected components is at most $3|E|-6$.

Consider a graph $T$, where vertices (further referred to as nodes for disambiguation) are triconnected components, and there is an edge between $a$ and $b$ in $T$, when $a$ and $b$ share a (copied) virtual edge.

Lemma 5**.**

[\citeauthoryearHopcroft and TarjanHopcroft and Tarjan1973b]** $T$ is a tree.

Example.

Figure 2 illustrates triconnected decomposition of a binconnected graph and intermediate steps towards it.

All triconnected components, and thus $T$, can be found in $O(N+|E|)$ steps [\citeauthoryearHopcroft and TarjanHopcroft and Tarjan1973b, \citeauthoryearGutwenger and MutzelGutwenger and Mutzel2001, \citeauthoryearVoVo1983]. Merging of two triconnected components is equivalent to contracting an edge in $T$ (VI on Figure 2). After all possible mergings, $G$ is recovered.

5.3 Inference via Dynamic Programming

Assume that there is a (small) number $C$ bounding the size of each nonplanar triconnected component. In the following, we present a polynomial time algorithm that computes $Z$ for a given (fixed) $C$.

First, one finds triconnected components of $G$ and $T$ in $O(N+|E|)$ steps. Choose a root node $d$ in $T$. For any node $a\neq d$ in $T$, let the next node $b$ (on a unique path from $a$ to $d$) be a parent of $a$, and $a$ be a child of $b$. Nodes, which do not have any children, are called leaves. For node $a$, let a subtree $T(a)$ denote a subgraph constructed from $a$, its children, grandchildren, and so on.

Our algorithm processes each node once. The node is only processed when all its children have been already processed, so a leaf is processed first and the root is processed last. Let $a=(V_{a},E_{a}),N_{a}=|V_{a}|$ be a currently processed node. Let $G^{T}_{a}=(V^{T}_{a},E^{T}_{a})$ be a graph obtained by merging all nodes in $T(a)$. If $a$ is a root, then $G^{T}_{a}=G$. Since the root is processed last, it outputs the desired PF, Z. Figure 3 provides a visualization of a node processing routine which is to be explained.

If $a$ is not a root, let $e_{\mathcal{V}}=\{p,t\}$ be a virtual edge shared between $a$ and its parent. The only virtual edge in $G^{T}_{a}$ is $e_{\mathcal{V}}$, and $G^{T}_{a}$ without $e_{\mathcal{V}}$ is a subgraph of $G$. Hence, pairwise interactions are defined for $E^{T}_{a}\setminus\{e_{\mathcal{V}}\}$. The result of node $a$’s processing is a quantity.

[TABLE]

where $x^{\prime},x^{\prime\prime}=\pm 1$. Notice that $\pi_{a}(+1,+1)=\pi_{a}(-1,-1)$, $\pi_{a}(+1,-1)=\pi_{a}(-1,+1)$, and hence $\pi_{a}(x^{\prime},x^{\prime\prime})=\pi_{a}(x^{\prime\prime},x^{\prime})$.

Processing nodes one by one we notice that the following cases are possible:

$a$** is a leaf**. Therefore, there is nothing to merge, and $a=G^{T}_{a}=(V_{a},E_{a})$. If $a$ is nonplanar, find $\pi_{a}(\pm 1,\pm 1)$ by brute force enumeration, completed in $O(1)$ steps. If $a$ is a multiple bond, $\pi_{a}(\pm 1,\pm 1)$ is found in $O(|E_{a}|)$ steps.

Assume now that node $a$ is (or corresponds to) a planar, normal graph. Define $J_{e_{\mathcal{V}}}=0$ and consider a ZFI model with the probability $\mathbb{P}_{a}(S_{a}=X_{a})$ defined over graph $a$ with $\{J_{e}\,|\,e\in E_{a}\}$ as pairwise interactions. Let $Z_{a}$ be the PF of the ZFI model. In the remaining part of this case we will only work with this induced ZFI model, so that one can assume that nodes in $V_{a}$ are ordered, $V_{a}=\{v_{1},...,v_{N_{a}}\}$, such that $v_{1}=p,v_{2}=t$. Then, one utilizes the notations $S_{a}=(s_{1},...,s_{N_{a}})$ and $X_{a}=(x_{1},...,x_{N_{a}})\in\{-1,+1\}^{N_{a}}$ and derives

[TABLE]

Next, one triangulates $a$ by adding enough edges with zero pairwise-interactions, similar to how it is done in Subsection 4.1. Assume that $a$ is triangulated, and observe that the right-hand side of Eq. (5) is not affected. Construct $G^{*}=(V^{*},E^{*})$, which is an expanded dual graph of $a$ with $E^{*}_{I},E^{*}_{C}$, and $g$ defined as in Subsection 4.1. Then, define mapping $M:\{-1,+1\}^{N_{a}}\to\text{PM}(G^{*})$, weights $c_{e^{*}}$, and the PF $Z^{*}$ as in 4.2. Denote $e^{*}_{\mathcal{V}}=g^{-1}(e_{\mathcal{V}})$.

According to the definition of $M$,

[TABLE]

Denote $G^{*}_{\mathcal{V}}=G^{*}(V^{*}\setminus e^{*}_{\mathcal{V}})$. We continue the chain of relations/equalities (6) observing that

[TABLE]

Then one arrives at

[TABLE]

where $Z^{*}_{\mathcal{V}}$ is a PF of the PM model over $G^{*}_{\mathcal{V}}$. Compute $Z^{*}$ and $Z_{a}$ in $O(N_{a}^{\frac{3}{2}})$ steps, as described in Section 4. Since $G^{*}_{\mathcal{V}}$ is planar of size $O(N_{a})$, $Z^{*}_{\mathcal{V}}$ can also be computed in $O(N_{a}^{\frac{3}{2}})$ steps, as Theorem 1 states. The following relations finalize computation of $\pi_{a}(\pm 1,\pm 1)$ in $O(N_{a}^{\frac{3}{2}})$ steps:

[TABLE] 2. 2.

$a$** is not a leaf, not a root**. Let $c_{1},...,c_{q}$ be $a$’s children, and $e^{i}_{\mathcal{V}}=\{p^{i},t^{i}\}$ be a virtual edge shared between $c_{i}$ and $a$, $1\leq i\leq q$. At this point, we already computed all $\pi_{c_{i}}(\pm 1,\pm 1)$. Each $\{p^{i},t^{i}\}$ is a separation pair in $G^{T}_{a}$ that splits it into $G^{T}_{c_{i}}$ and the rest of $G^{T}_{a}$, containing all $G^{T}_{c_{j}}$, $j\neq i$. Denote all virtual edges in $a$ as $E_{\mathcal{V}}$, and then the following relation holds:

[TABLE]

If $a$ is (or corresponds to) a multiple bond, (7) is computed trivially in $O(|E_{a}|)$ steps. Hence, one assumes next that $a$ is a normal graph.

Each $\pi_{c_{i}}(x^{\prime},x^{\prime\prime})$ is positive, and it essentially only depends on the product $x^{\prime}x^{\prime\prime}$, that is, there exist such $A_{i},B_{i}$ that $\log\pi_{c_{i}}(x^{\prime},x^{\prime\prime})=A_{i}+B_{i}x^{\prime}x^{\prime\prime}$. Using this relation, one rewrites (7) as

[TABLE]

Denote $J_{e_{\mathcal{V}}}=0$, $J_{e^{i}_{\mathcal{V}}}=B_{i}$ for each $1\leq i\leq q$. Then rewrite (8) as

[TABLE]

We compute (9) by brute force in $O(1)$ steps, if $a$ is nonplanar. If $a$ is normal planar, we once again consider a ZFI model with the probability $\mathbb{P}_{a}(S_{a}=X_{a})$, defined over $G_{a}$, where the pairwise weights are $\{J_{e}\,|\,e\in E_{a}\}$, and $Z_{a}$ is the respective PF. Then applying machinery from Case 1, one derives

[TABLE]

in $O(N_{a}^{\frac{3}{2}})$ steps. 3. 3.

$a$** is a root**. Once again, let $c_{1},...,c_{q}$ be children of $a$, $e^{i}_{\mathcal{V}}=\{p^{i},t^{i}\}$ be a virtual edge shared between $c_{i}$ and $a$, and $1\leq i\leq q$, $E_{\mathcal{V}}$ be the set of virtual edges in $E_{a}$ (which $a$ shares only with its children). Using considerations similar to those described while deriving Eq. (7), one arrives at

[TABLE]

Finally, one computes $Z$ similarly to how the $\pi$ values were derived in Case 2. It takes $O(|E_{a}|)$ steps if $a$ is a multiple bond. Otherwise, one constructs a ZFI model and finds the PF over the respective graphs in either $O(1)$ steps, if the graph is nonplanar, or in $O(N_{a}^{\frac{3}{2}})$ steps, if $a$ is normal planar.

5.4 Sampling via Dynamic Programming

The sampling algorithm, detailed below, follows naturally from the inference routine. Compute triconnected components of $G$ in $O(N+|E|)$ steps. If all the triconnected components of $G$ are multiple bonds, $G$ should be a multiple bond itself, but $G$ is normal. Therefore, there exists a component that is not a multiple bond; choose it as a root of $T$.

Use the inference routine (described in the previous Section) to compute $Z$. Now, do a backward pass through the tree, processing the root first, and then processing the node only when its parent has already been processed (Figure 4 visualizes the sampling algorithm).

Suppose $a$ is a root and it is processed by now. Since $a$ is not a multiple bond, it results in an Ising model, $\mathbb{P}_{a}(S_{a}=X_{a})$. Draw a spin configuration $X_{a}$ from this model. It will take $O(1)$ steps if $a$ is nonplanar or $O(N_{a}^{\frac{3}{2}})$ steps if $a$ is planar.

Suppose $a$ is not a root. If $a$ is a multiple bond, spin values were already assigned to its vertices (contained within the node/graph $a$). Otherwise, there exists a ZFI model $\mathbb{P}_{a}(S_{a}=X_{a})$ already constructed at the inference stage. Following the notation of Subsection 5.3, one has to sample from $\mathbb{P}_{a}(S_{a}=X_{a}|s_{p}=x_{p},s_{t}=x_{t})$, since spins $s_{p}$ and $s_{t}$ are shared with the parent model and have already been drawn as $x_{p}$ and $x_{t}$, respectively. If $x_{p}=x_{t}$, all valid $X_{a}$ are such that $e^{*}_{\mathcal{V}}\in M(X_{a})$, and the task is reduced to sampling PMs on $G^{*}_{\mathcal{V}}$. Otherwise, all valid $X_{a}$ are such that $e^{*}_{\mathcal{V}}\notin M(X_{a})$. Denote $\overline{G}^{*}_{\mathcal{V}}=(V^{*},E^{*}\setminus\{e^{*}_{\mathcal{V}}\})$ and notice that

[TABLE]

Therefore, the task is reduced to sampling PM over $\overline{G}^{*}_{\mathcal{V}}$.

6 $K_{33}$-free Topology

6.1 ZFI Model over $K_{33}$-free Graphs

Consider the ZFI model (1) over a normal connected graph $G$. Let $H$ be some graph. Then, $H$ is a minor of $G$, if it is isomorphic to $G$’s subgraph, in which some edges are contracted. (See [\citeauthoryearDiestelDiestel2006], Chapter 1.7, for a formal definition.)

$G$ is $K_{33}$-free, if $K_{33}$ is not a minor of $G$, that is, it cannot be derived from $G$’s subgraph by contraction of some edges.

Let a biconnected $G$ be decomposed into the tree of triconnected components. Then, the following lemma holds:

Lemma 6**.**

[\citeauthoryearHallHall1943]** Graph $G$ is $K_{33}$-free if and only if its nonplanar triconnected components are exactly $K_{5}$.

Therefore, if $G$ is $K_{33}$-free, it satisfies all the conditions needed for efficient inference and sampling, described in Section 5. According to the lemma, the graph in Fig. 2 is $K_{33}$-free. The next statement expresses the main contribution of this manuscript.

Theorem 2**.**

If $G$ is $K_{33}$-free, inference or sampling of (1) takes $O(N^{\frac{3}{2}})$ steps.

We point out that the family of models for which the algorithm from Section 5 applies is broader than just $K_{33}$-free models. However, we focus on $K_{33}$-free graphs because they have a fortunate characterization in terms of a missing minor.

6.2 Discussion: Genus of $K_{33}$-free Graphs

A remarkable feature of $K_{33}$-free models is related to considerations addressing the graph’s genus. Genus of a graph is a minimal genus (number of handles) of the orientable surface that the graph can be embedded into. Kasteleyn [\citeauthoryearKasteleynKasteleyn1963] has conjected that the complexity of evaluating the PF of a ZFI model embedded in a graph of genus $g$ is exponential in $g$. The result was proven and detailed in [\citeauthoryearRegge and ZecchinaRegge and Zecchina2000, \citeauthoryearGallucio and LoeblGallucio and Loebl1999, \citeauthoryearCimasoni and ReshetikhinCimasoni and Reshetikhin2007, \citeauthoryearCimasoni and ReshetikhinCimasoni and Reshetikhin2008]. One naturally asks what are genera of graphs over which the ZFI models are tractable. The following statement relates biconnectivity and graph topology (genus):

Theorem 3**.**

[\citeauthoryearBattle, Harary, and KodamaBattle et al.1962]** A graph’s genus is a sum of its biconnected component genera.

If a graph is not biconnected, its genus can be arbitrarily large, while inference and sampling may still be tractable in relation to the decomposition technique discussed in Subsection 5.1. Therefore, it becomes principally interesting to construct tractable biconnected models with large genus.

Lemma 7**.**

A biconnected $K_{33}$-free graph of size $5n$ can be of genus as big as $n$.

From this we conclude that $K_{33}$-free graphs can’t be tackled via the bounded-genus approach of [\citeauthoryearRegge and ZecchinaRegge and Zecchina2000, \citeauthoryearGallucio and LoeblGallucio and Loebl1999, \citeauthoryearCimasoni and ReshetikhinCimasoni and Reshetikhin2007, \citeauthoryearCimasoni and ReshetikhinCimasoni and Reshetikhin2008]. This justifies the novelty of our contribution.

7 Implementation and Tests

To test the correctness of inference, we generate random $K_{33}$-free models of a given size and then compare the value of PF computed in a brute force way (tractable for sufficiently small graphs) and by our algorithm. We simulate samples of sizes from $\{10,...,15\}$ ($1000$ samples per size) and verify that respective expressions coincide.

When testing sampling implementation, we take for granted that the produced samples do not correlate given that the sampling procedure (Section 5.4) accepts the Ising model as input and uses independent random number generation inside. The construction does not have any memory, therefore, it generates statistically independent samples. To test that the empirical distribution is approaching a theoretical one (in the limit of the infinite number of samples), we draw different numbers, $m$, of samples from a model of size $N$. Then we find Kullback-Leibler divergence between the probability distribution of the model (here we use our inference algorithm to compute the normalization, $Z$) and the empirical probability, obtained from samples. Fig. 5 shows that KL-divergence converges to zero as the sample size increases. Zero KL-divergence corresponds to equal distributions.

Finally, we simulate inference and sampling for random models of different size $N$ and observe that the computational time (efforts) scales as $O(N^{\frac{3}{2}})$ (Fig. 6)222Implementation of the algorithms is available at https://github.com/ValeryTyumen/planar_ising..

8 Conclusion

In this manuscript, we compiled results that were scattered over the literature on $O(N^{\frac{3}{2}})$ sampling and inference in the Ising model over planar graphs. To the best of our knowledge, we are the first to present a complete and mathematically accurate description of the tight asymptotic bounds.

We generalized the planar results to a new class of zero-field Ising models over graphs not containing $K_{33}$ as a minor. In this case, which is strictly more general than the planar case, we have shown that the complexity bounds for sampling and inference are the same as in the planar case. Along with the formal proof, we provided evidence of our algorithm’s correctness and complexity through simulations.

Acknowledgements

This work was supported by the U.S. Department of Energy through the Los Alamos National Laboratory as part of LDRD and the DOE Grid Modernization Laboratory Consortium (GMLC). Los Alamos National Laboratory is operated by Triad National Security, LLC, for the National Nuclear Security Administration of U.S. Department of Energy (Contract No. 89233218CNA000001).

Appendix A Technical Proofs

Lemma 1 proof. Let $E^{\prime}\in\text{PM}(G^{*})$. Call $e\in E$ saturated, if it intersects an edge from $E^{\prime}\cap E^{*}_{I}$. Each Fisher city is incident to an odd number of edges in $E^{\prime}\cap E^{*}_{I}$. Thus, each face of $G$ has an even number of unsaturated edges. This property is preserved, when two faces/cycles are merged into one by evaluating respective symmetric difference. Therefore, one gets that any cycle in $G$ has an even number of unsaturated edges.

For each $i$ define $x_{i}:=-1^{r_{i}}$, where $r_{i}$ is the number of unsaturated edges on the path connecting $v_{1}$ and $v_{i}$. The definition is consistent due to aforementioned cycle property. Now for each $e=\{v,w\}\in E$, $x_{v}=x_{w}$ if and only if $e$ is saturated. To conclude, we constructed $X$ such that $E^{\prime}=M(X)$. Such $X$ is unique, because parity of unsaturated edges on a path between $v_{1}$ and $v_{i}$ uniquely determines relationship between $x_{1}$ and $x_{i}$, and $x_{1}$ is always $+1$. ∎

Lemma 2 proof. Let $X^{\prime}=(x^{\prime}_{1},...,x^{\prime}_{N})\in\mathcal{C}_{+}$, $M(X^{\prime})=E^{\prime}$. The statement is justified by the following chain of transitions:

[TABLE]

Lemma 3 proof. The Algorithm 1 reduces sampling on $G$ to a series of samplings on $G_{1},...,G_{h}$.

Given the algorithm and inference formula in Lemma 3, the statement is obvious for $h=1$. Let $h=2$. Let $v$ be an articulation point shared by $G_{1}$ and $G_{2}$. Denote $G_{1}=(V_{1},E_{1})$, $G_{2}=(V_{2},E_{2})$. Without loss of generality assume that $v$ has index $1$ in $V,V_{1}$ and $V_{2}$. Let $\mathcal{C}^{i}_{+}=\{+1\}\times\{-1,+1\}^{|V_{i}|}$. Then one derives:

[TABLE]

where $Z_{i}$ is the PF of the ZFI model induced by $G_{i}$. As far as sampling is concerned, denote by $\mathbb{P}_{i}(S_{i}=X_{i})$ a probability distribution induced by the $i$-th ZFI model. Then, since $\mathbb{P}_{2}(s_{1}=x_{1})=\frac{1}{2}$:

[TABLE]

Assume that a method for sampling $S_{i}$ from $\mathbb{P}_{i}$ is available. Then, draw $X_{1}$ by sampling $S_{1}$ from $\mathbb{P}_{1}$. To sample $S_{2}$ conditional on $s_{1}=x_{1}$ from $\mathbb{P}_{2}$, draw $X^{\prime}_{2}=(x^{\prime}_{1},...)$ from $\mathbb{P}_{2}(S_{2}=X^{\prime}_{2})$. If $x^{\prime}_{1}=x_{1}$, then $X_{2}=X^{\prime}_{2}$, otherwise $X_{2}=-X^{\prime}_{2}$. This is consistent with Algorithm 1.

For graphs of $h>2$ the statement of lemma follows naturally by induction.

Theorem 2 proof. Since $G$ is normal and minor-free, it holds that $|E|=O(N)$ [\citeauthoryearThomasonThomason2001]. Find all biconnected components and for each construct a triconnected component tree in $O(N+|E|)=O(N)$.

As described above, the time (number of steps) of inference or sampling is a sum of inference or sampling times of each triconnected component of $G$. Let the set of all $G$’s triconnected components (that is, a union over all biconnected components) to consist of $k_{1}$ planar triconnected components of size $N_{1},...,N_{k_{1}}$ with $M^{p}_{1},...,M^{p}_{k_{1}}$ edges respectively, $k_{2}$ multiple bonds of $M^{b}_{1},...,M^{b}_{k_{2}}$ edges and $k_{3}$ $K_{5}$ graphs. Then the complexity of inference or sampling is $O(\sum_{i=1}^{k_{1}}N_{i}^{\frac{3}{2}}+\sum_{i=1}^{k_{2}}M^{b}_{i}+k_{3})$.

The edges of $G$ are partitioned among biconnected components. Inside each biconnected component apply second part of Lemma 4 to obtain that $\sum_{i=1}^{k_{1}}M^{p}_{i}+\sum_{i=1}^{k_{2}}M^{b}_{i}+10k_{3}=O(|E|)=O(N)$. This gives that $\sum_{i=1}^{k_{2}}M^{b}_{i}+k_{3}=O(N)$ and $\sum_{i=1}^{k_{1}}M^{p}_{i}=O(N)$. Since triconnected components are connected graphs, we get that $N_{i}=O(M^{p}_{i})$ for all $1\leq i\leq k_{1}$ and hence $\sum_{i=1}^{k_{1}}N_{i}=O(N)$. From convexity of $f(x)=x^{\frac{3}{2}}$ it follows that $\sum_{i=1}^{k_{1}}N_{i}^{\frac{3}{2}}=O(N^{\frac{3}{2}})$ and finally that $O(\sum_{i=1}^{k_{1}}N_{i}^{\frac{3}{2}}+\sum_{i=1}^{k_{2}}M^{b}_{i}+k_{3})=O(N^{\frac{3}{2}})$. ∎

Lemma 7 proof. A simple example illustrates that genus of a biconnected $K_{33}$-free graph can grow linearly with its size. First, notice that $K_{5}$ is a nonplanar graph, but it can be embedded in toroid (Fig. 7), therefore genus of the graph is unity. Consider a cycle of length $2n$, enumerate edges in the order of cycle traversal from $1$ to $2n$. Attach $K_{5}$ graph to each odd edge of the cycle (see Fig. 7). The resulting graph $G$ is of size $5n$, it is biconnected and $K_{33}$-free (see Figure 7). Remove an arbitrary even edge from the cycle. It results in a graph whose biconnected components are $n$ $K_{5}$ graphs and $n$ edges, so its genus is $n$. Since edge removal can only decrease genus, we conclude that $G$’s genus is at least $n$.

Appendix B Counting PMs of Planar $\hat{G}$ in $O(\hat{N}^{\frac{3}{2}})$ time

This section addresses inference part of Theorem 1.

B.1 Pfaffian Orientation

Let $\hat{G}$ be an oriented graph. Its cycle of even length (built on an even number of vertices) is said to be odd-oriented, if, when all edges along the cycle are traversed in any direction, an odd number of edges are directed along the traversal. An orientation of $\hat{G}$ is called Pfaffian, if all cycles $C$, such that $\text{PM}(\hat{G}(\hat{V}-C))\neq\varnothing$, are odd-oriented.

We will need $\hat{G}$ to contain a Pfaffian orientation, moreover the construction is easy.

Theorem 4**.**

Pfaffian orientation of $\hat{G}$ can be constructed in $O(\hat{N})$.

Proof.

This theorem is proven constructively, see e.g. [\citeauthoryearWilsonWilson1997, \citeauthoryearVaziraniVazirani1989], or [\citeauthoryearSchraudolph and KamenetskySchraudolph and Kamenetsky2009], where the latter construction is based on specifics of the expanded dual graph. ∎

Construct a skew-symmetric sparse matrix $K\in\mathbb{R}^{\hat{N}\times\hat{N}}$ ($\to$ denotes orientation of edges):

[TABLE]

The next result allows to compute PF $\hat{Z}$ of PM model on $\hat{G}$ in a polynomial time.

Theorem 5**.**

$\det K>0$, $\hat{Z}=\sqrt{\det K}$.

Proof.

See, e.g., [\citeauthoryearWilsonWilson1997] or [\citeauthoryearKasteleynKasteleyn1963]. ∎

B.2 Computing $\det K$

LU-decomposition of a matrix $A=LU$, found via Gaussian elimination, where $L$ is a lower-triangular matrix with unit diagonals and $U$ is an upper-triangular matrix, would be a standard way of computing $\det A$, which is then equal to a product of the diagonal elements of $U$. However, this standard way of constructing the LU decomposition applies only if all $A$’s leading principal submatrices are nonsingular (See e.g. [\citeauthoryearHorn and JohnsonHorn and Johnson2012], Section 3.5, for detailed discussions). And already the first, $1\times 1$, leading principal submatrix of $K$ is zero/singular.

Luckily, this difficulty can be resolved through the following construction. Take $\hat{G}$’s arbitrary perfect matching $E^{\prime}\in\text{PM}(\hat{G})$. In the case of a general planar graph $E^{\prime}$ can be found via e.g. Blum’s algorithm [\citeauthoryearBlumBlum1990] in $O(\sqrt{\hat{N}}|\hat{E}|)=O(\hat{N}^{\frac{3}{2}})$ time, while for graphs $G^{*},G^{*}_{v}$ and $\overline{G}^{*}_{v}$ appearing in this paper $E^{\prime}$ can be found in $O(N)$ from a spin configuration using $M$ mapping (e.g. $E^{\prime}=E^{*}_{I}=M(\{+1,...,+1\})\in\text{PM}(G^{*})$). Modify ordering of vertices, $\hat{V}=\{v_{1},v_{2},...,v_{\hat{N}}\}$, so that $E^{\prime}=\{\{v_{1},v_{2}\},...,\{v_{\hat{N}-1},v_{\hat{N}}\}\}$. Build $K$ according to the definition (10). Obtain $\overline{K}$ from $K$ by swapping column $1$ with column $2$, $3$ with $4$ and so on. This results in $\det K=|\det\overline{K}|$, where the new $\overline{K}$ is properly conditioned.

Lemma 8**.**

$\overline{K}$’s leading principal submatrices are nonsingular.

Proof.

The proof, presented in [\citeauthoryearWilsonWilson1997] for the case of unit weights $c_{e}$, generalizes to arbitrary positive $c_{e}$. ∎

Notice, that in the general case (of a matrix represented in terms of a general graph) complexity of the LU-decomposition is cubic in the size of the matrix. Fortunately, nested dissection technique, discussed in the following subsection, allows to reduce complexity of computing $\hat{Z}$ to $O(\hat{N}^{\frac{3}{2}})$.

B.3 Nested Dissection

The partition $P_{1},P_{2},P_{3}$ of set $\hat{V}$ is a separation of $\hat{G}$, if for any $v\in P_{1},w\in P_{2}$ it holds that $\{v,w\}\notin\hat{E}$. We refer to $P_{1},P_{2}$ as the parts, and to $P_{3}$ as the separator.

Lipton and Tarjan (LT) [\citeauthoryearLipton and TarjanLipton and Tarjan1979] found an $O(\hat{N})$ algorithm, which finds a separation $P_{1},P_{2},P_{3}$ such that $\max(|P_{1}|,|P_{2}|)\leq\frac{2}{3}\hat{N}$ and $|P_{3}|\leq 2^{\frac{3}{2}}\sqrt{\hat{N}}$. The LT algorithm can be used to construct the so called nested dissection ordering of $\hat{V}$. The ordering is built recursively, by first placing vertices of $P_{1}$, then $P_{2}$ and $P_{3}$, and finally permuting indices of $P_{1}$ and $P_{2}$ recursively according to the ordering of $\hat{G}(P_{1})$ and $\hat{G}(P_{2})$ (See [\citeauthoryearLipton, Rose, and TarjanLipton et al.1979] for accurate description of details, definitions and analysis of the nested dissection ordering). As shown in [\citeauthoryearLipton, Rose, and TarjanLipton et al.1979] the complexity of finding the nested dissection ordering is $O(\hat{N}\log\hat{N})$.

Let $A$ be a $\hat{N}\times\hat{N}$ matrix with a sparsity pattern of $\hat{G}$. That is, $A_{ij}$ can be nonzero only if $i=j$ or $\{v_{i},v_{j}\}\in\hat{E}$.

Theorem 6**.**

[\citeauthoryearLipton, Rose, and TarjanLipton et al.1979]** If $\hat{V}$ is ordered according to the nested dissection and $A$’s leading principal submatrices are nonsingular, computing the LU-decomposition of $A$ becomes a problem of the $O(N^{\frac{3}{2}})$ complexity.

Notice, however, that we cannot directly apply the Theorem to $\overline{K}$, because the sparsity pattern of $K$ is asymmetric and does not correspond, in general, to any graph.

Let $G^{**}=(V^{**},E^{**})$ be a planar graph, obtained from $\hat{G}$, by contracting each edge in $E^{\prime}$, $|V^{**}|=|E^{\prime}|=\frac{1}{2}\hat{N}$. Find and fix a nested dissection ordering over $V^{**}$ (it takes $O(\hat{N}\log\hat{N})$ steps) and let the $\{v_{1},v_{2}\},\dots,\{v_{\hat{N}-1},v_{\hat{N}}\}$ enumeration of $E^{\prime}$ correspond to this ordering. Split $K$ into $2\times 2$ cells and consider the sparsity pattern of the nonzero cells. One observes that the resulting sparsity pattern coincides with the sparsity patterns of $\overline{K}$ and $G^{**}$. Since LU-decomposition can be stated in the $2\times 2$ block elimination form, its complexity is reduced down to $O(\hat{N}^{\frac{3}{2}})$.

This concludes construction of an efficient inference (counting) algorithm for planar PM model.

Appendix C Sampling PMs of Planar $\hat{G}$ in $O(\hat{N}^{\frac{3}{2}})$ time (Wilson’s Algorithm)

This section addresses sampling part of Theorem 1. In this section we assume that degrees of $\hat{G}$’s vertices are upper-bounded by $3$. This is true for $G^{*}$, $G^{*}_{v}$ and $\overline{G}^{*}_{v}$ - the only PM models appearing in the paper. Any other constant substituting $3$ wouldn’t affect the analysis of complexity. Moreover, Wilson [\citeauthoryearWilsonWilson1997] shows that any PM model on a planar graph can be reduced to bounded-degree planar model without affecting $O(\hat{N}^{\frac{3}{2}})$ complexity.

C.1 Structure of the Algorithm

Denote a sampled PM as $M$, $\mathbb{P}(M)=\hat{Z}^{-1}\prod_{e\in M}c_{e}$. Wilson’s algorithm first applies LT algorithm of [\citeauthoryearLipton and TarjanLipton and Tarjan1979] to find a separation $P_{1},P_{2},P_{3}$ of $\hat{G}$ ($\max(|P_{1}|,|P_{2}|)\leq\frac{2}{3}\hat{N}$, $|P_{3}|\leq 2^{\frac{3}{2}}\sqrt{\hat{N}}$). Then it iterates over $v\in P_{3}$ and for each $v$ it draws an edge of $M$, saturating $v$. Then it appears that, given this intermediate result, drawing remaining edges of $M$ may be split into two independent drawings over $\hat{G}(P_{1})$ and $\hat{G}(P_{2})$, respectively, and then the process is repeated recursively.

It takes $O(\hat{N}^{\frac{3}{2}})$ steps to sample edges attached to $P_{3}$ at the first step of the recursion, therefore the overall complexity of the Wilson’s algorithm is also $O(\hat{N}^{\frac{3}{2}})$.

Subsection C.2 introduces probabilities required to draw the aforementioned PM samples. Subsections C.3 and C.4 describe how to sample edges attached to the separator, while Subsection C.5 focuses on describing the recursion.

C.2 Drawing Perfect Matchings

For some $Q\in\hat{E}$ consider the probability of getting $Q$ as a subset of $M$:

[TABLE]

Let $\hat{V}_{Q}=\cup_{e\in Q}e$ and $\hat{G}_{\setminus Q}=\hat{G}(\hat{V}\setminus\hat{V}_{Q})$. Then the set $\{M^{\prime}\setminus Q\,|\,M^{\prime}\in\text{PM}(\hat{G})\}$ coincides with $\text{PM}(\hat{G}_{\setminus Q})$. This yields the following expression

[TABLE]

where

[TABLE]

is a PF of the PM model on $\hat{G}_{\setminus Q}$ induced by the edge weights $c_{e}$.

For a square matrix $A$ let $A_{c_{1},...,c_{l}}^{r_{1},...,r_{l}}$ denote the matrix obtained by deleting rows $r_{1},...,r_{l}$ and columns $c_{1},...,c_{l}$ from $A$. Let $[A]_{c_{1},...,c_{l}}^{r_{1},...,r_{l}}$ be obtained by leaving only rows $r_{1},...,r_{l}$ and columns $c_{1},...,c_{l}$ of $A$ and placing them in this order.

Now let $\hat{V}_{Q}=\{v_{i_{1}},...,v_{i_{r}}\},i_{1}<...<i_{r}$. A simple check demonstrates that deleting vertex from a graph preserves the Pfaffian orientation. By induction this holds for any number of vertices deleted. From that it follows that $K_{i_{1},...,i_{r}}^{i_{1},...,i_{r}}$ is a Kasteleyn matrix for $\hat{G}_{\setminus Q}$ and then

[TABLE]

resulting in

[TABLE]

Linear algebra transformations, described in [\citeauthoryearWilsonWilson1997], suggest that if $A$ is non-singular, then

[TABLE]

This observation allows us to express probability (11) as

[TABLE]

Now we are in the position to describe the first step of the Wilson’s recursion.

C.3 Step 1: Computing Lower-Right Submatrix of $\overline{K}^{-1}$

Find a separation $P_{1},P_{2},P_{3}$ of $\hat{G}$. The goal is to sample an edge from every $v\in P_{3}$.

Let $T$ be a set of vertices from $P_{3}$ and their neighbors, then $|T|\leq 3|P_{3}|$ because each vertex in $\hat{G}$ is of degree at most $3$. Let $T^{**}\subseteq V^{**}$ be a set of the contracted edges (recall $G^{**}$ definition from Subsection B.3), containing at least one vertex from $T$, $|T^{**}|\leq|T|$. Then $T^{**}$ is a separator of $G^{**}$ such that

[TABLE]

where one uses that, $|V^{**}|=\frac{\hat{N}}{2}$. Find a nested dissection ordering (Subsection B.3) of $V^{**}$ with $T^{**}$ as a top-level separator. This is a correct nested dissection due to Eq. (12).

Utilizing this ordering, construct $\overline{K}$. Compute $L$ and $U$ - LU-decomposition of $\overline{K}$ ($O(\hat{N}^{\frac{3}{2}})$ time). Let $t=2|T^{**}|\leq 3\cdot 2^{\frac{5}{2}}\sqrt{\hat{N}}$ and let $\mathcal{I}$ be a shorthand notation for $(\hat{N}-t+1,...,\hat{N})$. Using $L$ and $U$, find $D=[\overline{K}^{-1}]_{\mathcal{I}}^{\mathcal{I}}$, which is a lower-right $\overline{K}^{-1}$’s submatrix of size $t\times t$.

It is straightforward to observe that the $i$-th column of $D$, $d_{i}$, satisfies

[TABLE]

where $e_{i}$ is a zero vector with unity at the $i$-th position. Therefore constructing $D$ is reduced to solving $2t$ triangular systems, each of size $t\times t$, resulting in $O(t^{3})=O(\hat{N}^{\frac{3}{2}})$ required steps.

C.4 Step 2: Sampling Edges in the Separator

Now, progressing iteratively, one finds $v\in P_{3}$ which is not yet paired and draw an edge emanating from it. Suppose that the edges, $e_{1}=\{v_{j_{1}},v_{j_{2}}\},...,e_{k}=\{v_{j_{2k-1}},v_{j_{2k}}\}$, are already sampled. We assume that by this point we have also computed LU-decomposition $A_{k}=[K^{-1}]_{j_{1},...,j_{2k}}^{j_{1},...,j_{2k}}=L_{k}U_{k}$ and we will update it to $A_{k+1}$ when the new edge is drawn. Then

[TABLE]

Next we choose $j_{2k+1}$ so that $v_{j_{2k+1}}$ is not saturated yet. We iterate over $v_{j_{2k+1}}$’s neighbors considered as candidates for becoming $v_{j_{2k+2}}$. Let $v_{j}$ to become the next candidate, denote $e_{k+1}=\{v_{j_{2k+1}},v_{j}\}$. For $n\in\mathbb{N}$ let $\alpha(n)=n+1$ if $n$ is odd and $\alpha(n)=n-1$ if $n$ is even. Then the identity

[TABLE]

follows from the definition of $\overline{K}$. One deduces from Eq. (14)

[TABLE]

Constructing $T^{**}$ one has $j_{1},...,j_{2k+1},j,\alpha(j_{1}),...,\alpha(j_{2k+1}),\alpha(j)>\hat{N}-t$. It means that $A_{k+1}$ is a submatrix of $D$ with permuted rows and columns, hence $A_{k+1}$ is known.

We further observe that

[TABLE]

Therefore to update $L_{k+1}$ and $U_{k+1}$, one just solves the triangular system of equations $RU_{k}=r$ and $L_{k}Y=y$, where $R^{\top},r^{\top},Y,y$ are of size $2k\times 2$ (this is done in $O(k^{2})$ steps), and then compute $z=d-RY$ which is of the size $2\times 2$, then set, $u=\det z$.

The probability to pair $v_{j_{2k+1}}$ and $v_{j}$ is

[TABLE]

Therefore maintaining $U_{k+1}$ allows us to compute the required probability and draw a new edge from $v_{j_{2k+1}}$. By construction of $\hat{G}$, $v_{j_{2k+1}}$ has only $3$ neighbors, therefore the complexity of this step is $O(\sum_{k=1}^{|P_{3}|}k^{2})=O(\hat{N}^{\frac{3}{2}})$ because $|P_{3}|\leq 2^{\frac{3}{2}}\sqrt{\hat{N}}$.

C.5 Step 3: Recursion

Let $M_{sep}=\{e_{1},e_{2},...\}$ be a set of edges drawn on the previous step, and $\hat{V}_{sep}$ be a set of vertices saturated by $M_{sep}$, $P_{3}\subseteq\hat{V}_{sep}$. Given $M_{sep}$, the task of sampling $M\in\text{PM}(\hat{G})$ such that $M_{sep}\subseteq M$ is reduced to sampling perfect matchings $M_{1}$ and $M_{2}$ over $\hat{G}(P_{1}\setminus V_{sep})$ and $\hat{G}(P_{2}\setminus V_{sep})$, respectively. Then $M=M_{1}\cup M_{2}\cup M_{sep}$ becomes the result of the perfect matching drawn from (2).

Even though only the first step of the Wilson’s recursion was discussed so far, any further step in the recursion is done in exactly the same way with the only exception that vertex degrees may become less than $3$, while in $\hat{G}$ they are exactly $3$. Obviously, this does not change the iterative procedure and it also does not affect the complexity analysis.

Appendix D Random Graph Generation

As our derivations cover the most general case of planar and $K_{33}$-free graphs, we want to test them on graphs which are as general as possible. Based on Lemma 6 (notice, that it provides necessary and sufficient conditions for a graph to be $K_{33}$-free) we implement a randomized construction of $K_{33}$-free graphs, which is assumed to cover most general $K_{33}$-free topologies.

Namely, one generates a set of $K_{5}$’s and random planar graphs, attaching them by edges to a tree-like structure. For simplicity, we slightly relax the condition that random planar components should be triconnected (because it is not clear how to generate such graphs efficiently) and simply require the components to be biconnected. This can be interpreted as constructing $T$, where some neighbor planar nodes are merged (merging planar graphs results in another planar graph). We refer to such non-unique decomposition $T^{\prime}$ as partially merged. Inference and sampling algorithm suggested in Section 5 is applied with no changes to the partially merged decomposition. Our generation process consists of the following two steps.

Planar graph generation. This step accepts $N\geq 3$ as an input and generates a normal biconnected planar graph of size $N$ along with its embedding on a plane. The details of the construction are as follows.

First, a random embedded tree is drawn iteratively. We start with a single vertex, on each iteration choose a random vertex of an already “grown” tree, and add a new vertex connected only to the chosen vertex. Items I-V in Fig. 8 illustrate this step.

Then we triangulate this tree by adding edges until the graph becomes biconnected and all faces are triangles, as in the Subsection 4.1 (VI in Figure 8). Next, to get a normal graph, we remove multiple edges possibly produced by triangulation (VII in Fig. 8). At this point the generation process is complete. 2. 2.

$K_{33}$-free graph generation. Here we take $N\geq 5$ as the input and generate a normal biconnected $K_{33}$-free graph $G$ in a form of its partially merged decomposition $T^{\prime}$. Namely, we generate a tree $T^{\prime}$ of graphs where each node is either a normal biconnected planar graph or $K_{5}$, and every two adjacent graphs share a virtual edge.

The construction is greedy and is essentially a tree generation process from Step 1. We start with $K_{5}$ root and then iteratively create and attach new nodes. Let $N^{\prime}<N$ be a size of the already generated graph, $N^{\prime}=5$ at first. Notice, that when a node of size $n$ is generated, it contributes $n-2$ new vertices to $G$.

An elementary step of iteration here is as follows. If $N-N^{\prime}\geq 3$, a coin is flipped and the type of new node is chosen - $K_{5}$ or planar. If $N-N^{\prime}<3$, $K_{5}$ cannot be added, so a planar type is chosen. If a planar node is added, its size is drawn uniformly in the range between $3$ and $N-N^{\prime}+2$ and then the graph itself is drawn as described in Step 1. Then we attach a new node to a randomly chosen free edge of a randomly chosen node of $T^{\prime}$. We repeat this process until $G$ is of the desired size $N$. Fig. 9 illustrates the algorithm.

To obtain an Ising model from $G$, we sample pairwise interactions for each edge of $G$ independently from $\mathcal{N}(0,0.1^{2})$.

Notice that the tractable Ising model generation procedure is designed in this section solely for the convenience of testing and it is not claimed to be sampling models of any particular practical interest (e.g. in statistical physics or computer science).

Appendix E Future Work

We conclude by discussing some future research directions:

• •

The class of models considered in the manuscript can be extended even further towards $K_{33}$-free generalizations of (a) the so-called outerplanar graphs, which can then be used for approximate inference and efficient learning in the spirit of [\citeauthoryearGloberson and JaakkolaGloberson and Jaakkola2007] and [\citeauthoryearJohnson, Oyen, Chertkov, and NetrapalliJohnson et al.2016] respectively; and (b) graphs embedded in the surfaces of $O(1)$ genus [\citeauthoryearRegge and ZecchinaRegge and Zecchina2000, \citeauthoryearGallucio and LoeblGallucio and Loebl1999, \citeauthoryearCimasoni and ReshetikhinCimasoni and Reshetikhin2007, \citeauthoryearCimasoni and ReshetikhinCimasoni and Reshetikhin2008].

• •

This manuscript was motivated by a larger task of using efficient inference and learning over the most general $K_{33}$-graphs for constructing more general (and thus, hopefully, more powerful) alternatives to traditional Neural Networks for efficient learning.