A note on parallel sampling in Markov graphs

Abstract

The paper proposes the use of parallel computing for Markov graphs as a subclass of exponential random graph models where the network statistics induce a conditional independence structure amongst the edges of the network. This conditional independence allows simulation of edges in parallel using multiple computing cores. Simulation in Markov models is helpful, since parameter estimation cannot be carried out analytically but requires simulation-based routines such as Markov chain Monte Carlo. In particular in large networks this can be computationally very demanding or even infeasible. Therefore, numerical enhancements are useful to accelerate computation.

Appendices

Appendix A: Proof of asymptotic Markovian structure in log changes

The Markovian independence assumption is violated in the example when using log changes instead of simple sums, as in the 2-star or triangle. We can however proof that we obtain an asymptotic independence structure which justifies parallel draws in large networks. We show that the conditional log odds ratio for $P(Y_{i,j}, Y_{l,m}|y_{rest})$ is of order $O_p(N^{-1})$ if $\{i,j\} \ne \{l,m\}$ and $O_p(1)$ otherwise where $y_{rest}$ refers to the network except of edges $Y_{i,j}$ and $Y_{l,m}$. In other words, we obtain asymptotic independence for large networks. To simplify notation we first look at $Y_{1,2}$ and $Y_{3,4}$. The log odds ratio for these ties which do not share a node results through

where

Note that $a_{i,j}(N)$ and $b_{i,j}(N)$ are of order $O_p(N)$ with N as number of nodes. Taylor series approximation yields

$$\begin{aligned} \log \left( 2 + a_{i,j}(N) \right) - \log \left( a_{i,j}(N) \right) = \frac{2}{ a_{i,j}(N)}. \end{aligned}$$

(14)

Similarly, the component multiplied with $\theta _3$ is of order $O\left( b_{i,j}^{-1}(N)\right)$ and hence decreases to zero. This implies asymptotic independence of $Y_{1,2}$ and $Y_{3,4}$ given the rest of the network.

Let us now look at the pair $Y_{1,2}$ and $Y_{1,3}$, that is we look at ties that do share a node. The log odds results through

$$\begin{aligned} \text {Log odds}= & {} \log \left( \frac{P\left( y_{1,2} = 1, y_{1,3}, y_{rest}\right) }{P\left( y_{1,2} = 0, y_{1,3}, y_{rest} \right) } \right) \\= & {} \frac{1}{2} \theta _1 \left( \log \left( y_{1,3} + \sum \limits _{\begin{array}{c} k \ne 2,3 \end{array}} y_{1,k} + \sum \limits _{\begin{array}{c} k \ne 1,2 \end{array}} y_{2,k} \right) \right. \\&\left. +\, \sum \limits _{\begin{array}{c} j>2 \end{array}} y_{1,j} \left[ \log \left( 1 + y_{1,3} + \underbrace{\sum \limits _{\begin{array}{c} k \ne 2,3,j \end{array}} y_{1,k} + \sum \limits _{\begin{array}{c} k \ne j,1 \end{array}} y_{k,j}}_{c_{1,j}(N)} \right) \right. \right. \\&\left. \left. -\log \left( y_{1,3} + \underbrace{\sum _{k \ne 2,3,j} y_{1,k} + \sum _{k \ne j,1} y_{k,j}}_{c_{1,j}(N)} \right) \right] \right) \\&+\,\frac{1}{3} \theta _2 \left( \log \left( y_{1,3} y_{2,3} + \sum _{k \ne 1,2,3} y_{1,k} y_{2,k} \right) \right. \\&\left. +\, \sum _{j> 2} y_{1,j} \log \left( 1 \cdot y_{j,2} + y_{1,3} y_{j,3} + \underbrace{\sum _{k \ne 1,j,2,3} y_{1,k} y_{j,k}}_{d_{1,j}(N)} \right) \right. \\&\left. - \,\sum _{j>2} y_{1,j} \log \left( y_{1,3} y_{j,3} + \underbrace{\sum _{k \ne 1,j,2,3} y_{1,k} y_{j,k}}_{d_{1,j}(N)} \right) \right) \\= & {} \frac{1}{2} \theta _1 \left( \log \left( y_{1,3} + \sum \limits _{\begin{array}{c} k \ne 2,3 \end{array}} y_{1,k} + \sum \limits _{\begin{array}{c} k \ne 1,2 \end{array}} y_{2,k} \right) \right. \\&\left. + \,\sum \limits _{\begin{array}{c} j>2 \end{array}} y_{1,j} \Big [ \log \left( 1 + y_{1,3} + c_{1,j}(N) \right) -\log \left( y_{1,3} + c_{1,j}(N) \right) \Big ] \right) \\&+\frac{1}{3} \theta _3 \left( \log \left( y_{1,3} y_{2,3} + \sum _{k \ne 1,2,3} y_{1,k} y_{2,k} \right) \right. \\&\left. + \,\sum _{j> 2} y_{1,j} \log \left( 1 \cdot y_{j,2} + y_{1,3} y_{j,3} + d_{1,j}(N) \right) \right. \\&\left. -\, \sum _{j>2} y_{1,j} \log \left( y_{1,3} y_{j,3} + d_{1,j}(N) \right) \right) \end{aligned}$$

Hence, the log odds ratio for ties sharing a node results through

$$\begin{aligned}&\text {Log odds (B)} \\&\quad = \log \left( \frac{P\left( y_{1,2} = 1, y_{1,3} = 1, y_{rest}\right) P\left( y_{1,2} = 0, y_{1,3} = 0, y_{rest} \right) }{ P\left( y_{1,2} = 0, y_{1,3} = 1, y_{rest} \right) P\left( y_{1,2} = 1, y_{1,3} = 0, y_{rest} \right) } \right) \\&\quad = \frac{1}{2} \theta _1 \left( \log \left( 1 + \sum _{k \ne 2,3} y_{1,k} + \sum _{k \ne 1,2} y_{2,k} \right) - \log \left( \sum _{k \ne 2,3} y_{1,k} + \sum _{k \ne 1, 2} y_{2,k} \right) \right. \\&\quad \quad \left. +\,\sum _{j>2} y_{1,j} \Big [ \log \left( 2 + c_{1,j}(N)\right) - \log \left( c_{1,j}(N) \right) \Big ] \right) \\&\quad \quad + \,\frac{1}{3} \theta _2 \left( \log \left( 1 \cdot y_{2,3} +\sum _{k \ne 1,2,3} y_{1,k} y_{2,k} \right) - \log \left( \sum _{k \ne 1,2,3}y_{1,k} y_{2,k}\right) \right. \\&\qquad \left. +\, \sum _{j>2} y_{1,j} \Big [ \log \left( 1 \cdot y_{j,2} + 1 \cdot y_{j,3} + d_{1,j}(N) \right) - \log \left( 1 \cdot y_{j,3} + d_{1,j}(N) \right) \right. \\&\qquad \left. -\, \log \left( 1 \cdot y_{j,2} + d_{1,j}(N) \right) + \log \left( d_{1,j}(N) \right) \Big ] \right) \end{aligned}$$

The first component in the bracketed term at $\theta _2$ is of order $O_p(N^{-1})$ since

$$\begin{aligned} \log \left( 1\cdot y_{2,3} +\sum _{k \ne 1,2,3} y_{1,k} y_{2,k} \right) = \log \left( \sum _{k \ne 1,2,3} y_{1,k} y_{2,k} \right) + \underbrace{\frac{1}{\sum _{k \ne 1,2,3} y_{1,k} y_{2,k} } y_{2,3}}_{O_p(N^{-1})} \end{aligned}$$

The second component in the bracketed term is however of order O(1), which follows since the components in the sum are of order $O_p(N^{-1})$ but the outer sum is of order $O_p(N)$ so that a component of order O(1) results. The bracketed term at $\theta _1$ shows the same behaviour.

Hence we can conclude with (A) and (B): The conditional log odds ratio results to

$$\begin{aligned}&\log \left( \frac{P\left( y_{i,j} = 1, y_{l,m} = 1 | y_{rest}\right) P\left( y_{i,j} = 0, y_{l,m} = 0 | y_{rest} \right) }{P\left( y_{i,j} = 0, y_{l,m} = 1| y_{rest} \right) P\left( y_{i,j} = 1, y_{l,m} = 0| y_{rest} \right) }\right) \\&\quad ={\left\{ \begin{array}{ll} O_p(N^{-1}) &{} \text {if}\;\; \{i,j\} \ne \{l,m\}\\ O_p(1) &{} \text {otherwise}. \end{array}\right. } \end{aligned}$$

In other words, we obtain asymptotically a Markovian independence structure, which allows to draw parallel draws.

Appendix B: Diagnostic plots

Different diagnostic plots show the quality of the proposed algorithms. We provide traceplots for two parameter constellations, resulting in a sparse and middle dense network with 2200 nodes. Figure 8 shows that the mixing in pergm is about the same as mixing for ergm with the random option of selecting ties. The TNT sampler of the ergm-Package reaches the target distribution faster. The contrary can be seen in Fig. 9 if the network is denser.

Fig. 8

Traceplot for a simulated network with edges, 2-stars and triangles parameter $\theta = (-0.1, -0.01, 0.02)$ for the different algorithms. These parameters correspond to a sparse network with an approximate density of 0.06

Fig. 9

Traceplot for a simulated network with edges, 2-stars and triangles parameter $\theta = (6, -0.01, 0.03)$ for the different algorithms. These parameters correspond to a middle dense network with an approximate density of 0.26

To assess adequate mixing, we simulate 1000 networks (iterations) per algorithm, each network consists of $N=2200$ nodes and $10 \cdot N^2$ update steps after a reasonable burnin. We provide the density, traceplot, running mean and autocorrelation of the resulting edge, 2-star and triangle counts (Figs. 10, 11, 12, 13).

Fig. 10

Density plots with parameter distribution of edge, 2-star and triangle counts of different algorithms

Fig. 11

Traceplots of edge, 2-star and triangle counts of different algorithms

Fig. 12

Running means per edge, 2-star and triangle counts and algorithms

Fig. 13

Autocorrelation per edge, 2-star and triangle counts and algorithms