Truncated Simulation and Inference in
Edge-Exchangeable Networks

A completely random measure (CRM) $\Theta$ on $\Psi$ is a random measure such that for any collection of $K\in\mathbb{N}$ disjoint measurable sets $A_{1},...,A_{K}\subset\Psi$, the values $\Theta(A_{1}),...,\Theta(A_{K})$ are independent random variables [15]. In this work, we focus on discrete CRMs taking the form $\Theta=\sum_{k}\theta_{k}\delta_{\psi_{k}}$, where $\delta_{x}$ is a Dirac measure on $\Psi$ at location $x\in\Psi$ (i.e., $\delta_{x}(A)=1$ if $x\in A$ and 0 otherwise), and $(\theta_{k},\psi_{k})_{k=1}^{\infty}$ are a sequence of rates $\theta_{k}$ and labels $\psi_{k}$ generated from a Poisson process on $\mathbb{R}_{+}\times\Psi$ with mean measure $\nu(\mathrm{d}\theta)\times L(\mathrm{d}\psi)$. Here $L$ is a diffuse probability measure, and $\nu$ is a $\sigma$-finite measure satisfying $\nu(\mathbb{R}_{+})=\infty$, which guarantees that the Poisson process has countably infinitely many points almost surely. The space $\Psi$ and distribution $L$ will not affect our analysis; thus as a shorthand, we write $\mathrm{CRM}(\nu)$ for the distribution of $\Theta$:

One can construct a multidimensional measure $\Theta^{(d)}$ on $\Psi^{d}$, $d\in\mathbb{N}$ from $\Theta$ defined in Eq. 2 by taking its self-product. In particular, we define

where $i$ is a $d$-dimensional multi-index, and $\mathbb{N}_{\neq}^{d}$ is the set of such indices with all distinct components. Note that $\Theta^{(d)}$ is no longer a CRM, as it does not satisfy the independence condition.

To simulate a realization $\Theta\sim\mathrm{CRM}(\nu)$—e.g., as a first step in the simulation of a self-product measure $\Theta^{(d)}$—the rates $\theta_{k}$ and labels $\psi_{k}$ may be generated in sequence using a series representation [39] of the CRM. In particular, we begin by simulating the jumps of a unit-rate homogeneous Poisson process $(\Gamma_{k})_{k=1}^{\infty}$ on $\mathbb{R}_{+}$ in increasing order. For a given distribution $g$ on $\mathbb{R}_{+}$ and nonnegative measurable function $\tau:\mathbb{R}_{+}\times\mathbb{R}_{+}\to\mathbb{R}_{+}$, we set

Depending on the particular choice of $\tau$ and $g$, one can construct several different series representations of a CRM [28]. For example, the inverse Lévy representation [35] has the form

In many cases, computing $\nu^{\leftarrow}(x)$ is intractable, making it hard to generate $\theta_{k}$ in this manner. Alternatively, we can generate a series of rates from $\mathrm{CRM}(\nu)$ with the rejection representation [34], which has the form

where $\mu$ is a measure on $\mathbb{R}_{+}$ chosen such that $\frac{\mathrm{d}\nu}{\mathrm{d}\mu}\leq 1$ uniformly and $\mu^{\leftarrow}(x)$ is easy to calculate in closed-form. While there are many other sequential representations of CRMs [28], the representations in Eqs. 4, 5 and 6 are broadly applicable and play a key role in our theoretical analysis.

Self-product measures $\Theta^{(d)}$ of the form Eq. 3 with $d=2$ have recently been used as priors in a number of probabilistic network models [3, 4].¹¹1There are also network models based on self-product measure priors that do not generate edge-exchangeable sequences [8, 11]; it is likely that many of the techniques in the present work would extend to these models, but we leave the study of this to future work. The focus of the present work are those models that associate each $\psi_{k}$ with a vertex, each tuple $\zeta_{i}=(\psi_{i_{1}},\dots,\psi_{i_{d}})$ with a (hyper)edge, and then build a growing sequence of networks by sequentially generating edges from $\Theta^{(d)}$ in rounds $n=1,\dots,N$. There are a number of choices for how to construct such a sequence. For example, in each round $n$ we may add multiple edges via an independent likelihood process $X_{n}\sim\mathrm{LP}(h,\Theta^{(d)})$ defined by

where $x_{ni}=k$ denotes that there were $k$ copies of edge $\zeta_{i}$ added at round $n$, and $h(\cdot|\vartheta)$ is a probability distribution on $\mathbb{N}\cup\{0\}$. We denote the mean $\mu(\vartheta):=\sum_{k=0}^{\infty}k\cdot h(k\,|\,\vartheta)$ and probability of 0 under $h$ to be $\pi(\vartheta):=h(0|\vartheta)$ for convenience. By the Slivnyak-Mecke theorem [40], if $h$ satisfies

then finitely many edges are added to the graph in each round. We make this assumption throughout this work. Alternatively, if

then $\Omega:=\Theta^{(d)}(\Psi^{d})<\infty$, and we may add only a single edge per round $n$ via a categorical likelihood process $X_{n}\sim\mathrm{Categorical}(\Theta^{(d)})$ defined by

This construction has appeared in [4], where $\Theta$ follows a Dirichlet process, which can be seen as a normalized gamma process [41]. Note that our definition of $\mathrm{Categorical}(\Theta^{(d)})$ involves simulating from the normalization; we use this definition to avoid introducing new notation for normalized processes.

Using either likelihood process, the edges in the network after $N$ rounds are

i.e., $x_{i}\in\mathbb{N}\cup\{0\}$ represents the count of edge $i$ after $N$ rounds.

There are three points to note about this formulation. First, since the atom locations $\zeta_{i}$ are not used, we can represent the network using only its array of edge counts

where $\mathcal{N}_{d}$ denotes the set of integer arrays indexed by $\mathbb{N}_{\neq}^{d}$ with finitely many nonzero entries. Note that $\mathcal{N}_{d}$ is a countable set. Second, by construction, the distribution of $E_{N}$ is invariant to reorderings of the arrival of edges, and thus the network is edge-exchangeable [3, 5, 6, 7]. Finally, note that the network $E_{N}$ as formulated in Eq. 12 is in general a directed multigraph with no self-loops (due to the restriction to indices $i\in\mathbb{N}_{\neq}^{d}$ rather than $\mathbb{N}^{d}$). Although the main theoretical results in this work are developed in this setting, we provide an additional result in Section 3.1 to translate to other common network structures (e.g. binary undirected networks).

We formulate the approximation error incurred by truncation as the total variation distance between the marginal network distributions, i.e., of $E_{N}$ and $E_{N,K}$. The first step in the analysis of truncation error—provided by Lemma 3.1—is to show that this is bounded above by the probability that there are edges in the full network $E_{N}$ involving vertices beyond the truncation level $K$. To this end, we denote the maximum vertex index of $E_{N}$ to be

and note that by definition, $I_{N}\leq K$ if and only if all edges in $E_{N}$ fall in the truncated region.

As mentioned in Section 2.3, the network $E_{N}$ is in general a directed multigraph with no self loops. However, Lemma 3.1—and the downstream truncation error bounds presented in Sections 3.2 and 3.3—also apply to any graph $E^{\prime}_{N}=(x^{\prime}_{i})_{i\in\mathbb{N}^{d}_{\neq}}$ that is a function of the original graph $E^{\prime}_{N}=f(E_{N})$ such that truncation commutes with the function, i.e., $E^{\prime}_{N,K}=f(E_{N,K})$. For example, to obtain a truncation error bound for the common setting of undirected binary graphs, we generate the directed multigraph $E_{N}$ as above and construct the undirected binary graph $E^{\prime}_{N}$ via

Corollary 3.2 provides the precise statement of the result; note that the bound is identical to that from Lemma 3.1.

We now specialize Lemma 3.1 to the setting where $\Theta$ is a CRM generated by a series representation of the form Eq. 4, and the network is generated via the independent likelihood process from Eq. 7. As a first step towards a bound on the truncation error for general hypergraphs with $d>1$ in Theorem 3.4, we present a simpler corollary in the case where $d=2$, which is of direct interest in analyzing the truncation error of popular edge-exchangeable networks.

The proof of this result (and Theorem 3.4 below) in Appendix B follows by representing the tail of the CRM as a unit-rate Poisson process on $[\Gamma_{K},\infty)$. Though perhaps complicated at first glance, an intuitive interpretation of the truncation error terms $B_{K,i}$ is provided by Fig. 1(b). $B_{K,1}$ corresponds to the truncation error arising from the upper right quadrant, where both vertices participating in the edge were in the discarded tail region. $B_{K,2}$ is the truncation error arising from the bottom right and upper left quadrants, where one of the two vertices participating in the edge was in the truncation, and the other was in the tail. Finally, $B_{K,3}$ represents the truncation error arising from edges in which one vertex was at the boundary of tail and truncation, and the other was in the tail.

Theorem 3.4 is the generalization of Corollary 3.3 from $d=2$ to the general setting of arbitrary hypergraphs with $d>1$. The bound is analogous to that in Corollary 3.3—with $B_{K}$ expressed as a sum of terms, each corresponding to whether vertices were in the tail, boundary, or truncation region—except that there are $d>1$ vertices participating in each edge, resulting in more terms in the sum. Note that Theorem 3.4 also guarantees that the bound is not vacuous, and indeed converges to 0 as the truncation level $K\rightarrow\infty$ as expected.

The same geometric intuition from the $d=2$-dimensional case applies to the general hypergraph truncation error in Eq. 26. $B_{K,1}$ corresponds to the error arising from the edges whose vertices all belong to the tail region $\Theta_{K+}$. Each term in the summation in $B_{K,2}$ corresponds to the error arising from edges that have $\ell$ out of $d$ vertices belonging to the truncation $\Theta_{K}$. Each term in the summation in $B_{K,3}$ corresponds to the error arising from the edges that have $\ell-1$ out of $d$ vertices belonging to the truncation $\Theta_{K}$ and have one vertex exactly on the boundary of the truncation. Note that we obtain Corollary 3.3 by taking $d=2$ in Eq. 26.

We may also specialize Lemma 3.1 to the setting where the network is generated via the single-edge-per-step categorical likelihood process in Eq. 10. However, truncation with the categorical likelihood process poses a few key challenges. From a practical angle, certain choices of series representation for generating $\Theta$ may be problematic. For instance, when using the rejection representation Eq. 6 in the typical case where $\mu\neq\nu$, there is a nonzero probability that

meaning there aren’t enough accepted vertices in the truncation to generate a single edge. In this case, the categorical likelihood process—which must generate exactly one edge per step—is ill-defined. An additional theoretical challenge arises from the normalization of the original and truncated networks in Eq. 10, which prevents the use of the usual theoretical tools for analyzing CRMs.

Fortunately, the inverse Lévy representation provides an avenue to address both issues. The rates $\theta_{k}$ are all guaranteed to be nonzero—meaning as long as $K\geq d$, the categorical likelihood process is well-defined—and are decreasing, which enables our theoretical analysis in Section B.1. However, as mentioned earlier, the inverse Lévy representation is well-known to be intractable to use in most practical settings.

Theorem 3.5 provides a solution: we use the rejection representation to simulate the rates $\theta_{k}$, but instead of simulating for iterations $k=1,\dots,K$, we simulate until we obtain $K$ nonzero rates. This is no longer a sample of a truncated rejection representation; but Theorem 3.5 shows that the first $K$ nonzero rates have the same distribution as simulating $K$ iterations of the inverse Lévy representation. Therefore, we can tailor the analysis of truncation error for the categorical likelihood process in Theorem 3.6 to the inverse Lévy representation, and simulate its truncation for any $K$ using the tractable rejection representation in practice.

We now apply the results of this section to binary undirected networks ($d=2$) constructed via Eq. 17 from common edge-exchangeable networks [3]. We derive the convergence rate of truncation error, and provide simulations of the expected number of edges and vertices under the infinite and truncated network. In each simulation we use the rejection representation to simulate $\Theta$, and run $N=10,000$ steps of network construction.

We begin by examining the density of the posterior distribution of the $K^{\text{th}}$ rate from the inverse Lévy representation $\theta_{K}$ for some fixed $K\in\mathbb{N}$, the unordered collection of rates $(\theta_{k})_{k=1}^{K-1}$ such that $\theta_{k}\geq\theta_{K}$, and the parameters $\sigma\in\mathbb{R}^{m}$ of the CRM rate measure $\nu_{\sigma}$ given the observed set of edges $E_{N}$. We denote $\nu_{\sigma}(\theta)$ to be the density of $\nu_{\sigma}(\mathrm{d}\theta)$ and $p(\sigma)$ to be the prior density of $\sigma$, both with respect to the Lebesgue measure. Given these definitions the posterior density can be expressed as

where $[K]^{d}_{\neq}$ is the subset of $\mathbb{N}^{d}_{\neq}$ such that $\max_{j\in[d]}i_{j}\leq K$, and $x_{1:N}^{K+}$ is shorthand for the set of $(x_{ni})_{n\in[N],i\in\mathbb{N}_{\neq}^{d}}$ such that $i\notin[K]^{d}_{\neq}$. All the factors on the first row correspond to the prior of $\sigma$ and $(\theta_{k})_{k=1}^{K}$, and the first factor on the second row corresponds to the likelihood of the portion of the network involving only vertices $1\dots K$; these are straightforward to evaluate, though $\nu_{\sigma}[\theta_{K},\infty)$ will occasionally require a 1-dimensional numerical integration. The last factor corresponds to the likelihood of the portion of the network involving vertices beyond $K$, and is not generally possible to evaluate exactly.

To handle this term, suppose that $K$ is large enough that $x_{1:N}^{K+}=0$. Then $p(x_{1:N}^{K+}|\theta_{1:K},\sigma)=\mathbb{P}\left(I_{N}\leq K\,|\,\Gamma_{1:K},\sigma\right)$, i.e., the probability that no portion of the network involves vertices of index $>K$. Theorem 4.1 provides upper and lower bounds on this probability akin to those of Theorem 3.4—indeed, Theorem 4.1 is an intermediate step in the proof of Theorem 3.4—that apply conditionally on the values of $U_{1:K}$, $\Gamma_{1:K}$ rather than marginally. This theorem also makes the dependence of the bound on the rate measure parameters $\sigma$ notationally explicit via parametrized series representation components $\tau_{\sigma}$ and $g_{\sigma}$ from Eq. 4. Finally, though Theorem 4.1 applies to general series representations, we require only the specific instantiation for the inverse Lévy representation in the present context.

Since $U_{1:K}$ is unused in the inverse Lévy representation, in the present context we use the notation $B(\Gamma_{1:K},\sigma)$ for brevity. The bound in Theorem 4.1 implies that as long as $K$ is set large enough such that both $x_{1:N}^{K+}=0$ and $B(\Gamma_{1:K},\sigma)\leq\epsilon/N$ for some $\epsilon>0$ then

Therefore as $K$ increases, this term should become $\approx 1$ and have a negligible effect on the posterior density. We use this intuition to propose a truncated Metropolis–Hastings algorithm that sets $K>I_{N}$, ignores the $p(x_{1:N}^{K+}|\theta_{1:K},\sigma)$ term in the acceptance ratio, and fixes $x_{N}^{K+}$ to 0. Theorem 4.2 provides a rigorous analysis of the error involved in using the truncated sampler.

Crucially, as long as one can obtain samples from the truncated posterior distribution, one can estimate the bound in Eq. 61 using sample estimates of the tail probability $\hat{\Pi}\left\{B(\Gamma_{1:K},\sigma)\leq\epsilon/N\right\}\geq 1-\eta$. Therefore, one can compute the bound Eq. 61 without needing to evaluate $p(x_{1:N}^{K+}|\theta_{1:K},\sigma)$ exactly. This suggests the following adaptive truncation procedure:

1. 1.

   Pick a value of $K>I_{N}$ and desired total variation error $\xi\in(0,1)$.

2. 2.

   Obtain samples from the truncated posterior.

3. 3.

   Minimize the bound in Eq. 61 over $\epsilon\in(0,1)$, using the samples to estimate $\eta=1-\hat{\Pi}\left\{B(\Gamma_{1:K},\sigma)\leq\epsilon/N\right\}$ for each value of $\epsilon$.

4. 4.

   If the minimum bound exceeds $\xi$, increase $K$ and return to 2.

5. 5.

   Otherwise, return the samples.

In this work, we start by initializing $K$ to $I_{N}+1$. In order to decide how much to increase $K$ by in each iteration, we use the sequential representation to extrapolate the total variation bound Eq. 61 to larger values of $K$ without actually performing MCMC sampling. In particular, for each posterior sample, we use its hyperparameters to generate additional rates from the sequential representation. We then use these extended samples to compute the total variation error guarantee as per step 3 above. We continue to generate additional rates (typically doubling the number each time) until the predicted total variation guarantee is below the desired threshold. Finally, we use linear interpolation between the last two predicted errors to find the next truncation level $K$ that matches the desired (log) error threshold. This fixes a new value of $K$; at this point, we return to step 2 above and iterate.

In this section, we examine the properties of the proposed adaptive truncated inference algorithm for the beta-independent Bernoulli network model in Eqs. 35 and 36 with discount $\alpha\in(0,1)$, concentration $\lambda>1$, mass $\gamma>0$, unordered collection of rates $\theta_{1:K-1}$, and $K^{\text{th}}$ rate from the sequential representation $\theta_{K}$. In order to simplify inference, we transform each of these parameters to an unconstrained version: