Variational approximations for stochastic dynamics on graphs

TL;DR

This paper explores mean-field-like variational approximations for stochastic graph dynamics, providing a unified framework that improves accuracy and is validated against Monte Carlo simulations in epidemic models.

Contribution

It introduces a systematic cluster-variational approach to improve mean-field approximations for stochastic processes on graphs.

Findings

The approach unifies various existing approximation schemes.

Systematic improvements lead to better accuracy over traditional methods.

Validated with Monte Carlo simulations on epidemic models.

Abstract

We investigate different mean-field-like approximations for stochastic dynamics on graphs, within the framework of a cluster-variational approach. In analogy with its equilibrium counterpart, this approach allows one to give a unified view of various (previously known) approximation schemes, and suggests quite a systematic way to improve the level of accuracy. We compare the different approximations with Monte Carlo simulations on a reversible (susceptible-infected-susceptible) discrete-time epidemic-spreading model on random graphs.

Equations

$$\wp (x^{(0)},\dots ,x^{(\tau )})\triangleq \mathbb{P}\{{\xi }^{(0)}=x^{(0)},\dots ,{\xi }^{(\tau )}=x^{(\tau )}\}.$$

$$\wp (x^{(0)},\dots ,x^{(\tau )})=p^{(0)}(x^{(0)})\prod _{t=0}^{\tau -1}{w}^{(t)}(x^{(t+1)}|x^{(t)}),$$

$$w^{(t)}(y|x)=\prod _i{w_i}^{(t)}(y_i|x_{i,\partial i}),$$

$$p^{(0)}(x)=\prod _i{p_i}^{(0)}(x_i).$$

$$\mathcal{F}[\wp ]=\sum _{x^{(0)},\dots ,x^{(\tau )}}\wp (x^{(0)},\dots ,x^{(\tau )})\ln \frac{\wp (x^{(0)},\dots ,x^{(\tau )})}{p^{(0)}(x^{(0)}){\prod }_{t=0}^{\tau -1}{w}^{(t)}(x^{(t+1)}|x^{(t)})},$$

$$\mathcal{S}[\wp ]=-\sum _{x^{(0)},\dots ,x^{(\tau )}}\wp (x^{(0)},\dots ,x^{(\tau )})\ln \wp (x^{(0)},\dots ,x^{(\tau )}),$$

$$\mathcal{F}[\wp ]=\sum _{x^{(0)},\dots ,x^{(\tau )}}\wp (x^{(0)},\dots ,x^{(\tau )})\epsilon (x^{(0)},\dots ,x^{(\tau )})-\mathcal{S}[\wp ],$$

$$\epsilon (x^{(0)},\dots ,x^{(\tau )})\triangleq -\ln p^{(0)}(x^{(0)})-\sum _{t=0}^{\tau -1}\ln w^{(t)}(x^{(t+1)}|x^{(t)}).$$

$$\epsilon (x^{(0)},\dots ,x^{(\tau )})=-\sum _i\ln {p_i}^{(0)}({x_i}^{(0)})-\sum _{t=0}^{\tau -1}\sum _i\ln {w_i}^{(t)}({x_i}^{(t+1)}|{x_{i,\partial i}}^{(t)}).$$

$${\mathcal{S}}_{\mathrm{CVM}}[\{{\wp }_{\alpha }{\}}_{\alpha \in \mathsf{R}}]=\sum _{\alpha \in \mathsf{R}}{a}_{\alpha }\mathcal{S}[{\wp }_{\alpha }],$$

$$\sum _{\beta \in \mathsf{R}|\beta \supseteq \alpha }{a}_{\beta }=1\forall \alpha \in \mathsf{R},$$

$${\mathcal{S}}_{\mathrm{P}}[P,Z,A]=\sum _{t=0}^{\tau -1}\left\{\sum _i\mathcal{S}[{P_i}^{(t)}]-\sum _{ij}\mathcal{S}[{Z_{ij}}^{(t)}]\right\}-\sum _{t=1}^{\tau -1}\sum _i\mathcal{S}[{A_i}^{(t)}],$$

$${\mathcal{S}}_{\mathrm{PQ}}[P,Q,T,Z,A]$$

$${\mathcal{S}}_{\mathrm{PQR}}[P,Q,R,S,T,U,V,Z]=\sum _{t=0}^{\tau -1}\{\sum _i\left[\mathcal{S}[{P_i}^{(t)}]+\mathcal{S}[{R_i}^{(t)}]+(d_i-1)\mathcal{S}[{V_i}^{(t)}]\right]+$$

$$+\sum _{ij}\mathcal{S}[{Q_{ij}}^{(t)}]-\sum _{i,j\in \partial i}\left[\mathcal{S}[{T_{i,ij}}^{(t)}]+\mathcal{S}[{U_{ij,i}}^{(t)}]\right]\}-\sum _{t=1}^{\tau -1}\left\{\sum _i\mathcal{S}[{S_i}^{(t)}]-\sum _{ij}\mathcal{S}[{Z_{ij}}^{(t)}]\right\},$$

$${\mathcal{S}}_{\mathrm{M}}[M,Q,S,Z]=\sum _{t=0}^{\tau -1}\left\{\sum _i\mathcal{S}[{M_i}^{(t)}]-\sum _{ij}\mathcal{S}[{Q_{ij}}^{(t)}]\right\}-\sum _{t=1}^{\tau -1}\left\{\sum _i\mathcal{S}[{S_i}^{(t)}]-\sum _{ij}\mathcal{S}[{Z_{ij}}^{(t)}]\right\}.$$

$$w_i(0|x_{i,\partial i})=x_i{\mu }_i+(1-x_i)\prod _{j\in \partial i}\left(1-x_j{\beta }_{j\to i}\right),$$

$$w_i(1|x_{i,\partial i})=1-w_i(0|x_{i,\partial i}).$$

$${p_i}^{(0)}(x_i)=x_i,$$

$${p_i}^{(0)}(x_i)=1-x_i.$$

$${p_i}^{(0)}(x_i)=0.01x_i+0.99(1-x_i).$$

$${{\rho }_i}^{(t)}\triangleq \mathbb{E}\{{{\xi }_i}^{(t)}\}=\mathbb{P}\{{{\xi }_i}^{(t)}=1\}={A_i}^{(t)}(1),$$

$${A_i}^{(t)}(x_i)=\sum _{x_{\partial i}}\sum _{y_i}{P_i}^{(t)}(y_i,x_{i,\partial i}).$$

$$\sum _{y_i}{P_i}^{(t)}(y_i,x_{i,\partial i})$$

$$\sum _{x_{\partial i\setminus j}}\sum _{y_i}{P_i}^{(t)}(y_i,x_{i,\partial i})$$

$$\sum _{x_{\partial i\setminus j}}{P_i}^{(t)}(y_i,x_{i,\partial i})$$

$$\sum _{x_{\partial i}}{P_i}^{(t)}(y_i,x_{i,\partial i})$$

$$\sum _{y_j}{Q_{ij}}^{(t)}(y_{i,j},x_{i,j})$$

$$\sum _{y_{i,j}}{Q_{ij}}^{(t)}(y_{i,j},x_{i,j})$$

$$\sum _{x_j}\sum _{y_j}{Q_{ij}}^{(t)}(y_{i,j},x_{i,j})$$

$$\sum _{x_j}{Q_{ij}}^{(t)}(y_{i,j},x_{i,j})$$

$$\sum _{x_{i,j}}{Q_{ij}}^{(t)}(y_{i,j},x_{i,j})$$

$$\sum _{y_{\partial i\setminus j}}{R_i}^{(t)}(y_{i,\partial i},x_i)$$

$$\sum _{y_{\partial i}}{R_i}^{(t)}(y_{i,\partial i},x_i)$$

$$\sum _{x_i}\sum _{y_{\partial i\setminus j}}{R_i}^{(t)}(y_{i,\partial i},x_i)$$

Full text

1 Introduction

In statistical mechanics, the cluster variational method (CVM), originally proposed by Kikuchi in 1951 [1], is a well-established technique, which provides a systematic way to improve the accuracy of the simplest mean-field approximation. The classical scope of application is that of equilibrium thermodynamics of lattice models with discrete degrees of freedom, such as Ising or Potts [2]. The basic idea is to cut off correlations among site variables beyond the range of certain maximal clusters (for instance, square “plaquettes” on a square lattice), so that the level of accuracy can be somehow “tailored” by means of a suitable choice of such clusters. The output of the calculation is the set of marginals of the Boltzmann distribution over the maximal clusters. In particular, choosing single sites or nearest-neighbor pairs as maximal clusters, the CVM recovers, respectively, the ordinary mean-field (Bragg-Williams) approximation and the Bethe-Peierls approximation [2]. More generally, one can show that a number of generalized mean-field approximations, originally developed by heuristic arguments, can be incorporated into the CVM framework [1, 3]. Let us also mention the fact that the method has gained a renewed interest in recent years, because of the formal analogy between equilibrium statistical mechanics and a number of inference problems of high technological relevance (error-correction decoding, pattern recognition, etc.), which ultimately require computing marginals of the Boltzmann distribution for (heterogeneous) statistical-mechanical models defined on (very large) graphs [4].

In the complementary field of non-equilibrium statistical mechanics, a rigorous foundation from first principles is still lacking, so that a considerable part of research work is devoted to model systems in which the (stochastic) dynamics is assumed to be known, and properly defined to describe some kind of phenomena. This class of models includes, for instance, epidemic-spreading models [5, 6, 7], kinetic Ising or Ising-like models [8] (including models of neural-network dynamics) [9, 10], models of opinion formation and social dynamics (Voter model and its numerous variants) [11, 12, 13], transport models (exclusion processes) [14] and so on. Increasing research efforts have also been oriented to investigating the aforementioned dynamics on complex networks, with special attention to the role played by the latter [6, 15]. Most of the stochastic processes involved in these models are of the Markovian (i.e., memoryless) type, fully characterized by a matrix of transition probabilities (stochastic matrix). The problem can be generally rephrased in a variational form through a thermodynamic formalism, fully analogous to that of equilibrium statistical mechanics, where the transition probabilities play the role of interactions, and, within this formalism, one can develop mean-field-like approximation strategies along the lines of the cluster variational method. Such a “dynamical analog” of the CVM was indeed put forward in the 1960s by Kikuchi himself [16], who named it path probability method (PPM). The latter method has not become so popular as its equilibrium counterpart, so that in fact most mean-field-like approaches to stochastic dynamics are still recently conceived on the basis of traditional heuristic reasoning, without resorting to the cluster-variational machinery [5, 6, 17].

At odds with the cited works, a path-probability approach has been recently followed by one of the authors [18], to develop a new type of approximation which, due to the peculiar “shape” of the maximal clusters, has been called the diamond approximation. The latter has been applied to kinetic Ising models (with Glauber dynamics, and either symmetric or asymmetric couplings), but it can be straightforwardly used for other Markov-chain models as well, provided the transition matrix has a suitable factorized form, namely that (the probability of) the next-time configuration of each site variable depends only on the current-time configuration of its neighbors. The method has been the subject of a detailed testing work by an independent group [19], still on various types of kinetic Ising-like models, demonstrating a remarkable accuracy for both the transient and the steady state, with just a slight breakdown in the notoriously difficult case of spin glasses at low temperature. Let us stress the fact that the derivation of the diamond approximation in reference [18] introduces a novelty item with respect to a “standard” PPM approach [20], in that it removes an unnecessary constraint in the choice of maximal clusters, potentially opening the way to new approximation schemes.

In this paper we analyze different alternative approximations that can be worked out in the framework of a path-probability approach along the lines of reference [18]. This analysis leads in particular to generalizing the diamond approximation in such a way that the transition probability of a given site variable may also depend on the current configuration of the site itself. This step forward is necessary for the method to be applicable to models with this feature, such as, for instance, epidemic models or the like. We also trace a number of connections with different approximation schemes, previously derived either by the conventional “heuristic” approach or by the “standard” PPM approach. In the second part of the paper, we evaluate the performance of these different approximations in terms of both accuracy and computational complexity, with respect to Monte Carlo simulations. As a test model, we consider a susceptible-infected-susceptible (SIS) epidemic model on large random graphs, possibly with heterogeneous infection parameters. It turns out that the generalized diamond approximation outperforms simpler, previously known approximation schemes, with only a slight increase in the computational complexity.

Let us note that our investigation is carried out under the assumption of a discrete-time dynamics, but all the approximations can be applied as well to describe continuous-time processes. Concerning this issue, we prove in particular that the (generalized) diamond approximation actually degenerates, in the continuous-time limit, into a simpler (pairwise) approximation. From the technical point of view, the latter result entails that, even though it is still possible to get a more accurate approximation for continuous-time models, one necessarily has to resort to even larger maximal clusters, which of course requires a larger computational effort.

The paper is organized as follows. Section 2 contains the main theoretical issues. Sub 2.1 defines the most general model we can deal with, which we denote as the reference model. In sub 2.2 we introduce the thermodynamic formalism and the basic concepts of the cluster variational method. In sub 2.3 we describe the specific approximations which are the subject of the present paper, in particular the one which generalizes the diamond approximation. Most analytical details, along with a scheme of the resulting computation procedures, are reported in A (generalized diamond approximation) and B (other approximations). C deals with the special case of continuous-time processes, mentioned above. Section 3 describes the numerical tests, whereas Section 4 contains a summary of the results and some concluding remarks.

2 Cluster-variational approach

2.1 The reference model

As mentioned above, we consider quite a generic model, made up of discrete random variables ${\xi_{i}}^{(t)}$, associated with the vertices (or sites) ${i=1,\dots,N}$ of a given undirected graph, and depending on a discrete time index ${t=0,\dots,\tau}$, where $\tau$ may possibly tend to infinity. This collection of random variables may be regarded as a multivariate stochastic process ${\xi}^{(t)}\equiv[{\xi_{1}}^{(t)},\dots,{\xi_{N}}^{(t)}]$, whereas a particular realization of such a process ${x^{(0)}\!,\dots,x^{(\tau)}}$ may be denoted as a path. Accordingly, we define the path probability as the probability of each given path, namely

[TABLE]

In the assumption of a Markovian dynamics, the path probability takes the general form

[TABLE]

where ${p}^{(0)}(x^{(0)})$ is the initial condition (i.e., the probability of the initial configuration $x^{(0)}$) and ${w}^{(t)}(x^{(t+1)}|\,x^{(t)})$ are the transition probabilities from $x^{(t)}$ to $x^{(t+1)}$ (i.e., the conditional probability of the latter configuration at time ${t+1}$, given the former at time ${t}$).

The graph structure defines which variables in the “current-time” configuration are actually relevant to determine the “next-time” configuration, and therefore the structure of the transition probabilities. In many physical models, such as those mentioned in the introduction, the next-time configuration of each given site depends only on the current-time configuration of its neighbors and also possibly of the same site. In formulae, this statement means that the transition probabilities can be written in the following factorized form

[TABLE]

where the product runs over all sites, $x$ and $y$ (respectively, $x_{i}$ and $y_{i}$) denote two generic configurations of the whole graph (respectively, of the $i$ site)222If possible, throughout the paper we shall denote a current-time configuration by a letter $x$ and a next-time configuration by a letter $y$., and $x_{i,\partial i}$ is a shorthand for a list of configuration variables $x_{i},\{x_{j}\}_{j\in\partial i}$ (i.e., those associated with the $i$ site and its neighbors). We shall also assume that the initial configurations of different sites are statistically independent, namely that

[TABLE]

The latter assumption could also be partially relaxed, but usually there is no specific interest in doing so.

2.2 Thermodynamic formalism and cluster variational method

Let us now consider the following variational functional

[TABLE]

where the sum runs over all possible paths, the path probability $\wp$ is regarded as a variational parameter, and the initial distribution ${p}^{(0)}$ and the transition probabilities ${w}^{(t)}$ are assumed to be “input parameters”, incorporating our knowledge of the stochastic dynamics. From the mathematical point of view, this functional is a so-called Kullback-Leibler divergence between the two sides of equation (2), and is known to have an absolute minimum (zero) where the two distributions are equal, that is when equation (2) is verified. Introducing the Shannon information entropy associated with the path probability, namely

[TABLE]

we can rewrite the variational functional as

[TABLE]

where we have also defined

[TABLE]

Equation (7) points out that the variational functional $\mathcal{F}[\wp]$ is analogous to a free energy for a system in which time plays the role of an extra dimension, and whose energy function $\varepsilon$ is defined by (8). More precisely, such a system can be viewed as a “stack” of ${\tau+1}$ “copies” of the original graph, each one labeled by a time index ${t=0,1,\dots,\tau}$, so that we shall denote it as a space-time thermodynamic system. The path probability minimizing $\mathcal{F}[\wp]$ is analogous to a Boltzmann distribution for the space-time thermodynamic system.

Up to this point, we have got just a reformulation of the problem, which is of no practical use in a fully general case. Nevertheless, our assumptions (3) and (4) about, respectively, transition probabilities and initial conditions allow us to rewrite the energy function as

[TABLE]

The latter expression points out that (the logarithm of) the elementary factors of the transition probabilities and of the initial conditions play respectively the role of interaction energies and external fields, and that both kinds of term are local in both space and time. As a consequence, since the energetic term of the free energy (7) is local, we expect that the (inherently nonlocal) entropic term could be effectively approximated by a truncation of its cumulant expansion, according to Kikuchi’s cluster-variational method [3, 21]. Each cumulant is associated with a given cluster (i.e., a given subset of lattice sites), so that different levels of approximation may be defined by the set of maximal clusters retained in the expansion. The approximate entropy can be finally recast into a linear combination of cluster entropies, i.e., information entropies of cluster distributions (the latter being obviously defined as joint probability distributions for the random variables associated with all sites in the given cluster). Thus, a generic CVM entropy takes the following form [3, 21]

[TABLE]

where the sum runs over the set $\mathsf{R}$ of all maximal clusters and their subclusters, ${{\wp}_{\alpha}}$ denotes the cluster distribution for cluster $\alpha$ (i.e., a particular marginal of $\wp$), and $\mathcal{S}[{{\wp}_{\alpha}}]$ is the corresponding information entropy. The coefficients $a_{\alpha}$ of the combination, usually denoted as Möbius numbers [3, 22], or (over)counting numbers [4, 22], depend only on the graph geometry and the choice of maximal clusters. They can be determined by a simple sum rule, which can be stated as follows [21]

[TABLE]

where $\alpha$ may be any maximal cluster or subcluster, and the sum runs over all (maximal clusters or subclusters) $\beta$ containing $\alpha$. The rule shows in particular that the coefficient of any maximal cluster is equal to $1$, because in that case the sum contains only one term.

Let us finally note that the choice of maximal clusters is generally driven by physical intuition, since it defines the range of correlations that a given approximation will be able to take into account. Therefore, in the framework of a space-time system like the one we are interested in, maximal clusters should necessarily extend over different time steps, in order to take into account time correlations. Kikuchi’s path probability method (PPM) can be regarded, in this framework, as a particular way of choosing the maximal clusters, namely, one first chooses maximal clusters for a one-time system, then every maximal cluster of the space-time system is the union of a maximal cluster at time $t$ with its translation at time ${t+1}$. As already noticed in reference [18], other criteria for this choice are possible, and may be more convenient, depending on the structure of the transition probabilities. In the following we shall see on test models that indeed one can obtain a considerably better tradeoff between accuracy and computational complexity, with respect to the conventional PPM choice.

2.3 Specific approximations

In this section we introduce the specific approximations (i.e., the different choices of maximal clusters), whose performance is analyzed in the following. First of all, let us define precisely the types of clusters (either maximal or not) that appear in at least one of the CVM entropies considered. Looking at table 1 we see that each cluster extends over either a single timestep or two consecutive ones. Of course, since we assume that the graph does not change with time, clusters at different timesteps are simply translations of one another. We introduce specific symbols (letters) to denote the different cluster types and the related probability distributions; such symbols are defined in table 1 as well.

As usual for the CVM, the simplest choice of maximal clusters is such that the latter coincide with the groups of variables coupled by the elementary interaction terms appearing in the energy function. Looking at equation (9), it is evident that these clusters coincides with the P-clusters, defined in table 1, for every ${t=0,\dots,\tau-1}$. We shall call this the P approximation. This approximation can be improved by including in the set of maximal clusters also the Q-clusters, which are expected to take into account the time-correlation of two (spatially-correlated) nearest neighbors. This will be denoted as the PQ approximation. A further improvement can be achieved by adding also the R-clusters, which are expected to better take into account the correlations among all the neighbors of a given site (at the same timestep), because all of them are influenced by the configuration of the “central site” at the previous timestep. The latter is the PQR approximation. Finally, looking at table 1, we recognize that every P, Q, and R cluster is contained in at least one M cluster, so that choosing the set of all M clusters as maximal clusters may lead to an even more accurate approximation (M approximation). Note that the M approximation is the only choice (among those considered here), which meets Kikuchi’s PPM criterion (as stated in the previous subsection), and that therefore can be regarded as an instance of PPM. Note also that it is the minimal PPM instance for the current reference model, because the M clusters are the minimal clusters, which simultaneously meet Kikuchi’s criterion and “contain” the elementary interaction terms.

Making use of the sum rule (11), one can determine the Möbius numbers, and therefore the approximate entropy functional, for each specific approximation. This is a rather standard calculation, whose results are summarized in table 2. Note that, as previously mentioned, all (though not only) the Möbius numbers associated with the maximal clusters are equal to unity. Note also that, in the M approximation, P, Q, and R clusters are not maximal, so their Möbius numbers do not need to be (and in fact they are not) equal to unity. In the remainder of this section we report the entropy functionals in formulae. As far as the P approximation is concerned, we have

[TABLE]

where $\sum_{i}$ and $\sum_{ij}$ denote sums over all sites and over all nearest neighbor pairs, respectively, and ${\mathcal{S}[\cdot]}$ denotes the Shannon entropy of each given probability distribution. The arguments ${P,Z,A}$ mean that the functional $\mathcal{S}_{\mathrm{P}}$ depend on all P-, Z-, and A-type cluster distributions, i.e., $P$ is a shorthand for $\{{P_{i}}^{(t)}:\forall t=0,1,\dots,\tau-1,\forall i\}$, and so on. With analogous notations, the PQ entropy approximation reads

[TABLE]

where $\sum_{i,j\in\partial i}$ denotes a double sum over all sites and over all neighbors of each given site. Furthermore, the PQR entropy reads

[TABLE]

whereas the M entropy reads

[TABLE]

For each approximation, the free-energy functional is obtained from equation (7) by replacing the entropic term $\mathcal{S}$ with one of the approximate entropies. As previously mentioned, the energetic term does not require any further approximation, because the maximal clusters are properly chosen to “contain” the elementary coupling terms of the energy function (8).

The minimization of the resulting functionals is far from being a straightforward problem. This is ultimately due to cluster overlapping, which entails that compatibility constraints between cluster distributions must be satisfied. In general, the constrained minimization problem for a generic CVM free energy can be solved iteratively by a message-passing procedure denoted as generalized belief propagation [3, 4] (which is not guaranteed to converge), or by more complex (provably convergent) methods [22]. Nevertheless, the peculiar structure of the couplings, arising in a space-time system describing a Markovian stochastic dynamics like that of our reference model, allows one to devise a much simpler procedure, which does not require iterative refinement, but only a “natural” sequence of iterations, starting from the initial conditions and following the (probabilistic) time evolution of the system. We report the details of such a procedure in A, only for the PQR approximation. The differences occurring in the other approximations (which are indeed very much analogous) are discussed in B. Furthermore, analogous procedures can be used to investigate continuous-time processes (in a discretized form), and, as mentioned in the introduction, in such a case one can argue that the PQR and PQ approximations become equivalent. The latter issues are discussed in C.

2.4 Related works

Before switching to a performance evaluation of the different approximations described above, we find it worth summarizing a number of connections with previous works, which emerge from our analysis. As mentioned in the introduction, in most cases equivalent approximations have been put forward without resorting to a variational scheme, but just on the basis of heuristic reasoning, and/or for specific model dynamics. In such cases, our approach retains some interest in that it offers a more general and unified view of the matter.

First of all, it is possible to recognize that the P and PQ schemes turn out to be respectively equivalent to the classical (non-variational) mean-field and pairwise approximations, which have been used, for instance, to investigate epidemic spreading phenomena on networks [5, 6, 17]. In such a context, the P approximation is customarily denoted as quenched mean field, in order to distinguish it from a naive mean-field approach, in which the infection probability is assumed to be a node-independent quantity [5], and from the popular heterogeneous mean-field approach by Pastor-Satorras and Vespignani [23], in which nodes are distinguished only on the basis of their degree. For epidemics on networks that are expected to be inherently homogeneous, such as random-regular graphs or regular lattices, the same approximations have been considered also by Petermann and De Los Rios [7], and other authors [24]. We point out these equivalences in B.1 along with some details about the respective iterative procedures. Secondly, the M approximation turns out to be an instance of Kikuchi’s (variational) path-probability method, as we have mentioned in the previous subsection, but it has also been derived independently (for the case of homogeneous epidemics) again by Petermann and De Los Rios [7], who denoted it as star approximation. The latter equivalence is detailed in B.2. Furthermore, let us note that, in the special case of a spatially one-dimensional system (i.e., when the graph is a linear chain), a recent analysis of nonlinear Voter models carried out by Schweitzer and Behera [12] makes use of two different approximation schemes (denoted as triplet and quintuplet approximations), which turn out to be special cases of the PQ and M approximations, respectively.

Another special case of our reference model, which nonetheless includes several models of physical relevance (such as the kinetic Ising model, or the ordinary Voter model), is that in which the elementary transition probability for a given site $i$ (namely, ${w_{i}}^{(t)}(y_{i}|x_{i,\partial i})$) does not depend on the current-time configuration of the site itself $x_{i}$ (so that we can denote it as ${w_{i}}^{(t)}(y_{i}|x_{\partial i})$). In other words, the next-time configuration of each site is conditioned only by the current-time configurations of its neighbors. For this reason, we shall denote this case as a purely neighbor-conditioned dynamics. As mentioned in the introduction, the star and diamond approximations of reference [18] refer to this kind of dynamics. Actually, it turns out that, for a generic purely-neighbor-conditioned dynamics, the PQ and P approximations degenerate into each other, and both coincide with the star approximation of reference [18]. Similarly, the M and PQR approximations degenerate into each other, and coincide with the diamond approximation. All these connections are detailed in A.4 and B. As a further consequence, we can state that in this case the diamond approximation is itself an instance of Kikuchi’s path probability method.

3 Test models

In this section we test our approximations against Monte Carlo simulations on specific instances of a SIS epidemic model on (sparse) graphs drawn from well known random ensembles, namely the Erdös-Rényi ensemble and the random-regular ensemble. In the former case, the degree distribution is poissonian, whereas in the latter case all vertices have the same degree, so that we also denote the two kinds of graphs as poissonian and uniform, respectively.

We consider a SIS model with parallel dynamics, whose definition is as follows. Each site variable ${\xi_{i}}^{(t)}$ can take two possible values, ${x_{i}=0,1}$, meaning that the $i$ site is susceptible or infected, respectively. An infected site $i$ has a probability $\mu_{i}$ to recover (i.e., to become susceptible) at the next timestep. Moreover, a susceptible site $i$ has a probability $\beta_{j\to i}$ to become infected because of any of its neighbors $j$, if the latter is infected at the current timestep, and [math] otherwise. In other words, the probability of an infection event from $j$ to $i$ is ${x_{j}\beta_{j\to i}}$. Since infection events by different neighbors are not mutually exclusive, the total probability of a given site $i$ to become infected can be more easily computed via its complementary probability (i.e., the probability that $i$ remains susceptible), which is equal to the probability that none of the neighbors transmits the contagion. Assuming independence among different infection events, the transition probability to a susceptible state can be finally written as

[TABLE]

whereas, by normalization, the transition probability to an infected state will be

[TABLE]

We consider two fixed graphs with ${N=1000}$ vertices, a poissonian one with average degree ${c=4}$ and a uniform one with degree ${c=4}$. As far as model parameters are concerned, in both cases we choose a uniform recovery probability, namely ${\mu_{i}=\mu=0.5}$, whereas the infection probability is chosen to be homogeneous in the former case (namely, ${\beta_{j\to i}=\beta=\mu/3}$) and heterogeneous in the latter (namely, $\beta_{j\to i}$ drawn at random with uniform distribution between [math] and $\mu$). Note that, as a result, both systems turn out to be heterogeneous: in the Erdös-Rényi system, heterogeneity arises from the random graph structure, whereas, in the random-regular system, heterogeneity arises from the randomness of infection probabilities. The initial conditions are chosen in such a way that both systems have on average 1% of infected sites, but with a slight difference. In the poissonian case, we use “quenched” initial conditions, that is, 1% of (randomly drawn) sites are assigned the initial distribution

[TABLE]

whereas the remaining 99% are assigned the distribution

[TABLE]

In the uniform case we use “annealed” initial conditions, that is, all sites are assigned the distribution

[TABLE]

In the simulation, quenched conditions correspond to all realizations initialized with exactly the same configuration (with 1% of infected sites) drawn at random at the beginning, whereas annealed conditions correspond to a different initial configuration (still with 1% of infected sites) drawn at random for each realization. Note that there is no specific reason for associating a particular type of initial condition with a particular type of graph. We simply observe that different initial conditions give rise to a slightly different (probabilistic) short-time behavior of the system, appearing independently of the graph type, and we are meant to analyze whether such a behavior is more or less accurately reproduced by the approximations.

For each approximation, each graph and parameters choice, and each site ${i=1,\dots,N}$, we can compute the probability of the infected state (infection probability) as a function of time, namely

[TABLE]

where the single-site (A-cluster) distribution ${A_{i}}^{(t)}(x_{i})$ can be determined as a marginal of a maximal-cluster distribution (for instance the P cluster):

[TABLE]

The maximal cluster distributions are determined in turn by the minimization procedure of the free-energy functional, as mentioned in the previous section (and detailed in A). In the simulations, the infection probability is evaluated as an empirical average of ${\xi_{i}}^{(t)}$ over ${10000}$ realizations (let us denote it as ${\bar{\rho}_{i}}\vphantom{\rho_{i}}^{(t)}$). The steady-state is easily identified in the approximations, since we observe that, as $t$ increases, ${\rho_{i}}^{(t)}$ tends to a well-defined plateau value $\rho_{i}$. Conversely, in the simulations ${\bar{\rho}_{i}}\vphantom{\rho_{i}}^{(t)}$ exhibits a fluctuating behavior for large $t$, so we compute the steady-state value $\bar{\rho}_{i}$ as an average over a large number of timesteps, which allows us to get a suitably small statistical uncertainty. We then compute the relative errors between theory and simulation (in the steady state) as ${|\rho_{i}/\bar{\rho}_{i}-1|}$. We display all these data in figure LABEL:fig:ss, for each approximation, as a function of the corresponding data for the PQR approximation. Note that, as a consequence, PQR data fall on a straight line. More importantly, this way of showing data allows us to argue at a glance that the M approximation is systematically (i.e., for all sites) better than the PQR approximation, which is in turn systematically better than PQ. Let us also observe that, in the poissonian graph case (left panel), the PQ and P data are so separated that also PQ clearly turns out to be systematically better than P. For the uniform graph case (right panel), the boundary between the two data “clusters” appears a little more fuzzy, so that one might suspect that, at least for some sites, the P approximation could give a better approximation than PQ. Nevertheless, plotting as a function of the PQ data (right panel inset) shows that this is not the case, i.e., that PQ is still systematically better than P. From figure LABEL:fig:ss, one can also assess that the accuracy improvements that one can obtain by “upgrading” the approximation level from P, to PQ, to PQR are, in relative terms, of the same order of magnitude, whereas a further upgrade to the M approximation achieves a minor improvement.

In order to assess also the transient behavior, we display in figure LABEL:fig:wc the whole time evolution of the infection probability ${\rho_{i}}^{(t)}$ for a specific site $i$. In particular, we choose the site $i$ for which the PQR approximation gives the largest relative error (worst case) in the steady state. We can observe that the systematic accuracy improvement, already noticed for the steady state, also takes place in the short-time evolution. Such a behavior is indeed common for all sites. We also note that the PQR and M results are almost indistinguishable at this scale, whereas they provide a noticeable improvement over the PQ (pair) approximation. Moreover, we can observe that, in the case of quenched initial conditions (left panel), the infection probability remains sharply zero up to a given timestep (the $i$ site being one of those initially susceptible). In such a case, even the simplest P approximation (mean field) correctly reproduces the very short time behavior (including the precise timestep at which the $i$ site first exhibits a nonzero infection probability), but otherwise it rapidly deviates from the trend predicted by the simulations.

One more test case we consider is that of a uniform graph (still with ${c=4}$), with homogeneous recovery probability (${\mu_{i}=\mu=0.5}$) and, as well, homogeneous infection probability (${\beta_{j\to i}=\beta}$). Initial conditions are chosen to be of the annealed type (as defined above). In such a system we expect a fully homogeneous infection probability (i.e., ${{\rho_{i}}^{(t)}=\rho^{(t)}}$ for all $i$), so that in the simulation we can evaluate the latter by an extra average over the sites, thus obtaining much less noisy results. In order to enhance the “noise-reduction” effect, we also choose a larger graph size, namely ${N=100000}$, and we verify that, with an average over just $100$ realizations, the statistical uncertainty turns out to be negligible at the scale of the graphs. We display the results in figure LABEL:fig:uu for the specific cases ${\beta=\mu/2}$ and ${\beta=\mu/3}$. In the former case (left panel), we can observe that almost all the approximations (except possibly P) give a rather accurate prediction of the steady-state infection probability. This may be partially ascribed to the fact that in this setting the steady state is expected to be an equilibrium state, at odds with the previously considered cases. It is remarkable that, even though the M, PQR, and PQ approximations are all very accurate in the steady state, the PQ approximation shows quite a relevant discrepancy during the transient. Furthermore, upon decreasing $\beta$ (right panel) one can observe that the various approximations exhibit increasingly different results also in the steady state. This is related to the fact that the model approaches its critical point, at which the active (endemic) steady state disappears, and the long-time infection probability becomes rigorously zero. In analogy with equilibrium phase transitions, different approximations give different predictions for the critical $\beta$ value, as already pointed out in the continuous-time case [6, 17].

Let us finally compare the different approximations in terms of computational complexity. In figure LABEL:fig:tempi (left panel) we display the CPU time (Intel Core 2 Duo) needed to perform $1000$ evolution timesteps of a uniform graph with ${N=1000}$ vertices, as a function of the (fixed) degree $c$. Of course, greater accuracy corresponds to larger computational effort, though not with a simple proportionality. In particular, we note that the M approximation entails a very large increase of the complexity, in spite of the fact that its accuracy is just slightly better than that of PQR (at least for the specific epidemic model considered here). The scaling behavior with respect to $N$ is obviously linear, so that we do not need to report detailed results, whereas it turns out to be exponential with respect to $c$. In principle, the complexity scales with the number of configurations of the largest cluster retained in the CVM expansion, so that the P, PQ and PQR approximations are all exponential in $c$, whereas the M approximation is exponential in ${2c}$. From the theoretical point of view of computational complexity, this is an irrelevant difference, but from the practical point of view this means that the computational cost of the M approximation may rapidly become similar to that of a simulation, upon increasing the degree, as shown in figure LABEL:fig:tempi. On the other hand, the P, PQ and PQR approximations present an extra computational advantage, which can be stated as follows. If the transition probabilities can be written as a product of terms depending each on a single neighbor (or, as well, as a sum of a few terms of this kind), simple dynamic programming precautions allow one to reduce the complexity scaling to be linear in the degree. Details about this issue are beyond the scope of the current paper, but we only mention the fact that the transition probabilities of the SIS model satisfy the above prescription. Execution times of the optimized programs are reported in the right panel of figure LABEL:fig:tempi.

4 Conclusions

In this paper we have performed a detailed investigation of different types of mean-field-like approximations for stochastic dynamics on graphs, within the unifying framework of a cluster-variational approach. In analogy with its equilibrium counterpart, this approach naturally opens the way to improving the approximation, by a suitable choice of the maximal clusters, which basically define the range of correlations that a given approximation scheme is able to take into account. On the one hand, we have pointed out that the simplest choices of maximal clusters (namely, the P and PQ approximations, according to our nomenclature) coincide respectively with the ordinary mean-field and the pairwise approximations, previously derived in different contexts without resorting to a variational scheme. On the other hand, we have proposed a slightly more complex cluster choice (PQR in our nomenclature), which provides a generalization of the so-called diamond approximation, to dynamics that are not purely neighbor-conditioned (i.e., such that the transition probability of a given site variable may also depend on the current configuration of the site itself). Furthermore, we have considered an even larger cluster choice (M in our nomenclature), which turns out to be an instance of Kikuchi’s path-probability method. We have shown explicitly that, in the case of a purely neighbor-conditioned dynamics, the P (mean-field) and PQ (pairwise) approximations become equivalent, and similarly the M approximation (path-probability method) degenerates into the PQR (generalized diamond) approximation.

We have tested the different approximations on a (discrete-time) epidemic model of the SIS type (i.e., with recurrent dynamics) on random graphs, possibly with heterogeneous infection parameters. We have compared the results with Monte Carlo simulations, and demonstrated that increasing the maximal cluster size leads to a systematic improvement of the level of accuracy. Such an improvement can obviously be obtained at the cost of a higher computational complexity, but not with a simple proportionality. In particular, the PQR approximation yields nearly the same accuracy as the M approximation, with a considerably lower complexity. We have also confirmed that the PQR approximation yields very good accuracy even in the transient regime, as previously reported for the diamond approximation in the case of kinetic Ising models [19]. Let us observe that our choice of a discrete-time dynamics is a bit artificial, as continuous time is the most typical setting of epidemic models (even though with some exceptions [25]). In our comparative analysis, this choice is meant to point out the difference between the PQR and PQ approximations, which we indeed prove to become equivalent in the continuous-time limit.

Let us finally mention an important alternative class of approximate methods for the same kind of problems, which go under the name of dynamic message-passing (DMP) or the like. Basically, these methods rely on the so-called dynamic cavity equations [26, 27, 28], which are exact on tree graphs (and therefore likely accurate on treelike graphs), but require taking into account the whole time evolution of each given site variable. Different approximations over the time dimension can then be worked out to make the problem computationally tractable. Such alternative approximations include for instance the early “dynamic cavity” approximation by Aurell and Mahmoudi [29, 30, 31], the Markovian closure scheme by Del Ferraro and Aurell [27, 28], and the very refined (and computationally demanding) matrix-product algorithm by Barthel and coworkers [32]. An analogous (but apparently different) DMP approach, directly formulated for the continuous-time SIS model, has also been put forward recently [33]. Interestingly, the latter method has been reported to perform slightly better than the pairwise (PQ) approximation, but mainly in the stationary state, whereas it seems to be generally less accurate in transient states. On the other hand, the aforementioned Markovian-closure scheme [27] (which can in principle be applied to any kind of dynamics falling within our reference model) has been thoroughly tested [19] on various (Glauber) kinetic Ising models, and reported to perform (slightly but quite systematically) worse than the diamond (PQR) approximation. Note that, as in this case the dynamics is purely neighbor-conditioned, the pair approximation is no longer a good term of comparison, because it degenerates into the ordinary mean field, as we show in B.1. We have not included in our paper a detailed comparative analysis with respect to dynamic-message-passing approaches, because, as previously mentioned, the main focus of our work is on the unifying framework of the cluster variational method. Nevertheless, such an investigation might be in fact worth being performed, so that we shall consider it as a subject for future work.

Appendix A Details of the PQR approximation

In this appendix we report in detail the method we have used to determine the minimum of the cluster-variational functional for the case of the PQR approximation. As we shall see, this method naturally leads to an “iterative” procedure (the one we have implemented numerically), which, starting from the initial conditions, builds up the cluster distributions at subsequent timesteps, up to the eventual steady state. Let us first discuss the issue of compatibility constraints, which will be of use in the following.

A.1 Compatibility constraints

The probability distributions for maximal clusters and subclusters (which for the PQR approximations have been denoted as $P,Q,R,S,T,U,V,Z$) are of course not independent variational parameters, because they have to satisfy (besides normalization) a number of compatibility constraints. By compatibility we mean that the marginalization of probability distributions of overlapping maximal clusters must give, for the same subcluster, the same distribution. Writing down these conditions for all subclusters of all maximal clusters, we obtain a sufficient (indeed redundant) set of constraint equations. Let us group together these equations depending on the type of maximal cluster involved. We understand that all the equations must hold for ${t=0,1,\dots,\tau-1}$.

P-type maximal clusters:

[TABLE]

Q-type maximal clusters:

[TABLE]

R-type maximal clusters:

[TABLE]

This set of conditions also entails direct compatibility equations between subcluster distributions (for overlapping subclusters). Let us write explicitly some of these equations below.

From equations (23a) and (23b), or equivalently (25c) and (25d), we have

[TABLE]

From equations (23b) and (23c), or equivalently (24a) and (24b), we have

[TABLE]

From equations (23c) and (23d), or equivalently (24a) and (24c), we have

[TABLE]

From equations (24c) and (24d), or equivalently (25a) and (25b), we have

[TABLE]

From equations (24d) and (24e), or equivalently (25a) and (25c), we have

[TABLE]

Let us finally remark that the largely redundant form, in which we have written the compatibility constraints, is not an end in itself. Indeed, it will enable us to get a very compact proof of the fact that our method actually produces compatible solutions. This will be shown, along with the computation algorithm, in A.3.

A.2 The variational functional

As mentioned in section 2.3, we derive the PQR free-energy functional (denoted as $\mathcal{F}_{\mathrm{PQR}}$) from equation (7), replacing the entropic term $\mathcal{S}$ with $\mathcal{S}_{\mathrm{PQR}}$ (equation (14)). Taking into account that we are interested in solutions respecting the constraint equations, specifically (23a), (24a), (24b), (25a), (25b), it is possible to rewrite the functional in such a way that the subcluster distributions appear only in the logarithmic terms, that is in the following form

[TABLE]

where the inner sums run over all possible configurations of P-, Q-, and R-clusters, respectively. Recall that the terms ${w_{i}}^{(t)}(y_{i}|x_{i,\partial i})$ are transition probabilities, so that by construction we assume they satisfy the following normalization conditions

[TABLE]

Note also that, given the Möbius numbers of the PQR approximation (see table 2 and equation (14)), the probability distributions ${S_{i}}^{(0)}$ (i.e. the S-cluster distributions at time ${t=0}$) and ${Z_{ij}}^{(0)}$ (i.e. the Z-cluster distributions at time ${t=0}$) actually do not belong to the entropic term. We have nonetheless included them in the free-energy expression, to take into account the part of energetic term that represents initial conditions, assuming therefore

[TABLE]

Now, a crucial observation is that the variational functional in the form (31) is made up of a sum of Kullback-Leibler divergence terms, and therefore that the absolute minimum conditions (in which every term equals [math]) are, for ${t=0,1,\dots,\tau-1}$,

[TABLE]

Actually, the minimum condition of a Kullback-Leibler divergence holds provided its two arguments (that is in our case the two sides of each equation (34)) are “good” (i.e. normalized) probability distributions. Assuming of course the normalization of all cluster distributions, it is easy to see that the right-hand sides of equations (34a), (34b) and (34c) are actually good distributions, if, respectively, the normalization condition (32) and the compatibility conditions (27) and (29) are verified.

Let us finally observe that equations (34) are somehow “weak”, meaning that they do not define a point in the space of variational parameters, but rather a manifold. In order to determine completely the variational parameters, that is the cluster distributions ${P,Q,R,S,T,U,V,Z}$, we have to intersect the aforementioned manifold with that defined by the constraints, that is to solve simultaneously the minimum equations (34) with the compatibility equations (23), (24), (25) (whence also (26), (27), (28), (29), (30)) and the normalizations. Let us stress the fact that it is not necessary to resort to the Lagrange multiplier method to solve the constrained minimum problem, just because it turns out that the intersection between the two manifolds is nonempty. In general, the latter is not a necessary condition for the existence of a constrained minimum. In our specific problem, such a condition arises from the fact that the exponentials (Boltzmann factors) of the coupling terms are actually (transition) probabilities, and therefore that they enjoy the normalization property (32).

A.3 Computation algorithm

In this section we present an iterative procedure, which solves the problem described above. In particular we see that, starting with a set of $S$ and $Z$ distributions satisfying (26) at a given timestep $t$ (initially ${t=0}$, with the distributions defined by (33)), one can determine all the other distributions $P,Q,R$ and $T,U,V$ at timestep $t$, as well as $S$ and $Z$ at timestep ${t+1}$, satisfying all the required equations at time $t$ and (26) at time $t+1$. Moreover, all the computed distributions turn out to be correctly normalized.333In the numerical implementation, it turns out to be necessary to restore normalization explicitly. The steps of the procedure are described in detail below, and summarized as a flow chart in figure LABEL:fig:flow-chart.

Let us begin, at a given timestep $t$ (initially ${t=0}$) with a set of (normalized) S-cluster distributions ${S_{i}}^{(t)}(x_{i,\partial i})$ ($\forall i$) and Z-cluster distributions ${Z_{ij}}^{(t)}(x_{i,j})$ ($\forall ij$), that we assume to be compatible, i.e. satisfying equations (26) (which is obviously verified by the initial distributions (33)). 2. 2.

We determine the P-cluster distributions at time $t$ by means of equations (34a). Thanks to (32), we see that such distributions satisfy equations (23a). Moreover, since (23a) and (26) hold simultaneously, also (23b) turns out to be satisfied. We then determine the following marginals.

• –

T-cluster distributions ${T_{i,ij}}^{(t)}(y_{i},x_{i,j})$ ($\forall i,\forall j\in\partial i$) by equations (23c). Since (23c) and (23b) hold simultaneously, such distributions also satisfy (27).

• –

V-cluster distributions ${V_{i}}^{(t)}(y_{i},x_{i})$ ($\forall i$) by equations (23d). Since (23d) and (23c) hold simultaneously, also (28) turns out to be satisfied. 3. 3.

We determine the Q-cluster distributions at time $t$ by means of equations (34b). Thanks to (27), we see that such distributions satisfy equations (24a) and (24b). Moreover, since (24a) and (28) hold simultaneously, also (24c) turns out to be satisfied. We then determine the following marginals.

• –

U-cluster distributions ${U_{ij,i}}^{(t)}(y_{i,j},x_{i})$ ($\forall i,\forall j\in\partial i$) by equations (24d). Since (24d) and (24c) hold simultaneously, such distributions also satisfy (29).

• –

Z-cluster distributions (at the next timestep) ${Z_{ij}}^{(t+1)}(y_{i,j})$ ($\forall ij$) by equations (24e). Since (24e) and (24d) hold simultaneously, also (30) turns out to be satisfied. 4. 4.

We determine the R-cluster distributions at time $t$ by means of equations (34c). Thanks to (29), we see that such distributions satisfy equations (25a) and (25b). Moreover, since (25a) and (30) hold simultaneously, also (25c) turns out to be satisfied. We then determine the following marginals.

• –

S-cluster distributions (at the next timestep) ${S_{i}}^{(t+1)}(y_{i,\partial i})$ ($\forall i$) by equations (25d). Since (25d) and (25c) hold simultaneously, such distributions also satisfy (26) at time ${t+1}$. 5. 5.

We repeat the procedure from step (2) for the distributions at time ${t+1}$.

Let us stress the fact that, in the above procedure, all the compatibility conditions stated in A.1 have been either used explicitly to compute the marginals or proved to be satisfied by the computed distributions.

A.4 The diamond approximation

In this section we show that, as mentioned in section 2.4, in the case of a pure neighbor-conditioned dynamics the PQR approximation coincides with the diamond approximation, previously proposed by one of us in reference [18].

First of all it is useful to define some auxiliary clusters, namely, the B-clusters, made up of two neighbor sites at subsequent timesteps $t$ and ${t+1}$, for ${t=0,1,\dots,\tau-1}$, and the C-clusters, made up of all the neighbors of a given site at a given timestep $t$, for ${t=0,1,\dots,\tau}$. The related distributions are denoted as

[TABLE]

These distributions can be written as marginals of maximal cluster distributions, so that they inherit all the required compatibility conditions. We shall use in particular the following relationship

[TABLE]

which, together with

[TABLE]

also imply

[TABLE]

Let us now hypothesize that the $S$ and $Z$ distributions at a given timestep $t$ factor in the following fashion

[TABLE]

obviously assuming also that $C$ and $A$ distributions are compatible, that is they satisfy

[TABLE]

Retracing the steps of the algorithm described in the previous section, let us show that the hypothesized factorizations are stable along the time evolution (i.e., our hypotheses imply the same factorizations at timestep $t+1$), and the evolution equations are identical to those of reference [18].

We easily see that the hypotheses (40) and (41) satisfy the compatibility (26) between $S$ and $Z$ distributions at time $t$. 2. 2.

Equation (34a) with the factorization (40a) reads

[TABLE]

Therefore, equation (23c) together with (37) gives

[TABLE]

where

[TABLE]

whereas equation (23d) together with (38) gives

[TABLE]

where

[TABLE] 3. 3.

Equation (34b) with the factorizations (43) and (40b) reads

[TABLE]

Therefore, equation (24d) together with (39) gives

[TABLE]

whereas equation (24e) together with (39) gives

[TABLE]

We thus argue that the $Z$ distributions at time ${t+1}$ take a factorized form analogous to (40b). 4. 4.

Equation (34c) with the factorizations (48) and (45) reads

[TABLE]

Therefore, equation (25d) proves that the $S$ distributions at time ${t+1}$ take a factorized form analogous to (40a), namely

[TABLE]

where

[TABLE]

Note that the “diamond clusters” of reference [18] never appear explicitly in the above discussion. Nevertheless, we can see that, plugging equation (52) into (44) and (46), one obtains evolution equations for B- and A-cluster distributions, which exactly coincide with those one could obtain, in reference [18], replacing the diamond cluster distribution (26) into the marginalization equations (25) and (23).

Appendix B Other approximations

In this appendix we discuss the other approximations we have introduced in the text, namely the PQ and P approximations, and the M approximation. We do not go into the same level of detail, as done for the PQR approximation, because the derivations are in fact very much analogous. We only give the final results in terms of the respective computation procedures.

B.1 PQ and P approximations

The PQ approximation can be regarded as a simplification of the PQR approximation, in which the S-cluster distributions are factorized in terms of Z-cluster (nearest-neighbor pair) and A-cluster (single site) distributions, according to

[TABLE]

As a consequence, one can replace step (4) of the PQR procedure with the previous equation (evaluated at timestep ${t+1}$), where the Z-cluster distributions have been determined at step (3), and the A-cluster distributions, compatible with the former, i.e. satisfying

[TABLE]

can be written by the obvious marginalization (38). In conclusion, as previously mentioned, the PQ approximation turns out to be equivalent to the pairwise approximation, considered by different authors [5, 6, 7]. Such an equivalence might indeed be not fully evident, because in the cited papers the approximations are immediately specialized to a particular model (namely, the SIS model) and to the continuous-time limit, whereas our treatment is more general, at least under the assumptions of a discrete-time dynamics and a transition probability factorized according to (3).

The P approximation can be regarded as a further simplification, in which also the pair distributions are factorized, as in equation (40b), so that the S-cluster distributions (53) read

[TABLE]

As a consequence, one can skip also step (3) of the PQR procedure, because there is no need to compute pair distributions at the next timestep, to be able to close the equations. Indeed, replacing equation (55) into (34a) and then into (38), one obtains a unique time-evolution equation for the single-site distributions, namely

[TABLE]

which coincides with the (quenched) mean-field approximation [5, 6, 7].

This argument also proves that, for a pure neighbor-conditioned dynamics, for which (40b) holds, the PQ approximation degenerates into the P approximation, and therefore that both coincide with the star approximation of reference [18], as previously stated in section 2.4.

B.2 M approximation

The iterative procedure for the M approximation turns out to coincide with that of the PQR approximation until step (3). Step (4) is modified in that the R-cluster distributions are no longer computed by equation (34c), but rather as marginals of the M-cluster distributions, namely

[TABLE]

the latter being previously computed as

[TABLE]

In the simpler case of a pure neighbor-conditioned dynamics, we have shown for the PQR approximation that the factorizations (40b), (42), and (43) hold. These equations remain valid also for the M approximation, because they stay within step (3) of the procedure, so that we can plug them into (58) and then into (57). Taking into account also equation (46) (still valid), we obtain exactly equation (50). This argument proves that, for a pure neighbor-conditioned dynamics, the M approximation degenerates into the PQR approximation, and therefore that both coincide with the diamond approximation of reference [18], as previously stated in section 2.4.

Appendix C Continuous-time limit

The approximation schemes we have presented in this paper can be used as well to analyze continuous-time processes (in a discretized form), or equivalently discrete-time processes with sequential dynamics. This can be done by defining the following transition probabilities

[TABLE]

where it is understood that ${\tau\to 0}$, whereas $\delta(\cdot,\cdot)$ denotes a Kronecker delta and $\tilde{w}{\vphantom{w}{}_{i}}^{(t)}(y_{i}|x_{i,\partial i})$ are suitable conditional probabilities, normalized in the proper way

[TABLE]

Equation (59) can be interpreted as follows. With a large probability (${1-\tau}$) each vertex variable $\xi_{i}^{(t)}$ stays in its current value, and with a small probability ($\tau$) it is allowed to undergo a transition, governed by the probabilities $\tilde{w}{\vphantom{w}{}_{i}}^{(t)}$. For the case of sequential dynamics, ${\tau=1/N}$ is the probability of choosing one specific vertex out of $N$. Alternatively, for the case of continuous-time processes, one can think of $\tau$ as a small time interval, and of $\tilde{w}{\vphantom{w}{}_{i}}^{(t)}(y_{i}|x_{i,\partial i})$ as the rate of the transition ${x_{i}\to y_{i}}$ (for $y_{i}\neq x_{i}$), which is assumed to depend on the neighborhood configuration $x_{\partial i}$.

As mentioned in the text, we can show that, in the limit ${\tau\to 0}$ (continuous-time limit), the PQR and PQ approximations become equivalent. Making use of equations (34a) and (23c) (step (2) of the computation algorithm) with the transition probabilities (59), and taking into account the compatibility equations (26), the T-cluster distributions can be written as

[TABLE]

where we have defined444Comparing equations (61a) and (62), one can argue that $\tilde{T}$ is the T-cluster distribution one would get if the transition probabilities were barely $\tilde{w}$. In the following we use the same notation as well for other cluster types.

[TABLE]

Making use of equations (34b) and (24e) (step (3) of the computation algorithm), still with the transition probabilities (59), the Z-cluster distributions at the next timestep (now denoted as ${t+\tau}$ rather than ${t+1}$) can be written as

[TABLE]

where we have defined

[TABLE]

Moreover, by (34a) and (38) we get

[TABLE]

where

[TABLE]

In principle, taking the limit ${\tau\to 0}$ of equations (63) and (65), one can obtain a system of ordinary differential equations for the pair and site distributions. Indeed, equations (63) and (65) can be viewed as discretized forms thereof, which are directly suitable for numerical implementation.

Let us observe that equations (63) and (65) hold for both the PQ and PQR approximations, since steps (2) and (3) of the computational procedure are identical, as mentioned in B.1. For the PQ approximation, these equations are immediately closed, because the $S$ distribution (appearing in the expressions of $\tilde{A}$ and $\tilde{B}$) directly depends on $Z$ and $A$, according to equation (53). As we shall see below, the reason why the PQ and PQR approximations turn out to be equivalent is that the difference between the $S$ distributions for the two cases is $\mathcal{O}(\tau)$ for ${\tau\to 0}$, so that, being multiplied by $\tau$ in both equations (63) and (65), its contribution becomes irrelevant. Let us first consider the PQ approximation. Equation (53) at time ${t+\tau}$ reads

[TABLE]

Plugging equations (63) and (65) into the latter, we easily get

[TABLE]

As far as the PQR approximation is concerned, the $S$ distribution is computed by equation (25d) as a marginal of the $R$ distribution, which in turn depends on the $U$ and $V$ distributions via equation (34c) (step (4) of the numerical procedure). Using equations (24d) and (34b) (step (3) of the procedure) with the $T$ distribution expressed by (61b), we can write

[TABLE]

hence, from equations (29) and (54),

[TABLE]

Plugging equations (69) and (70) into (34c), we get

[TABLE]

The latter expression, marginalized via (25d), yields again (68), which is equivalent to the PQ result.

We warmly express our thanks to Gino Del Ferraro, Eduardo Dominguez and Federico Ricci-Tersenghi for useful discussions on the subject.